SIEMENSFile

SIEMENSFile is a Python package for previewing, reading and processing Siemens .dat raw MRI files, and for reconstructing MR images. It extracts the raw k-space data and the scan metadata, and performs image reconstruction with centered inverse Fourier transforms. The Cartesian reconstruction pipeline supports 2D multi-slice, 2D multi-stack (e.g. sag+cor+tra localizers) and 3D acquisitions (e.g. MPRAGE), including Partial Fourier handling. Non-Cartesian reconstruction (radial/spiral via NUFFT) is planned.

Features

• Reads Siemens VD/VE raw data files (.dat, twix format) through the bundled twixtools.
• Cartesian reconstruction of:
  • 3D acquisitions: each partition (cPar) is treated as an individual slice — a 3D MPRAGE volume yields one image per partition.
  • 2D multi-stack acquisitions: stacks sharing the same slice counter are separated by spatial position from the acquisition header.
  • Partial Fourier 6/8: phase lines are mapped onto the full acquisition matrix using CenterLin; missing lines are zero-filled.
  • Partial Fourier 1/2 (interleaved): odd phase lines are interpolated from their even neighbours to remove the Nyquist ghost produced by zero-filling.
  • Readout centering: the echo peak (k-space center) is re-centered using CenterCol with an exact circular roll.
• Root-sum-of-squares (RMS) coil combination over all channels.
• Outputs per acquisition: k-space preview (PNG), reconstruction mosaic (PNG), one standard-compliant DICOM MR file per slice (MRImageStorage), and full metadata (JSON).
• Robust metadata extraction from the Siemens Phoenix protocol: TR/TE (ms), FOV, pixel spacing, slice thickness/spacing, 3D flag.

Installation

Install from PyPI:

pip install siemensfile

For a development setup with conda (see requirements.txt for pinned versions):

conda create -n siemensfile python=3.12 -y --override-channels -c conda-forge
conda activate siemensfile
pip install -r requirements.txt
pip install -e .

Requires Python 3.12+.

Usage

from siemensfile import siemensfile

metadata, kspace = siemensfile(r"path/to/meas_MID00068_FID09111_t1_mprage_tra.dat",
                               reconstruction="Cartesian")

kspace.shape   # [line, channel, column, slice] (complex k-space)

reconstruction accepts "Cartesian" (implemented) or "NonCartesian" (not yet implemented — raises NotImplementedError).

API change in 0.2.0: the keyword argument and its values were renamed from reconstruccion="Cartesiana"/"NoCartesiana" to reconstruction="Cartesian"/"NonCartesian", and all module-level functions were renamed to English (e.g. lectura_twix → read_twix_pipeline, extraer_metadata_recursivamente → extract_metadata_recursively).

Outputs

Results are written to an ordered, timestamped layout next to the input .dat file (each run creates its own job folder — nothing is deleted, so previous runs are preserved):

output/
└── 2026-09-13/                                  # date
    └── 214831_t1_mprage_tra/                    # HHMMSS_<base_name> (job folder)
        ├── metadata.json                        # flattened twix headers
        ├── <base>_kspace.png                    # k-space preview (log scale)
        ├── <base>_reconstruction.png            # reconstruction mosaic (all slices)
        └── Cartesian/                           # one folder per reconstruction method
            ├── png/magnitude/
            │   ├── slice_001.png                # magnitude slice 1 (windowed)
            │   └── slice_NNN.png
            ├── dcm/magnitude/
            │   ├── slice_001.dcm                # DICOM MR (MRImageStorage)
            │   └── slice_NNN.dcm
            └── numpy/magnitude/
                └── image_magnitude.npy          # magnitude volume [slice, line, column]

All slice images are magnitude images (root-sum-of-squares over coils). The PNG and DICOM sets share a consistent per-volume percentile window so slices can be compared.

Example: 3D MPRAGE

metadata, kspace = siemensfile(r"meas_MID00068_FID09111_t1_mprage_tra.dat")
kspace.shape    # (224, 15, 352, 102) -> 102 axial partitions reconstructed

Example: multi-stack localizer

metadata, kspace = siemensfile(r"meas_MID00062_FID09105_localizer_sag+cor+tra.dat")
kspace.shape    # (288, 15, 512, 26) -> 26 slices (8 cor + 8 sag + 10 tra), ghost-free

Mathematics & physics of the reconstruction

This section documents every equation the pipeline applies, so you know exactly what the package computes. Symbols: ρ = proton-density (spin) magnetization, x = spatial position, k = k-space position, γ = gyromagnetic ratio, G = gradient, Nc = number of coils, s = slice/partition index.

1. Physics: what the scanner measures

Gradient-encoding MRI samples the spatial Fourier transform of the object. During readout, the received signal is

S(t) = ∫ ρ(x) · exp( −i·2π·k(t)·x ) dx         with  k(t) = (γ/2π) ∫₀ᵗ G(τ) dτ

so k-space is the Fourier transform of the image: reconstructing means computing an inverse Fourier transform. Each line of the raw file is one phase-encoded acquisition S(cLin, cPar)[channel, kx].

2. Acquisition → k-space matrix mapping (io.twix_reader)

The physical phase line is placed in the full matrix using the header's k-space center (CenterLin):

row = cLin + ( N/2 − CenterLin )

N = lPhaseEncodingLines                      (from the Phoenix protocol), or
N = 2 · max( CenterLin, cLin_max + 1 − CenterLin )   (inferred)

