From d3e6603c32d91a338dd1d505eadcbcdde6df4faf Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?C=C3=A9cile?= Date: Thu, 26 Feb 2026 15:10:26 +0100 Subject: [PATCH 01/24] Start adding functions for generating T1 maps and computing concentration maps --- pyproject.toml | 4 +- src/mritk/concentration/__init__.py | 8 + src/mritk/concentration/concentration.py | 56 +++++++ src/mritk/masking/__init__.py | 8 + src/mritk/masking/masks.py | 80 ++++++++++ src/mritk/masking/utils.py | 15 ++ src/mritk/t1_maps/__init__.py | 6 + src/mritk/t1_maps/dicom_to_nifti.py | 138 +++++++++++++++++ src/mritk/t1_maps/t1_maps.py | 184 +++++++++++++++++++++++ src/mritk/t1_maps/t1_to_r1.py | 56 +++++++ src/mritk/t1_maps/utils.py | 165 ++++++++++++++++++++ 11 files changed, 719 insertions(+), 1 deletion(-) create mode 100644 src/mritk/concentration/__init__.py create mode 100644 src/mritk/concentration/concentration.py create mode 100644 src/mritk/masking/__init__.py create mode 100644 src/mritk/masking/masks.py create mode 100644 src/mritk/masking/utils.py create mode 100644 src/mritk/t1_maps/__init__.py create mode 100644 src/mritk/t1_maps/dicom_to_nifti.py create mode 100644 src/mritk/t1_maps/t1_maps.py create mode 100644 src/mritk/t1_maps/t1_to_r1.py create mode 100644 src/mritk/t1_maps/utils.py diff --git a/pyproject.toml b/pyproject.toml index 1c57521..f8f354b 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -23,7 +23,9 @@ dependencies = [ "nibabel", "pandas", "scipy", -] + "scikit-image", + "pydicom", + "dcm2niix"] [project.optional-dependencies] show = [ diff --git a/src/mritk/concentration/__init__.py b/src/mritk/concentration/__init__.py new file mode 100644 index 0000000..8db514f --- /dev/null +++ b/src/mritk/concentration/__init__.py @@ -0,0 +1,8 @@ +# Copyright (C) 2026 Jørgen Riseth (jnriseth@gmail.com) +# Copyright (C) 2026 Cécile Daversin-Catty (cecile@simula.no) +# Copyright (C) 2026 Simula Research Laboratory + + +from . import concentration + +__all__ = ["concentration"] \ No newline at end of file diff --git a/src/mritk/concentration/concentration.py b/src/mritk/concentration/concentration.py new file mode 100644 index 0000000..89c1101 --- /dev/null +++ b/src/mritk/concentration/concentration.py @@ -0,0 +1,56 @@ +"""Concentration maps module + +Copyright (C) 2026 Jørgen Riseth (jnriseth@gmail.com) +Copyright (C) 2026 Cécile Daversin-Catty (cecile@simula.no) +Copyright (C) 2026 Simula Research Laboratory +""" + +from pathlib import Path +from typing import Optional + +import numpy as np +from ..data.base import MRIData +from ..data.io import load_mri_data, save_mri_data +from ..data.orientation import assert_same_space + + +def concentration_from_T1(T1: np.ndarray, T1_0: np.ndarray, r1: float) -> np.ndarray: + C = 1 / r1 * (1 / T1 - 1 / T1_0) + return C + + +def concentration_from_R1(R1: np.ndarray, R1_0: np.ndarray, r1: float) -> np.ndarray: + C = 1 / r1 * (R1 - R1_0) + return C + + +def concentration( + input: Path, + reference: Path, + output: Optional[Path] = None, + r1: float = 0.0045, + mask: Optional[Path] = None, +)-> MRIData: + T1_mri = load_mri_data(input, np.single) + T10_mri = load_mri_data(reference, np.single) + assert_same_space(T1_mri, T10_mri) + + if mask is not None: + mask_mri = load_mri_data(mask, bool) + assert_same_space(mask_mri, T10_mri) + mask_data = mask_mri.data * (T10_mri.data > 1e-10) * (T1_mri.data > 1e-10) + T1_mri.data *= mask_data + T10_mri.data *= mask_data + else: + mask_data = (T10_mri.data > 1e-10) * (T1_mri.data > 1e-10) + T1_mri.data[~mask_data] = np.nan + T10_mri.data[~mask_data] = np.nan + + concentrations = np.nan * np.zeros_like(T10_mri.data) + concentrations[mask_data] = concentration_from_T1( + T1=T1_mri.data[mask_data], T1_0=T10_mri.data[mask_data], r1=r1 + ) + mri_data = MRIData(data=concentrations, affine=T10_mri.affine) + if output is not None: + save_mri_data(mri_data, output, np.single) + return mri_data diff --git a/src/mritk/masking/__init__.py b/src/mritk/masking/__init__.py new file mode 100644 index 0000000..08e02c3 --- /dev/null +++ b/src/mritk/masking/__init__.py @@ -0,0 +1,8 @@ +# Copyright (C) 2026 Jørgen Riseth (jnriseth@gmail.com) +# Copyright (C) 2026 Cécile Daversin-Catty (cecile@simula.no) +# Copyright (C) 2026 Simula Research Laboratory + + +from . import masks + +__all__ = ["masks", "utils"] \ No newline at end of file diff --git a/src/mritk/masking/masks.py b/src/mritk/masking/masks.py new file mode 100644 index 0000000..f4689c5 --- /dev/null +++ b/src/mritk/masking/masks.py @@ -0,0 +1,80 @@ +"""Intracranial and CSF masks generation module + +Copyright (C) 2026 Jørgen Riseth (jnriseth@gmail.com) +Copyright (C) 2026 Cécile Daversin-Catty (cecile@simula.no) +Copyright (C) 2026 Simula Research Laboratory +""" + +import numpy as np +import skimage +from typing import Optional +from pathlib import Path + +from ..data.base import MRIData +from ..data.io import load_mri_data, save_mri_data +from ..data.orientation import assert_same_space +from .utils import largest_island + + +def create_csf_mask( + vol: np.ndarray, + connectivity: Optional[int] = 2, + use_li: bool = False, +) -> np.ndarray: + connectivity = connectivity or vol.ndim + if use_li: + thresh = skimage.filters.threshold_li(vol) + binary = vol > thresh + binary = largest_island(binary, connectivity=connectivity) + else: + (hist, bins) = np.histogram( + vol[(vol > 0) * (vol < np.quantile(vol, 0.999))], bins=512 + ) + thresh = skimage.filters.threshold_yen(hist=(hist, bins)) + binary = vol > thresh + binary = largest_island(binary, connectivity=connectivity) + return binary + + +def csf_mask( + input: Path, + connectivity: Optional[int] = 2, + use_li: bool = False, + output: Path = None, +) -> MRIData: + input_vol = load_mri_data(input, dtype=np.single) + mask = create_csf_mask(input_vol.data, connectivity, use_li) + assert np.max(mask) > 0, "Masking failed, no voxels in mask" + mri_data = MRIData(data=mask, affine=input_vol.affine) + if output is not None: + save_mri_data(mri_data, output, dtype=np.uint8) + return mri_data + + + +def create_intracranial_mask( + csf_mask: MRIData, + segmentation: MRIData +) -> np.ndarray: + assert_same_space(csf_mask, segmentation) + combined_mask = csf_mask.data + segmentation.data.astype(bool) + background_mask = largest_island(~combined_mask, connectivity=1) + opened = skimage.morphology.binary_opening( + background_mask, skimage.morphology.ball(3) + ) + return ~opened + #return MRIData(data=~opened, affine=segmentation.affine) + + +def intracranial_mask( + csf_mask: Path, + segmentation: Path, + output: Optional[Path] = None, +) -> MRIData: + input_csf_mask = load_mri_data(csf_mask, dtype=bool) + segmentation_data = load_mri_data(segmentation, dtype=bool) + mask_data = create_intracranial_mask(input_csf_mask, segmentation_data) + mri_data = MRIData(data=mask_data, affine=segmentation_data.affine) + if output is not None: + save_mri_data(mri_data, output, dtype=np.uint8) + return mri_data \ No newline at end of file diff --git a/src/mritk/masking/utils.py b/src/mritk/masking/utils.py new file mode 100644 index 0000000..cab7c8c --- /dev/null +++ b/src/mritk/masking/utils.py @@ -0,0 +1,15 @@ +"""Masking utils + +Copyright (C) 2026 Jørgen Riseth (jnriseth@gmail.com) +Copyright (C) 2026 Cécile Daversin-Catty (cecile@simula.no) +Copyright (C) 2026 Simula Research Laboratory +""" + +import numpy as np +import skimage + +def largest_island(mask: np.ndarray, connectivity: int = 1) -> np.ndarray: + newmask = skimage.measure.label(mask, connectivity=connectivity) + regions = skimage.measure.regionprops(newmask) + regions.sort(key=lambda x: x.num_pixels, reverse=True) + return newmask == regions[0].label \ No newline at end of file diff --git a/src/mritk/t1_maps/__init__.py b/src/mritk/t1_maps/__init__.py new file mode 100644 index 0000000..6bf694c --- /dev/null +++ b/src/mritk/t1_maps/__init__.py @@ -0,0 +1,6 @@ +""" +Copyright (C) 2026 Jørgen Riseth (jnriseth@gmail.com) +Copyright (C) 2026 Cécile Daversin-Catty (cecile@simula.no) +Copyright (C) 2026 Simula Research Laboratory +""" +__all__ = ['dicom_to_nifti', 't1_maps', 't1_to_r1', 'utils'] \ No newline at end of file diff --git a/src/mritk/t1_maps/dicom_to_nifti.py b/src/mritk/t1_maps/dicom_to_nifti.py new file mode 100644 index 0000000..6636815 --- /dev/null +++ b/src/mritk/t1_maps/dicom_to_nifti.py @@ -0,0 +1,138 @@ +"""MRI DICOM to NIfTI conversion Module + +Copyright (C) 2026 Jørgen Riseth (jnriseth@gmail.com) +Copyright (C) 2026 Cécile Daversin-Catty (cecile@simula.no) +Copyright (C) 2026 Simula Research Laboratory +""" + +import shutil +import subprocess +import tempfile +from pathlib import Path +from typing import Optional +import nibabel +import json + +import numpy as np + +from ..data.io import load_mri_data, save_mri_data +from ..dicom.utils import VOLUME_LABELS, read_dicom_trigger_times + + +def extract_mixed_dicom(dcmpath: Path, subvolumes: list[str]): + import pydicom + + dcm = pydicom.dcmread(dcmpath) + frames_total = int(dcm.NumberOfFrames) + frames_per_volume = dcm[0x2001, 0x1018].value # [Number of Slices MR] + num_volumes = frames_total // frames_per_volume + assert num_volumes * frames_per_volume == frames_total, ( + "Subvolume dimensions do not match" + ) + + D = dcm.pixel_array.astype(np.single) + frame_fg_sequence = dcm.PerFrameFunctionalGroupsSequence + + vols_out = [] + for volname in subvolumes: + vol_idx = VOLUME_LABELS.index(volname) + + # Find volume slices representing current subvolume + subvol_idx_start = vol_idx * frames_per_volume + subvol_idx_end = (vol_idx + 1) * frames_per_volume + frame_fg = frame_fg_sequence[subvol_idx_start] + logger.info( + ( + f"Converting volume {vol_idx + 1}/{len(VOLUME_LABELS)}: {volname} between indices" + + f"{subvol_idx_start, subvol_idx_end} / {frames_total}." + ) + ) + mri = extract_single_volume(D[subvol_idx_start:subvol_idx_end], frame_fg) + + nii_oriented = nibabel.nifti1.Nifti1Image(mri.data, mri.affine) + nii_oriented.set_sform(nii_oriented.affine, "scanner") + nii_oriented.set_qform(nii_oriented.affine, "scanner") + + # Include meta-data + description = { + "TR": float( + frame_fg.MRTimingAndRelatedParametersSequence[0].RepetitionTime + ), + "TE": float(frame_fg.MREchoSequence[0].EffectiveEchoTime), + } + if hasattr(frame_fg.MRModifierSequence[0], "InversionTimes"): + description["TI"] = frame_fg.MRModifierSequence[0].InversionTimes[0] + if hasattr(frame_fg.MRTimingAndRelatedParametersSequence[0], "EchoTrainLength"): + description["ETL"] = frame_fg.MRTimingAndRelatedParametersSequence[ + 0 + ].EchoTrainLength + vols_out.append({"nifti": nii_oriented, "descrip": description}) + return vols_out + + +def dicom_to_looklocker( + dicomfile: Path, + outpath: Path +): + outdir, form = outpath.parent, outpath.stem + outdir.mkdir(exist_ok=True, parents=True) + times = read_dicom_trigger_times(dicomfile) + np.savetxt(f"{outdir}/{form}" + "_trigger_times.txt", times) + + with tempfile.TemporaryDirectory(prefix=outpath.stem) as tmpdir: + tmppath = Path(tmpdir) + cmd = f"dcm2niix -f {form} -z y --ignore_trigger_times -o '{tmppath}' '{dicomfile}' > /tmp/dcm2niix.txt" + subprocess.run(cmd, shell=True, check=True) + shutil.copy( + tmppath / f"{form}.json", + outpath.with_suffix(".json"), + ) + mri = load_mri_data( + tmppath / f"{form}.nii.gz", dtype=np.double + ) + save_mri_data( + mri, outpath.with_suffix(".nii.gz"), dtype=np.single, intent_code=2001 + ) + + +def dicom_to_mixed( + dcmpath: Path, + outpath: Path, + subvolumes: Optional[list[str]] = None, +): + subvolumes = subvolumes or VOLUME_LABELS + assert all([volname in VOLUME_LABELS for volname in subvolumes]), ( + f"Invalid subvolume name in {subvolumes}, not in {VOLUME_LABELS}" + ) + outdir, form = outpath.parent, outpath.stem + outdir.mkdir(exist_ok=True, parents=True) + + vols = extract_mixed_dicom(dcmpath, subvolumes) + meta = {} + for vol, volname in zip(vols, subvolumes): + output = outpath.with_name(outpath.stem + "_" + volname + ".nii.gz") + + nii = vol["nifti"] + descrip = vol["descrip"] + nibabel.nifti1.save(nii, output) + try: + if volname == "SE-modulus": + meta["TR_SE"] = descrip["TR"] + meta["TE"] = descrip["TE"] + meta["ETL"] = descrip["ETL"] + elif volname == "IR-corrected-real": + meta["TR_IR"] = descrip["TR"] + meta["TI"] = descrip["TI"] + except KeyError as e: + print(volname, descrip) + raise e + + with open(outpath.parent / f"{form}_meta.json", "w") as f: + json.dump(meta, f) + + try: + cmd = f"dcm2niix -w 0 --terse -b o -f '{form}' -o '{outdir}' '{dcmpath}' >> /tmp/dcm2niix.txt " + subprocess.run(cmd, shell=True).check_returncode() + except (ValueError, subprocess.CalledProcessError) as e: + print(str(e)) + pass \ No newline at end of file diff --git a/src/mritk/t1_maps/t1_maps.py b/src/mritk/t1_maps/t1_maps.py new file mode 100644 index 0000000..5ffc148 --- /dev/null +++ b/src/mritk/t1_maps/t1_maps.py @@ -0,0 +1,184 @@ +"""T1 Maps generation module + +Copyright (C) 2026 Jørgen Riseth (jnriseth@gmail.com) +Copyright (C) 2026 Cécile Daversin-Catty (cecile@simula.no) +Copyright (C) 2026 Simula Research Laboratory +""" + +import numpy as np +import tqdm +from functools import partial +import skimage +from typing import Optional +from pathlib import Path +import nibabel +import json +import scipy + +from ..data.base import MRIData +from ..data.io import load_mri_data +from ..masking.masks import create_csf_mask +from .utils import ( + mri_facemask, + fit_voxel, + nan_filter_gaussian, + T1_lookup_table, +) + + +def looklocker_t1map( + t_data: np.ndarray, + D: np.ndarray, + affine: np.ndarray +) -> MRIData: + + T1_ROOF = 10000 + assert len(D.shape) >= 4, ( + f"data should be 4-dimensional, got data with shape {D.shape}" + ) + mask = mri_facemask(D[..., 0]) + valid_voxels = (np.nanmax(D, axis=-1) > 0) * mask + + D_normalized = np.nan * np.zeros_like(D) + D_normalized[valid_voxels] = ( + D[valid_voxels] / np.nanmax(D, axis=-1)[valid_voxels, np.newaxis] + ) + voxel_mask = np.array(np.where(valid_voxels)).T + Dmasked = np.array([D_normalized[i, j, k] for (i, j, k) in voxel_mask]) + + with tqdm.tqdm(total=len(Dmasked)) as pbar: + voxel_fitter = partial(fit_voxel, t_data, pbar) + vfunc = np.vectorize(voxel_fitter, signature="(n) -> (3)") + fitted_coefficients = vfunc(Dmasked) + + x1, x2, x3 = ( + fitted_coefficients[:, 0], + fitted_coefficients[:, 1], + fitted_coefficients[:, 2], + ) + + I, J, K = voxel_mask.T + T1map = np.nan * np.zeros_like(D[..., 0]) + T1map[I, J, K] = (x2 / x3) ** 2 * 1000.0 # convert to ms + T1map = np.minimum(T1map, T1_ROOF) + return MRIData(T1map.astype(np.single), affine) + + +def looklocker_t1map_postprocessing( + T1map_mri: MRIData, + T1_lo: float, + T1_hi: float, + radius: int = 10, + erode_dilate_factor: float = 1.3, + mask: Optional[np.ndarray] = None, +) -> MRIData: + T1map = T1map_mri.data.copy() + if mask is None: + # Create mask for largest island. + mask = skimage.measure.label(np.isfinite(T1map)) + regions = skimage.measure.regionprops(mask) + regions.sort(key=lambda x: x.num_pixels, reverse=True) + mask = mask == regions[0].label + skimage.morphology.remove_small_holes( + mask, 10 ** (mask.ndim), connectivity=2, out=mask + ) + skimage.morphology.binary_dilation( + mask, skimage.morphology.ball(radius), out=mask + ) + skimage.morphology.binary_erosion( + mask, skimage.morphology.ball(erode_dilate_factor * radius), out=mask + ) + + # Remove non-zero artifacts outside of the mask. + surface_vox = np.isfinite(T1map) * (~mask) + print(f"Removing {surface_vox.sum()} voxels outside of the head mask") + T1map[~mask] = np.nan + + # Remove outliers within the mask. + outliers = np.logical_or(T1map < T1_lo, T1_hi < T1map) + print("Removing", outliers.sum(), f"voxels outside the range ({T1_lo}, {T1_hi}).") + T1map[outliers] = np.nan + if np.isfinite(T1map).sum() / T1map.size < 0.01: + raise RuntimeError( + "After outlier removal, less than 1% of the image is left. Check image units." + ) + + # Fill internallly missing values + fill_mask = np.isnan(T1map) * mask + while fill_mask.sum() > 0: + print(f"Filling in {fill_mask.sum()} voxels within the mask.") + T1map[fill_mask] = nan_filter_gaussian(T1map, 1.0)[fill_mask] + fill_mask = np.isnan(T1map) * mask + return MRIData(T1map, T1map_mri.affine) + + +def mixed_t1map( + SE_nii_path: Path, + IR_nii_path: Path, + meta_path: Path, + T1_low: float, + T1_hi: float, +) -> nibabel.nifti1.Nifti1Image: + SE = load_mri_data(SE_nii_path, dtype=np.single) + IR = load_mri_data(IR_nii_path, dtype=np.single) + with open(meta_path, "r") as f: + meta = json.load(f) + + nonzero_mask = SE.data != 0 + F_data = np.nan * np.zeros_like(IR.data) + F_data[nonzero_mask] = IR.data[nonzero_mask] / SE.data[nonzero_mask] + + TR_se, TI, TE, ETL = meta["TR_SE"], meta["TI"], meta["TE"], meta["ETL"] + F, T1_grid = T1_lookup_table(TR_se, TI, TE, ETL, T1_low, T1_hi) + interpolator = scipy.interpolate.interp1d( + F, T1_grid, kind="nearest", bounds_error=False, fill_value=np.nan + ) + T1_volume = interpolator(F_data).astype(np.single) + nii = nibabel.nifti1.Nifti1Image(T1_volume, IR.affine) + nii.set_sform(nii.affine, "scanner") + nii.set_qform(nii.affine, "scanner") + return nii + + +def mixed_t1map_postprocessing( + se: Path, + t1: Path, + output: Path +): + T1map_nii = nibabel.nifti1.load(t1) + + SE_mri = load_mri_data(se, np.single) + mask = create_csf_mask(SE_mri.data, use_li=True) + mask = skimage.morphology.binary_erosion(mask) + + masked_T1map = T1map_nii.get_fdata(dtype=np.single) + masked_T1map[~mask] = np.nan + masked_T1map_nii = nibabel.nifti1.Nifti1Image( + masked_T1map, T1map_nii.affine, T1map_nii.header + ) + nibabel.nifti1.save(masked_T1map_nii, output) + + +def hybrid_t1_map( + ll_path: Path, + mixed_path: Path, + csf_mask_path: Path, + threshold: float, + erode: int = 0, +) -> nibabel.nifti1.Nifti1Image: + mixed_mri = nibabel.nifti1.load(mixed_path) + mixed = mixed_mri.get_fdata() + + ll_mri = nibabel.nifti1.load(ll_path) + ll = ll_mri.get_fdata() + csf_mask_mri = nibabel.nifti1.load(csf_mask_path) + csf_mask = csf_mask_mri.get_fdata().astype(bool) + if erode > 0: + csf_mask = skimage.morphology.binary_erosion( + csf_mask, skimage.morphology.ball(erode) + ) + + hybrid = ll + newmask = csf_mask * (ll > threshold) * (mixed > threshold) + hybrid[newmask] = mixed[newmask] + return nibabel.nifti1.Nifti1Image(hybrid, affine=ll_mri.affine, header=ll_mri.header) \ No newline at end of file diff --git a/src/mritk/t1_maps/t1_to_r1.py b/src/mritk/t1_maps/t1_to_r1.py new file mode 100644 index 0000000..a385207 --- /dev/null +++ b/src/mritk/t1_maps/t1_to_r1.py @@ -0,0 +1,56 @@ +"""T1 to R1 module + +Copyright (C) 2026 Jørgen Riseth (jnriseth@gmail.com) +Copyright (C) 2026 Cécile Daversin-Catty (cecile@simula.no) +Copyright (C) 2026 Simula Research Laboratory +""" + +import numpy as np +from pathlib import Path +from typing import Union + +from ..data.base import MRIData +from ..data.io import load_mri_data, save_mri_data + + +def convert_T1_to_R1( + T1map_mri: MRIData, + scale: float = 1000, + t1_low: float = 1, + t1_high: float = float("inf"), +) -> MRIData: + valid_t1 = (t1_low <= T1map_mri.data) * (T1map_mri.data <= t1_high) + R1map = np.nan * np.zeros_like(T1map_mri.data) + R1map[valid_t1] = scale / np.minimum( + t1_high, np.maximum(t1_low, T1map_mri.data[valid_t1]) + ) + + R1map_mri = MRIData(data=R1map, affine=T1map_mri.affine) + return R1map_mri + + +def T1_to_R1( + input_mri: Union[Path, MRIData], + output: Path = None, + scale: float = 1000, + t1_low: float = 1, + t1_high: float = float("inf"), +) -> MRIData: + if isinstance(input_mri, Path): + T1map_mri = load_mri_data(input_mri, dtype=np.single) + elif isinstance(input_mri, MRIData): + T1map_mri = input_mri + else: + raise ValueError(f"Input should be a Path or MRIData, got {type(input_mri)}") + + valid_t1 = (t1_low <= T1map_mri.data) * (T1map_mri.data <= t1_high) + R1map = np.nan * np.zeros_like(T1map_mri.data) + R1map[valid_t1] = scale / np.minimum( + t1_high, np.maximum(t1_low, T1map_mri.data[valid_t1]) + ) + + R1map_mri = MRIData(data=R1map, affine=T1map_mri.affine) + + if output is not None: + save_mri_data(R1map_mri, output, dtype=np.single) + return R1map_mri \ No newline at end of file diff --git a/src/mritk/t1_maps/utils.py b/src/mritk/t1_maps/utils.py new file mode 100644 index 0000000..fd2c0bc --- /dev/null +++ b/src/mritk/t1_maps/utils.py @@ -0,0 +1,165 @@ +"""MRI DICOM to NIfTI conversion - utils + +Copyright (C) 2026 Jørgen Riseth (jnriseth@gmail.com) +Copyright (C) 2026 Cécile Daversin-Catty (cecile@simula.no) +Copyright (C) 2026 Simula Research Laboratory +""" + +from pathlib import Path +import numpy as np +import nibabel +import logging +import scipy as sp +import skimage +import warnings +from scipy.optimize import OptimizeWarning + +from ..data.orientation import data_reorientation, change_of_coordinates_map +from ..data.base import MRIData + +logger = logging.getLogger(__name__) + +VOLUME_LABELS = [ + "IR-modulus", + "IR-real", + "IR-corrected-real", + "SE-modulus", + "SE-real", + "T1map-scanner", +] + + +def read_dicom_trigger_times(dicomfile): + import pydicom + + dcm = pydicom.dcmread(dicomfile) + all_frame_times = [ + f.CardiacSynchronizationSequence[0].NominalCardiacTriggerDelayTime + for f in dcm.PerFrameFunctionalGroupsSequence + ] + return np.unique(all_frame_times) + + +def dicom_standard_affine( + frame_fg, +) -> np.ndarray: + # Get the original data shape + df = float(frame_fg.PixelMeasuresSequence[0].SpacingBetweenSlices) + dr, dc = (float(x) for x in frame_fg.PixelMeasuresSequence[0].PixelSpacing) + plane_orientation = frame_fg.PlaneOrientationSequence[0] + orientation = np.array(plane_orientation.ImageOrientationPatient) + + # Find orientation of data array relative to LPS-coordinate system. + row_cosine = orientation[:3] + col_cosine = orientation[3:] + frame_cosine = np.cross(row_cosine, col_cosine) + + # Create DICOM-definition affine map to LPS. + T_1 = np.array(frame_fg.PlanePositionSequence[0].ImagePositionPatient) + + # Create DICOM-definition affine map to LPS. + M_dcm = np.zeros((4, 4)) + M_dcm[:3, 0] = row_cosine * dc + M_dcm[:3, 1] = col_cosine * dr + M_dcm[:3, 2] = frame_cosine * df + M_dcm[:3, 3] = T_1 + M_dcm[3, 3] = 1.0 + + # Reorder from "natural index order" to DICOM affine map definition order. + N_order = np.eye(4)[[2, 1, 0, 3]] + return M_dcm @ N_order + + +def extract_single_volume( + D: np.ndarray, + frame_fg, +) -> MRIData: + # Find scaling values (should potentially be inside scaling loop) + pixel_value_transform = frame_fg.PixelValueTransformationSequence[0] + slope = float(pixel_value_transform.RescaleSlope) + intercept = float(pixel_value_transform.RescaleIntercept) + private = frame_fg[0x2005, 0x140F][0] + scale_slope = private[0x2005, 0x100E].value + + # Loop over and scale values. + volume = np.zeros_like(D, dtype=np.single) + for idx in range(D.shape[0]): + volume[idx] = (intercept + slope * D[idx]) / (scale_slope * slope) + + A_dcm = dicom_standard_affine(frame_fg) + C = change_of_coordinates_map("LPS", "RAS") + mri = data_reorientation(MRIData(volume, C @ A_dcm)) + + return mri + + +def mri_facemask(vol: np.ndarray, smoothing_level=5): + thresh = skimage.filters.threshold_triangle(vol) + binary = vol > thresh + binary = sp.ndimage.binary_fill_holes(binary) + binary = skimage.filters.gaussian(binary, sigma=smoothing_level) + binary = binary > skimage.filters.threshold_isodata(binary) + return binary + + +def curve_fit_wrapper(f, t, y, p0): + """Raises error instead of catching numpy warnings, such that + these cases may be treated.""" + with warnings.catch_warnings(): + warnings.simplefilter("error", OptimizeWarning) + popt, _ = sp.optimize.curve_fit(f, xdata=t, ydata=y, p0=p0, maxfev=1000) + return popt + + +def fit_voxel(time_s: np.ndarray, pbar, m: np.ndarray) -> np.ndarray: + if pbar is not None: + pbar.update(1) + x1 = 1.0 + x2 = np.sqrt(1.25) + T1 = time_s[np.argmin(m)] / np.log(1 + x2**2) + x3 = np.sqrt(1 / T1) + p0 = np.array((x1, x2, x3)) + if not np.all(np.isfinite(m)): + return np.nan * np.zeros_like(p0) + try: + popt = curve_fit_wrapper(f, time_s, m, p0) + except (OptimizeWarning, FloatingPointError): + return np.nan * np.zeros_like(p0) + except RuntimeError as e: + if "maxfev" in str(e): + return np.nan * np.zeros_like(p0) + raise e + return popt + + +def nan_filter_gaussian( + U: np.ndarray, sigma: float, truncate: float = 4.0 +) -> np.ndarray: + V = U.copy() + V[np.isnan(U)] = 0 + VV = sp.ndimage.gaussian_filter(V, sigma=sigma, truncate=truncate) + + W = np.ones_like(U) + W[np.isnan(U)] = 0 + WW = sp.ndimage.gaussian_filter(W, sigma=sigma, truncate=truncate) + mask = ~((WW == 0) * (VV == 0)) + out = np.nan * np.zeros_like(U) + out[mask] = VV[mask] / WW[mask] + return out + + +def estimate_se_free_relaxation_time(TRse, TE, ETL): + """Compute free relaxation time following spin echo image from effective echo + time TE and echo train length ETL, corrected for 20 dummy echoes.""" + return TRse - TE * (1 + 0.5 * (ETL - 1) / (0.5 * (ETL + 1) + 20)) + + +def T1_lookup_table( + TRse: float, TI: float, TE: float, ETL: int, T1_low: float, T1_hi: float +) -> tuple[np.ndarray, np.ndarray]: + TRfree = estimate_se_free_relaxation_time(TRse, TE, ETL) + T1_grid = np.arange(int(T1_low), int(T1_hi + 1)) + Sse = 1 - np.exp(-TRfree / T1_grid) + Sir = 1 - (1 + Sse) * np.exp(-TI / T1_grid) + fractionCurve = Sir / Sse + return fractionCurve, T1_grid \ No newline at end of file From 1fedb19a9cd60c8b338e55eaf70587bd900047d0 Mon Sep 17 00:00:00 2001 From: "pre-commit-ci[bot]" <66853113+pre-commit-ci[bot]@users.noreply.github.com> Date: Tue, 3 Mar 2026 18:55:56 +0000 Subject: [PATCH 02/24] [pre-commit.ci] auto fixes from pre-commit.com hooks for more information, see https://pre-commit.ci --- src/mritk/concentration/__init__.py | 2 +- src/mritk/concentration/concentration.py | 6 +-- src/mritk/masking/__init__.py | 2 +- src/mritk/masking/masks.py | 18 +++------ src/mritk/masking/utils.py | 3 +- src/mritk/t1_maps/__init__.py | 3 +- src/mritk/t1_maps/dicom_to_nifti.py | 27 ++++--------- src/mritk/t1_maps/t1_maps.py | 51 ++++++------------------ src/mritk/t1_maps/t1_to_r1.py | 10 ++--- src/mritk/t1_maps/utils.py | 17 +++----- 10 files changed, 40 insertions(+), 99 deletions(-) diff --git a/src/mritk/concentration/__init__.py b/src/mritk/concentration/__init__.py index 8db514f..9e93ffa 100644 --- a/src/mritk/concentration/__init__.py +++ b/src/mritk/concentration/__init__.py @@ -5,4 +5,4 @@ from . import concentration -__all__ = ["concentration"] \ No newline at end of file +__all__ = ["concentration"] diff --git a/src/mritk/concentration/concentration.py b/src/mritk/concentration/concentration.py index 89c1101..08c0a67 100644 --- a/src/mritk/concentration/concentration.py +++ b/src/mritk/concentration/concentration.py @@ -30,7 +30,7 @@ def concentration( output: Optional[Path] = None, r1: float = 0.0045, mask: Optional[Path] = None, -)-> MRIData: +) -> MRIData: T1_mri = load_mri_data(input, np.single) T10_mri = load_mri_data(reference, np.single) assert_same_space(T1_mri, T10_mri) @@ -47,9 +47,7 @@ def concentration( T10_mri.data[~mask_data] = np.nan concentrations = np.nan * np.zeros_like(T10_mri.data) - concentrations[mask_data] = concentration_from_T1( - T1=T1_mri.data[mask_data], T1_0=T10_mri.data[mask_data], r1=r1 - ) + concentrations[mask_data] = concentration_from_T1(T1=T1_mri.data[mask_data], T1_0=T10_mri.data[mask_data], r1=r1) mri_data = MRIData(data=concentrations, affine=T10_mri.affine) if output is not None: save_mri_data(mri_data, output, np.single) diff --git a/src/mritk/masking/__init__.py b/src/mritk/masking/__init__.py index 08e02c3..4bd6155 100644 --- a/src/mritk/masking/__init__.py +++ b/src/mritk/masking/__init__.py @@ -5,4 +5,4 @@ from . import masks -__all__ = ["masks", "utils"] \ No newline at end of file +__all__ = ["masks", "utils"] diff --git a/src/mritk/masking/masks.py b/src/mritk/masking/masks.py index f4689c5..15df91f 100644 --- a/src/mritk/masking/masks.py +++ b/src/mritk/masking/masks.py @@ -27,9 +27,7 @@ def create_csf_mask( binary = vol > thresh binary = largest_island(binary, connectivity=connectivity) else: - (hist, bins) = np.histogram( - vol[(vol > 0) * (vol < np.quantile(vol, 0.999))], bins=512 - ) + (hist, bins) = np.histogram(vol[(vol > 0) * (vol < np.quantile(vol, 0.999))], bins=512) thresh = skimage.filters.threshold_yen(hist=(hist, bins)) binary = vol > thresh binary = largest_island(binary, connectivity=connectivity) @@ -51,19 +49,13 @@ def csf_mask( return mri_data - -def create_intracranial_mask( - csf_mask: MRIData, - segmentation: MRIData -) -> np.ndarray: +def create_intracranial_mask(csf_mask: MRIData, segmentation: MRIData) -> np.ndarray: assert_same_space(csf_mask, segmentation) combined_mask = csf_mask.data + segmentation.data.astype(bool) background_mask = largest_island(~combined_mask, connectivity=1) - opened = skimage.morphology.binary_opening( - background_mask, skimage.morphology.ball(3) - ) + opened = skimage.morphology.binary_opening(background_mask, skimage.morphology.ball(3)) return ~opened - #return MRIData(data=~opened, affine=segmentation.affine) + # return MRIData(data=~opened, affine=segmentation.affine) def intracranial_mask( @@ -77,4 +69,4 @@ def intracranial_mask( mri_data = MRIData(data=mask_data, affine=segmentation_data.affine) if output is not None: save_mri_data(mri_data, output, dtype=np.uint8) - return mri_data \ No newline at end of file + return mri_data diff --git a/src/mritk/masking/utils.py b/src/mritk/masking/utils.py index cab7c8c..2858a5c 100644 --- a/src/mritk/masking/utils.py +++ b/src/mritk/masking/utils.py @@ -8,8 +8,9 @@ import numpy as np import skimage + def largest_island(mask: np.ndarray, connectivity: int = 1) -> np.ndarray: newmask = skimage.measure.label(mask, connectivity=connectivity) regions = skimage.measure.regionprops(newmask) regions.sort(key=lambda x: x.num_pixels, reverse=True) - return newmask == regions[0].label \ No newline at end of file + return newmask == regions[0].label diff --git a/src/mritk/t1_maps/__init__.py b/src/mritk/t1_maps/__init__.py index 6bf694c..1b67101 100644 --- a/src/mritk/t1_maps/__init__.py +++ b/src/mritk/t1_maps/__init__.py @@ -3,4 +3,5 @@ Copyright (C) 2026 Cécile Daversin-Catty (cecile@simula.no) Copyright (C) 2026 Simula Research Laboratory """ -__all__ = ['dicom_to_nifti', 't1_maps', 't1_to_r1', 'utils'] \ No newline at end of file + +__all__ = ["dicom_to_nifti", "t1_maps", "t1_to_r1", "utils"] diff --git a/src/mritk/t1_maps/dicom_to_nifti.py b/src/mritk/t1_maps/dicom_to_nifti.py index 6636815..5ced6d4 100644 --- a/src/mritk/t1_maps/dicom_to_nifti.py +++ b/src/mritk/t1_maps/dicom_to_nifti.py @@ -26,9 +26,7 @@ def extract_mixed_dicom(dcmpath: Path, subvolumes: list[str]): frames_total = int(dcm.NumberOfFrames) frames_per_volume = dcm[0x2001, 0x1018].value # [Number of Slices MR] num_volumes = frames_total // frames_per_volume - assert num_volumes * frames_per_volume == frames_total, ( - "Subvolume dimensions do not match" - ) + assert num_volumes * frames_per_volume == frames_total, "Subvolume dimensions do not match" D = dcm.pixel_array.astype(np.single) frame_fg_sequence = dcm.PerFrameFunctionalGroupsSequence @@ -55,25 +53,18 @@ def extract_mixed_dicom(dcmpath: Path, subvolumes: list[str]): # Include meta-data description = { - "TR": float( - frame_fg.MRTimingAndRelatedParametersSequence[0].RepetitionTime - ), + "TR": float(frame_fg.MRTimingAndRelatedParametersSequence[0].RepetitionTime), "TE": float(frame_fg.MREchoSequence[0].EffectiveEchoTime), } if hasattr(frame_fg.MRModifierSequence[0], "InversionTimes"): description["TI"] = frame_fg.MRModifierSequence[0].InversionTimes[0] if hasattr(frame_fg.MRTimingAndRelatedParametersSequence[0], "EchoTrainLength"): - description["ETL"] = frame_fg.MRTimingAndRelatedParametersSequence[ - 0 - ].EchoTrainLength + description["ETL"] = frame_fg.MRTimingAndRelatedParametersSequence[0].EchoTrainLength vols_out.append({"nifti": nii_oriented, "descrip": description}) return vols_out -def dicom_to_looklocker( - dicomfile: Path, - outpath: Path -): +def dicom_to_looklocker(dicomfile: Path, outpath: Path): outdir, form = outpath.parent, outpath.stem outdir.mkdir(exist_ok=True, parents=True) times = read_dicom_trigger_times(dicomfile) @@ -87,12 +78,8 @@ def dicom_to_looklocker( tmppath / f"{form}.json", outpath.with_suffix(".json"), ) - mri = load_mri_data( - tmppath / f"{form}.nii.gz", dtype=np.double - ) - save_mri_data( - mri, outpath.with_suffix(".nii.gz"), dtype=np.single, intent_code=2001 - ) + mri = load_mri_data(tmppath / f"{form}.nii.gz", dtype=np.double) + save_mri_data(mri, outpath.with_suffix(".nii.gz"), dtype=np.single, intent_code=2001) def dicom_to_mixed( @@ -135,4 +122,4 @@ def dicom_to_mixed( subprocess.run(cmd, shell=True).check_returncode() except (ValueError, subprocess.CalledProcessError) as e: print(str(e)) - pass \ No newline at end of file + pass diff --git a/src/mritk/t1_maps/t1_maps.py b/src/mritk/t1_maps/t1_maps.py index 5ffc148..44091bc 100644 --- a/src/mritk/t1_maps/t1_maps.py +++ b/src/mritk/t1_maps/t1_maps.py @@ -26,23 +26,14 @@ ) -def looklocker_t1map( - t_data: np.ndarray, - D: np.ndarray, - affine: np.ndarray -) -> MRIData: - +def looklocker_t1map(t_data: np.ndarray, D: np.ndarray, affine: np.ndarray) -> MRIData: T1_ROOF = 10000 - assert len(D.shape) >= 4, ( - f"data should be 4-dimensional, got data with shape {D.shape}" - ) + assert len(D.shape) >= 4, f"data should be 4-dimensional, got data with shape {D.shape}" mask = mri_facemask(D[..., 0]) valid_voxels = (np.nanmax(D, axis=-1) > 0) * mask D_normalized = np.nan * np.zeros_like(D) - D_normalized[valid_voxels] = ( - D[valid_voxels] / np.nanmax(D, axis=-1)[valid_voxels, np.newaxis] - ) + D_normalized[valid_voxels] = D[valid_voxels] / np.nanmax(D, axis=-1)[valid_voxels, np.newaxis] voxel_mask = np.array(np.where(valid_voxels)).T Dmasked = np.array([D_normalized[i, j, k] for (i, j, k) in voxel_mask]) @@ -79,15 +70,9 @@ def looklocker_t1map_postprocessing( regions = skimage.measure.regionprops(mask) regions.sort(key=lambda x: x.num_pixels, reverse=True) mask = mask == regions[0].label - skimage.morphology.remove_small_holes( - mask, 10 ** (mask.ndim), connectivity=2, out=mask - ) - skimage.morphology.binary_dilation( - mask, skimage.morphology.ball(radius), out=mask - ) - skimage.morphology.binary_erosion( - mask, skimage.morphology.ball(erode_dilate_factor * radius), out=mask - ) + skimage.morphology.remove_small_holes(mask, 10 ** (mask.ndim), connectivity=2, out=mask) + skimage.morphology.binary_dilation(mask, skimage.morphology.ball(radius), out=mask) + skimage.morphology.binary_erosion(mask, skimage.morphology.ball(erode_dilate_factor * radius), out=mask) # Remove non-zero artifacts outside of the mask. surface_vox = np.isfinite(T1map) * (~mask) @@ -99,9 +84,7 @@ def looklocker_t1map_postprocessing( print("Removing", outliers.sum(), f"voxels outside the range ({T1_lo}, {T1_hi}).") T1map[outliers] = np.nan if np.isfinite(T1map).sum() / T1map.size < 0.01: - raise RuntimeError( - "After outlier removal, less than 1% of the image is left. Check image units." - ) + raise RuntimeError("After outlier removal, less than 1% of the image is left. Check image units.") # Fill internallly missing values fill_mask = np.isnan(T1map) * mask @@ -130,9 +113,7 @@ def mixed_t1map( TR_se, TI, TE, ETL = meta["TR_SE"], meta["TI"], meta["TE"], meta["ETL"] F, T1_grid = T1_lookup_table(TR_se, TI, TE, ETL, T1_low, T1_hi) - interpolator = scipy.interpolate.interp1d( - F, T1_grid, kind="nearest", bounds_error=False, fill_value=np.nan - ) + interpolator = scipy.interpolate.interp1d(F, T1_grid, kind="nearest", bounds_error=False, fill_value=np.nan) T1_volume = interpolator(F_data).astype(np.single) nii = nibabel.nifti1.Nifti1Image(T1_volume, IR.affine) nii.set_sform(nii.affine, "scanner") @@ -140,11 +121,7 @@ def mixed_t1map( return nii -def mixed_t1map_postprocessing( - se: Path, - t1: Path, - output: Path -): +def mixed_t1map_postprocessing(se: Path, t1: Path, output: Path): T1map_nii = nibabel.nifti1.load(t1) SE_mri = load_mri_data(se, np.single) @@ -153,9 +130,7 @@ def mixed_t1map_postprocessing( masked_T1map = T1map_nii.get_fdata(dtype=np.single) masked_T1map[~mask] = np.nan - masked_T1map_nii = nibabel.nifti1.Nifti1Image( - masked_T1map, T1map_nii.affine, T1map_nii.header - ) + masked_T1map_nii = nibabel.nifti1.Nifti1Image(masked_T1map, T1map_nii.affine, T1map_nii.header) nibabel.nifti1.save(masked_T1map_nii, output) @@ -174,11 +149,9 @@ def hybrid_t1_map( csf_mask_mri = nibabel.nifti1.load(csf_mask_path) csf_mask = csf_mask_mri.get_fdata().astype(bool) if erode > 0: - csf_mask = skimage.morphology.binary_erosion( - csf_mask, skimage.morphology.ball(erode) - ) + csf_mask = skimage.morphology.binary_erosion(csf_mask, skimage.morphology.ball(erode)) hybrid = ll newmask = csf_mask * (ll > threshold) * (mixed > threshold) hybrid[newmask] = mixed[newmask] - return nibabel.nifti1.Nifti1Image(hybrid, affine=ll_mri.affine, header=ll_mri.header) \ No newline at end of file + return nibabel.nifti1.Nifti1Image(hybrid, affine=ll_mri.affine, header=ll_mri.header) diff --git a/src/mritk/t1_maps/t1_to_r1.py b/src/mritk/t1_maps/t1_to_r1.py index a385207..a783928 100644 --- a/src/mritk/t1_maps/t1_to_r1.py +++ b/src/mritk/t1_maps/t1_to_r1.py @@ -21,9 +21,7 @@ def convert_T1_to_R1( ) -> MRIData: valid_t1 = (t1_low <= T1map_mri.data) * (T1map_mri.data <= t1_high) R1map = np.nan * np.zeros_like(T1map_mri.data) - R1map[valid_t1] = scale / np.minimum( - t1_high, np.maximum(t1_low, T1map_mri.data[valid_t1]) - ) + R1map[valid_t1] = scale / np.minimum(t1_high, np.maximum(t1_low, T1map_mri.data[valid_t1])) R1map_mri = MRIData(data=R1map, affine=T1map_mri.affine) return R1map_mri @@ -45,12 +43,10 @@ def T1_to_R1( valid_t1 = (t1_low <= T1map_mri.data) * (T1map_mri.data <= t1_high) R1map = np.nan * np.zeros_like(T1map_mri.data) - R1map[valid_t1] = scale / np.minimum( - t1_high, np.maximum(t1_low, T1map_mri.data[valid_t1]) - ) + R1map[valid_t1] = scale / np.minimum(t1_high, np.maximum(t1_low, T1map_mri.data[valid_t1])) R1map_mri = MRIData(data=R1map, affine=T1map_mri.affine) if output is not None: save_mri_data(R1map_mri, output, dtype=np.single) - return R1map_mri \ No newline at end of file + return R1map_mri diff --git a/src/mritk/t1_maps/utils.py b/src/mritk/t1_maps/utils.py index fd2c0bc..2ccf14c 100644 --- a/src/mritk/t1_maps/utils.py +++ b/src/mritk/t1_maps/utils.py @@ -5,9 +5,7 @@ Copyright (C) 2026 Simula Research Laboratory """ -from pathlib import Path import numpy as np -import nibabel import logging import scipy as sp import skimage @@ -34,8 +32,7 @@ def read_dicom_trigger_times(dicomfile): dcm = pydicom.dcmread(dicomfile) all_frame_times = [ - f.CardiacSynchronizationSequence[0].NominalCardiacTriggerDelayTime - for f in dcm.PerFrameFunctionalGroupsSequence + f.CardiacSynchronizationSequence[0].NominalCardiacTriggerDelayTime for f in dcm.PerFrameFunctionalGroupsSequence ] return np.unique(all_frame_times) @@ -109,7 +106,7 @@ def curve_fit_wrapper(f, t, y, p0): warnings.simplefilter("error", OptimizeWarning) popt, _ = sp.optimize.curve_fit(f, xdata=t, ydata=y, p0=p0, maxfev=1000) return popt - + def fit_voxel(time_s: np.ndarray, pbar, m: np.ndarray) -> np.ndarray: if pbar is not None: @@ -132,9 +129,7 @@ def fit_voxel(time_s: np.ndarray, pbar, m: np.ndarray) -> np.ndarray: return popt -def nan_filter_gaussian( - U: np.ndarray, sigma: float, truncate: float = 4.0 -) -> np.ndarray: +def nan_filter_gaussian(U: np.ndarray, sigma: float, truncate: float = 4.0) -> np.ndarray: V = U.copy() V[np.isnan(U)] = 0 VV = sp.ndimage.gaussian_filter(V, sigma=sigma, truncate=truncate) @@ -154,12 +149,10 @@ def estimate_se_free_relaxation_time(TRse, TE, ETL): return TRse - TE * (1 + 0.5 * (ETL - 1) / (0.5 * (ETL + 1) + 20)) -def T1_lookup_table( - TRse: float, TI: float, TE: float, ETL: int, T1_low: float, T1_hi: float -) -> tuple[np.ndarray, np.ndarray]: +def T1_lookup_table(TRse: float, TI: float, TE: float, ETL: int, T1_low: float, T1_hi: float) -> tuple[np.ndarray, np.ndarray]: TRfree = estimate_se_free_relaxation_time(TRse, TE, ETL) T1_grid = np.arange(int(T1_low), int(T1_hi + 1)) Sse = 1 - np.exp(-TRfree / T1_grid) Sir = 1 - (1 + Sse) * np.exp(-TI / T1_grid) fractionCurve = Sir / Sse - return fractionCurve, T1_grid \ No newline at end of file + return fractionCurve, T1_grid From 91e1e9a8d49dca4b3a4eb29807ce0126c4cc5a3b Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?C=C3=A9cile?= Date: Thu, 5 Mar 2026 09:28:57 +0100 Subject: [PATCH 03/24] Fix module import + start updating functions in t1 map --- src/mritk/data/orientation.py | 1 + src/mritk/t1_maps/dicom_to_nifti.py | 3 ++- src/mritk/t1_maps/t1_maps.py | 22 ++++++++++++++++++---- 3 files changed, 21 insertions(+), 5 deletions(-) diff --git a/src/mritk/data/orientation.py b/src/mritk/data/orientation.py index 67be168..03ae413 100644 --- a/src/mritk/data/orientation.py +++ b/src/mritk/data/orientation.py @@ -210,3 +210,4 @@ def assert_same_space(mri1: MRIData, mri2: MRIData, rtol: float = 1e-5): ) raise ValueError(msg) + diff --git a/src/mritk/t1_maps/dicom_to_nifti.py b/src/mritk/t1_maps/dicom_to_nifti.py index 5ced6d4..3587c8f 100644 --- a/src/mritk/t1_maps/dicom_to_nifti.py +++ b/src/mritk/t1_maps/dicom_to_nifti.py @@ -16,7 +16,8 @@ import numpy as np from ..data.io import load_mri_data, save_mri_data -from ..dicom.utils import VOLUME_LABELS, read_dicom_trigger_times +from ..t1_maps.utils import VOLUME_LABELS, read_dicom_trigger_times +from .utils import extract_single_volume, logger def extract_mixed_dicom(dcmpath: Path, subvolumes: list[str]): diff --git a/src/mritk/t1_maps/t1_maps.py b/src/mritk/t1_maps/t1_maps.py index 44091bc..a0e19fb 100644 --- a/src/mritk/t1_maps/t1_maps.py +++ b/src/mritk/t1_maps/t1_maps.py @@ -16,7 +16,7 @@ import scipy from ..data.base import MRIData -from ..data.io import load_mri_data +from ..data.io import load_mri_data, save_mri_data from ..masking.masks import create_csf_mask from .utils import ( mri_facemask, @@ -26,7 +26,17 @@ ) -def looklocker_t1map(t_data: np.ndarray, D: np.ndarray, affine: np.ndarray) -> MRIData: +def looklocker_t1map( + looklocker_input: Path, + timestamps: Path, + output: Path = None +) -> MRIData: + + LL_mri = load_mri_data(looklocker_input, dtype=np.single) + D = LL_mri.data + affine = LL_mri.affine + t_data = np.loadtxt(timestamps) / 1000 + T1_ROOF = 10000 assert len(D.shape) >= 4, f"data should be 4-dimensional, got data with shape {D.shape}" mask = mri_facemask(D[..., 0]) @@ -52,7 +62,11 @@ def looklocker_t1map(t_data: np.ndarray, D: np.ndarray, affine: np.ndarray) -> M T1map = np.nan * np.zeros_like(D[..., 0]) T1map[I, J, K] = (x2 / x3) ** 2 * 1000.0 # convert to ms T1map = np.minimum(T1map, T1_ROOF) - return MRIData(T1map.astype(np.single), affine) + T1map_mri = MRIData(T1map.astype(np.single), affine) + if output is not None: + save_mri_data(T1map_mri, output, dtype=np.single) + + return T1map_mri def looklocker_t1map_postprocessing( @@ -134,7 +148,7 @@ def mixed_t1map_postprocessing(se: Path, t1: Path, output: Path): nibabel.nifti1.save(masked_T1map_nii, output) -def hybrid_t1_map( +def hybrid_t1map( ll_path: Path, mixed_path: Path, csf_mask_path: Path, From 2acf187b912043a0de0f78f5873fa7c7a91ac1e6 Mon Sep 17 00:00:00 2001 From: "pre-commit-ci[bot]" <66853113+pre-commit-ci[bot]@users.noreply.github.com> Date: Thu, 5 Mar 2026 08:30:54 +0000 Subject: [PATCH 04/24] [pre-commit.ci] auto fixes from pre-commit.com hooks for more information, see https://pre-commit.ci --- src/mritk/data/orientation.py | 1 - src/mritk/t1_maps/t1_maps.py | 7 +------ 2 files changed, 1 insertion(+), 7 deletions(-) diff --git a/src/mritk/data/orientation.py b/src/mritk/data/orientation.py index 03ae413..67be168 100644 --- a/src/mritk/data/orientation.py +++ b/src/mritk/data/orientation.py @@ -210,4 +210,3 @@ def assert_same_space(mri1: MRIData, mri2: MRIData, rtol: float = 1e-5): ) raise ValueError(msg) - diff --git a/src/mritk/t1_maps/t1_maps.py b/src/mritk/t1_maps/t1_maps.py index a0e19fb..2924b8c 100644 --- a/src/mritk/t1_maps/t1_maps.py +++ b/src/mritk/t1_maps/t1_maps.py @@ -26,12 +26,7 @@ ) -def looklocker_t1map( - looklocker_input: Path, - timestamps: Path, - output: Path = None -) -> MRIData: - +def looklocker_t1map(looklocker_input: Path, timestamps: Path, output: Path = None) -> MRIData: LL_mri = load_mri_data(looklocker_input, dtype=np.single) D = LL_mri.data affine = LL_mri.affine From 056ddee6c7fccf4735d51af4a38137b7f0fc4671 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?C=C3=A9cile?= Date: Thu, 5 Mar 2026 12:47:37 +0100 Subject: [PATCH 05/24] Fix downloading test data command --- CONTRIBUTING.md | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/CONTRIBUTING.md b/CONTRIBUTING.md index 6f4c966..dd573cc 100644 --- a/CONTRIBUTING.md +++ b/CONTRIBUTING.md @@ -90,7 +90,7 @@ We use **pytest** for testing. Before submitting a PR, ensure all tests pass loc Some tests require specific data. You can download this using the CLI included in the toolkit. ```bash # Downloads data to the 'test_data' folder (or your preferred location) - python -m mritk download-test-data test_data + python -m mritk datasets download test-data -o test_data ``` *Note: You may need to set the `MRITK_TEST_DATA_FOLDER` environment variable if you download the data to a custom location.* From 923600834040b6d403fccf6cda101a8ced4e546e Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?C=C3=A9cile?= Date: Thu, 5 Mar 2026 13:15:58 +0100 Subject: [PATCH 06/24] Add mri-dataset (from Gonzo) to the test-data dataset --- src/mritk/datasets.py | 1 + 1 file changed, 1 insertion(+) diff --git a/src/mritk/datasets.py b/src/mritk/datasets.py index e44834a..1ccd553 100644 --- a/src/mritk/datasets.py +++ b/src/mritk/datasets.py @@ -33,6 +33,7 @@ def get_datasets() -> dict[str, Dataset]: license="CC-BY-4.0", links={ "mri-processed.zip": "https://zenodo.org/records/14266867/files/mri-processed.zip?download=1", + "mri-dataset.zip": "https://zenodo.org/records/14266867/files/mri-dataset.zip?download=1", "timetable.tsv": "https://github.com/jorgenriseth/gonzo/blob/main/mri_dataset/timetable.tsv?raw=true", }, ), From f12d58fa0eb50aaf48ec44783d2d4b5ef781823d Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?C=C3=A9cile?= Date: Fri, 6 Mar 2026 09:43:17 +0100 Subject: [PATCH 07/24] Updates in T1maps functions and add test for T1maps (mixed, LL, hybrid) --- src/mritk/t1_maps/t1_maps.py | 62 ++++++++++++++++++++------- src/mritk/t1_maps/utils.py | 48 ++++++++++++++++++++- test/test_mri_t1_maps.py | 82 ++++++++++++++++++++++++++++++++++++ 3 files changed, 176 insertions(+), 16 deletions(-) create mode 100644 test/test_mri_t1_maps.py diff --git a/src/mritk/t1_maps/t1_maps.py b/src/mritk/t1_maps/t1_maps.py index 2924b8c..fb1a0c8 100644 --- a/src/mritk/t1_maps/t1_maps.py +++ b/src/mritk/t1_maps/t1_maps.py @@ -26,7 +26,11 @@ ) -def looklocker_t1map(looklocker_input: Path, timestamps: Path, output: Path = None) -> MRIData: +def looklocker_t1map( + looklocker_input: Path, + timestamps: Path, + output: Path = None +) -> MRIData: LL_mri = load_mri_data(looklocker_input, dtype=np.single) D = LL_mri.data affine = LL_mri.affine @@ -66,11 +70,12 @@ def looklocker_t1map(looklocker_input: Path, timestamps: Path, output: Path = No def looklocker_t1map_postprocessing( T1map_mri: MRIData, - T1_lo: float, - T1_hi: float, + T1_low: float, + T1_high: float, radius: int = 10, erode_dilate_factor: float = 1.3, mask: Optional[np.ndarray] = None, + output: Path = None ) -> MRIData: T1map = T1map_mri.data.copy() if mask is None: @@ -89,8 +94,8 @@ def looklocker_t1map_postprocessing( T1map[~mask] = np.nan # Remove outliers within the mask. - outliers = np.logical_or(T1map < T1_lo, T1_hi < T1map) - print("Removing", outliers.sum(), f"voxels outside the range ({T1_lo}, {T1_hi}).") + outliers = np.logical_or(T1map < T1_low, T1_high < T1map) + print("Removing", outliers.sum(), f"voxels outside the range ({T1_low}, {T1_high}).") T1map[outliers] = np.nan if np.isfinite(T1map).sum() / T1map.size < 0.01: raise RuntimeError("After outlier removal, less than 1% of the image is left. Check image units.") @@ -101,7 +106,12 @@ def looklocker_t1map_postprocessing( print(f"Filling in {fill_mask.sum()} voxels within the mask.") T1map[fill_mask] = nan_filter_gaussian(T1map, 1.0)[fill_mask] fill_mask = np.isnan(T1map) * mask - return MRIData(T1map, T1map_mri.affine) + + processed_T1map = MRIData(T1map, T1map_mri.affine) + if output is not None: + save_mri_data(processed_T1map, output, dtype=np.single) + + return processed_T1map def mixed_t1map( @@ -109,7 +119,8 @@ def mixed_t1map( IR_nii_path: Path, meta_path: Path, T1_low: float, - T1_hi: float, + T1_high: float, + output: Path = None ) -> nibabel.nifti1.Nifti1Image: SE = load_mri_data(SE_nii_path, dtype=np.single) IR = load_mri_data(IR_nii_path, dtype=np.single) @@ -121,39 +132,55 @@ def mixed_t1map( F_data[nonzero_mask] = IR.data[nonzero_mask] / SE.data[nonzero_mask] TR_se, TI, TE, ETL = meta["TR_SE"], meta["TI"], meta["TE"], meta["ETL"] - F, T1_grid = T1_lookup_table(TR_se, TI, TE, ETL, T1_low, T1_hi) + F, T1_grid = T1_lookup_table(TR_se, TI, TE, ETL, T1_low, T1_high) interpolator = scipy.interpolate.interp1d(F, T1_grid, kind="nearest", bounds_error=False, fill_value=np.nan) T1_volume = interpolator(F_data).astype(np.single) nii = nibabel.nifti1.Nifti1Image(T1_volume, IR.affine) nii.set_sform(nii.affine, "scanner") nii.set_qform(nii.affine, "scanner") + + # T1map_mri = MRIData(T1_volume, nii.affine) + if output is not None: + nibabel.nifti1.save(nii, output) + # save_mri_data(T1map_mri, output, dtype=np.single) + return nii + # return T1map_mri -def mixed_t1map_postprocessing(se: Path, t1: Path, output: Path): - T1map_nii = nibabel.nifti1.load(t1) +def mixed_t1map_postprocessing( + SE_nii_path: Path, + T1_path: Path, + output: Path = None +) -> nibabel.nifti1.Nifti1Image: + T1map_nii = nibabel.nifti1.load(T1_path) - SE_mri = load_mri_data(se, np.single) + SE_mri = load_mri_data(SE_nii_path, np.single) mask = create_csf_mask(SE_mri.data, use_li=True) mask = skimage.morphology.binary_erosion(mask) masked_T1map = T1map_nii.get_fdata(dtype=np.single) masked_T1map[~mask] = np.nan masked_T1map_nii = nibabel.nifti1.Nifti1Image(masked_T1map, T1map_nii.affine, T1map_nii.header) - nibabel.nifti1.save(masked_T1map_nii, output) + + if output is not None: + nibabel.nifti1.save(masked_T1map_nii, output) + + return masked_T1map_nii def hybrid_t1map( - ll_path: Path, + LL_path: Path, mixed_path: Path, csf_mask_path: Path, threshold: float, erode: int = 0, + output: Path = None ) -> nibabel.nifti1.Nifti1Image: mixed_mri = nibabel.nifti1.load(mixed_path) mixed = mixed_mri.get_fdata() - ll_mri = nibabel.nifti1.load(ll_path) + ll_mri = nibabel.nifti1.load(LL_path) ll = ll_mri.get_fdata() csf_mask_mri = nibabel.nifti1.load(csf_mask_path) csf_mask = csf_mask_mri.get_fdata().astype(bool) @@ -163,4 +190,9 @@ def hybrid_t1map( hybrid = ll newmask = csf_mask * (ll > threshold) * (mixed > threshold) hybrid[newmask] = mixed[newmask] - return nibabel.nifti1.Nifti1Image(hybrid, affine=ll_mri.affine, header=ll_mri.header) + + hybrid_nii = nibabel.nifti1.Nifti1Image(hybrid, affine=ll_mri.affine, header=ll_mri.header) + if output is not None: + nibabel.nifti1.save(hybrid_nii, output) + + return hybrid_nii diff --git a/src/mritk/t1_maps/utils.py b/src/mritk/t1_maps/utils.py index 2ccf14c..547ce05 100644 --- a/src/mritk/t1_maps/utils.py +++ b/src/mritk/t1_maps/utils.py @@ -11,6 +11,8 @@ import skimage import warnings from scipy.optimize import OptimizeWarning +import os +import nibabel from ..data.orientation import data_reorientation, change_of_coordinates_map from ..data.base import MRIData @@ -99,6 +101,10 @@ def mri_facemask(vol: np.ndarray, smoothing_level=5): return binary +def voxel_fit_function(t, x1, x2, x3): + return np.abs(x1 * (1.0 - (1 + x2**2) * np.exp(-(x3**2) * t))) + +@np.errstate(divide="raise", invalid="raise", over="raise") def curve_fit_wrapper(f, t, y, p0): """Raises error instead of catching numpy warnings, such that these cases may be treated.""" @@ -119,7 +125,7 @@ def fit_voxel(time_s: np.ndarray, pbar, m: np.ndarray) -> np.ndarray: if not np.all(np.isfinite(m)): return np.nan * np.zeros_like(p0) try: - popt = curve_fit_wrapper(f, time_s, m, p0) + popt = curve_fit_wrapper(voxel_fit_function, time_s, m, p0) except (OptimizeWarning, FloatingPointError): return np.nan * np.zeros_like(p0) except RuntimeError as e: @@ -156,3 +162,43 @@ def T1_lookup_table(TRse: float, TI: float, TE: float, ETL: int, T1_low: float, Sir = 1 - (1 + Sse) * np.exp(-TI / T1_grid) fractionCurve = Sir / Sse return fractionCurve, T1_grid + +def compare_nifti_images(img_path1, img_path2, data_tolerance=0.0): + """ + Compares two NIfTI images for equality of data, affine, and header. + + Args: + img_path1 (str): Path to the first NIfTI file. + img_path2 (str): Path to the second NIfTI file. + data_tolerance (float): Tolerance for data comparison (use 0.0 for exact equality). + + Returns: + bool: True if images are considered the same, False otherwise. + list: A list of differences found. + """ + if not os.path.exists(img_path1): + return False, [f"File not found: {img_path1}"] + if not os.path.exists(img_path2): + return False, [f"File not found: {img_path2}"] + + img1 = nibabel.load(img_path1) + img2 = nibabel.load(img_path2) + differences = [] + + # 1. Compare Image Data + data1 = img1.get_fdata() + data2 = img2.get_fdata() + # Use np.allclose for data comparison with tolerance, which is often needed + # for floating-point data, or np.array_equal for exact comparison. + if data_tolerance > 0: + data_equal = np.allclose(data1, data2, atol=data_tolerance) + else: + data_equal = np.array_equal(data1, data2) + + deviation = np.mean(np.abs(data1 - data2)) + assert data_equal, f"Data mismatch (mean absolute deviation: {deviation:.4f})" + + # Overall result + is_same = not differences + assert False + return is_same, differences \ No newline at end of file diff --git a/test/test_mri_t1_maps.py b/test/test_mri_t1_maps.py new file mode 100644 index 0000000..fbb8f2b --- /dev/null +++ b/test/test_mri_t1_maps.py @@ -0,0 +1,82 @@ +"""MRI T1 maps - Tests + +Copyright (C) 2026 Jørgen Riseth (jnriseth@gmail.com) +Copyright (C) 2026 Cécile Daversin-Catty (cecile@simula.no) +Copyright (C) 2026 Simula Research Laboratory +""" + +from pathlib import Path + + +from mritk.t1_maps.dicom_to_nifti import dicom_to_looklocker, dicom_to_mixed +from mritk.t1_maps.t1_maps import ( + looklocker_t1map, + looklocker_t1map_postprocessing, + mixed_t1map, + mixed_t1map_postprocessing, + hybrid_t1map +) +import mritk.cli as cli + +from mritk.t1_maps.utils import compare_nifti_images + + +def test_looklocker_t1map(tmp_path, mri_data_dir: Path): + + LL_path = mri_data_dir / "mri-dataset/mri_dataset/sub-01" / "ses-01/anat/sub-01_ses-01_acq-looklocker_IRT1.nii.gz" + timestamps = mri_data_dir / "mri-dataset/mri_dataset/sub-01" / "ses-01/anat/sub-01_ses-01_acq-looklocker_IRT1_trigger_times.txt" + ## CHECK : In snakemake scripts, output is T1raw="mri_dataset/derivatives/{subject}/{session}/{subject}_{session}_acq-looklocker_T1map_raw.nii.gz" + ## but this file is not in the dataset ? + ref_output = mri_data_dir / "mri-dataset/mri_dataset/derivatives/sub-01" / "ses-01/sub-01_ses-01_acq-looklocker_T1map_nICE.nii.gz" + test_output = tmp_path / "output_acq-looklocker_T1map_raw.nii.gz" + + looklocker_t1map(looklocker_input=LL_path, timestamps=timestamps, output=test_output) + compare_nifti_images(test_output, ref_output, data_tolerance=1e-4) + + +def test_looklocker_t1map_postprocessing(tmp_path, mri_data_dir: Path): + T1_path = tmp_path / "output_acq-looklocker_T1map_raw.nii.gz" + ref_output = mri_data_dir / "mri-dataset/mri_dataset/derivatives/sub-01" / "ses-01/sub-01_ses-01_acq-looklocker_T1map.nii.gz", + test_output = tmp_path / "output_acq-looklocker_T1map.nii.gz" + + T1_low=100 + T1_high=6000 + + looklocker_t1map_postprocessing(T1map_mri=T1_path, T1_low=T1_low, T1_high=T1_high, output=test_output) + compare_nifti_images(test_output, ref_output, data_tolerance=1e-4) + + +def test_mixed_t1map(tmp_path, mri_data_dir: Path): + SE_path = mri_data_dir / "mri-dataset/mri_dataset/sub-01" / "ses-01/mixed/sub-01_ses-01_acq-mixed_SE-modulus.nii.gz" + IR_path = mri_data_dir / "mri-dataset/mri_dataset/sub-01" / "ses-01/mixed/sub-01_ses-01_acq-mixed_IR-corrected-real.nii.gz" + meta_path = mri_data_dir / "mri-dataset/mri_dataset/sub-01" / "ses-01/mixed/sub-01_ses-01_acq-mixed_meta.json" + ref_output = mri_data_dir / "mri-dataset/mri_dataset/derivatives/sub-01" / "ses-01/sub-01_ses-01_acq-mixed_T1map_raw.nii.gz", + test_output = tmp_path / "output_acq-mixed_T1map_raw.nii.gz" + + T1_low=100 + T1_high=10000 + + mixed_t1map(SE_nii_path=SE_path, IR_nii_path=IR_path, meta_path=meta_path, T1_low=T1_low, T1_high=T1_high, output=test_output) + compare_nifti_images(test_output, ref_output, data_tolerance=1e-4) + + +def test_mixed_t1map_postprocessing(tmp_path, mri_data_dir: Path): + SE_path = mri_data_dir / "mri-dataset/mri_dataset/sub-01" / "ses-01/mixed/sub-01_ses-01_acq-mixed_SE-modulus.nii.gz" + T1_path = tmp_path / "output_acq-mixed_T1map_raw.nii.gz" + ref_output = mri_data_dir / "mri-dataset/mri_dataset/derivatives/sub-01" / "ses-01/sub-01_ses-01_acq-mixed_T1map.nii.gz", + test_output = tmp_path / "output_acq-mixed_T1map.nii.gz" + + mixed_t1map_postprocessing(SE_nii_path=SE_path, T1_path=T1_path, output=test_output) + compare_nifti_images(test_output, ref_output, data_tolerance=1e-4) + +def test_hybrid_t1map(tmp_path, mri_data_dir: Path): + LL_path = mri_data_dir / "mri-processed/mri_processed_data/sub-01/registered/sub-01_ses-01_acq-looklocker_T1map_registered.nii.gz" + mixed_path = mri_data_dir / "mri-processed/mri_processed_data/sub-01/registered/sub-01_ses-01_acq-mixed_T1map_registered.nii.gz" + csf_mask_path = mri_data_dir / "mri-processed/mri_processed_data/sub-01/segmentations/sub-01_seg-csf_binary.nii.gz" + test_output = tmp_path / "output_T1map_hybrid.nii.gz" + ref_output = mri_data_dir / "mri-processed/mri_processed_data/sub-01/T1maps/sub-01_ses-01_T1map_hybrid.nii.gz" + threshold = 1500 + erode = 1 + + hybrid_t1map(LL_path=LL_path, mixed_path=mixed_path, csf_mask_path=csf_mask_path, threshold=threshold, erode=erode, output=test_output) + compare_nifti_images(test_output, ref_output, data_tolerance=1e-4) \ No newline at end of file From 73415d66c0c65160294a41a72f86b229a96a4fc8 Mon Sep 17 00:00:00 2001 From: "pre-commit-ci[bot]" <66853113+pre-commit-ci[bot]@users.noreply.github.com> Date: Fri, 6 Mar 2026 08:43:41 +0000 Subject: [PATCH 08/24] [pre-commit.ci] auto fixes from pre-commit.com hooks for more information, see https://pre-commit.ci --- src/mritk/t1_maps/t1_maps.py | 28 ++++-------------- src/mritk/t1_maps/utils.py | 6 ++-- test/test_mri_t1_maps.py | 56 ++++++++++++++++++++++-------------- 3 files changed, 43 insertions(+), 47 deletions(-) diff --git a/src/mritk/t1_maps/t1_maps.py b/src/mritk/t1_maps/t1_maps.py index fb1a0c8..07aa516 100644 --- a/src/mritk/t1_maps/t1_maps.py +++ b/src/mritk/t1_maps/t1_maps.py @@ -26,11 +26,7 @@ ) -def looklocker_t1map( - looklocker_input: Path, - timestamps: Path, - output: Path = None -) -> MRIData: +def looklocker_t1map(looklocker_input: Path, timestamps: Path, output: Path = None) -> MRIData: LL_mri = load_mri_data(looklocker_input, dtype=np.single) D = LL_mri.data affine = LL_mri.affine @@ -75,7 +71,7 @@ def looklocker_t1map_postprocessing( radius: int = 10, erode_dilate_factor: float = 1.3, mask: Optional[np.ndarray] = None, - output: Path = None + output: Path = None, ) -> MRIData: T1map = T1map_mri.data.copy() if mask is None: @@ -115,12 +111,7 @@ def looklocker_t1map_postprocessing( def mixed_t1map( - SE_nii_path: Path, - IR_nii_path: Path, - meta_path: Path, - T1_low: float, - T1_high: float, - output: Path = None + SE_nii_path: Path, IR_nii_path: Path, meta_path: Path, T1_low: float, T1_high: float, output: Path = None ) -> nibabel.nifti1.Nifti1Image: SE = load_mri_data(SE_nii_path, dtype=np.single) IR = load_mri_data(IR_nii_path, dtype=np.single) @@ -148,11 +139,7 @@ def mixed_t1map( # return T1map_mri -def mixed_t1map_postprocessing( - SE_nii_path: Path, - T1_path: Path, - output: Path = None -) -> nibabel.nifti1.Nifti1Image: +def mixed_t1map_postprocessing(SE_nii_path: Path, T1_path: Path, output: Path = None) -> nibabel.nifti1.Nifti1Image: T1map_nii = nibabel.nifti1.load(T1_path) SE_mri = load_mri_data(SE_nii_path, np.single) @@ -170,12 +157,7 @@ def mixed_t1map_postprocessing( def hybrid_t1map( - LL_path: Path, - mixed_path: Path, - csf_mask_path: Path, - threshold: float, - erode: int = 0, - output: Path = None + LL_path: Path, mixed_path: Path, csf_mask_path: Path, threshold: float, erode: int = 0, output: Path = None ) -> nibabel.nifti1.Nifti1Image: mixed_mri = nibabel.nifti1.load(mixed_path) mixed = mixed_mri.get_fdata() diff --git a/src/mritk/t1_maps/utils.py b/src/mritk/t1_maps/utils.py index 547ce05..c413bb4 100644 --- a/src/mritk/t1_maps/utils.py +++ b/src/mritk/t1_maps/utils.py @@ -104,6 +104,7 @@ def mri_facemask(vol: np.ndarray, smoothing_level=5): def voxel_fit_function(t, x1, x2, x3): return np.abs(x1 * (1.0 - (1 + x2**2) * np.exp(-(x3**2) * t))) + @np.errstate(divide="raise", invalid="raise", over="raise") def curve_fit_wrapper(f, t, y, p0): """Raises error instead of catching numpy warnings, such that @@ -163,6 +164,7 @@ def T1_lookup_table(TRse: float, TI: float, TE: float, ETL: int, T1_low: float, fractionCurve = Sir / Sse return fractionCurve, T1_grid + def compare_nifti_images(img_path1, img_path2, data_tolerance=0.0): """ Compares two NIfTI images for equality of data, affine, and header. @@ -188,7 +190,7 @@ def compare_nifti_images(img_path1, img_path2, data_tolerance=0.0): # 1. Compare Image Data data1 = img1.get_fdata() data2 = img2.get_fdata() - # Use np.allclose for data comparison with tolerance, which is often needed + # Use np.allclose for data comparison with tolerance, which is often needed # for floating-point data, or np.array_equal for exact comparison. if data_tolerance > 0: data_equal = np.allclose(data1, data2, atol=data_tolerance) @@ -201,4 +203,4 @@ def compare_nifti_images(img_path1, img_path2, data_tolerance=0.0): # Overall result is_same = not differences assert False - return is_same, differences \ No newline at end of file + return is_same, differences diff --git a/test/test_mri_t1_maps.py b/test/test_mri_t1_maps.py index fbb8f2b..48d85a7 100644 --- a/test/test_mri_t1_maps.py +++ b/test/test_mri_t1_maps.py @@ -8,26 +8,27 @@ from pathlib import Path -from mritk.t1_maps.dicom_to_nifti import dicom_to_looklocker, dicom_to_mixed from mritk.t1_maps.t1_maps import ( - looklocker_t1map, - looklocker_t1map_postprocessing, - mixed_t1map, - mixed_t1map_postprocessing, - hybrid_t1map + looklocker_t1map, + looklocker_t1map_postprocessing, + mixed_t1map, + mixed_t1map_postprocessing, + hybrid_t1map, ) -import mritk.cli as cli from mritk.t1_maps.utils import compare_nifti_images def test_looklocker_t1map(tmp_path, mri_data_dir: Path): - LL_path = mri_data_dir / "mri-dataset/mri_dataset/sub-01" / "ses-01/anat/sub-01_ses-01_acq-looklocker_IRT1.nii.gz" - timestamps = mri_data_dir / "mri-dataset/mri_dataset/sub-01" / "ses-01/anat/sub-01_ses-01_acq-looklocker_IRT1_trigger_times.txt" + timestamps = ( + mri_data_dir / "mri-dataset/mri_dataset/sub-01" / "ses-01/anat/sub-01_ses-01_acq-looklocker_IRT1_trigger_times.txt" + ) ## CHECK : In snakemake scripts, output is T1raw="mri_dataset/derivatives/{subject}/{session}/{subject}_{session}_acq-looklocker_T1map_raw.nii.gz" ## but this file is not in the dataset ? - ref_output = mri_data_dir / "mri-dataset/mri_dataset/derivatives/sub-01" / "ses-01/sub-01_ses-01_acq-looklocker_T1map_nICE.nii.gz" + ref_output = ( + mri_data_dir / "mri-dataset/mri_dataset/derivatives/sub-01" / "ses-01/sub-01_ses-01_acq-looklocker_T1map_nICE.nii.gz" + ) test_output = tmp_path / "output_acq-looklocker_T1map_raw.nii.gz" looklocker_t1map(looklocker_input=LL_path, timestamps=timestamps, output=test_output) @@ -35,12 +36,14 @@ def test_looklocker_t1map(tmp_path, mri_data_dir: Path): def test_looklocker_t1map_postprocessing(tmp_path, mri_data_dir: Path): - T1_path = tmp_path / "output_acq-looklocker_T1map_raw.nii.gz" - ref_output = mri_data_dir / "mri-dataset/mri_dataset/derivatives/sub-01" / "ses-01/sub-01_ses-01_acq-looklocker_T1map.nii.gz", + T1_path = tmp_path / "output_acq-looklocker_T1map_raw.nii.gz" + ref_output = ( + mri_data_dir / "mri-dataset/mri_dataset/derivatives/sub-01" / "ses-01/sub-01_ses-01_acq-looklocker_T1map.nii.gz", + ) test_output = tmp_path / "output_acq-looklocker_T1map.nii.gz" - T1_low=100 - T1_high=6000 + T1_low = 100 + T1_high = 6000 looklocker_t1map_postprocessing(T1map_mri=T1_path, T1_low=T1_low, T1_high=T1_high, output=test_output) compare_nifti_images(test_output, ref_output, data_tolerance=1e-4) @@ -50,11 +53,13 @@ def test_mixed_t1map(tmp_path, mri_data_dir: Path): SE_path = mri_data_dir / "mri-dataset/mri_dataset/sub-01" / "ses-01/mixed/sub-01_ses-01_acq-mixed_SE-modulus.nii.gz" IR_path = mri_data_dir / "mri-dataset/mri_dataset/sub-01" / "ses-01/mixed/sub-01_ses-01_acq-mixed_IR-corrected-real.nii.gz" meta_path = mri_data_dir / "mri-dataset/mri_dataset/sub-01" / "ses-01/mixed/sub-01_ses-01_acq-mixed_meta.json" - ref_output = mri_data_dir / "mri-dataset/mri_dataset/derivatives/sub-01" / "ses-01/sub-01_ses-01_acq-mixed_T1map_raw.nii.gz", + ref_output = ( + mri_data_dir / "mri-dataset/mri_dataset/derivatives/sub-01" / "ses-01/sub-01_ses-01_acq-mixed_T1map_raw.nii.gz", + ) test_output = tmp_path / "output_acq-mixed_T1map_raw.nii.gz" - T1_low=100 - T1_high=10000 + T1_low = 100 + T1_high = 10000 mixed_t1map(SE_nii_path=SE_path, IR_nii_path=IR_path, meta_path=meta_path, T1_low=T1_low, T1_high=T1_high, output=test_output) compare_nifti_images(test_output, ref_output, data_tolerance=1e-4) @@ -63,20 +68,27 @@ def test_mixed_t1map(tmp_path, mri_data_dir: Path): def test_mixed_t1map_postprocessing(tmp_path, mri_data_dir: Path): SE_path = mri_data_dir / "mri-dataset/mri_dataset/sub-01" / "ses-01/mixed/sub-01_ses-01_acq-mixed_SE-modulus.nii.gz" T1_path = tmp_path / "output_acq-mixed_T1map_raw.nii.gz" - ref_output = mri_data_dir / "mri-dataset/mri_dataset/derivatives/sub-01" / "ses-01/sub-01_ses-01_acq-mixed_T1map.nii.gz", + ref_output = (mri_data_dir / "mri-dataset/mri_dataset/derivatives/sub-01" / "ses-01/sub-01_ses-01_acq-mixed_T1map.nii.gz",) test_output = tmp_path / "output_acq-mixed_T1map.nii.gz" mixed_t1map_postprocessing(SE_nii_path=SE_path, T1_path=T1_path, output=test_output) compare_nifti_images(test_output, ref_output, data_tolerance=1e-4) + def test_hybrid_t1map(tmp_path, mri_data_dir: Path): - LL_path = mri_data_dir / "mri-processed/mri_processed_data/sub-01/registered/sub-01_ses-01_acq-looklocker_T1map_registered.nii.gz" - mixed_path = mri_data_dir / "mri-processed/mri_processed_data/sub-01/registered/sub-01_ses-01_acq-mixed_T1map_registered.nii.gz" + LL_path = ( + mri_data_dir / "mri-processed/mri_processed_data/sub-01/registered/sub-01_ses-01_acq-looklocker_T1map_registered.nii.gz" + ) + mixed_path = ( + mri_data_dir / "mri-processed/mri_processed_data/sub-01/registered/sub-01_ses-01_acq-mixed_T1map_registered.nii.gz" + ) csf_mask_path = mri_data_dir / "mri-processed/mri_processed_data/sub-01/segmentations/sub-01_seg-csf_binary.nii.gz" test_output = tmp_path / "output_T1map_hybrid.nii.gz" ref_output = mri_data_dir / "mri-processed/mri_processed_data/sub-01/T1maps/sub-01_ses-01_T1map_hybrid.nii.gz" threshold = 1500 erode = 1 - hybrid_t1map(LL_path=LL_path, mixed_path=mixed_path, csf_mask_path=csf_mask_path, threshold=threshold, erode=erode, output=test_output) - compare_nifti_images(test_output, ref_output, data_tolerance=1e-4) \ No newline at end of file + hybrid_t1map( + LL_path=LL_path, mixed_path=mixed_path, csf_mask_path=csf_mask_path, threshold=threshold, erode=erode, output=test_output + ) + compare_nifti_images(test_output, ref_output, data_tolerance=1e-4) From 6c978b551932941a28fc08e90ba7d003f86d1efc Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?C=C3=A9cile?= Date: Fri, 6 Mar 2026 09:47:30 +0100 Subject: [PATCH 09/24] Fix compare_nifti_images --- src/mritk/t1_maps/utils.py | 6 +----- 1 file changed, 1 insertion(+), 5 deletions(-) diff --git a/src/mritk/t1_maps/utils.py b/src/mritk/t1_maps/utils.py index 547ce05..8d45611 100644 --- a/src/mritk/t1_maps/utils.py +++ b/src/mritk/t1_maps/utils.py @@ -174,7 +174,6 @@ def compare_nifti_images(img_path1, img_path2, data_tolerance=0.0): Returns: bool: True if images are considered the same, False otherwise. - list: A list of differences found. """ if not os.path.exists(img_path1): return False, [f"File not found: {img_path1}"] @@ -183,7 +182,6 @@ def compare_nifti_images(img_path1, img_path2, data_tolerance=0.0): img1 = nibabel.load(img_path1) img2 = nibabel.load(img_path2) - differences = [] # 1. Compare Image Data data1 = img1.get_fdata() @@ -199,6 +197,4 @@ def compare_nifti_images(img_path1, img_path2, data_tolerance=0.0): assert data_equal, f"Data mismatch (mean absolute deviation: {deviation:.4f})" # Overall result - is_same = not differences - assert False - return is_same, differences \ No newline at end of file + return data_equal \ No newline at end of file From 3c945a8683fe6edcf1464392b0f906a2eecf831d Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?C=C3=A9cile?= Date: Fri, 6 Mar 2026 09:56:32 +0100 Subject: [PATCH 10/24] Replace binary_erosion (deprecated) by erosion --- src/mritk/t1_maps/t1_maps.py | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/src/mritk/t1_maps/t1_maps.py b/src/mritk/t1_maps/t1_maps.py index 07aa516..936c269 100644 --- a/src/mritk/t1_maps/t1_maps.py +++ b/src/mritk/t1_maps/t1_maps.py @@ -82,7 +82,7 @@ def looklocker_t1map_postprocessing( mask = mask == regions[0].label skimage.morphology.remove_small_holes(mask, 10 ** (mask.ndim), connectivity=2, out=mask) skimage.morphology.binary_dilation(mask, skimage.morphology.ball(radius), out=mask) - skimage.morphology.binary_erosion(mask, skimage.morphology.ball(erode_dilate_factor * radius), out=mask) + skimage.morphology.erosion(mask, skimage.morphology.ball(erode_dilate_factor * radius), out=mask) # Remove non-zero artifacts outside of the mask. surface_vox = np.isfinite(T1map) * (~mask) @@ -144,7 +144,7 @@ def mixed_t1map_postprocessing(SE_nii_path: Path, T1_path: Path, output: Path = SE_mri = load_mri_data(SE_nii_path, np.single) mask = create_csf_mask(SE_mri.data, use_li=True) - mask = skimage.morphology.binary_erosion(mask) + mask = skimage.morphology.erosion(mask) masked_T1map = T1map_nii.get_fdata(dtype=np.single) masked_T1map[~mask] = np.nan @@ -167,7 +167,7 @@ def hybrid_t1map( csf_mask_mri = nibabel.nifti1.load(csf_mask_path) csf_mask = csf_mask_mri.get_fdata().astype(bool) if erode > 0: - csf_mask = skimage.morphology.binary_erosion(csf_mask, skimage.morphology.ball(erode)) + csf_mask = skimage.morphology.erosion(csf_mask, skimage.morphology.ball(erode)) hybrid = ll newmask = csf_mask * (ll > threshold) * (mixed > threshold) From f2d7e9c58dbbe0e1b6f38fce70f6260d4a4243bd Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?C=C3=A9cile?= Date: Fri, 6 Mar 2026 10:43:09 +0100 Subject: [PATCH 11/24] Skipping looklocker T1map test (too long) - Mixed and Hybrid T1maps tests now passing --- src/mritk/t1_maps/t1_maps.py | 18 +++++++----- test/test_mri_t1_maps.py | 53 ++++++++++++------------------------ 2 files changed, 28 insertions(+), 43 deletions(-) diff --git a/src/mritk/t1_maps/t1_maps.py b/src/mritk/t1_maps/t1_maps.py index 936c269..111a027 100644 --- a/src/mritk/t1_maps/t1_maps.py +++ b/src/mritk/t1_maps/t1_maps.py @@ -65,7 +65,7 @@ def looklocker_t1map(looklocker_input: Path, timestamps: Path, output: Path = No def looklocker_t1map_postprocessing( - T1map_mri: MRIData, + T1map: Path, T1_low: float, T1_high: float, radius: int = 10, @@ -73,6 +73,8 @@ def looklocker_t1map_postprocessing( mask: Optional[np.ndarray] = None, output: Path = None, ) -> MRIData: + + T1map_mri = load_mri_data(T1map, dtype=np.single) T1map = T1map_mri.data.copy() if mask is None: # Create mask for largest island. @@ -80,8 +82,8 @@ def looklocker_t1map_postprocessing( regions = skimage.measure.regionprops(mask) regions.sort(key=lambda x: x.num_pixels, reverse=True) mask = mask == regions[0].label - skimage.morphology.remove_small_holes(mask, 10 ** (mask.ndim), connectivity=2, out=mask) - skimage.morphology.binary_dilation(mask, skimage.morphology.ball(radius), out=mask) + skimage.morphology.remove_small_holes(mask, max_size=10 ** (mask.ndim), connectivity=2, out=mask) + skimage.morphology.dilation(mask, skimage.morphology.ball(radius), out=mask) skimage.morphology.erosion(mask, skimage.morphology.ball(erode_dilate_factor * radius), out=mask) # Remove non-zero artifacts outside of the mask. @@ -111,7 +113,12 @@ def looklocker_t1map_postprocessing( def mixed_t1map( - SE_nii_path: Path, IR_nii_path: Path, meta_path: Path, T1_low: float, T1_high: float, output: Path = None + SE_nii_path: Path, + IR_nii_path: Path, + meta_path: Path, + T1_low: float, + T1_high: float, + output: Path = None ) -> nibabel.nifti1.Nifti1Image: SE = load_mri_data(SE_nii_path, dtype=np.single) IR = load_mri_data(IR_nii_path, dtype=np.single) @@ -130,13 +137,10 @@ def mixed_t1map( nii.set_sform(nii.affine, "scanner") nii.set_qform(nii.affine, "scanner") - # T1map_mri = MRIData(T1_volume, nii.affine) if output is not None: nibabel.nifti1.save(nii, output) - # save_mri_data(T1map_mri, output, dtype=np.single) return nii - # return T1map_mri def mixed_t1map_postprocessing(SE_nii_path: Path, T1_path: Path, output: Path = None) -> nibabel.nifti1.Nifti1Image: diff --git a/test/test_mri_t1_maps.py b/test/test_mri_t1_maps.py index 48d85a7..c733d1e 100644 --- a/test/test_mri_t1_maps.py +++ b/test/test_mri_t1_maps.py @@ -6,7 +6,7 @@ """ from pathlib import Path - +import pytest from mritk.t1_maps.t1_maps import ( looklocker_t1map, @@ -18,61 +18,42 @@ from mritk.t1_maps.utils import compare_nifti_images - +@pytest.mark.skip(reason="Takes too long") def test_looklocker_t1map(tmp_path, mri_data_dir: Path): LL_path = mri_data_dir / "mri-dataset/mri_dataset/sub-01" / "ses-01/anat/sub-01_ses-01_acq-looklocker_IRT1.nii.gz" timestamps = ( mri_data_dir / "mri-dataset/mri_dataset/sub-01" / "ses-01/anat/sub-01_ses-01_acq-looklocker_IRT1_trigger_times.txt" ) - ## CHECK : In snakemake scripts, output is T1raw="mri_dataset/derivatives/{subject}/{session}/{subject}_{session}_acq-looklocker_T1map_raw.nii.gz" - ## but this file is not in the dataset ? - ref_output = ( - mri_data_dir / "mri-dataset/mri_dataset/derivatives/sub-01" / "ses-01/sub-01_ses-01_acq-looklocker_T1map_nICE.nii.gz" - ) - test_output = tmp_path / "output_acq-looklocker_T1map_raw.nii.gz" - - looklocker_t1map(looklocker_input=LL_path, timestamps=timestamps, output=test_output) - compare_nifti_images(test_output, ref_output, data_tolerance=1e-4) - + T1_low = 100 + T1_high = 6000 -def test_looklocker_t1map_postprocessing(tmp_path, mri_data_dir: Path): - T1_path = tmp_path / "output_acq-looklocker_T1map_raw.nii.gz" ref_output = ( - mri_data_dir / "mri-dataset/mri_dataset/derivatives/sub-01" / "ses-01/sub-01_ses-01_acq-looklocker_T1map.nii.gz", + mri_data_dir / "mri-dataset/mri_dataset/derivatives/sub-01" / "ses-01/sub-01_ses-01_acq-looklocker_T1map.nii.gz" ) + test_output_raw = tmp_path / "output_acq-looklocker_T1map_raw.nii.gz" test_output = tmp_path / "output_acq-looklocker_T1map.nii.gz" - T1_low = 100 - T1_high = 6000 - - looklocker_t1map_postprocessing(T1map_mri=T1_path, T1_low=T1_low, T1_high=T1_high, output=test_output) - compare_nifti_images(test_output, ref_output, data_tolerance=1e-4) + looklocker_t1map(looklocker_input=LL_path, timestamps=timestamps, output=test_output_raw) + looklocker_t1map_postprocessing(T1map=test_output_raw, T1_low=T1_low, T1_high=T1_high, output=test_output) + compare_nifti_images(test_output, ref_output, data_tolerance=1e-12) def test_mixed_t1map(tmp_path, mri_data_dir: Path): SE_path = mri_data_dir / "mri-dataset/mri_dataset/sub-01" / "ses-01/mixed/sub-01_ses-01_acq-mixed_SE-modulus.nii.gz" IR_path = mri_data_dir / "mri-dataset/mri_dataset/sub-01" / "ses-01/mixed/sub-01_ses-01_acq-mixed_IR-corrected-real.nii.gz" meta_path = mri_data_dir / "mri-dataset/mri_dataset/sub-01" / "ses-01/mixed/sub-01_ses-01_acq-mixed_meta.json" - ref_output = ( - mri_data_dir / "mri-dataset/mri_dataset/derivatives/sub-01" / "ses-01/sub-01_ses-01_acq-mixed_T1map_raw.nii.gz", - ) - test_output = tmp_path / "output_acq-mixed_T1map_raw.nii.gz" + + ref_output = mri_data_dir / "mri-dataset/mri_dataset/derivatives/sub-01" / "ses-01/sub-01_ses-01_acq-mixed_T1map.nii.gz" + test_output_raw = tmp_path / "output_acq-mixed_T1map_raw.nii.gz" + test_output = tmp_path / "output_acq-mixed_T1map.nii.gz" T1_low = 100 T1_high = 10000 - mixed_t1map(SE_nii_path=SE_path, IR_nii_path=IR_path, meta_path=meta_path, T1_low=T1_low, T1_high=T1_high, output=test_output) - compare_nifti_images(test_output, ref_output, data_tolerance=1e-4) - - -def test_mixed_t1map_postprocessing(tmp_path, mri_data_dir: Path): - SE_path = mri_data_dir / "mri-dataset/mri_dataset/sub-01" / "ses-01/mixed/sub-01_ses-01_acq-mixed_SE-modulus.nii.gz" - T1_path = tmp_path / "output_acq-mixed_T1map_raw.nii.gz" - ref_output = (mri_data_dir / "mri-dataset/mri_dataset/derivatives/sub-01" / "ses-01/sub-01_ses-01_acq-mixed_T1map.nii.gz",) - test_output = tmp_path / "output_acq-mixed_T1map.nii.gz" + mixed_t1map(SE_nii_path=SE_path, IR_nii_path=IR_path, meta_path=meta_path, T1_low=T1_low, T1_high=T1_high, output=test_output_raw) + mixed_t1map_postprocessing(SE_nii_path=SE_path, T1_path=test_output_raw, output=test_output) - mixed_t1map_postprocessing(SE_nii_path=SE_path, T1_path=T1_path, output=test_output) - compare_nifti_images(test_output, ref_output, data_tolerance=1e-4) + compare_nifti_images(test_output, ref_output, data_tolerance=1e-12) def test_hybrid_t1map(tmp_path, mri_data_dir: Path): @@ -91,4 +72,4 @@ def test_hybrid_t1map(tmp_path, mri_data_dir: Path): hybrid_t1map( LL_path=LL_path, mixed_path=mixed_path, csf_mask_path=csf_mask_path, threshold=threshold, erode=erode, output=test_output ) - compare_nifti_images(test_output, ref_output, data_tolerance=1e-4) + compare_nifti_images(test_output, ref_output, data_tolerance=1e-12) From e3dcd0f2f08912632f89731b78a71e7c38f8c20f Mon Sep 17 00:00:00 2001 From: "pre-commit-ci[bot]" <66853113+pre-commit-ci[bot]@users.noreply.github.com> Date: Fri, 6 Mar 2026 09:44:40 +0000 Subject: [PATCH 12/24] [pre-commit.ci] auto fixes from pre-commit.com hooks for more information, see https://pre-commit.ci --- src/mritk/t1_maps/t1_maps.py | 8 +------- test/test_mri_t1_maps.py | 11 ++++++----- 2 files changed, 7 insertions(+), 12 deletions(-) diff --git a/src/mritk/t1_maps/t1_maps.py b/src/mritk/t1_maps/t1_maps.py index 111a027..88cced6 100644 --- a/src/mritk/t1_maps/t1_maps.py +++ b/src/mritk/t1_maps/t1_maps.py @@ -73,7 +73,6 @@ def looklocker_t1map_postprocessing( mask: Optional[np.ndarray] = None, output: Path = None, ) -> MRIData: - T1map_mri = load_mri_data(T1map, dtype=np.single) T1map = T1map_mri.data.copy() if mask is None: @@ -113,12 +112,7 @@ def looklocker_t1map_postprocessing( def mixed_t1map( - SE_nii_path: Path, - IR_nii_path: Path, - meta_path: Path, - T1_low: float, - T1_high: float, - output: Path = None + SE_nii_path: Path, IR_nii_path: Path, meta_path: Path, T1_low: float, T1_high: float, output: Path = None ) -> nibabel.nifti1.Nifti1Image: SE = load_mri_data(SE_nii_path, dtype=np.single) IR = load_mri_data(IR_nii_path, dtype=np.single) diff --git a/test/test_mri_t1_maps.py b/test/test_mri_t1_maps.py index c733d1e..d40f934 100644 --- a/test/test_mri_t1_maps.py +++ b/test/test_mri_t1_maps.py @@ -18,6 +18,7 @@ from mritk.t1_maps.utils import compare_nifti_images + @pytest.mark.skip(reason="Takes too long") def test_looklocker_t1map(tmp_path, mri_data_dir: Path): LL_path = mri_data_dir / "mri-dataset/mri_dataset/sub-01" / "ses-01/anat/sub-01_ses-01_acq-looklocker_IRT1.nii.gz" @@ -27,9 +28,7 @@ def test_looklocker_t1map(tmp_path, mri_data_dir: Path): T1_low = 100 T1_high = 6000 - ref_output = ( - mri_data_dir / "mri-dataset/mri_dataset/derivatives/sub-01" / "ses-01/sub-01_ses-01_acq-looklocker_T1map.nii.gz" - ) + ref_output = mri_data_dir / "mri-dataset/mri_dataset/derivatives/sub-01" / "ses-01/sub-01_ses-01_acq-looklocker_T1map.nii.gz" test_output_raw = tmp_path / "output_acq-looklocker_T1map_raw.nii.gz" test_output = tmp_path / "output_acq-looklocker_T1map.nii.gz" @@ -42,7 +41,7 @@ def test_mixed_t1map(tmp_path, mri_data_dir: Path): SE_path = mri_data_dir / "mri-dataset/mri_dataset/sub-01" / "ses-01/mixed/sub-01_ses-01_acq-mixed_SE-modulus.nii.gz" IR_path = mri_data_dir / "mri-dataset/mri_dataset/sub-01" / "ses-01/mixed/sub-01_ses-01_acq-mixed_IR-corrected-real.nii.gz" meta_path = mri_data_dir / "mri-dataset/mri_dataset/sub-01" / "ses-01/mixed/sub-01_ses-01_acq-mixed_meta.json" - + ref_output = mri_data_dir / "mri-dataset/mri_dataset/derivatives/sub-01" / "ses-01/sub-01_ses-01_acq-mixed_T1map.nii.gz" test_output_raw = tmp_path / "output_acq-mixed_T1map_raw.nii.gz" test_output = tmp_path / "output_acq-mixed_T1map.nii.gz" @@ -50,7 +49,9 @@ def test_mixed_t1map(tmp_path, mri_data_dir: Path): T1_low = 100 T1_high = 10000 - mixed_t1map(SE_nii_path=SE_path, IR_nii_path=IR_path, meta_path=meta_path, T1_low=T1_low, T1_high=T1_high, output=test_output_raw) + mixed_t1map( + SE_nii_path=SE_path, IR_nii_path=IR_path, meta_path=meta_path, T1_low=T1_low, T1_high=T1_high, output=test_output_raw + ) mixed_t1map_postprocessing(SE_nii_path=SE_path, T1_path=test_output_raw, output=test_output) compare_nifti_images(test_output, ref_output, data_tolerance=1e-12) From b833bffd37d6edca7f7521070f0cfa8a82dde6ed Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?C=C3=A9cile?= Date: Fri, 6 Mar 2026 13:23:58 +0100 Subject: [PATCH 13/24] Add test for concentration map - passing --- src/mritk/t1_maps/utils.py | 8 ++++++++ test/test_concentration.py | 36 ++++++++++++++++++++++++++++++++++++ 2 files changed, 44 insertions(+) create mode 100644 test/test_concentration.py diff --git a/src/mritk/t1_maps/utils.py b/src/mritk/t1_maps/utils.py index a4767d8..ef71181 100644 --- a/src/mritk/t1_maps/utils.py +++ b/src/mritk/t1_maps/utils.py @@ -188,6 +188,14 @@ def compare_nifti_images(img_path1, img_path2, data_tolerance=0.0): # 1. Compare Image Data data1 = img1.get_fdata() data2 = img2.get_fdata() + + # Convert NaN to zero (can have NaNs in concentration maps) + data1 = np.nan_to_num(data1, nan=0.0) + data2 = np.nan_to_num(data2, nan=0.0) + + print("data 1= ", np.unique(data1)) + print("data 1= ", np.unique(data2)) + # Use np.allclose for data comparison with tolerance, which is often needed # for floating-point data, or np.array_equal for exact comparison. if data_tolerance > 0: diff --git a/test/test_concentration.py b/test/test_concentration.py new file mode 100644 index 0000000..1daee67 --- /dev/null +++ b/test/test_concentration.py @@ -0,0 +1,36 @@ +"""MRI Concentration maps - Tests + +Copyright (C) 2026 Jørgen Riseth (jnriseth@gmail.com) +Copyright (C) 2026 Cécile Daversin-Catty (cecile@simula.no) +Copyright (C) 2026 Simula Research Laboratory +""" + +from pathlib import Path +import pytest + +from mritk.concentration.concentration import ( + concentration +) + +from mritk.t1_maps.utils import compare_nifti_images + +def test_intracranial_concentration(tmp_path, mri_data_dir: Path): + baseline_path = mri_data_dir / "mri-processed/mri_processed_data/sub-01/T1maps/sub-01_ses-01_T1map_hybrid.nii.gz" + sessions = range(2, 5) + + images_path = [ + mri_data_dir / f"mri-processed/mri_processed_data/sub-01/T1maps/sub-01_ses-0{i}_T1map_hybrid.nii.gz" for i in sessions + ] + mask_path = mri_data_dir / "mri-processed/mri_processed_data/sub-01/segmentations/sub-01_seg-intracranial_binary.nii.gz" + r1 = 0.0032 + + ref_outputs = [ + mri_data_dir / f"mri-processed/mri_processed_data/sub-01/concentrations/sub-01_ses-0{i}_concentration.nii.gz" for i in sessions + ] + test_outputs = [ + tmp_path / f"output_ses-0{i}_concentration.nii.gz" for i in sessions + ] + + for i, s in enumerate(sessions): + concentration(input=images_path[i], reference=baseline_path, output=test_outputs[i], r1=r1, mask=mask_path) + compare_nifti_images(test_outputs[i], ref_outputs[i], data_tolerance=1e-4) \ No newline at end of file From 07291a500a45a6cc2cd34e6f5a8f5dc6575d3e4e Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?C=C3=A9cile?= Date: Fri, 6 Mar 2026 13:32:32 +0100 Subject: [PATCH 14/24] Minor - Change tolerance for test_concentration to 1e-12 --- test/test_concentration.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/test/test_concentration.py b/test/test_concentration.py index 1daee67..3070660 100644 --- a/test/test_concentration.py +++ b/test/test_concentration.py @@ -33,4 +33,4 @@ def test_intracranial_concentration(tmp_path, mri_data_dir: Path): for i, s in enumerate(sessions): concentration(input=images_path[i], reference=baseline_path, output=test_outputs[i], r1=r1, mask=mask_path) - compare_nifti_images(test_outputs[i], ref_outputs[i], data_tolerance=1e-4) \ No newline at end of file + compare_nifti_images(test_outputs[i], ref_outputs[i], data_tolerance=1e-12) \ No newline at end of file From b8254e7c06244828c9232748373eca01a1f7326a Mon Sep 17 00:00:00 2001 From: "pre-commit-ci[bot]" <66853113+pre-commit-ci[bot]@users.noreply.github.com> Date: Fri, 6 Mar 2026 12:36:10 +0000 Subject: [PATCH 15/24] [pre-commit.ci] auto fixes from pre-commit.com hooks for more information, see https://pre-commit.ci --- src/mritk/t1_maps/utils.py | 2 +- test/test_concentration.py | 17 +++++++---------- 2 files changed, 8 insertions(+), 11 deletions(-) diff --git a/src/mritk/t1_maps/utils.py b/src/mritk/t1_maps/utils.py index ef71181..d2a3258 100644 --- a/src/mritk/t1_maps/utils.py +++ b/src/mritk/t1_maps/utils.py @@ -190,7 +190,7 @@ def compare_nifti_images(img_path1, img_path2, data_tolerance=0.0): data2 = img2.get_fdata() # Convert NaN to zero (can have NaNs in concentration maps) - data1 = np.nan_to_num(data1, nan=0.0) + data1 = np.nan_to_num(data1, nan=0.0) data2 = np.nan_to_num(data2, nan=0.0) print("data 1= ", np.unique(data1)) diff --git a/test/test_concentration.py b/test/test_concentration.py index 3070660..d0d7b39 100644 --- a/test/test_concentration.py +++ b/test/test_concentration.py @@ -6,14 +6,12 @@ """ from pathlib import Path -import pytest - -from mritk.concentration.concentration import ( - concentration -) + +from mritk.concentration.concentration import concentration from mritk.t1_maps.utils import compare_nifti_images + def test_intracranial_concentration(tmp_path, mri_data_dir: Path): baseline_path = mri_data_dir / "mri-processed/mri_processed_data/sub-01/T1maps/sub-01_ses-01_T1map_hybrid.nii.gz" sessions = range(2, 5) @@ -25,12 +23,11 @@ def test_intracranial_concentration(tmp_path, mri_data_dir: Path): r1 = 0.0032 ref_outputs = [ - mri_data_dir / f"mri-processed/mri_processed_data/sub-01/concentrations/sub-01_ses-0{i}_concentration.nii.gz" for i in sessions - ] - test_outputs = [ - tmp_path / f"output_ses-0{i}_concentration.nii.gz" for i in sessions + mri_data_dir / f"mri-processed/mri_processed_data/sub-01/concentrations/sub-01_ses-0{i}_concentration.nii.gz" + for i in sessions ] + test_outputs = [tmp_path / f"output_ses-0{i}_concentration.nii.gz" for i in sessions] for i, s in enumerate(sessions): concentration(input=images_path[i], reference=baseline_path, output=test_outputs[i], r1=r1, mask=mask_path) - compare_nifti_images(test_outputs[i], ref_outputs[i], data_tolerance=1e-12) \ No newline at end of file + compare_nifti_images(test_outputs[i], ref_outputs[i], data_tolerance=1e-12) From fab2c2dd1028574ea53cbbac9277bb6dcfaa9c22 Mon Sep 17 00:00:00 2001 From: Henrik Finsberg Date: Fri, 6 Mar 2026 14:13:24 +0100 Subject: [PATCH 16/24] Add new test data --- .github/workflows/main.yml | 2 +- .github/workflows/setup-data.yml | 2 +- src/mritk/datasets.py | 12 ++++++------ test/test_concentration.py | 19 ++++++++----------- 4 files changed, 16 insertions(+), 19 deletions(-) diff --git a/.github/workflows/main.yml b/.github/workflows/main.yml index b8c0068..47e4d06 100644 --- a/.github/workflows/main.yml +++ b/.github/workflows/main.yml @@ -12,7 +12,7 @@ jobs: build: needs: prepare-data env: - MRITK_TEST_DATA_FOLDER: ./test_data + MRITK_TEST_DATA_FOLDER: ./test_data/mritk-test-data/mritk-test-data name: Test on ${{ matrix.os }} with Python ${{ matrix.python-version }} runs-on: ${{ matrix.os }} diff --git a/.github/workflows/setup-data.yml b/.github/workflows/setup-data.yml index 98869b9..51defe4 100644 --- a/.github/workflows/setup-data.yml +++ b/.github/workflows/setup-data.yml @@ -22,7 +22,7 @@ jobs: path: data/ # The folder you want to cache # The key determines if we have a match. # Change 'v1' to 'v2' manually to force a re-download in the future. - key: test-data-v2 + key: test-data-v3 # 2. DOWNLOAD ONLY IF CACHE MISS diff --git a/src/mritk/datasets.py b/src/mritk/datasets.py index 1ccd553..1fadd84 100644 --- a/src/mritk/datasets.py +++ b/src/mritk/datasets.py @@ -15,6 +15,11 @@ logger = logging.getLogger(__name__) +def download_link_google_drive(file_id: str) -> str: + # https://gist.github.com/tanaikech/f0f2d122e05bf5f971611258c22c110f + return f"https://drive.usercontent.google.com/download?id={file_id}&confirm=xxx" + + @dataclass class Dataset: name: str @@ -29,13 +34,8 @@ def get_datasets() -> dict[str, Dataset]: "test-data": Dataset( name="Test Data", description="A small test dataset for testing functionality (based on the Gonzo dataset).", - doi="10.5281/zenodo.14266867", license="CC-BY-4.0", - links={ - "mri-processed.zip": "https://zenodo.org/records/14266867/files/mri-processed.zip?download=1", - "mri-dataset.zip": "https://zenodo.org/records/14266867/files/mri-dataset.zip?download=1", - "timetable.tsv": "https://github.com/jorgenriseth/gonzo/blob/main/mri_dataset/timetable.tsv?raw=true", - }, + links={"mritk-test-data.zip": download_link_google_drive("1Rx-DnT3fBq-3S-0DXgr0ajkANYk8v_dZ")}, ), "gonzo": Dataset( name="The Gonzo Dataset", diff --git a/test/test_concentration.py b/test/test_concentration.py index 1daee67..5258415 100644 --- a/test/test_concentration.py +++ b/test/test_concentration.py @@ -6,17 +6,15 @@ """ from pathlib import Path -import pytest - -from mritk.concentration.concentration import ( - concentration -) + +from mritk.concentration.concentration import concentration from mritk.t1_maps.utils import compare_nifti_images + def test_intracranial_concentration(tmp_path, mri_data_dir: Path): baseline_path = mri_data_dir / "mri-processed/mri_processed_data/sub-01/T1maps/sub-01_ses-01_T1map_hybrid.nii.gz" - sessions = range(2, 5) + sessions = [1, 2] images_path = [ mri_data_dir / f"mri-processed/mri_processed_data/sub-01/T1maps/sub-01_ses-0{i}_T1map_hybrid.nii.gz" for i in sessions @@ -25,12 +23,11 @@ def test_intracranial_concentration(tmp_path, mri_data_dir: Path): r1 = 0.0032 ref_outputs = [ - mri_data_dir / f"mri-processed/mri_processed_data/sub-01/concentrations/sub-01_ses-0{i}_concentration.nii.gz" for i in sessions - ] - test_outputs = [ - tmp_path / f"output_ses-0{i}_concentration.nii.gz" for i in sessions + mri_data_dir / f"mri-processed/mri_processed_data/sub-01/concentrations/sub-01_ses-0{i}_concentration.nii.gz" + for i in sessions ] + test_outputs = [tmp_path / f"output_ses-0{i}_concentration.nii.gz" for i in sessions] for i, s in enumerate(sessions): concentration(input=images_path[i], reference=baseline_path, output=test_outputs[i], r1=r1, mask=mask_path) - compare_nifti_images(test_outputs[i], ref_outputs[i], data_tolerance=1e-4) \ No newline at end of file + compare_nifti_images(test_outputs[i], ref_outputs[i], data_tolerance=1e-4) From 6eaeb73cae20b904d4a1febbb73d76329876ccc2 Mon Sep 17 00:00:00 2001 From: Henrik Finsberg Date: Fri, 6 Mar 2026 14:14:18 +0100 Subject: [PATCH 17/24] Add new test data --- src/mritk/t1_maps/utils.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/mritk/t1_maps/utils.py b/src/mritk/t1_maps/utils.py index ef71181..d2a3258 100644 --- a/src/mritk/t1_maps/utils.py +++ b/src/mritk/t1_maps/utils.py @@ -190,7 +190,7 @@ def compare_nifti_images(img_path1, img_path2, data_tolerance=0.0): data2 = img2.get_fdata() # Convert NaN to zero (can have NaNs in concentration maps) - data1 = np.nan_to_num(data1, nan=0.0) + data1 = np.nan_to_num(data1, nan=0.0) data2 = np.nan_to_num(data2, nan=0.0) print("data 1= ", np.unique(data1)) From ac0c02921b536d94d6bb6f089c1ffafb88e0e63e Mon Sep 17 00:00:00 2001 From: Henrik Finsberg Date: Fri, 6 Mar 2026 14:17:48 +0100 Subject: [PATCH 18/24] Fix typing issues --- src/mritk/masking/masks.py | 2 +- src/mritk/t1_maps/t1_maps.py | 34 +++++++++++++++++----------------- src/mritk/t1_maps/t1_to_r1.py | 2 +- 3 files changed, 19 insertions(+), 19 deletions(-) diff --git a/src/mritk/masking/masks.py b/src/mritk/masking/masks.py index 15df91f..1ead7ae 100644 --- a/src/mritk/masking/masks.py +++ b/src/mritk/masking/masks.py @@ -38,7 +38,7 @@ def csf_mask( input: Path, connectivity: Optional[int] = 2, use_li: bool = False, - output: Path = None, + output: Path | None = None, ) -> MRIData: input_vol = load_mri_data(input, dtype=np.single) mask = create_csf_mask(input_vol.data, connectivity, use_li) diff --git a/src/mritk/t1_maps/t1_maps.py b/src/mritk/t1_maps/t1_maps.py index 88cced6..a6e4184 100644 --- a/src/mritk/t1_maps/t1_maps.py +++ b/src/mritk/t1_maps/t1_maps.py @@ -26,7 +26,7 @@ ) -def looklocker_t1map(looklocker_input: Path, timestamps: Path, output: Path = None) -> MRIData: +def looklocker_t1map(looklocker_input: Path, timestamps: Path, output: Path | None = None) -> MRIData: LL_mri = load_mri_data(looklocker_input, dtype=np.single) D = LL_mri.data affine = LL_mri.affine @@ -47,7 +47,7 @@ def looklocker_t1map(looklocker_input: Path, timestamps: Path, output: Path = No vfunc = np.vectorize(voxel_fitter, signature="(n) -> (3)") fitted_coefficients = vfunc(Dmasked) - x1, x2, x3 = ( + _, x2, x3 = ( fitted_coefficients[:, 0], fitted_coefficients[:, 1], fitted_coefficients[:, 2], @@ -71,13 +71,13 @@ def looklocker_t1map_postprocessing( radius: int = 10, erode_dilate_factor: float = 1.3, mask: Optional[np.ndarray] = None, - output: Path = None, + output: Path | None = None, ) -> MRIData: T1map_mri = load_mri_data(T1map, dtype=np.single) - T1map = T1map_mri.data.copy() + T1map_data = T1map_mri.data.copy() if mask is None: # Create mask for largest island. - mask = skimage.measure.label(np.isfinite(T1map)) + mask = skimage.measure.label(np.isfinite(T1map_data)) regions = skimage.measure.regionprops(mask) regions.sort(key=lambda x: x.num_pixels, reverse=True) mask = mask == regions[0].label @@ -86,25 +86,25 @@ def looklocker_t1map_postprocessing( skimage.morphology.erosion(mask, skimage.morphology.ball(erode_dilate_factor * radius), out=mask) # Remove non-zero artifacts outside of the mask. - surface_vox = np.isfinite(T1map) * (~mask) + surface_vox = np.isfinite(T1map_data) * (~mask) print(f"Removing {surface_vox.sum()} voxels outside of the head mask") - T1map[~mask] = np.nan + T1map_data[~mask] = np.nan # Remove outliers within the mask. - outliers = np.logical_or(T1map < T1_low, T1_high < T1map) + outliers = np.logical_or(T1map_data < T1_low, T1_high < T1map_data) print("Removing", outliers.sum(), f"voxels outside the range ({T1_low}, {T1_high}).") - T1map[outliers] = np.nan - if np.isfinite(T1map).sum() / T1map.size < 0.01: + T1map_data[outliers] = np.nan + if np.isfinite(T1map_data).sum() / T1map_data.size < 0.01: raise RuntimeError("After outlier removal, less than 1% of the image is left. Check image units.") # Fill internallly missing values - fill_mask = np.isnan(T1map) * mask + fill_mask = np.isnan(T1map_data) * mask while fill_mask.sum() > 0: print(f"Filling in {fill_mask.sum()} voxels within the mask.") - T1map[fill_mask] = nan_filter_gaussian(T1map, 1.0)[fill_mask] - fill_mask = np.isnan(T1map) * mask + T1map_data[fill_mask] = nan_filter_gaussian(T1map_data, 1.0)[fill_mask] + fill_mask = np.isnan(T1map_data) * mask - processed_T1map = MRIData(T1map, T1map_mri.affine) + processed_T1map = MRIData(T1map_data, T1map_mri.affine) if output is not None: save_mri_data(processed_T1map, output, dtype=np.single) @@ -112,7 +112,7 @@ def looklocker_t1map_postprocessing( def mixed_t1map( - SE_nii_path: Path, IR_nii_path: Path, meta_path: Path, T1_low: float, T1_high: float, output: Path = None + SE_nii_path: Path, IR_nii_path: Path, meta_path: Path, T1_low: float, T1_high: float, output: Path | None = None ) -> nibabel.nifti1.Nifti1Image: SE = load_mri_data(SE_nii_path, dtype=np.single) IR = load_mri_data(IR_nii_path, dtype=np.single) @@ -137,7 +137,7 @@ def mixed_t1map( return nii -def mixed_t1map_postprocessing(SE_nii_path: Path, T1_path: Path, output: Path = None) -> nibabel.nifti1.Nifti1Image: +def mixed_t1map_postprocessing(SE_nii_path: Path, T1_path: Path, output: Path | None = None) -> nibabel.nifti1.Nifti1Image: T1map_nii = nibabel.nifti1.load(T1_path) SE_mri = load_mri_data(SE_nii_path, np.single) @@ -155,7 +155,7 @@ def mixed_t1map_postprocessing(SE_nii_path: Path, T1_path: Path, output: Path = def hybrid_t1map( - LL_path: Path, mixed_path: Path, csf_mask_path: Path, threshold: float, erode: int = 0, output: Path = None + LL_path: Path, mixed_path: Path, csf_mask_path: Path, threshold: float, erode: int = 0, output: Path | None = None ) -> nibabel.nifti1.Nifti1Image: mixed_mri = nibabel.nifti1.load(mixed_path) mixed = mixed_mri.get_fdata() diff --git a/src/mritk/t1_maps/t1_to_r1.py b/src/mritk/t1_maps/t1_to_r1.py index a783928..af4290d 100644 --- a/src/mritk/t1_maps/t1_to_r1.py +++ b/src/mritk/t1_maps/t1_to_r1.py @@ -29,7 +29,7 @@ def convert_T1_to_R1( def T1_to_R1( input_mri: Union[Path, MRIData], - output: Path = None, + output: Path | None = None, scale: float = 1000, t1_low: float = 1, t1_high: float = float("inf"), From f9e900c89f6592b2eaad054f9bc87acfeb87a08e Mon Sep 17 00:00:00 2001 From: Henrik Finsberg Date: Fri, 6 Mar 2026 14:24:00 +0100 Subject: [PATCH 19/24] Update test for datasets --- test/test_datasets.py | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/test/test_datasets.py b/test/test_datasets.py index 16bcbe3..e4493c7 100644 --- a/test/test_datasets.py +++ b/test/test_datasets.py @@ -31,9 +31,9 @@ def mock_datasets(): def test_get_datasets_structure(): """Ensure get_datasets returns a dict of Dataset objects.""" datasets = mritk.datasets.get_datasets() - assert "test-data" in datasets - assert isinstance(datasets["test-data"], Dataset) - assert datasets["test-data"].doi == "10.5281/zenodo.14266867" + assert "gonzo" in datasets + assert isinstance(datasets["gonzo"], Dataset) + assert datasets["gonzo"].doi == "10.5281/zenodo.14266867" def test_add_arguments(): From 6f151082d15f5f553d1cc535b686cdd75c41c513 Mon Sep 17 00:00:00 2001 From: Henrik Finsberg Date: Fri, 6 Mar 2026 14:26:20 +0100 Subject: [PATCH 20/24] List files in test data folder --- .github/workflows/main.yml | 7 ++++++- 1 file changed, 6 insertions(+), 1 deletion(-) diff --git a/.github/workflows/main.yml b/.github/workflows/main.yml index 47e4d06..ff03e5e 100644 --- a/.github/workflows/main.yml +++ b/.github/workflows/main.yml @@ -12,7 +12,7 @@ jobs: build: needs: prepare-data env: - MRITK_TEST_DATA_FOLDER: ./test_data/mritk-test-data/mritk-test-data + MRITK_TEST_DATA_FOLDER: ./test_data/mritk-test-data name: Test on ${{ matrix.os }} with Python ${{ matrix.python-version }} runs-on: ${{ matrix.os }} @@ -33,6 +33,11 @@ jobs: name: shared-test-data path: ${{ env.MRITK_TEST_DATA_FOLDER }} # Where you want the data to appear in this runner + - name: List downloaded files + run: | + echo "Files in ${{ env.MRITK_TEST_DATA_FOLDER }}:" + ls -R ${{ env.MRITK_TEST_DATA_FOLDER }} + - name: Set up Python ${{ matrix.python-version }} uses: actions/setup-python@v6 with: From 6209bf284d8878417f48bd26b19bd59352b84d52 Mon Sep 17 00:00:00 2001 From: Henrik Finsberg Date: Fri, 6 Mar 2026 14:27:46 +0100 Subject: [PATCH 21/24] Fix path to test data --- .github/workflows/main.yml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/.github/workflows/main.yml b/.github/workflows/main.yml index ff03e5e..2a593d4 100644 --- a/.github/workflows/main.yml +++ b/.github/workflows/main.yml @@ -12,7 +12,7 @@ jobs: build: needs: prepare-data env: - MRITK_TEST_DATA_FOLDER: ./test_data/mritk-test-data + MRITK_TEST_DATA_FOLDER: ./test_data/mritk-test-data/mritk-test-data/mritk-test-data/ name: Test on ${{ matrix.os }} with Python ${{ matrix.python-version }} runs-on: ${{ matrix.os }} From 1c1e103ec24fc54af029c6d51ee34806f8b00d17 Mon Sep 17 00:00:00 2001 From: Henrik Finsberg Date: Fri, 6 Mar 2026 14:34:05 +0100 Subject: [PATCH 22/24] Add new linkt to data --- .github/workflows/main.yml | 2 +- .github/workflows/setup-data.yml | 2 +- src/mritk/datasets.py | 2 +- 3 files changed, 3 insertions(+), 3 deletions(-) diff --git a/.github/workflows/main.yml b/.github/workflows/main.yml index 2a593d4..875ba03 100644 --- a/.github/workflows/main.yml +++ b/.github/workflows/main.yml @@ -12,7 +12,7 @@ jobs: build: needs: prepare-data env: - MRITK_TEST_DATA_FOLDER: ./test_data/mritk-test-data/mritk-test-data/mritk-test-data/ + MRITK_TEST_DATA_FOLDER: ./test_data/mritk-test-data/mritk-test-data name: Test on ${{ matrix.os }} with Python ${{ matrix.python-version }} runs-on: ${{ matrix.os }} diff --git a/.github/workflows/setup-data.yml b/.github/workflows/setup-data.yml index 51defe4..33dcf74 100644 --- a/.github/workflows/setup-data.yml +++ b/.github/workflows/setup-data.yml @@ -22,7 +22,7 @@ jobs: path: data/ # The folder you want to cache # The key determines if we have a match. # Change 'v1' to 'v2' manually to force a re-download in the future. - key: test-data-v3 + key: test-data-v4 # 2. DOWNLOAD ONLY IF CACHE MISS diff --git a/src/mritk/datasets.py b/src/mritk/datasets.py index 1fadd84..bf7a23e 100644 --- a/src/mritk/datasets.py +++ b/src/mritk/datasets.py @@ -35,7 +35,7 @@ def get_datasets() -> dict[str, Dataset]: name="Test Data", description="A small test dataset for testing functionality (based on the Gonzo dataset).", license="CC-BY-4.0", - links={"mritk-test-data.zip": download_link_google_drive("1Rx-DnT3fBq-3S-0DXgr0ajkANYk8v_dZ")}, + links={"mritk-test-data.zip": download_link_google_drive("1CSj3CHd4ztcU4Aqdlw9K2OWjPi5u75bd")}, ), "gonzo": Dataset( name="The Gonzo Dataset", From 3053e9cd3d9985b69ed4588f761aa77cad0d8fc4 Mon Sep 17 00:00:00 2001 From: Henrik Finsberg Date: Fri, 6 Mar 2026 14:38:20 +0100 Subject: [PATCH 23/24] Try different test folder --- .github/workflows/main.yml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/.github/workflows/main.yml b/.github/workflows/main.yml index 875ba03..1545d40 100644 --- a/.github/workflows/main.yml +++ b/.github/workflows/main.yml @@ -12,7 +12,7 @@ jobs: build: needs: prepare-data env: - MRITK_TEST_DATA_FOLDER: ./test_data/mritk-test-data/mritk-test-data + MRITK_TEST_DATA_FOLDER: ./test_data/mritk-test-data/mritk-test-data/mritk-test-data name: Test on ${{ matrix.os }} with Python ${{ matrix.python-version }} runs-on: ${{ matrix.os }} From d559c8cac017b9cd23661b743909082f4c331fde Mon Sep 17 00:00:00 2001 From: Henrik Finsberg Date: Fri, 6 Mar 2026 14:45:57 +0100 Subject: [PATCH 24/24] Separate between test folder and root test folder --- .github/workflows/main.yml | 9 +++++---- 1 file changed, 5 insertions(+), 4 deletions(-) diff --git a/.github/workflows/main.yml b/.github/workflows/main.yml index 1545d40..6b84d91 100644 --- a/.github/workflows/main.yml +++ b/.github/workflows/main.yml @@ -12,7 +12,8 @@ jobs: build: needs: prepare-data env: - MRITK_TEST_DATA_FOLDER: ./test_data/mritk-test-data/mritk-test-data/mritk-test-data + MRITK_TEST_DATA_FOLDER_ROOT: ./test_data + MRITK_TEST_DATA_FOLDER: ./test_data/mritk-test-data name: Test on ${{ matrix.os }} with Python ${{ matrix.python-version }} runs-on: ${{ matrix.os }} @@ -31,12 +32,12 @@ jobs: uses: actions/download-artifact@v8 with: name: shared-test-data - path: ${{ env.MRITK_TEST_DATA_FOLDER }} # Where you want the data to appear in this runner + path: ${{ env.MRITK_TEST_DATA_FOLDER_ROOT }} # Where you want the data to appear in this runner - name: List downloaded files run: | - echo "Files in ${{ env.MRITK_TEST_DATA_FOLDER }}:" - ls -R ${{ env.MRITK_TEST_DATA_FOLDER }} + echo "Files in ${{ env.MRITK_TEST_DATA_FOLDER_ROOT }}:" + ls -R ${{ env.MRITK_TEST_DATA_FOLDER_ROOT }} - name: Set up Python ${{ matrix.python-version }} uses: actions/setup-python@v6