Metadata-Version: 2.4
Name: siemensfile
Version: 0.3.1
Summary: Python package to read Siemens MRI raw data files (.dat) and reconstruct MR images.
Home-page: https://github.com/cenarius1985/SIEMENSFile
Author: Fernando Jose Ramirez
Author-email: fernando.ramirez.sarmiento@gmail.com
Classifier: Programming Language :: Python :: 3
Classifier: License :: OSI Approved :: MIT License
Classifier: Operating System :: OS Independent
Requires-Python: >=3.12.4
Description-Content-Type: text/markdown
License-File: LICENSE
Requires-Dist: numpy
Requires-Dist: scipy
Requires-Dist: pandas
Requires-Dist: matplotlib
Requires-Dist: mri-nufft
Requires-Dist: pydicom
Requires-Dist: openpyxl
Requires-Dist: h5py
Dynamic: author
Dynamic: author-email
Dynamic: classifier
Dynamic: description
Dynamic: description-content-type
Dynamic: home-page
Dynamic: license-file
Dynamic: requires-dist
Dynamic: requires-python
Dynamic: summary

# 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:

```bash
pip install siemensfile
```

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

```bash
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

```python
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

```python
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

```python
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

```bash
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](https://github.com/pehses/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.
