From f7ccfc53146048a26570befb2de2b409c47dcb43 Mon Sep 17 00:00:00 2001 From: Oscar Esteban Date: Thu, 26 Nov 2020 15:10:06 +0100 Subject: [PATCH 1/2] FIX: Convert SEI fieldmaps given in rad/s into Hz Adds a lightweight node (`run_without_submitting=True`) wrapping an interface that just checks the units of the input image and divides by 2pi when units are rad/s. Unfortunatelly, we don't currently have data to test these fieldmaps in some integration/smoke test. Resolves: #124. --- sdcflows/interfaces/fmap.py | 32 +++++++++++++++++++++++++- sdcflows/interfaces/tests/test_fmap.py | 24 +++++++++++++++++++ sdcflows/workflows/fit/fieldmap.py | 14 ++++++++--- 3 files changed, 66 insertions(+), 4 deletions(-) create mode 100644 sdcflows/interfaces/tests/test_fmap.py diff --git a/sdcflows/interfaces/fmap.py b/sdcflows/interfaces/fmap.py index 3dc9b3b961..a011aaec5d 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,31 @@ 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 From 42bc9caeeb11ba4f59bbf5c799ee776740245701 Mon Sep 17 00:00:00 2001 From: Oscar Esteban Date: Thu, 26 Nov 2020 15:14:42 +0100 Subject: [PATCH 2/2] sty: black these files --- sdcflows/interfaces/fmap.py | 7 +++++-- 1 file changed, 5 insertions(+), 2 deletions(-) diff --git a/sdcflows/interfaces/fmap.py b/sdcflows/interfaces/fmap.py index a011aaec5d..4d8ec6529e 100644 --- a/sdcflows/interfaces/fmap.py +++ b/sdcflows/interfaces/fmap.py @@ -124,8 +124,11 @@ def _run_interface(self, runtime): return runtime self._results["out_file"] = fname_presuffix( - self.inputs.in_file, suffix="_Hz", newpath=runtime.cwd) + 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"]) + img.__class__(data, img.affine, img.header).to_filename( + self._results["out_file"] + ) return runtime