enh: increase test coverage of utilities

sdcflows/interfaces/fmap.py

@@ -9,7 +9,6 @@
     File,
     traits,
     SimpleInterface,
-    InputMultiObject,
 )
 
 LOGGER = logging.getLogger("nipype.interface")
@@ -94,70 +93,32 @@ class Phasediff2Fieldmap(SimpleInterface):
     output_spec = _Phasediff2FieldmapOutputSpec
 
     def _run_interface(self, runtime):
-        from ..utils.phasemanip import phdiff2fmap
+        from ..utils.phasemanip import phdiff2fmap, delta_te as _delta_te
 
         self._results["out_file"] = phdiff2fmap(
             self.inputs.in_file, _delta_te(self.inputs.metadata), newpath=runtime.cwd
         )
         return runtime
 
 
-class _Coefficients2WarpInputSpec(BaseInterfaceInputSpec):
-    in_target = File(exist=True, mandatory=True, desc="input EPI data to be corrected")
-    in_coeff = InputMultiObject(File(exists=True), mandatory=True,
-                                desc="input coefficients, after alignment to the EPI data")
-    ro_time = traits.Float(mandatory=True, desc="EPI readout time")
+# Code stub for #119 (B-Spline smoothing)
+# class _Coefficients2WarpInputSpec(BaseInterfaceInputSpec):
+#     in_target = File(exist=True, mandatory=True, desc="input EPI data to be corrected")
+#     in_coeff = InputMultiObject(File(exists=True), mandatory=True,
+#                                 desc="input coefficients, after alignment to the EPI data")
+#     ro_time = traits.Float(mandatory=True, desc="EPI readout time")
 
 
-class _Coefficients2WarpOutputSpec(TraitedSpec):
-    out_warp = File(exists=True)
+# class _Coefficients2WarpOutputSpec(TraitedSpec):
+#     out_warp = File(exists=True)
 
 
-class Coefficients2Warp(SimpleInterface):
-    """Convert coefficients to a full displacements map."""
+# class Coefficients2Warp(SimpleInterface):
+#     """Convert coefficients to a full displacements map."""
 
-    input_spec = _Coefficients2WarpInputSpec
-    output_spec = _Coefficients2WarpOutputSpec
+#     input_spec = _Coefficients2WarpInputSpec
+#     output_spec = _Coefficients2WarpOutputSpec
 
-    def _run_interface(self, runtime):
-        raise NotImplementedError
-        return runtime
-
-
-def _delta_te(in_values, te1=None, te2=None):
-    r"""Read :math:`\Delta_\text{TE}` from BIDS metadata dict."""
-    if isinstance(in_values, float):
-        te2 = in_values
-        te1 = 0.0
-
-    if isinstance(in_values, dict):
-        te1 = in_values.get("EchoTime1")
-        te2 = in_values.get("EchoTime2")
-
-        if not all((te1, te2)):
-            te2 = in_values.get("EchoTimeDifference")
-            te1 = 0
-
-    if isinstance(in_values, list):
-        te2, te1 = in_values
-        if isinstance(te1, list):
-            te1 = te1[1]
-        if isinstance(te2, list):
-            te2 = te2[1]
-
-    # For convienience if both are missing we should give one error about them
-    if te1 is None and te2 is None:
-        raise RuntimeError(
-            "EchoTime1 and EchoTime2 metadata fields not found. "
-            "Please consult the BIDS specification."
-        )
-    if te1 is None:
-        raise RuntimeError(
-            "EchoTime1 metadata field not found. Please consult the BIDS specification."
-        )
-    if te2 is None:
-        raise RuntimeError(
-            "EchoTime2 metadata field not found. Please consult the BIDS specification."
-        )
-
-    return abs(float(te2) - float(te1))
+#     def _run_interface(self, runtime):
+#         raise NotImplementedError
+#         return runtime

sdcflows/utils/phasemanip.py

@@ -21,7 +21,7 @@ def au2rads(in_file, newpath=None):
 
     hdr.set_data_dtype(np.float32)
     hdr.set_xyzt_units("mm")
-    out_file = fname_presuffix(in_file, suffix="_rads", newpath=newpath)
+    out_file = fname_presuffix(str(in_file), suffix="_rads", newpath=newpath)
     nb.Nifti1Image(data, None, hdr).to_filename(out_file)
     return out_file
 
