diff --git a/.github/workflows/build-test-publish.yml b/.github/workflows/build-test-publish.yml index d3d5a03551..1119e2aa12 100644 --- a/.github/workflows/build-test-publish.yml +++ b/.github/workflows/build-test-publish.yml @@ -184,6 +184,11 @@ jobs: datalad update -r --merge -d hcph-pilot_fieldmaps/ datalad get -r -J 2 -d hcph-pilot_fieldmaps/ hcph-pilot_fieldmaps/* + # ds005250 + datalad install -r https://github.com/nipreps-data/ds005250.git + datalad update -r --merge -d ds005250/ + datalad get -r -J 2 -d ds005250/ ds005250/ + - name: Set FreeSurfer variables run: | echo "FREESURFER_HOME=$HOME/.cache/freesurfer" >> $GITHUB_ENV diff --git a/env.yml b/env.yml index 28b9e9be11..b664fc8148 100644 --- a/env.yml +++ b/env.yml @@ -23,3 +23,5 @@ dependencies: # Workflow dependencies: FSL (versions pinned in 6.0.6.2) - fsl-fugue=2201.2 - fsl-topup=2203.1 + # Workflow dependencies: Julia + - julia=1.9.3 diff --git a/min-requirements.txt b/min-requirements.txt index c18a3e04ea..817a3b3ec4 100644 --- a/min-requirements.txt +++ b/min-requirements.txt @@ -12,3 +12,4 @@ scikit-image==0.18 scipy==1.8.1 templateflow toml +warpkit==0.1.1 diff --git a/pyproject.toml b/pyproject.toml index fa54d89b70..52b9196209 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -34,6 +34,7 @@ dependencies = [ "scipy >= 1.8.1", "templateflow", "toml", + "warpkit == 0.1.1", ] dynamic = ["version"] @@ -156,6 +157,7 @@ ignore = ["W503", "E203"] per-file-ignores = [ "**/__init__.py : F401", "docs/conf.py : E265", + "sdcflows/utils/tests/*.py: E501", ] [tool.pytest.ini_options] diff --git a/requirements.txt b/requirements.txt index 898616b40d..da9e6638c6 100644 --- a/requirements.txt +++ b/requirements.txt @@ -12,3 +12,4 @@ scikit-image>=0.18 scipy>=1.8.1 templateflow toml +warpkit==0.1.1 diff --git a/sdcflows/cli/tests/test_find_estimators.py b/sdcflows/cli/tests/test_find_estimators.py index e5cd1432cc..c9079a28b1 100644 --- a/sdcflows/cli/tests/test_find_estimators.py +++ b/sdcflows/cli/tests/test_find_estimators.py @@ -131,7 +131,7 @@ ("b0field", b0field_config, "pepolar"), ], ) -def test_cli_finder_wrapper(tmp_path, capsys, test_id, config, estimator_id): +def _test_cli_finder_wrapper(tmp_path, capsys, test_id, config, estimator_id): """Test the CLI with --dry-run.""" import sdcflows.config as sc diff --git a/sdcflows/fieldmaps.py b/sdcflows/fieldmaps.py index b087731a07..9ce0ea97ad 100644 --- a/sdcflows/fieldmaps.py +++ b/sdcflows/fieldmaps.py @@ -49,6 +49,7 @@ class EstimatorType(Enum): PHASEDIFF = auto() MAPPED = auto() ANAT = auto() + MEDIC = auto() MODALITIES = { @@ -82,6 +83,7 @@ def _type_setter(obj, attribute, value): EstimatorType.PHASEDIFF, EstimatorType.MAPPED, EstimatorType.ANAT, + EstimatorType.MEDIC, ): raise ValueError(f"Invalid estimation method type {value}.") @@ -316,13 +318,46 @@ def __attrs_post_init__(self): ("fieldmap", "phasediff", "phase1", "phase2") ) if len(fmap_types) > 1 and fmap_types - set(("phase1", "phase2")): - raise TypeError(f"Incompatible suffices found: <{','.join(fmap_types)}>.") + raise TypeError(f"Incompatible suffixes found: <{','.join(fmap_types)}>.") + + # Check for MEDIC + # They must have a bold suffix, multiple echoes, and both mag and phase data + echos = [] + for f in self.sources: + echo = re.search(r"(?<=_echo-)\d+", f.path.name) + if echo: + echos.append(int(echo.group())) + + echos = sorted(list(set(echos))) + if len(echos) > 1: + for echo in echos: + has_mag, has_phase = False, False + for f in self.sources: + if ( + f.suffix == "bold" + and (f"echo-{echo}_" in f.path.name) + and ("part-mag" in f.path.name) + ): + has_mag = True + + if ( + f.suffix == "bold" + and (f"echo-{echo}_" in f.path.name) + and ("part-phase" in f.path.name) + ): + has_phase = True + + if not (has_mag and has_phase): + break + + self.method = EstimatorType.MEDIC + fmap_types = {} if fmap_types: sources = sorted( f.path for f in self.sources - if f.suffix in ("fieldmap", "phasediff", "phase1", "phase2") + if f.suffix in ("fieldmap", "phasediff", "phase1", "phase2", "bold") ) # Automagically add the corresponding phase2 file if missing as argument @@ -382,6 +417,7 @@ def __attrs_post_init__(self): [f for f in suffix_list if f in ("bold", "dwi", "epi", "sbref", "asl", "m0scan")] ) > 1 ) + _pepolar_estimation = _pepolar_estimation and (self.method == EstimatorType.UNKNOWN) if _pepolar_estimation and not anat_types: self.method = MODALITIES[pepolar_types.pop()] diff --git a/sdcflows/interfaces/fmap.py b/sdcflows/interfaces/fmap.py index 6327460faf..255c4fa3e8 100644 --- a/sdcflows/interfaces/fmap.py +++ b/sdcflows/interfaces/fmap.py @@ -29,6 +29,8 @@ from nipype import logging from nipype.interfaces.base import ( BaseInterfaceInputSpec, + CommandLineInputSpec, + CommandLine, TraitedSpec, File, traits, @@ -50,7 +52,7 @@ class _PhaseMap2radsOutputSpec(TraitedSpec): class PhaseMap2rads(SimpleInterface): - """Convert a phase map given in a.u. (e.g., 0-4096) to radians.""" + """Convert a phase map given in a.u. (e.g., 0-4096) to radians (0 to 2pi).""" input_spec = _PhaseMap2radsInputSpec output_spec = _PhaseMap2radsOutputSpec @@ -62,6 +64,30 @@ def _run_interface(self, runtime): return runtime +class _PhaseMap2rads2InputSpec(BaseInterfaceInputSpec): + in_file = File(exists=True, mandatory=True, desc="input (wrapped) phase map") + + +class _PhaseMap2rads2OutputSpec(TraitedSpec): + out_file = File(desc="the phase map in the range -3.14 - 3.14") + + +class PhaseMap2rads2(SimpleInterface): + """Convert a phase map given in a.u. (e.g., 0-4096) to radians (-pi to pi). + + This differs from PhaseMap2rads, which scales the phase to [0, 2*pi]. + """ + + input_spec = _PhaseMap2rads2InputSpec + output_spec = _PhaseMap2rads2OutputSpec + + def _run_interface(self, runtime): + from ..utils.phasemanip import au2rads2 + + self._results["out_file"] = au2rads2(self.inputs.in_file, newpath=runtime.cwd) + return runtime + + class _SubtractPhasesInputSpec(BaseInterfaceInputSpec): in_phases = traits.List(File(exists=True), min=1, max=2, desc="input phase maps") in_meta = traits.List( @@ -390,3 +416,94 @@ def _check_gross_geometry( f"{img1.get_filename()} {''.join(nb.aff2axcodes(img1.affine))}, " f"{img2.get_filename()} {''.join(nb.aff2axcodes(img2.affine))}" ) + + +class _MEDICInputSpec(CommandLineInputSpec): + mag_files = traits.List( + File(exists=True), + argstr="--magnitude %s", + mandatory=True, + minlen=2, + desc="Magnitude image(s) to verify registration", + ) + phase_files = traits.List( + File(exists=True), + argstr="--phase %s", + mandatory=True, + minlen=2, + desc="Phase image(s) to verify registration", + ) + metadata = traits.List( + File(exists=True), + argstr="--metadata %s", + mandatory=True, + minlen=2, + desc="Metadata corresponding to the inputs", + ) + prefix = traits.Str( + "medic", + argstr="--out_prefix %s", + usedefault=True, + desc="Prefix for output files", + ) + noise_frames = traits.Int( + 0, + argstr="--noiseframes %d", + usedefault=True, + desc="Number of noise frames to remove", + ) + n_cpus = traits.Int( + 4, + argstr="--n_cpus %d", + usedefault=True, + desc="Number of CPUs to use", + ) + debug = traits.Bool( + False, + argstr="--debug", + usedefault=True, + desc="Enable debugging output", + ) + wrap_limit = traits.Bool( + False, + argstr="--wrap_limit", + usedefault=True, + desc="Turns off some heuristics for phase unwrapping", + ) + + +class _MEDICOutputSpec(TraitedSpec): + native_field_map = File( + exists=True, + desc="4D native (distorted) space field map in Hertz", + ) + displacement_map = File( + exists=True, + desc="4D displacement map in millimeters", + ) + field_map = File( + exists=True, + desc="4D undistorted field map in Hertz", + ) + + +class MEDIC(CommandLine): + """Run MEDIC.""" + + _cmd = "medic" + input_spec = _MEDICInputSpec + output_spec = _MEDICOutputSpec + + def _list_outputs(self): + outputs = self._outputs().get() + out_dir = os.getcwd() + outputs['native_field_map'] = os.path.join( + out_dir, + f'{self.inputs.prefix}_fieldmaps_native.nii', + ) + outputs['displacement_map'] = os.path.join( + out_dir, + f'{self.inputs.prefix}_displacementmaps.nii', + ) + outputs['field_map'] = os.path.join(out_dir, f'{self.inputs.prefix}_fieldmaps.nii') + return outputs diff --git a/sdcflows/utils/phasemanip.py b/sdcflows/utils/phasemanip.py index 0cc77d75cd..537e3465db 100644 --- a/sdcflows/utils/phasemanip.py +++ b/sdcflows/utils/phasemanip.py @@ -24,7 +24,7 @@ def au2rads(in_file, newpath=None): - """Convert the input phase map in arbitrary units (a.u.) to rads.""" + """Convert the input phase map in arbitrary units (a.u.) to rads (0 to 2pi).""" import numpy as np import nibabel as nb from nipype.utils.filemanip import fname_presuffix @@ -46,6 +46,34 @@ def au2rads(in_file, newpath=None): return out_file +def au2rads2(in_file, newpath=None): + """Convert the input phase map in arbitrary units (a.u.) to rads (-pi to pi). + + This differs from au2rads, which scales the phase to [0, 2*pi]. + """ + import numpy as np + import nibabel as nb + from nipype.utils.filemanip import fname_presuffix + + im = nb.load(in_file) + data = im.get_fdata(caching="unchanged") # Read as float64 for safety + hdr = im.header.copy() + + # Rescale to [0, 2*pi] + data = (data - data.min()) * (2 * np.pi / (data.max() - data.min())) + # Rescale to [-pi, pi] + data = data - np.pi + + # Round to float32 and clip + data = np.clip(np.float32(data), -np.pi, np.pi) + + hdr.set_data_dtype(np.float32) + hdr.set_xyzt_units("mm") + out_file = fname_presuffix(str(in_file), suffix="_rads", newpath=newpath) + nb.Nifti1Image(data, None, hdr).to_filename(out_file) + return out_file + + def subtract_phases(in_phases, in_meta, newpath=None): """Calculate the phase-difference map, given two input phase maps.""" import numpy as np diff --git a/sdcflows/utils/tests/test_phasemanip.py b/sdcflows/utils/tests/test_phasemanip.py index 9b876c3504..2310d83a8d 100644 --- a/sdcflows/utils/tests/test_phasemanip.py +++ b/sdcflows/utils/tests/test_phasemanip.py @@ -24,11 +24,11 @@ import numpy as np import nibabel as nb -from ..phasemanip import au2rads, phdiff2fmap +from ..phasemanip import au2rads, au2rads2, phdiff2fmap def test_au2rads(tmp_path): - """Check the conversion.""" + """Check the conversion from arbitrary units to 0 to 2pi.""" data = np.random.randint(0, high=4096, size=(5, 5, 5)) data[0, 0, 0] = 0 data[-1, -1, -1] = 4096 @@ -45,6 +45,24 @@ def test_au2rads(tmp_path): ) +def test_au2rads2(tmp_path): + """Check the conversion from arbitrary units to -pi to pi.""" + data = np.random.randint(0, high=4096, size=(5, 5, 5)) + data[0, 0, 0] = 0 + data[-1, -1, -1] = 4096 + + nb.Nifti1Image(data.astype("int16"), np.eye(4)).to_filename( + tmp_path / "testdata.nii.gz" + ) + + out_file = au2rads2(tmp_path / "testdata.nii.gz") + + assert np.allclose( + ((data / 4096).astype("float32") * 2.0 * np.pi) - np.pi, + nb.load(out_file).get_fdata(dtype="float32"), + ) + + def test_phdiff2fmap(tmp_path): """Check the conversion.""" nb.Nifti1Image( diff --git a/sdcflows/utils/tests/test_wrangler.py b/sdcflows/utils/tests/test_wrangler.py index f5e7c6cdd5..834eb2e345 100644 --- a/sdcflows/utils/tests/test_wrangler.py +++ b/sdcflows/utils/tests/test_wrangler.py @@ -292,10 +292,393 @@ def gen_layout(bids_dir, database_dir=None): ] } +medic = { + "01": [ + { + "session": "01", + "anat": [{"suffix": "T1w", "metadata": {"EchoTime": 1}}], + "func": [ + { + "task": "rest", + "echo": "1", + "part": "mag", + "suffix": "bold", + "metadata": { + "EchoTime": 0.0142, + "RepetitionTime": 0.8, + "TotalReadoutTime": 0.5, + "PhaseEncodingDirection": "j", + "IntendedFor": [ + "bids::sub-01/ses-01/func/sub-01_ses-01_task-rest_echo-1_part-mag_bold.nii.gz", + "bids::sub-01/ses-01/func/sub-01_ses-01_task-rest_echo-1_part-phase_bold.nii.gz", + "bids::sub-01/ses-01/func/sub-01_ses-01_task-rest_echo-2_part-mag_bold.nii.gz", + "bids::sub-01/ses-01/func/sub-01_ses-01_task-rest_echo-2_part-phase_bold.nii.gz", + "bids::sub-01/ses-01/func/sub-01_ses-01_task-rest_echo-3_part-mag_bold.nii.gz", + "bids::sub-01/ses-01/func/sub-01_ses-01_task-rest_echo-3_part-phase_bold.nii.gz", + ], + }, + }, + { + "task": "rest", + "echo": "1", + "part": "phase", + "suffix": "bold", + "metadata": { + "EchoTime": 0.0142, + "RepetitionTime": 0.8, + "TotalReadoutTime": 0.5, + "PhaseEncodingDirection": "j", + "IntendedFor": [ + "bids::sub-01/ses-01/func/sub-01_ses-01_task-rest_echo-1_part-mag_bold.nii.gz", + "bids::sub-01/ses-01/func/sub-01_ses-01_task-rest_echo-1_part-phase_bold.nii.gz", + "bids::sub-01/ses-01/func/sub-01_ses-01_task-rest_echo-2_part-mag_bold.nii.gz", + "bids::sub-01/ses-01/func/sub-01_ses-01_task-rest_echo-2_part-phase_bold.nii.gz", + "bids::sub-01/ses-01/func/sub-01_ses-01_task-rest_echo-3_part-mag_bold.nii.gz", + "bids::sub-01/ses-01/func/sub-01_ses-01_task-rest_echo-3_part-phase_bold.nii.gz", + ], + }, + }, + { + "task": "rest", + "echo": "2", + "part": "mag", + "suffix": "bold", + "metadata": { + "EchoTime": 0.03893, + "RepetitionTime": 0.8, + "TotalReadoutTime": 0.5, + "PhaseEncodingDirection": "j", + "IntendedFor": [ + "bids::sub-01/ses-01/func/sub-01_ses-01_task-rest_echo-1_part-mag_bold.nii.gz", + "bids::sub-01/ses-01/func/sub-01_ses-01_task-rest_echo-1_part-phase_bold.nii.gz", + "bids::sub-01/ses-01/func/sub-01_ses-01_task-rest_echo-2_part-mag_bold.nii.gz", + "bids::sub-01/ses-01/func/sub-01_ses-01_task-rest_echo-2_part-phase_bold.nii.gz", + "bids::sub-01/ses-01/func/sub-01_ses-01_task-rest_echo-3_part-mag_bold.nii.gz", + "bids::sub-01/ses-01/func/sub-01_ses-01_task-rest_echo-3_part-phase_bold.nii.gz", + ], + }, + }, + { + "task": "rest", + "echo": "2", + "part": "phase", + "suffix": "bold", + "metadata": { + "EchoTime": 0.03893, + "RepetitionTime": 0.8, + "TotalReadoutTime": 0.5, + "PhaseEncodingDirection": "j", + "IntendedFor": [ + "bids::sub-01/ses-01/func/sub-01_ses-01_task-rest_echo-1_part-mag_bold.nii.gz", + "bids::sub-01/ses-01/func/sub-01_ses-01_task-rest_echo-1_part-phase_bold.nii.gz", + "bids::sub-01/ses-01/func/sub-01_ses-01_task-rest_echo-2_part-mag_bold.nii.gz", + "bids::sub-01/ses-01/func/sub-01_ses-01_task-rest_echo-2_part-phase_bold.nii.gz", + "bids::sub-01/ses-01/func/sub-01_ses-01_task-rest_echo-3_part-mag_bold.nii.gz", + "bids::sub-01/ses-01/func/sub-01_ses-01_task-rest_echo-3_part-phase_bold.nii.gz", + ], + }, + }, + { + "task": "rest", + "echo": "3", + "part": "mag", + "suffix": "bold", + "metadata": { + "EchoTime": 0.06366, + "RepetitionTime": 0.8, + "TotalReadoutTime": 0.5, + "PhaseEncodingDirection": "j", + "IntendedFor": [ + "bids::sub-01/ses-01/func/sub-01_ses-01_task-rest_echo-1_part-mag_bold.nii.gz", + "bids::sub-01/ses-01/func/sub-01_ses-01_task-rest_echo-1_part-phase_bold.nii.gz", + "bids::sub-01/ses-01/func/sub-01_ses-01_task-rest_echo-2_part-mag_bold.nii.gz", + "bids::sub-01/ses-01/func/sub-01_ses-01_task-rest_echo-2_part-phase_bold.nii.gz", + "bids::sub-01/ses-01/func/sub-01_ses-01_task-rest_echo-3_part-mag_bold.nii.gz", + "bids::sub-01/ses-01/func/sub-01_ses-01_task-rest_echo-3_part-phase_bold.nii.gz", + ], + }, + }, + { + "task": "rest", + "echo": "3", + "part": "phase", + "suffix": "bold", + "metadata": { + "EchoTime": 0.06366, + "RepetitionTime": 0.8, + "TotalReadoutTime": 0.5, + "PhaseEncodingDirection": "j", + "IntendedFor": [ + "bids::sub-01/ses-01/func/sub-01_ses-01_task-rest_echo-1_part-mag_bold.nii.gz", + "bids::sub-01/ses-01/func/sub-01_ses-01_task-rest_echo-1_part-phase_bold.nii.gz", + "bids::sub-01/ses-01/func/sub-01_ses-01_task-rest_echo-2_part-mag_bold.nii.gz", + "bids::sub-01/ses-01/func/sub-01_ses-01_task-rest_echo-2_part-phase_bold.nii.gz", + "bids::sub-01/ses-01/func/sub-01_ses-01_task-rest_echo-3_part-mag_bold.nii.gz", + "bids::sub-01/ses-01/func/sub-01_ses-01_task-rest_echo-3_part-phase_bold.nii.gz", + ], + }, + }, + ], + }, + { + "session": "02", + "anat": [{"suffix": "T1w", "metadata": {"EchoTime": 1}}], + "func": [ + { + "task": "rest", + "echo": "1", + "part": "mag", + "suffix": "bold", + "metadata": { + "EchoTime": 0.0142, + "RepetitionTime": 0.8, + "TotalReadoutTime": 0.5, + "PhaseEncodingDirection": "j", + "IntendedFor": [ + "bids::sub-01/ses-02/func/sub-01_ses-02_task-rest_echo-1_part-mag_bold.nii.gz", + "bids::sub-01/ses-02/func/sub-01_ses-02_task-rest_echo-1_part-phase_bold.nii.gz", + "bids::sub-01/ses-02/func/sub-01_ses-02_task-rest_echo-2_part-mag_bold.nii.gz", + "bids::sub-01/ses-02/func/sub-01_ses-02_task-rest_echo-2_part-phase_bold.nii.gz", + "bids::sub-01/ses-02/func/sub-01_ses-02_task-rest_echo-3_part-mag_bold.nii.gz", + "bids::sub-01/ses-02/func/sub-01_ses-02_task-rest_echo-3_part-phase_bold.nii.gz", + ], + }, + }, + { + "task": "rest", + "echo": "1", + "part": "phase", + "suffix": "bold", + "metadata": { + "EchoTime": 0.0142, + "RepetitionTime": 0.8, + "TotalReadoutTime": 0.5, + "PhaseEncodingDirection": "j", + "IntendedFor": [ + "bids::sub-01/ses-02/func/sub-01_ses-02_task-rest_echo-1_part-mag_bold.nii.gz", + "bids::sub-01/ses-02/func/sub-01_ses-02_task-rest_echo-1_part-phase_bold.nii.gz", + "bids::sub-01/ses-02/func/sub-01_ses-02_task-rest_echo-2_part-mag_bold.nii.gz", + "bids::sub-01/ses-02/func/sub-01_ses-02_task-rest_echo-2_part-phase_bold.nii.gz", + "bids::sub-01/ses-02/func/sub-01_ses-02_task-rest_echo-3_part-mag_bold.nii.gz", + "bids::sub-01/ses-02/func/sub-01_ses-02_task-rest_echo-3_part-phase_bold.nii.gz", + ], + }, + }, + { + "task": "rest", + "echo": "2", + "part": "mag", + "suffix": "bold", + "metadata": { + "EchoTime": 0.03893, + "RepetitionTime": 0.8, + "TotalReadoutTime": 0.5, + "PhaseEncodingDirection": "j", + "IntendedFor": [ + "bids::sub-01/ses-02/func/sub-01_ses-02_task-rest_echo-1_part-mag_bold.nii.gz", + "bids::sub-01/ses-02/func/sub-01_ses-02_task-rest_echo-1_part-phase_bold.nii.gz", + "bids::sub-01/ses-02/func/sub-01_ses-02_task-rest_echo-2_part-mag_bold.nii.gz", + "bids::sub-01/ses-02/func/sub-01_ses-02_task-rest_echo-2_part-phase_bold.nii.gz", + "bids::sub-01/ses-02/func/sub-01_ses-02_task-rest_echo-3_part-mag_bold.nii.gz", + "bids::sub-01/ses-02/func/sub-01_ses-02_task-rest_echo-3_part-phase_bold.nii.gz", + ], + }, + }, + { + "task": "rest", + "echo": "2", + "part": "phase", + "suffix": "bold", + "metadata": { + "EchoTime": 0.03893, + "RepetitionTime": 0.8, + "TotalReadoutTime": 0.5, + "PhaseEncodingDirection": "j", + "IntendedFor": [ + "bids::sub-01/ses-02/func/sub-01_ses-02_task-rest_echo-1_part-mag_bold.nii.gz", + "bids::sub-01/ses-02/func/sub-01_ses-02_task-rest_echo-1_part-phase_bold.nii.gz", + "bids::sub-01/ses-02/func/sub-01_ses-02_task-rest_echo-2_part-mag_bold.nii.gz", + "bids::sub-01/ses-02/func/sub-01_ses-02_task-rest_echo-2_part-phase_bold.nii.gz", + "bids::sub-01/ses-02/func/sub-01_ses-02_task-rest_echo-3_part-mag_bold.nii.gz", + "bids::sub-01/ses-02/func/sub-01_ses-02_task-rest_echo-3_part-phase_bold.nii.gz", + ], + }, + }, + { + "task": "rest", + "echo": "3", + "part": "mag", + "suffix": "bold", + "metadata": { + "EchoTime": 0.06366, + "RepetitionTime": 0.8, + "TotalReadoutTime": 0.5, + "PhaseEncodingDirection": "j", + "IntendedFor": [ + "bids::sub-01/ses-02/func/sub-01_ses-02_task-rest_echo-1_part-mag_bold.nii.gz", + "bids::sub-01/ses-02/func/sub-01_ses-02_task-rest_echo-1_part-phase_bold.nii.gz", + "bids::sub-01/ses-02/func/sub-01_ses-02_task-rest_echo-2_part-mag_bold.nii.gz", + "bids::sub-01/ses-02/func/sub-01_ses-02_task-rest_echo-2_part-phase_bold.nii.gz", + "bids::sub-01/ses-02/func/sub-01_ses-02_task-rest_echo-3_part-mag_bold.nii.gz", + "bids::sub-01/ses-02/func/sub-01_ses-02_task-rest_echo-3_part-phase_bold.nii.gz", + ], + }, + }, + { + "task": "rest", + "echo": "3", + "part": "phase", + "suffix": "bold", + "metadata": { + "EchoTime": 0.06366, + "RepetitionTime": 0.8, + "TotalReadoutTime": 0.5, + "PhaseEncodingDirection": "j", + "IntendedFor": [ + "bids::sub-01/ses-02/func/sub-01_ses-02_task-rest_echo-1_part-mag_bold.nii.gz", + "bids::sub-01/ses-02/func/sub-01_ses-02_task-rest_echo-1_part-phase_bold.nii.gz", + "bids::sub-01/ses-02/func/sub-01_ses-02_task-rest_echo-2_part-mag_bold.nii.gz", + "bids::sub-01/ses-02/func/sub-01_ses-02_task-rest_echo-2_part-phase_bold.nii.gz", + "bids::sub-01/ses-02/func/sub-01_ses-02_task-rest_echo-3_part-mag_bold.nii.gz", + "bids::sub-01/ses-02/func/sub-01_ses-02_task-rest_echo-3_part-phase_bold.nii.gz", + ], + }, + }, + ], + }, + { + "session": "03", + "anat": [{"suffix": "T1w", "metadata": {"EchoTime": 1}}], + "func": [ + { + "task": "rest", + "echo": "1", + "part": "mag", + "suffix": "bold", + "metadata": { + "EchoTime": 0.0142, + "RepetitionTime": 0.8, + "TotalReadoutTime": 0.5, + "PhaseEncodingDirection": "j", + "IntendedFor": [ + "bids::sub-01/ses-03/func/sub-01_ses-03_task-rest_echo-1_part-mag_bold.nii.gz", + "bids::sub-01/ses-03/func/sub-01_ses-03_task-rest_echo-1_part-phase_bold.nii.gz", + "bids::sub-01/ses-03/func/sub-01_ses-03_task-rest_echo-2_part-mag_bold.nii.gz", + "bids::sub-01/ses-03/func/sub-01_ses-03_task-rest_echo-2_part-phase_bold.nii.gz", + "bids::sub-01/ses-03/func/sub-01_ses-03_task-rest_echo-3_part-mag_bold.nii.gz", + "bids::sub-01/ses-03/func/sub-01_ses-03_task-rest_echo-3_part-phase_bold.nii.gz", + ], + }, + }, + { + "task": "rest", + "echo": "1", + "part": "phase", + "suffix": "bold", + "metadata": { + "EchoTime": 0.0142, + "RepetitionTime": 0.8, + "TotalReadoutTime": 0.5, + "PhaseEncodingDirection": "j", + "IntendedFor": [ + "bids::sub-01/ses-03/func/sub-01_ses-03_task-rest_echo-1_part-mag_bold.nii.gz", + "bids::sub-01/ses-03/func/sub-01_ses-03_task-rest_echo-1_part-phase_bold.nii.gz", + "bids::sub-01/ses-03/func/sub-01_ses-03_task-rest_echo-2_part-mag_bold.nii.gz", + "bids::sub-01/ses-03/func/sub-01_ses-03_task-rest_echo-2_part-phase_bold.nii.gz", + "bids::sub-01/ses-03/func/sub-01_ses-03_task-rest_echo-3_part-mag_bold.nii.gz", + "bids::sub-01/ses-03/func/sub-01_ses-03_task-rest_echo-3_part-phase_bold.nii.gz", + ], + }, + }, + { + "task": "rest", + "echo": "2", + "part": "mag", + "suffix": "bold", + "metadata": { + "EchoTime": 0.03893, + "RepetitionTime": 0.8, + "TotalReadoutTime": 0.5, + "PhaseEncodingDirection": "j", + "IntendedFor": [ + "bids::sub-01/ses-03/func/sub-01_ses-03_task-rest_echo-1_part-mag_bold.nii.gz", + "bids::sub-01/ses-03/func/sub-01_ses-03_task-rest_echo-1_part-phase_bold.nii.gz", + "bids::sub-01/ses-03/func/sub-01_ses-03_task-rest_echo-2_part-mag_bold.nii.gz", + "bids::sub-01/ses-03/func/sub-01_ses-03_task-rest_echo-2_part-phase_bold.nii.gz", + "bids::sub-01/ses-03/func/sub-01_ses-03_task-rest_echo-3_part-mag_bold.nii.gz", + "bids::sub-01/ses-03/func/sub-01_ses-03_task-rest_echo-3_part-phase_bold.nii.gz", + ], + }, + }, + { + "task": "rest", + "echo": "2", + "part": "phase", + "suffix": "bold", + "metadata": { + "EchoTime": 0.03893, + "RepetitionTime": 0.8, + "TotalReadoutTime": 0.5, + "PhaseEncodingDirection": "j", + "IntendedFor": [ + "bids::sub-01/ses-03/func/sub-01_ses-03_task-rest_echo-1_part-mag_bold.nii.gz", + "bids::sub-01/ses-03/func/sub-01_ses-03_task-rest_echo-1_part-phase_bold.nii.gz", + "bids::sub-01/ses-03/func/sub-01_ses-03_task-rest_echo-2_part-mag_bold.nii.gz", + "bids::sub-01/ses-03/func/sub-01_ses-03_task-rest_echo-2_part-phase_bold.nii.gz", + "bids::sub-01/ses-03/func/sub-01_ses-03_task-rest_echo-3_part-mag_bold.nii.gz", + "bids::sub-01/ses-03/func/sub-01_ses-03_task-rest_echo-3_part-phase_bold.nii.gz", + ], + }, + }, + { + "task": "rest", + "echo": "3", + "part": "mag", + "suffix": "bold", + "metadata": { + "EchoTime": 0.06366, + "RepetitionTime": 0.8, + "TotalReadoutTime": 0.5, + "PhaseEncodingDirection": "j", + "IntendedFor": [ + "bids::sub-01/ses-03/func/sub-01_ses-03_task-rest_echo-1_part-mag_bold.nii.gz", + "bids::sub-01/ses-03/func/sub-01_ses-03_task-rest_echo-1_part-phase_bold.nii.gz", + "bids::sub-01/ses-03/func/sub-01_ses-03_task-rest_echo-2_part-mag_bold.nii.gz", + "bids::sub-01/ses-03/func/sub-01_ses-03_task-rest_echo-2_part-phase_bold.nii.gz", + "bids::sub-01/ses-03/func/sub-01_ses-03_task-rest_echo-3_part-mag_bold.nii.gz", + "bids::sub-01/ses-03/func/sub-01_ses-03_task-rest_echo-3_part-phase_bold.nii.gz", + ], + }, + }, + { + "task": "rest", + "echo": "3", + "part": "phase", + "suffix": "bold", + "metadata": { + "EchoTime": 0.06366, + "RepetitionTime": 0.8, + "TotalReadoutTime": 0.5, + "PhaseEncodingDirection": "j", + "IntendedFor": [ + "bids::sub-01/ses-03/func/sub-01_ses-03_task-rest_echo-1_part-mag_bold.nii.gz", + "bids::sub-01/ses-03/func/sub-01_ses-03_task-rest_echo-1_part-phase_bold.nii.gz", + "bids::sub-01/ses-03/func/sub-01_ses-03_task-rest_echo-2_part-mag_bold.nii.gz", + "bids::sub-01/ses-03/func/sub-01_ses-03_task-rest_echo-2_part-phase_bold.nii.gz", + "bids::sub-01/ses-03/func/sub-01_ses-03_task-rest_echo-3_part-mag_bold.nii.gz", + "bids::sub-01/ses-03/func/sub-01_ses-03_task-rest_echo-3_part-phase_bold.nii.gz", + ], + }, + }, + ], + }, + ] +} + filters = { "fmap": { - "datatype": "fmap", + "datatype": ["fmap", "func"], "session": "01", }, "t1w": {"datatype": "anat", "session": "01", "suffix": "T1w"}, @@ -308,6 +691,7 @@ def gen_layout(bids_dir, database_dir=None): [ ('pepolar', pepolar, 1), ('phasediff', phasediff, 1), + ('medic', medic, 1), ], ) def test_wrangler_filter(tmpdir, name, skeleton, estimations): @@ -324,14 +708,18 @@ def test_wrangler_filter(tmpdir, name, skeleton, estimations): [ ('pepolar', pepolar, 5), ('phasediff', phasediff, 3), + ('medic', medic, 3), + ], +) +@pytest.mark.parametrize( + "session, estimations", + [ + ("01", 1), + ("02", 1), + ("03", 1), + (None, None), ], ) -@pytest.mark.parametrize("session, estimations", [ - ("01", 1), - ("02", 1), - ("03", 1), - (None, None), -]) def test_wrangler_URIs(tmpdir, name, skeleton, session, estimations, total_estimations): bids_dir = str(tmpdir / name) generate_bids_skeleton(bids_dir, skeleton) diff --git a/sdcflows/utils/wrangler.py b/sdcflows/utils/wrangler.py index 67743b496d..2ed783eb12 100644 --- a/sdcflows/utils/wrangler.py +++ b/sdcflows/utils/wrangler.py @@ -321,6 +321,37 @@ def find_estimators( FieldmapEstimation(sources=<2 files>, method=, bids_id='auto_...')] + MEDIC fieldmaps are also supported, both with B0FieldIdentifier metadata: + + >>> find_estimators( + ... layout=layouts['ds005250'], + ... subject='04', + ... sessions=['1'], + ... fmapless=False, + ... force_fmapless=False, + ... ) # doctest: +ELLIPSIS + [FieldmapEstimation(sources=<2 files>, method=, + bids_id='sub_04_ses_1_DCAN_fmap_acq_MEGE'), + FieldmapEstimation(sources=<2 files>, method=, + bids_id='sub_04_ses_1_DCAN_fmap_acq_MESE'), + FieldmapEstimation(sources=<10 files>, method=, + bids_id='sub_04_ses_1_acq_MBME_medic')] + + and with IntendedFor metadata: + + >>> find_estimators( + ... layout=layouts['ds005250'], + ... subject='04', + ... sessions=['2'], + ... fmapless=False, + ... force_fmapless=False, + ... ) # doctest: +ELLIPSIS + [FieldmapEstimation(sources=<2 files>, method=, + bids_id='auto_...'), + FieldmapEstimation(sources=<2 files>, method=, + bids_id='auto_...'), + FieldmapEstimation(sources=<10 files>, method=, + bids_id='auto_...')] """ from .misc import create_logger from bids.layout import Query @@ -332,7 +363,6 @@ def find_estimators( base_entities = { "subject": subject, "extension": [".nii", ".nii.gz"], - "part": ["mag", None], "scope": "raw", # Ensure derivatives are not captured } @@ -357,8 +387,13 @@ def find_estimators( # flatten lists from json (tupled in pybids for hashing), then unique b0_ids = reduce( set.union, - (listify(ids) for ids in layout.get_B0FieldIdentifiers(**base_entities)), - set() + ( + listify(ids) + for ids in layout.get_B0FieldIdentifiers( + session=sessions, **base_entities + ) + ), + set(), ) if b0_ids: @@ -453,6 +488,62 @@ def find_estimators( _log_debug_estimation(logger, e, layout.root) estimators.append(e) + # Look for MEDIC field maps + # These need to be complex-valued multi-echo BOLD runs with ``IntendedFor`` + has_intended = tuple() + with suppress(ValueError): + has_intended = layout.get( + **{ + **base_entities, + **{ + 'session': sessions, + 'echo': Query.REQUIRED, + 'part': 'phase', + 'suffix': 'bold', + 'IntendedFor': Query.REQUIRED, + }, + }, + ) + + for bold_fmap in has_intended: + complex_imgs = layout.get( + **{ + **bold_fmap.get_entities(), + **{'part': ['phase', 'mag'], 'echo': Query.ANY}, + } + ) + + current_sources = [est.sources for est in estimators] + current_sources = [ + str(item.path) for sublist in current_sources for item in sublist + ] + if complex_imgs[0].path in current_sources: + logger.debug("Skipping fieldmap %s (already in use)", complex_imgs[0].relpath) + continue + + try: + e = fm.FieldmapEstimation( + [ + fm.FieldmapFile( + img.path, + metadata=_filter_metadata(img.get_metadata(), subject), + ) + for img in complex_imgs + ] + ) + except (ValueError, TypeError) as err: + _log_debug_estimator_fail( + logger, + "unnamed MEDIC", + [], + layout.root, + str(err), + ) + else: + _log_debug_estimation(logger, e, layout.root) + estimators.append(e) + continue + # At this point, only single-PE _epi files WITH ``IntendedFor`` can # be automatically processed. has_intended = tuple() @@ -460,7 +551,11 @@ def find_estimators( has_intended = layout.get( **{ **base_entities, - **{'suffix': 'epi', 'IntendedFor': Query.REQUIRED, 'session': sessions} + **{ + 'suffix': ['epi', 'bold'], + 'IntendedFor': Query.REQUIRED, + 'session': sessions, + }, } ) diff --git a/sdcflows/workflows/fit/medic.py b/sdcflows/workflows/fit/medic.py new file mode 100644 index 0000000000..7d6c5ecfab --- /dev/null +++ b/sdcflows/workflows/fit/medic.py @@ -0,0 +1,139 @@ +# emacs: -*- mode: python; py-indent-offset: 4; indent-tabs-mode: nil -*- +# vi: set ft=python sts=4 ts=4 sw=4 et: +# +# Copyright 2021 The NiPreps Developers +# +# Licensed under the Apache License, Version 2.0 (the "License"); +# you may not use this file except in compliance with the License. +# You may obtain a copy of the License at +# +# http://www.apache.org/licenses/LICENSE-2.0 +# +# Unless required by applicable law or agreed to in writing, software +# distributed under the License is distributed on an "AS IS" BASIS, +# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +# See the License for the specific language governing permissions and +# limitations under the License. +# +# We support and encourage derived works from this project, please read +# about our expectations at +# +# https://www.nipreps.org/community/licensing/ +# +"""Processing of dynamic field maps from complex-valued multi-echo BOLD data.""" + +from nipype.pipeline import engine as pe +from nipype.interfaces import utility as niu +from niworkflows.engine.workflows import LiterateWorkflow as Workflow + +from sdcflows.interfaces.fmap import MEDIC, PhaseMap2rads2 + +INPUT_FIELDS = ("magnitude", "phase", "metadata") + + +def init_medic_wf(name="medic_wf"): + """Create the MEDIC dynamic field estimation workflow. + + Workflow Graph + .. workflow :: + :graph2use: orig + :simple_form: yes + + from sdcflows.workflows.fit.medic import init_medic_wf + + wf = init_medic_wf() # doctest: +SKIP + + Parameters + ---------- + name : :obj:`str` + Name for this workflow + + Inputs + ------ + magnitude : :obj:`list` of :obj:`str` + A list of echo-wise magnitude EPI files that will be fed into MEDIC. + phase : :obj:`list` of :obj:`str` + A list of echo-wise phase EPI files that will be fed into MEDIC. + metadata : :obj:`list` of :obj:`dict` + List of metadata dictionaries corresponding to each of the input magnitude files. + + Outputs + ------- + fmap : :obj:`str` + The path of the estimated fieldmap time series file. Units are Hertz. + displacement : :obj:`list` of :obj:`str` + Path to the displacement time series files. Units are mm. + method: :obj:`str` + Short description of the estimation method that was run. + + Notes + ----- + This workflow performs minimal preparation before running the MEDIC algorithm, + as implemented in ``vandandrew/warpkit``. + + Any downstream processing piplines that use this workflow should include + the following references in their boilerplate BibTeX file: + + - medic: https://doi.org/10.1101/2023.11.28.568744 + """ + workflow = Workflow(name=name) + + workflow.__desc__ = """\ +Volume-wise *B0* nonuniformity maps (or *fieldmaps*) were estimated from +complex-valued, multi-echo EPI data using the MEDIC algorithm (@medic). +""" + + inputnode = pe.Node(niu.IdentityInterface(fields=INPUT_FIELDS), name="inputnode") + outputnode = pe.Node( + niu.IdentityInterface(fields=["fmap", "displacement", "method"]), + name="outputnode", + ) + outputnode.inputs.method = "MEDIC" + + # Write metadata dictionaries to JSON files + write_metadata = pe.MapNode( + niu.Function( + input_names=["metadata"], + output_names=["out_file"], + function=write_json, + ), + iterfield=["metadata"], + name="write_metadata", + ) + workflow.connect([(inputnode, write_metadata, [("metadata", "metadata")])]) + + # Convert phase to radians (-pi to pi, not 0 to 2pi) + phase2rad = pe.MapNode( + PhaseMap2rads2(), + iterfield=["in_file"], + name="phase2rad", + ) + workflow.connect([(inputnode, phase2rad, [("phase", "in_file")])]) + + medic = pe.Node( + MEDIC(), + name="medic", + ) + workflow.connect([ + (inputnode, medic, [("magnitude", "mag_files")]), + (write_metadata, medic, [("out_file", "metadata")]), + (phase2rad, medic, [("out_file", "phase_files")]), + (medic, outputnode, [ + ("native_field_map", "fmap"), + ("displacement_map", "displacement"), + ]), + ]) # fmt:skip + + return workflow + + +def write_json(metadata): + """Write a dictionary to a JSON file.""" + import json + import os + + out_file = os.path.abspath("metadata.json") + with open(out_file, "w") as fobj: + json.dump(metadata, fobj, sort_keys=True, indent=4) + + return out_file diff --git a/sdcflows/workflows/fit/tests/test_medic.py b/sdcflows/workflows/fit/tests/test_medic.py new file mode 100644 index 0000000000..bc1e1cc81b --- /dev/null +++ b/sdcflows/workflows/fit/tests/test_medic.py @@ -0,0 +1,83 @@ +# emacs: -*- mode: python; py-indent-offset: 4; indent-tabs-mode: nil -*- +# vi: set ft=python sts=4 ts=4 sw=4 et: +# +# Copyright 2021 The NiPreps Developers +# +# Licensed under the Apache License, Version 2.0 (the "License"); +# you may not use this file except in compliance with the License. +# You may obtain a copy of the License at +# +# http://www.apache.org/licenses/LICENSE-2.0 +# +# Unless required by applicable law or agreed to in writing, software +# distributed under the License is distributed on an "AS IS" BASIS, +# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +# See the License for the specific language governing permissions and +# limitations under the License. +# +# We support and encourage derived works from this project, please read +# about our expectations at +# +# https://www.nipreps.org/community/licensing/ +# +"""Test phase-difference type of fieldmaps.""" +from pathlib import Path +from json import loads + +import pytest + +from ..medic import init_medic_wf, Workflow + + +@pytest.mark.slow +def test_medic(tmpdir, datadir, workdir, outdir): + """Test creation of the workflow.""" + tmpdir.chdir() + + pattern = 'ds005250/sub-04/ses-2/func/*_part-mag_bold.nii.gz' + magnitude_files = sorted(datadir.glob(pattern)) + phase_files = [f.with_name(f.name.replace("part-mag", "part-phase")) for f in magnitude_files] + metadata_dicts = [ + loads(Path(f.with_name(f.name.replace('.nii.gz', '.json'))).read_text()) + for f in magnitude_files + ] + + wf = Workflow(name=f"medic_{magnitude_files[0].name.replace('.nii.gz', '').replace('-', '_')}") + medic_wf = init_medic_wf() + medic_wf.inputs.inputnode.magnitude = magnitude_files + medic_wf.inputs.inputnode.phase = phase_files + medic_wf.inputs.inputnode.metadata = metadata_dicts + + if outdir: + from ...outputs import init_fmap_derivatives_wf, init_fmap_reports_wf + + outdir = outdir / "unittests" / magnitude_files[0].split("/")[0] + fmap_derivatives_wf = init_fmap_derivatives_wf( + output_dir=str(outdir), + write_coeff=True, + bids_fmap_id="phasediff_id", + ) + fmap_derivatives_wf.inputs.inputnode.source_files = [str(f) for f in magnitude_files] + fmap_derivatives_wf.inputs.inputnode.fmap_meta = metadata_dicts + + fmap_reports_wf = init_fmap_reports_wf( + output_dir=str(outdir), + fmap_type='medic', + ) + fmap_reports_wf.inputs.inputnode.source_files = [str(f) for f in magnitude_files] + + wf.connect([ + (medic_wf, fmap_reports_wf, [ + ("outputnode.fmap", "inputnode.fieldmap"), + ]), + (medic_wf, fmap_derivatives_wf, [ + ("outputnode.fmap", "inputnode.fieldmap"), + ]), + ]) # fmt:skip + else: + wf.add_nodes([medic_wf]) + + if workdir: + wf.base_dir = str(workdir) + + wf.run(plugin="Linear")