diff --git a/sdcflows/interfaces/fmap.py b/sdcflows/interfaces/fmap.py index 3dc9b3b961..4d8ec6529e 100644 --- a/sdcflows/interfaces/fmap.py +++ b/sdcflows/interfaces/fmap.py @@ -1,7 +1,9 @@ # emacs: -*- mode: python; py-indent-offset: 4; indent-tabs-mode: nil -*- # vi: set ft=python sts=4 ts=4 sw=4 et: """Interfaces to deal with the various types of fieldmap sources.""" - +import numpy as np +import nibabel as nb +from nipype.utils.filemanip import fname_presuffix from nipype import logging from nipype.interfaces.base import ( BaseInterfaceInputSpec, @@ -99,3 +101,34 @@ def _run_interface(self, runtime): self.inputs.in_file, _delta_te(self.inputs.metadata), newpath=runtime.cwd ) return runtime + + +class _CheckB0UnitsInputSpec(BaseInterfaceInputSpec): + in_file = File(exists=True, mandatory=True, desc="input fieldmap") + units = traits.Enum("Hz", "rad/s", mandatory=True, desc="fieldmap units") + + +class _CheckB0UnitsOutputSpec(TraitedSpec): + out_file = File(exists=True, desc="output fieldmap in Hz") + + +class CheckB0Units(SimpleInterface): + """Ensure the input fieldmap is given in Hz.""" + + input_spec = _CheckB0UnitsInputSpec + output_spec = _CheckB0UnitsOutputSpec + + def _run_interface(self, runtime): + if self.inputs.units == "Hz": + self._results["out_file"] = self.inputs.in_file + return runtime + + self._results["out_file"] = fname_presuffix( + self.inputs.in_file, suffix="_Hz", newpath=runtime.cwd + ) + img = nb.load(self.inputs.in_file) + data = np.asanyarray(img.dataobj) / (2.0 * np.pi) + img.__class__(data, img.affine, img.header).to_filename( + self._results["out_file"] + ) + return runtime diff --git a/sdcflows/interfaces/tests/test_fmap.py b/sdcflows/interfaces/tests/test_fmap.py new file mode 100644 index 0000000000..5919f5cbce --- /dev/null +++ b/sdcflows/interfaces/tests/test_fmap.py @@ -0,0 +1,24 @@ +"""Test fieldmap interfaces.""" +import numpy as np +import nibabel as nb +import pytest + +from ..fmap import CheckB0Units + + +@pytest.mark.parametrize("units", ("rad/s", "Hz")) +def test_units(tmpdir, units): + """Check the conversion of units.""" + tmpdir.chdir() + hz = np.ones((5, 5, 5), dtype="float32") * 100 + data = hz.copy() + + if units == "rad/s": + data *= 2.0 * np.pi + + nb.Nifti1Image(data, np.eye(4), None).to_filename("data.nii.gz") + out_data = nb.load( + CheckB0Units(units=units, in_file="data.nii.gz").run().outputs.out_file + ).get_fdata(dtype="float32") + + assert np.allclose(hz, out_data) diff --git a/sdcflows/workflows/fit/fieldmap.py b/sdcflows/workflows/fit/fieldmap.py index 0030169a21..7c99670e99 100644 --- a/sdcflows/workflows/fit/fieldmap.py +++ b/sdcflows/workflows/fit/fieldmap.py @@ -148,6 +148,9 @@ def init_fmap_wf(omp_nthreads=1, debug=False, mode="phasediff", name="fmap_wf"): Path to the corresponding magnitude image for anatomical reference. fieldmap : :obj:`str` Path to the fieldmap acquisition (``*_fieldmap.nii[.gz]`` of BIDS). + units : :obj:`str` + Units (`"Hz"` or `"rad/s"`) of the fieldmap (only direct :math:`B_0` + acquisitions with :abbr:`SEI (Spiral-Echo Imaging)` fieldmaps). Outputs ------- @@ -223,19 +226,24 @@ def init_fmap_wf(omp_nthreads=1, debug=False, mode="phasediff", name="fmap_wf"): # fmt: on else: from niworkflows.interfaces.images import IntraModalMerge + from ...interfaces.fmap import CheckB0Units workflow.__desc__ = """\ A *B0* nonuniformity map (or *fieldmap*) was directly measured with -an MRI scheme designed with that purpose (e.g., a spiral pulse sequence). +an MRI scheme designed with that purpose such as SEI (Spiral-Echo Imaging). """ - # Merge input fieldmap images + # Merge input fieldmap images (assumes all are given in the same units!) fmapmrg = pe.Node( IntraModalMerge(zero_based_avg=False, hmc=False), name="fmapmrg" ) + units = pe.Node(CheckB0Units(), name="units", run_without_submitting=True) + # fmt: off workflow.connect([ + (inputnode, units, [("units", "units")]), (inputnode, fmapmrg, [("fieldmap", "in_files")]), - (fmapmrg, bs_filter, [("out_avg", "in_data")]), + (fmapmrg, units, [("out_avg", "in_file")]), + (units, bs_filter, [("out_file", "in_data")]), ]) # fmt: on