Lines not acquired (Partial Fourier) remain zero (zero-filled reconstruction — a low-pass approximation).

3. Slice identification (io.twix_reader.get_slice_key)

Acquisitions are grouped into slices by their spatial position and counters:

key = ( Sag, Cor, Tra, cSli, cPar )        # SliceData.SlicePos + Counter

• 2D multi-slice: slices differ in cSli (position differs too).
• 2D multi-stack (localizer sag+cor+tra): cSli repeats per stack; the position separates stacks.
• 3D (MPRAGE): all partitions share the slab position; each partition is distinguished by cPar and reconstructed as its own slice.

4. Averages (domain.kspace)

Repeated acquisitions of the same (slice, line) are summed:

S[key, row] = Σₐ  data_a           (Nₐ averages → SNR gain ∝ √Nₐ)

5. Partial Fourier 6/8 (zero-filled)

Only the lower 6/8 of the phase lines were acquired (e.g. 196 of 224). With the CenterLin mapping the acquired band sits in its true place and the missing outer band stays zero:

S_full = [ 0 … 0 | S_acquired | 0 ]          →  slightly blurred magnitude

(Homodyne reconstruction would sharpen this; not implemented yet.)

6. Partial Fourier 1/2 — odd-line interpolation (domain.kspace)

Localizers acquire only the even phase lines. Zero-filling that checkerboard creates a Nyquist ghost (duplicated anatomy); SIEMENSFile fills the odd lines by linear interpolation:

S[2j+1] = ( S[2j] + S[2j+2] ) / 2            S[last odd] = S[last even]

7. Readout centering (domain.kspace)

The echo (k-space center in the readout direction) sits at sample CenterCol, not necessarily at N_col/2. Because the readout is fully sampled, a circular roll is exact:

S ← roll( S, N_col/2 − CenterCol, axis = kx )

8. Centered inverse FFT — the reconstruction itself (domain.strategies)

For every slice, the image is the centered 2D inverse DFT over the phase (line) and readout (column) axes, with ortho normalization (Parseval: energy is preserved between k-space and image):

ρ(x, y) = FFTSHIFT{ IFFT₂⁰ [ IFFTSHIFT( S(ky, kx) ) ] }        (norm = "ortho")

9. 3D acquisitions (MPRAGE)

A 3D volume is a stack of partitions; each partition is an independent 2D image through (phase ky, readout kx):

ρₛ(x, y) = IFFT₂{ S[:, :, s] }          s = 1 … N_partitions

10. Coil combination — root sum of squares (domain.strategies)

Each coil channel c is reconstructed separately and combined by root-sum-of-squares (RMS):

ρ(x, y) = √( (1/Nc) · Σ𝒸 |ρ𝒸(x, y)|² )

This is phase-insensitive (no coil phases needed) but keeps the coil sensitivity profile — see Known limitations.

11. DICOM percentile windowing (io.dicom_writer)

Each slice is scaled to uint16 with a percentile window (P0.5–P99.5 of that slice); plain min–max would crush the contrast in the presence of a few hyperintense pixels:

DICOM = clip( (ρ − P₀.₅) / (P₉₉.₅ − P₀.₅), 0, 1 ) · 65535

12. Geometry from the protocol (domain.metadata)

Δx = FOV_readout / N_columns          Δy = FOV_phase / N_lines
slice thickness (2D)  = dThickness
slice thickness (3D)  = dThickness(slab) / lImagesPerSlab
TR, TE = alTR[0], alTE[0]             (µs → ms: ÷ 1000)

Tests

python -m pytest tests/ -q

Architecture

Since 0.3.0 the package is organized in layers with a strict dependency rule (api → pipeline → {io, domain, visualization}; the domain layer is pure numpy/stdlib):

siemensfile/
├── siemensfile.py     # public facade: siemensfile(), read_twix_pipeline()
├── pipeline.py        # orchestration + logging
├── domain/            # pure numeric core (no I/O)
│   ├── data.py        # dataclasses: Acquisition, KSpaceData, SeriesMetadata, ReconstructionResult
│   ├── kspace.py      # k-space construction, Partial Fourier, readout centering
│   ├── strategies.py  # Strategy pattern + registry (Cartesian; NonCartesian stub)
│   └── metadata.py    # SeriesMetadata from the Phoenix protocol
├── io/                # adapters: twix_reader, dicom_writer, exporters, ismrmrd_export
├── visualization/     # PNG previews
└── utils.py, reconstruction.py, ismrmrd_format.py   # deprecated 0.2.x compatibility shims

Adding a new reconstruction method means implementing a strategy class and registering it in domain.strategies.RECONSTRUCTION_REGISTRY — the pipeline stays untouched.

Known limitations

• No coil intensity inhomogeneity correction: sum-of-squares images keep the coil sensitivity profile (brighter periphery than the vendor reconstruction).
• Partial Fourier 6/8 is zero-filled (slight blur); homodyne reconstruction is not implemented.
• No parallel imaging (GRAPPA/R) support.
• No non-Cartesian (radial/spiral) reconstruction yet.
• ismrmrd_formato.py (ISMRMRD HDF5 export) is experimental and untested.

Credits and acknowledgements

This project builds on twixtools by Philipp Ehses, bundled under src/twixtools, which provides the core Siemens .dat reading functionality. This package extends it with image reconstruction. If you need a more complete tool for reading/writing Siemens raw data, check twixtools directly.