diff --git a/CHANGELOG.md b/CHANGELOG.md index cee741f..b22b007 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -1,3 +1,7 @@ +- 1.3.0 + - use the official [`bifti`](https://pypi.org/project/bifti/) package for BIfTI phantoms instead of the internal implementation (#201). `NiftiPhantom`, `NiftiTissue` and `ResliceConfig` are now aliases of `bifti.BiftiPhantom`, `bifti.BiftiTissue` and `bifti.ResliceTo`; `TissueDict.load` uses the `bifti` loader (density-weighted resampling, `patient` position support, complex B1 maps, `func` mappings are parsed instead of `eval`ed) + - `TissueDict.save` uses `bifti.NumpyPhantom.save`: NIfTIs are named `_.nii.gz` (e.g. `phantom_B1plus_0.nii.gz`) and stored as float64 / complex128; the `.json` no longer contains `reslice_to` + - requires Python >= 3.12 (required by `bifti`) - 1.2.0 - fix: phantoms with zero size along an axis (2D, single voxel) no longer crash the prepass with `Encountered a NaN while sorting + dists by mag` when simulating normalized-gradient sequences; such axes are not scaled - add `fid_graph`: builds a graph with only `z0` and one `+` state per repetition, as an alternative to `compute_graph` without running the prepass diff --git a/Cargo.lock b/Cargo.lock index 4597482..4bde265 100644 --- a/Cargo.lock +++ b/Cargo.lock @@ -67,7 +67,7 @@ dependencies = [ [[package]] name = "mrzero_core" -version = "1.2.0" +version = "1.3.0" dependencies = [ "num-complex", "pulseq-rs", diff --git a/Cargo.toml b/Cargo.toml index c5e18ff..46a2a11 100644 --- a/Cargo.toml +++ b/Cargo.toml @@ -1,6 +1,6 @@ [package] name = "mrzero_core" -version = "1.2.0" +version = "1.3.0" edition = "2021" # See more keys and their definitions at https://doc.rust-lang.org/cargo/reference/manifest.html diff --git a/documentation/src/api-nifti.md b/documentation/src/api-nifti.md index 265e494..806b1f4 100644 --- a/documentation/src/api-nifti.md +++ b/documentation/src/api-nifti.md @@ -6,6 +6,10 @@ You can load a complete phantom directly with [`mr0.TissueDict.load`](api-phanto *But you can also* load only the configuration with [`mr0.NiftiPhantom.load`](#niftiphantom). This allows you to modify it as needed before loading the phantom or saving the changes. +These classes are provided by the official [`bifti`](https://github.com/mrx-org/bifti-phantoms) package and re-exported by MRzero. +`NiftiPhantom`, `NiftiTissue` and `ResliceConfig` are aliases of `bifti.BiftiPhantom`, `bifti.BiftiTissue` and `bifti.ResliceTo`. +See the [bifti documentation](https://github.com/mrx-org/bifti-phantoms/blob/main/python/bifti/README.md) for the full reference. + All classes are [dataclasses](https://docs.python.org/3/library/dataclasses.html), augmented with additional functions. ## `PhantomUnits` @@ -118,7 +122,7 @@ where `Prop = float | NiftiRef | NiftiMapping` Represents the configuration given by a phantom `.json` file. ```python -NiftiPhantom(units, system, tissues) +NiftiPhantom(units, system, tissues, reslice_to=None, schema=NiftiPhantom.DEFAULT_SCHEMA, patient=None) ``` | Parameter | Description | @@ -126,8 +130,9 @@ NiftiPhantom(units, system, tissues) | `units` | [`PhantomUnits`](#phantomunits) object | | `system` | [`PhantomSystem`](#phantomsystem) object | | `tissues` | Dictionary of tissues; key specifies name, values are [`NiftiTissue`](#niftitissue) objects | - -The class attribute `schema_uri` is always `"https://github.com/mrx-org/bifti-phantoms/bifti-phantom-v1"` and is not a constructor parameter. +| `reslice_to` | Optional `ResliceConfig` (`affine`, `resolution`): target grid the tissues are resampled to on load | +| `schema` | Value of the `$schema` field; defaults to the schema hosted in the bifti-phantoms repository | +| `patient` | Optional `Patient` (`position`, e.g. `"HFS"`); omitted means `FFS` (no transformation) | - `NiftiPhantom.default(gyro=42.5764, B0=3)`: return a `NiftiPhantom` without tissues - `NiftiPhantom.load(path)`: load from `.json` file, given by the path diff --git a/documentation/src/api-phantom.md b/documentation/src/api-phantom.md index 8d4215a..bd784d3 100644 --- a/documentation/src/api-phantom.md +++ b/documentation/src/api-phantom.md @@ -103,7 +103,7 @@ TissueDict(dict[str, VoxelGridPhantom]) ``` - `TissueDict.load(path, config=None)`: either load from `.json` or from a dir + [`NiftiPhantom`](api-nifti.md#niftiphantom) -- `TissueDict.save(path, gyro, B0)`: save as [NIfTI phantom](nifti.md) +- `TissueDict.save(path, gyro, B0)`: save as [NIfTI phantom](nifti.md) (via `bifti.NumpyPhantom.save`; NIfTIs are named after the `.json`, nearly constant maps are stored as values) - `TissueDict.interpolate()`: apply [`VoxelGridPhantom.interpolate()] to all tissues - `TissueDict.slices()`: apply [`VoxelGridPhantom.slices()] to all tissues - `TissueDict.combine()`: merge all tissues into a single [`VoxelGridPhantom`](#voxelgridphantom) diff --git a/documentation/src/api.md b/documentation/src/api.md index 8ba37aa..15c1c76 100644 --- a/documentation/src/api.md +++ b/documentation/src/api.md @@ -27,6 +27,7 @@ import MRzeroCore as mr0 - [`mr0.NiftiMapping`](api-nifti.md#niftimapping) - [`mr0.NiftiTissue`](api-nifti.md#niftitissue) - [`mr0.NiftiPhantom`](api-nifti.md#niftiphantom) + - [`mr0.ResliceConfig`](api-nifti.md#niftiphantom) - [Simulation](api-simulation.md) - [`mr0.isochromat_sim`](api-simulation.md#isochromat_sim) - [`mr0.compute_graph`](api-simulation.md#compute_graph) diff --git a/pyproject.toml b/pyproject.toml index 5191343..0b0e7a4 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -4,7 +4,7 @@ build-backend = "maturin" [project] name = "mrzerocore" -version = "1.2.0" +version = "1.3.0" description = "Core functionality of MRzero" authors = [ {name = "Jonathan Endres", email = "jonathan.endres@uk-erlangen.de"}, @@ -16,8 +16,9 @@ classifiers = [ "Programming Language :: Python :: Implementation :: PyPy", "License :: OSI Approved :: GNU Affero General Public License v3", ] -requires-python = ">=3.10" +requires-python = ">=3.12" dependencies = [ + "bifti>=0.3.0", "ismrmrd", "matplotlib>=3.5", "pydisseqt>=0.1.13", diff --git a/python/MRzeroCore/__init__.py b/python/MRzeroCore/__init__.py index 6a5c97b..8302991 100644 --- a/python/MRzeroCore/__init__.py +++ b/python/MRzeroCore/__init__.py @@ -19,7 +19,10 @@ from .phantom.custom_voxel_phantom import CustomVoxelPhantom from .phantom.sim_data import SimData from .phantom.brainweb import generate_brainweb_phantoms -from .phantom.nifti_phantom import NiftiPhantom, NiftiTissue, NiftiRef, NiftiMapping, ResliceConfig +from .phantom.nifti_phantom import ( + NiftiPhantom, NiftiTissue, NiftiRef, NiftiMapping, ResliceConfig, + BiftiPhantom, BiftiTissue, ResliceTo, Patient, PhantomUnits, PhantomSystem, +) from .phantom.tissue_dict import TissueDict from .simulation.isochromat_sim import isochromat_sim from .simulation.pre_pass import compute_graph, compute_graph_ext, merge_graphs, Graph diff --git a/python/MRzeroCore/phantom/nifti_phantom.py b/python/MRzeroCore/phantom/nifti_phantom.py index 8e4f6a7..827bf36 100644 --- a/python/MRzeroCore/phantom/nifti_phantom.py +++ b/python/MRzeroCore/phantom/nifti_phantom.py @@ -1,233 +1,17 @@ -from dataclasses import dataclass, field -from typing import Literal, Any -from pathlib import Path - - -@dataclass -class PhantomUnits: - gyro: Literal["MHz/T"] - B0: Literal["T"] - T1: Literal["s"] - T2: Literal["s"] - T2dash: Literal["s"] - ADC: Literal["10^-3 mm^2/s"] - dB0: Literal["Hz"] - B1_tx: Literal["rel"] - B1_rx: Literal["rel"] - - @classmethod - def default(cls): - return cls( - gyro="MHz/T", - B0="T", - T1="s", - T2="s", - T2dash="s", - ADC="10^-3 mm^2/s", - dB0="Hz", - B1_tx="rel", - B1_rx="rel", - ) - - @classmethod - def from_dict(cls, config: dict[str, str]): - # Currently this is only for documentation and no other units than the - # default units are supported. This definetly can change in the future - # but the implementation has no priority right now - default = cls.default() - assert default.to_dict() == config, "Only default units are supported for now" - return default - - def to_dict(self) -> dict[str, str]: - return { - "gyro": self.gyro, - "B0": self.B0, - "T1": self.T1, - "T2": self.T2, - "T2'": self.T2dash, - "ADC": self.ADC, - "dB0": self.dB0, - "B1+": self.B1_tx, - "B1-": self.B1_rx, - } - - -@dataclass -class PhantomSystem: - gyro: float - B0: float - - @classmethod - def from_dict(cls, config: dict[str, float]): - return cls(**config) - - def to_dict(self) -> dict[str, float]: - return {"gyro": self.gyro, "B0": self.B0} - - -@dataclass -class ResliceConfig: - """Target resampling grid declared inside the phantom JSON.""" - resolution: list # [nx, ny, nz] — output voxel count - affine: list # 3×4 NIfTI-style affine in mm (rotation×voxel_size | origin) - - @classmethod - def from_dict(cls, config: dict): - return cls(resolution=config["resolution"], affine=config["affine"]) - - def to_dict(self) -> dict: - return {"resolution": self.resolution, "affine": self.affine} - - -@dataclass -class NiftiRef: - file_name: Path - tissue_index: int - - @classmethod - def parse(cls, config: str): - import re - - regex = re.compile(r"(?P.+?)\[(?P\d+)\]$") - m = regex.match(config) - if not m: - raise ValueError("Invalid file_ref", m) - return cls(file_name=Path(m.group("file")), tissue_index=int(m.group("idx"))) - def to_str(self) -> str: - return f"{self.file_name}[{self.tissue_index}]" - - -@dataclass -class NiftiMapping: - file: NiftiRef - func: str - - @classmethod - def parse(cls, config: dict[str, Any]): - return cls(file=NiftiRef.parse(config["file"]), func=config["func"]) - - def to_dict(self) -> dict[str, Any]: - return {"file": self.file.to_str(), "func": self.func} - - -@dataclass -class NiftiTissue: - density: NiftiRef - T1: float | NiftiRef | NiftiMapping - T2: float | NiftiRef | NiftiMapping - T2dash: float | NiftiRef | NiftiMapping - ADC: float | NiftiRef | NiftiMapping - dB0: float | NiftiRef | NiftiMapping - B1_tx: list[float | NiftiRef | NiftiMapping] - B1_rx: list[float | NiftiRef | NiftiMapping] - - @classmethod - def default(cls, density: NiftiRef): - return cls.from_dict({"density": density}) - - @classmethod - def from_dict(cls, config: dict[str, Any]): - def parse_prop(prop): - if isinstance(prop, (float, int)): - return float(prop) - elif isinstance(prop, str): - return NiftiRef.parse(prop) - else: - return NiftiMapping.parse(prop) - - return cls( - density=NiftiRef.parse(config["density"]), - T1=parse_prop(config.get("T1", float("inf"))), - T2=parse_prop(config.get("T2", float("inf"))), - T2dash=parse_prop(config.get("T2'", float("inf"))), - ADC=parse_prop(config.get("ADC", 0.0)), - dB0=parse_prop(config.get("dB0", 0.0)), - B1_tx=[parse_prop(ch) for ch in config.get("B1+", [1.0])], - B1_rx=[parse_prop(ch) for ch in config.get("B1-", [1.0])], - ) - - def to_dict(self) -> dict: - def serialize_prop(prop): - if isinstance(prop, (float, int)): - return prop - elif isinstance(prop, NiftiRef): - return prop.to_str() - elif isinstance(prop, NiftiMapping): - return prop.to_dict() - else: - raise ValueError("Unsupported property type", type(prop)) - - return { - "density": self.density.to_str(), - "T1": serialize_prop(self.T1), - "T2": serialize_prop(self.T2), - "T2'": serialize_prop(self.T2dash), - "ADC": serialize_prop(self.ADC), - "dB0": serialize_prop(self.dB0), - "B1+": [serialize_prop(ch) for ch in self.B1_tx], - "B1-": [serialize_prop(ch) for ch in self.B1_rx], - } - - -@dataclass -class NiftiPhantom: - schema_uri = "https://github.com/mrx-org/bifti-phantoms/bifti-phantom-v1" - units: PhantomUnits - system: PhantomSystem - tissues: dict[str, NiftiTissue] - reslice_to: ResliceConfig | None = field(default=None) - - @classmethod - def default(cls, gyro=42.5764, B0=3.0): - return cls(PhantomUnits.default(), PhantomSystem(gyro, B0), {}) - - @classmethod - def load(cls, path: Path | str): - import json - - with open(path, "r") as f: - config = json.load(f) - return cls.from_dict(config) - - def save(self, path: Path | str): - import json - import re - import os - path = Path(path) - - os.makedirs(path.parent, exist_ok=True) - text = json.dumps(self.to_dict(), indent=2) - # Compact arrays of numbers onto a single line - text = re.sub( - r'\[(\n\s+[-\d.eE+]+,?)+\n\s*\]', - lambda m: '[' + ', '.join(re.findall(r'[-\d.eE+]+', m.group(0))) + ']', - text - ) - with open(path, "w") as f: - f.write(text) - - @classmethod - def from_dict(cls, config: dict): - units = PhantomUnits.from_dict(config["units"]) - system = PhantomSystem.from_dict(config["system"]) - tissues = { - name: NiftiTissue.from_dict(tissue) - for name, tissue in config["tissues"].items() - } - reslice_to = (ResliceConfig.from_dict(config["reslice_to"]) - if "reslice_to" in config else None) - - return cls(units, system, tissues, reslice_to) - - def to_dict(self) -> dict: - d = { - "$schema": self.schema_uri, - "units": self.units.to_dict(), - "system": self.system.to_dict(), - "tissues": { - name: tissue.to_dict() for name, tissue in self.tissues.items() - }, - } - if self.reslice_to is not None: - d["reslice_to"] = self.reslice_to.to_dict() - return d +# The BIfTI phantom config is implemented by the official `bifti` package. +# This module only re-exports it under the names MRzero used before, so that +# existing code like `mr0.NiftiPhantom.load(...)` keeps working. +from bifti import ( + BiftiPhantom, + BiftiTissue, + NiftiMapping, + NiftiRef, + Patient, + PhantomSystem, + PhantomUnits, + ResliceTo, +) + +NiftiPhantom = BiftiPhantom +NiftiTissue = BiftiTissue +ResliceConfig = ResliceTo diff --git a/python/MRzeroCore/phantom/tissue_dict.py b/python/MRzeroCore/phantom/tissue_dict.py index b5e0298..428b4e4 100644 --- a/python/MRzeroCore/phantom/tissue_dict.py +++ b/python/MRzeroCore/phantom/tissue_dict.py @@ -1,12 +1,11 @@ from .voxel_grid_phantom import VoxelGridPhantom from .sim_data import SimData -from .nifti_phantom import NiftiPhantom, NiftiTissue, NiftiRef, NiftiMapping, ResliceConfig +from .nifti_phantom import NiftiPhantom, NiftiTissue, NiftiRef, ResliceConfig, Patient from pathlib import Path import torch import numpy as np import matplotlib.pyplot as plt from typing import Literal, Self -from functools import lru_cache class TissueDict(dict[str, VoxelGridPhantom]): @@ -41,103 +40,50 @@ def load(cls, path: Path | str, config: NiftiPhantom | None = None): # units are supported (conversion factor 1); this might change in the future return TissueDict({ - name: load_tissue(tissue, base_dir, reslice=config.reslice_to) + name: load_tissue(tissue, base_dir, config.reslice_to, config.patient) for name, tissue in config.tissues.items() }) def save(self, path_to_json: str | Path, gyro=42.5764, B0=3.0): - from pathlib import Path - import os - import nibabel as nib + """Save as BIfTI phantom: the `.json` and NIfTIs named after it. - path_to_json = Path(path_to_json) - base_name = path_to_json.stem - base_dir = path_to_json.parent - os.makedirs(base_dir, exist_ok=True) - - density = [] - T1 = [] - T2 = [] - T2dash = [] - ADC = [] - dB0 = [] - B1_tx = [] - B1_rx = [] - - def save_tissue(tissue: VoxelGridPhantom): - config = {} - - def save_map(name, map, nifti): - ext = f"_{name}" if name != "density" else "" - # Multi-channel data when setting the same property multiple times - def set(value): - if name in config: - config[name].append(value) - else: - config[name] = value - - if map.std() < 1e-5: - set(float(map.mean())) - else: - # Check if map is shared with other tissues - for idx, nifti_map in enumerate(nifti): - if torch.equal(map, nifti_map): - set(f"{base_name}{ext}.nii.gz[{idx}]") - return - # Not shared, write new map - set(f"{base_name}{ext}.nii.gz[{len(nifti)}]") - nifti.append(map) - - save_map("density", tissue.PD, density) - save_map("T1", tissue.T1, T1) - save_map("T2", tissue.T2, T2) - save_map("T2'", tissue.T2dash, T2dash) - save_map("ADC", tissue.D, ADC) - save_map("dB0", tissue.B0, dB0) - config["B1+"] = [] - for channel in tissue.B1: - save_map("B1+", channel, B1_tx) - config["B1-"] = [] - for channel in tissue.coil_sens: - save_map("B1-", channel, B1_rx) - - return NiftiTissue.from_dict(config) - - # Generate all tissues (and fill the prop maps) - tissues = {tissue: save_tissue(self[tissue]) for tissue in self.keys()} - - # Write the NIfTIs - affine = np.asarray(next(iter(self.values())).affine) - if affine.shape == (3, 4): - affine = np.vstack([affine, [0, 0, 0, 1]], dtype=np.float32) - - def save_nifti(prop, name): - if len(prop) > 0: - ext = f"_{name}" if name != "density" else "" - file_name = base_dir / f"{base_name}{ext}.nii.gz" - data = np.stack(prop, -1) - - print(f"Storing '{file_name}' - {data.shape}") - nib.save(nib.nifti1.Nifti1Image(data, affine), file_name) - - save_nifti(density, "density") - save_nifti(T1, "T1") - save_nifti(T2, "T2") - save_nifti(T2dash, "T2'") - save_nifti(ADC, "ADC") - save_nifti(dB0, "dB0") - save_nifti(B1_tx, "B1+") - save_nifti(B1_rx, "B1-") - - reslice_to = { - "resolution": list(next(iter(self.values())).PD.shape), - "affine": affine[:3, :4].tolist() - } + Maps that are (nearly) constant are stored as values in the `.json`. + """ + import bifti + path_to_json = Path(path_to_json) config = NiftiPhantom.default(gyro, B0) - config.tissues = tissues - config.reslice_to = ResliceConfig.from_dict(reslice_to) - config.save(path_to_json) + tissues = {} + + def np_map(t: torch.Tensor): + return t.detach().cpu().numpy() + + for name, tissue in self.items(): + # Placeholder references, replaced by the files written by bifti + ref = NiftiRef(Path(), 0) + config.tissues[name] = NiftiTissue( + density=ref, T1=ref, T2=ref, T2dash=ref, ADC=ref, dB0=ref, + B1_tx=[ref] * tissue.B1.shape[0], + B1_rx=[ref] * tissue.coil_sens.shape[0], + ) + tissues[name] = bifti.NumpyTissue( + density=np_map(tissue.PD), + T1=np_map(tissue.T1), + T2=np_map(tissue.T2), + T2dash=np_map(tissue.T2dash), + ADC=np_map(tissue.D), + dB0=np_map(tissue.B0), + B1_tx=np_map(tissue.B1), + B1_rx=np_map(tissue.coil_sens), + resliced=ResliceConfig( + affine=tissue.affine[:3, :].tolist(), + resolution=list(tissue.PD.shape), + ), + ) + + bifti.NumpyPhantom(config, tissues).save( + path_to_json.parent, path_to_json.name, constant_threshold=1e-5 + ) def interpolate(self, x: int, y: int, z: int): return TissueDict({ @@ -298,174 +244,25 @@ def plot(self, tissue="all", plot_masks=False, plot_slice="center", time_unit='s # ============================ def load_tissue(config: NiftiTissue, base_dir: Path, - reslice: ResliceConfig | None = None) -> VoxelGridPhantom: - density, nifti_affine = load_file_ref(base_dir, config.density) + reslice: ResliceConfig | None = None, + patient: Patient | None = None) -> VoxelGridPhantom: + """Load a tissue with the `bifti` reference loader into a VoxelGridPhantom. - def lp(cfg): - return load_property(cfg, base_dir, density, nifti_affine) - - T1 = lp(config.T1) - T2 = lp(config.T2) - T2dash = lp(config.T2dash) - ADC = lp(config.ADC) - B0 = lp(config.dB0) - B1 = [lp(cfg) for cfg in config.B1_tx] - coil = [lp(cfg) for cfg in config.B1_rx] + The voxel-to-scanner affine includes the phantom's patient position (if any). + """ + import bifti.loader - if reslice is None: - target_shape = density.shape - aff_mm = nifti_affine[:3, :] - else: - target_shape = tuple(reslice.resolution) - aff_mm = np.array(reslice.affine, dtype=float) - - def rs(arr): - return _resample_nifti(arr, nifti_affine, target_shape, aff_mm) - density = rs(density) - T1, T2, T2dash, ADC, B0 = rs(T1), rs(T2), rs(T2dash), rs(ADC), rs(B0) - B1 = [rs(b) for b in B1] - coil = [rs(c) for c in coil] + tissue = bifti.loader.load_tissue(config, base_dir, reslice) + affine = bifti.to_scanner_affine(tissue.affine, patient)[:3, :] return VoxelGridPhantom( - PD=torch.as_tensor(density), - T1=torch.as_tensor(T1), - T2=torch.as_tensor(T2), - T2dash=torch.as_tensor(T2dash), - D=torch.as_tensor(ADC), - B0=torch.as_tensor(B0), - B1=torch.stack([torch.as_tensor(b) for b in B1], 0), - coil_sens=torch.stack([torch.as_tensor(c) for c in coil], 0), - affine=torch.as_tensor(aff_mm), + PD=torch.as_tensor(tissue.density), + T1=torch.as_tensor(tissue.T1), + T2=torch.as_tensor(tissue.T2), + T2dash=torch.as_tensor(tissue.T2dash), + D=torch.as_tensor(tissue.ADC), + B0=torch.as_tensor(tissue.dB0), + B1=torch.as_tensor(tissue.B1_tx), + coil_sens=torch.as_tensor(tissue.B1_rx), + affine=torch.as_tensor(affine), ) - - -def load_property(config: float | NiftiRef | NiftiMapping, - base_dir: Path, density_mat: np.ndarray, target_affine: np.ndarray - ) -> np.ndarray: - - if isinstance(config, float): - return np.full_like(density_mat, config) - - if isinstance(config, NiftiRef): - data, affine = load_file_ref(base_dir, config) - assert np.all(affine == target_affine) - return data - - if isinstance(config, NiftiMapping): - data, affine = load_mapping(base_dir, config) - assert np.all(affine == target_affine) - return data - - raise TypeError("Config must be a float, file_ref or mapping", type(config)) - -def load_mapping(base_dir: Path, file_mapping: NiftiMapping) -> tuple[np.ndarray, np.ndarray]: - data, affine = load_file_ref(base_dir, file_mapping.file) - - # TODO - SAFETY: Don't use eval but a custom (imported?) expression parser. - print(f"Executing mapping function: '{file_mapping.func}'") - return eval( - file_mapping.func, - {"__builtins__": None}, - { - "x": data, - "x_min": data.min(), - "x_max": data.max(), - "x_mean": data.mean(), - "x_std": data.std() - } - ), affine - - -def load_file_ref(base_dir: Path, file_ref: NiftiRef) -> tuple[np.ndarray, np.ndarray]: - print(f"Loading NIfTI (file={file_ref.file_name}, index={file_ref.tissue_index})") - file = file_ref.file_name - index = file_ref.tissue_index - if not file.is_absolute(): - file = (base_dir / file).resolve() - - data, affine = _load_cached(str(file)) - return data[:, :, :, index], affine - - -# Use a small cache to avoid reloading NIfTIs every time -@lru_cache(maxsize=20) -def _load_cached(file_name): - import nibabel - img = nibabel.loadsave.load(file_name) - return np.asarray(img.dataobj), img.get_sform() - - -def _resample_nifti(data: np.ndarray, nifti_affine: np.ndarray, - target_shape: tuple, - target_affine_mm: np.ndarray) -> np.ndarray: - """Resample a 3D array onto a target grid via trilinear interpolation, - averaging over the slice thickness using the source voxel size as step. - - Parameters - ---------- - data: - Source 3D numpy array (native NIfTI voxel space). - nifti_affine: - 4×4 sform affine of the source NIfTI (mm units). - target_shape: - Output shape ``(nx, ny, nz)``. - target_affine_mm: - 3×4 or 4×4 NIfTI-style affine of the target grid in mm. - Maps target voxel ``[i, j, k]`` to physical coordinates in mm. - """ - from scipy.ndimage import affine_transform - - # A constant field stays constant under any resampling. Short-circuit this - # case so that constant defaults (e.g. T2 == inf for "no decay") survive: - # feeding non-finite values to ``affine_transform`` would yield NaNs. - flat = np.asarray(data).reshape(-1) - if flat.size > 0 and np.all(flat == flat[0]): - return np.full(tuple(target_shape), flat[0], dtype=float) - - A = np.array(target_affine_mm, dtype=float) - A_rot = A[:3, :3] - A_trans = A[:3, 3] - - A_nifti_inv = np.linalg.inv(nifti_affine[:3, :3]) - - # Voxel size of source phantom (column norms of affine rotation) - src_voxel_mm = np.linalg.norm(nifti_affine[:3, :3], axis=0) - # Slice thickness and normal from Z-column of target affine - slice_vec_mm = A_rot[:, 2] - slice_thickness = np.linalg.norm(slice_vec_mm) - slice_unit = slice_vec_mm / slice_thickness - # Step size = source voxel projected onto slice normal - step_mm = float(np.dot(src_voxel_mm, np.abs(slice_unit))) - n_samples = max(int(round(slice_thickness / step_mm)), 1) - - offsets_mm = np.arange(n_samples, dtype=float) * step_mm - offsets_mm -= offsets_mm.mean() - - M = A_nifti_inv @ A_rot - o0 = A_nifti_inv @ (A_trans - nifti_affine[:3, 3]) - o_step = A_nifti_inv @ slice_unit # offset change per mm along slice normal - - kwargs = dict(output_shape=tuple(target_shape), order=1, - mode='constant', cval=0.0, prefilter=False) - - # Convert once outside the loop - if np.iscomplexobj(data): - data_r = data.real.copy() - data_i = data.imag.copy() - else: - data_r = data.astype(float, copy=False) - data_i = None - - accumulator = np.zeros(target_shape, - dtype=np.complex128 if data_i is not None else float) - - for delta_mm in offsets_mm: - o = o0 + delta_mm * o_step - if data_i is not None: - accumulator += (affine_transform(data_r, M, offset=o, **kwargs) + - 1j * affine_transform(data_i, M, offset=o, **kwargs)) - else: - accumulator += affine_transform(data_r, M, offset=o, **kwargs) - - return accumulator / n_samples -