@@ -70,15 +70,67 @@ def subtract_phases(in_phases, in_meta, newpath=None):
 
 def phdiff2fmap(in_file, delta_te, newpath=None):
     """Convert the input phase-difference map into a *fieldmap* in Hz."""
-    import math
     import numpy as np
     import nibabel as nb
     from nipype.utils.filemanip import fname_presuffix
 
-    out_file = fname_presuffix(in_file, suffix="_fmap", newpath=newpath)
+    out_file = fname_presuffix(str(in_file), suffix="_fmap", newpath=newpath)
     image = nb.load(in_file)
-    data = image.get_fdata(dtype="float32") / (2.0 * math.pi * delta_te)
+    data = image.get_fdata(dtype="float32") / (2.0 * np.pi * delta_te)
     nii = nb.Nifti1Image(data, image.affine, image.header)
     nii.set_data_dtype(np.float32)
     nii.to_filename(out_file)
     return out_file
+
+
+def delta_te(in_values):
+    r"""
+    Read :math:`\Delta_\text{TE}` from BIDS metadata dict.
+
+    Examples
+    --------
+    >>> t = delta_te({"EchoTime1": 0.00522, "EchoTime2": 0.00768})
+    >>> f"{t:.5f}"
+    '0.00246'
+
+    >>> t = delta_te({'EchoTimeDifference': 0.00246})
+    >>> f"{t:.5f}"
+    '0.00246'
+
+    >>> delta_te({"EchoTime1": "a"})  # doctest: +IGNORE_EXCEPTION_DETAIL
+    Traceback (most recent call last):
+    ValueError:
+
+    >>> delta_te({"EchoTime2": "a"})  # doctest: +IGNORE_EXCEPTION_DETAIL
+    Traceback (most recent call last):
+    ValueError:
+
+    >>> delta_te({})  # doctest: +IGNORE_EXCEPTION_DETAIL
+    Traceback (most recent call last):
+    ValueError:
+
+    """
+    te2 = in_values.get("EchoTime2")
+    te1 = in_values.get("EchoTime1")
+
+    if te2 is None and te1 is None:
+        try:
+            te2 = float(in_values.get("EchoTimeDifference"))
+            return abs(te2)
+        except TypeError:
+            raise ValueError(
+                "Phase/phase-difference fieldmaps: no echo-times information."
+            )
+        except ValueError:
+            raise ValueError(
+                f"Could not interpret metadata <EchoTimeDifference={te2}>."
+            )
+    try:
+        te2 = float(te2 or "unknown")
+        te1 = float(te1 or "unknown")
+    except ValueError:
+        raise ValueError(
+            f"Could not interpret metadata <EchoTime(1,2)={(te1, te2)}>."
+        )
+
+    return abs(te2 - te1)

sdcflows/utils/tests/test_phasemanip.py

@@ -0,0 +1,39 @@
+"""Test phase manipulation routines."""
+import numpy as np
+import nibabel as nb
+
+from ..phasemanip import au2rads, phdiff2fmap
+
+
+def test_au2rads(tmp_path):
+    """Check the conversion."""
+    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 = au2rads(tmp_path / "testdata.nii.gz")
+
+    assert np.allclose(
+        (data / 4096).astype("float32") * 2.0 * np.pi,
+        nb.load(out_file).get_fdata(dtype="float32")
+    )
+
+
+def test_phdiff2fmap(tmp_path):
+    """Check the conversion."""
+    nb.Nifti1Image(
+        np.ones((5, 5, 5), dtype="float32") * 2.0 * np.pi * 2.46e-3,
+        np.eye(4)
+    ).to_filename(tmp_path / "testdata.nii.gz")
+
+    out_file = phdiff2fmap(tmp_path / "testdata.nii.gz", 2.46e-3)
+
+    assert np.allclose(
+        np.ones((5, 5, 5)),
+        nb.load(out_file).get_fdata(dtype="float32")
+    )

sdcflows/workflows/base.py

@@ -1,13 +1,8 @@
 # emacs: -*- mode: python; py-indent-offset: 4; indent-tabs-mode: nil -*-
 # vi: set ft=python sts=4 ts=4 sw=4 et:
 """Estimate fieldmaps for :abbr:`SDC (susceptibility distortion correction)`."""
-from nipype.pipeline import engine as pe
-from nipype.interfaces import utility as niu
 from nipype import logging
 
-from niworkflows.engine.workflows import LiterateWorkflow as Workflow
-
-
 LOGGER = logging.getLogger('nipype.workflow')
 DEFAULT_MEMORY_MIN_GB = 0.01