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

@@ -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

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)

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 *B<sub>0</sub>* 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