Draft MEDIC dynamic distortion correction method (second attempt)

.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://gin.g-node.org/tsalo/ds005250-sdcflows.git
+        datalad update -r --merge -d ds005250-sdcflows/
+        datalad get -r -J 2 -d ds005250-sdcflows/ ds005250-sdcflows/sub-04/ses-2/
+
     - name: Set FreeSurfer variables
       run: |
         echo "FREESURFER_HOME=$HOME/.cache/freesurfer" >> $GITHUB_ENV

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

min-requirements.txt

@@ -12,3 +12,4 @@ scikit-image==0.18
 scipy==1.8.1
 templateflow
 toml
+warpkit==0.1.1

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]

requirements.txt

@@ -12,3 +12,4 @@ scikit-image>=0.18
 scipy>=1.8.1
 templateflow
 toml
+warpkit==0.1.1

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
 

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()]

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 ative (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

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

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(