Draft MEDIC dynamic distortion correction method

Dockerfile

@@ -162,6 +162,13 @@ ENV PATH="/opt/afni-latest:$PATH" \
     AFNI_IMSAVE_WARNINGS="NO" \
     AFNI_PLUGINPATH="/opt/afni-latest"
 
+# CompileMRI 4.0.6
+RUN mkdir /opt/CompileMRI && \
+curl -fsSL --retry 5 https://github.com/korbinian90/CompileMRI.jl/releases/download/v4.0.6/mritools_ubuntu-22.04_4.0.6.tar.gz \
+| tar -xz -C /opt/CompileMRI --strip-components 1
+
+ENV PATH="/opt/CompileMRI/bin:$PATH"
+
 # Create a shared $HOME directory
 RUN useradd -m -s /bin/bash -G users sdcflows
 WORKDIR /home/sdcflows

pyproject.toml

@@ -31,7 +31,7 @@ dependencies = [
     "numpy >= 1.21.0",
     "pybids >= 0.16.4",
     "scikit-image >= 0.18",
-    "scipy >= 1.8.1",
+    "scipy >= 1.8.1,<=1.12.0",
     "templateflow",
     "toml",
 ]

sdcflows/fieldmaps.py

@@ -49,6 +49,7 @@ class EstimatorType(Enum):
     PHASEDIFF = auto()
     MAPPED = auto()
     ANAT = auto()
+    MEDIC = auto()
 
 
 MODALITIES = {

sdcflows/interfaces/fmap.py

@@ -29,6 +29,8 @@
 from nipype import logging
 from nipype.interfaces.base import (
     BaseInterfaceInputSpec,
+    CommandLineInputSpec,
+    CommandLine,
     TraitedSpec,
     File,
     traits,
@@ -62,6 +64,27 @@ 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."""
+
+    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 +413,326 @@ def _check_gross_geometry(
             f"{img1.get_filename()} {''.join(nb.aff2axcodes(img1.affine))}, "
             f"{img2.get_filename()} {''.join(nb.aff2axcodes(img2.affine))}"
         )
+
+
+class _ROMEOInputSpec(CommandLineInputSpec):
+    """Input specification for ApplyAffine."""
+
+    phase_file = File(
+        exists=True,
+        argstr="--phase %s",
+        desc="The phase image that should be unwrapped",
+    )
+    mag_file = File(
+        exists=True,
+        argstr="--magnitude %s",
+        desc="The magnitude image (better unwrapping if specified)",
+    )
+    out_file = File(
+        "unwrapped.nii",
+        argstr="--output %s",
+        usedefault=True,
+        desc="The output path or filename (default: unwrapped.nii)",
+    )
+    echo_times = traits.List(
+        traits.Float,
+        argstr="--echo-times [%s]",
+        desc=(
+            "The echo times required for temporal unwrapping specified in array or range syntax "
+            "(e.g., '[1.5,3.0]' or '3.5:3.5:14'). "
+            "For identical echo times, '-t epi' can be used with the possibility to specify the "
+            "echo time as e.g. '-t epi 5.3' (for B0 calculation)."
+        ),
+    )
+    mask = traits.Either(
+        File(exists=True),
+        traits.Enum(
+            "nomask",
+            "robustmask",
+        ),
+        argstr="--mask %s",
+        desc=(
+            "nomask | qualitymask <threshold> | robustmask "
+            "| <mask_file>. <threshold>=0.1 for qualitymask "
+            "in [0;1] (default: ['robustmask']). "
+            "qualitymask <threshold> isn't supported in this interface."
+        ),
+    )
+    mask_unwrapped = traits.Bool(
+        argstr="--mask-unwrapped",
+        desc=(
+            "Apply the mask on the unwrapped result. "
+            "If mask is 'nomask', sets it to 'robustmask'."
+        ),
+    )
+    weights = traits.Enum(
+        "romeo",
+        "romeo2",
+        "romeo3",
+        "romeo4",
+        "romeo6",
+        "bestpath",
+        argstr="--weights %s",
+        desc=(
+            "romeo | romeo2 | romeo3 | romeo4 | romeo6 | "
+            "bestpath | <4d-weights-file> | <flags>. "
+            "<flags> are up to 6 bits to activate individual weights (eg. '1010'). "
+            "The weights are (1)phasecoherence (2)phasegradientcoherence "
+            "(3)phaselinearity (4)magcoherence (5)magweight (6)magweight2 "
+            "(default: 'romeo')."
+            "4d-weights-file and flags aren't supported in this interface."
+        ),
+    )
+    # TODO: Figure out what the output file would be and populate outputs.
+    calculate_b0 = traits.Bool(
+        argstr="--compute-B0",
+        desc=(
+            "Calculate combined B0 map in [Hz]. "
+            "This activates MCPC3Ds phase offset correction (monopolar) for multi-echo data."
+        ),
+    )
+    phase_offset_correction = traits.Enum(
+        "on",
+        "off",
+        "bipolar",
+        argstr="--phase-offset-correction %s",
+        desc=(
+            "on | off | bipolar. "
+            "Applies the MCPC3Ds method to perform phase offset determination and removal "
+            "(for multi-echo). "
+            "'bipolar' removes eddy current artefacts (requires >= 3 echoes). "
+            "(default: 'off', without arg: 'on')"
+        ),
+    )
+    phase_offset_smoothing_sigma_mm = traits.List(
+        [7, 7, 7],
+        traits.Float,
+        minlen=3,
+        maxlen=3,
+        argstr="--phase-offset-smoothing-sigma-mm %s",
+        usedefault=True,
+        desc=(
+            "default: [7,7,7] "
+            "Only applied if phase-offset-correction is activated. "
+            "The given sigma size is divided by the voxel size from the nifti phase file "
+            "to obtain a smoothing size in voxels. "
+            "A value of [0,0,0] deactivates phase offset smoothing (not recommended)."
+        ),
+    )
+    # TODO: Figure out what the output file would be and populate outputs.
+    write_phase_offsets = traits.Bool(
+        argstr="--write-phase-offsets",
+        desc="Saves the estimated phase offsets to the output folder",
+    )
+    individual_unwrapping = traits.Bool(
+        argstr="--individual-unwrapping",
+        desc=(
+            "Unwraps the echoes individually (not temporal). "
+            "This might be necessary if there is large movement (timeseries) or "
+            "phase-offset-correction is not applicable."
+        ),
+    )
+    template_echo = traits.Int(
+        argstr="--template %d",
+        default_value=1,
+        usedefault=True,
+        desc=(
+            "Template echo that is spatially unwrapped and used for temporal unwrapping "
+            "(type: Int64, default: 1)"
+        ),
+    )
+    no_mmap = traits.Bool(
+        argstr="--no-mmap",
+        desc="Deactivate memory mapping. Memory mapping might cause problems on network storage",
+    )
+    no_rescale = traits.Bool(
+        argstr="--no-rescale",
+        desc=(
+            "Deactivate rescaling of input images. "
+            "By default the input phase is rescaled to the range [-π;π]. "
+            "This option allows inputting already unwrapped phase images without "
+            "manually wrapping them first."
+        ),
+    )
+    threshold = traits.Float(
+        argstr="--threshold %f",
+        desc=(
+            "<maximum number of wraps>. "
+            "Threshold the unwrapped phase to the maximum number of wraps and sets exceeding "
+            "values to 0 (type: Float64, default: Inf)"
+        ),
+    )
+    verbose = traits.Bool(
+        argstr="--verbose",
+        desc="verbose output messages",
+    )
+    correct_global = traits.Bool(
+        argstr="--correct-global",
+        desc=(
+            "Phase is corrected to remove global n2π phase offset. "
+            "The median of phase values (inside mask if given) is used to calculate the "
+            "correction term"
+        ),
+    )
+    # TODO: Figure out what the output file would be and populate outputs.
+    write_quality = traits.Bool(
+        argstr="--write-quality",
+        desc="Writes out the ROMEO quality map as a 3D image with one value per voxel",
+    )
+    # TODO: Figure out what the output files would be and populate outputs.
+    write_quality_all = traits.Bool(
+        argstr="--write-quality-all",
+        desc="Writes out an individual quality map for each of the ROMEO weights.",
+    )
+    max_seeds = traits.Int(
+        argstr="--max-seeds %d",
+        default_value=1,
+        usedefault=True,
+        desc=(
+            "EXPERIMENTAL! "
+            "Sets the maximum number of seeds for unwrapping. "
+            "Higher values allow more separated regions. "
+            "(type: Int64, default: 1)"
+        ),
+    )
+    merge_regions = traits.Bool(
+        argstr="--merge-regions",
+        desc="EXPERIMENTAL! Spatially merges neighboring regions after unwrapping.",
+    )
+    correct_regions = traits.Bool(
+        argstr="--correct-regions",
+        desc=(
+            "EXPERIMENTAL! "
+            "Performed after merging. "
+            "Brings the median of each region closest to 0 (mod 2π)."
+        ),
+    )
+    wrap_addition = traits.Float(
+        argstr="--wrap-addition %f",
+        desc=(
+            "[0;π] "
+            "EXPERIMENTAL! "
+            "Usually the true phase difference of neighboring voxels cannot exceed π "
+            "to be able to unwrap them. "
+            "This setting increases the limit and uses 'linear unwrapping' of 3 voxels in a line. "
+            "Neighbors can have (π + wrap-addition) phase difference. "
+            "(type: Float64, default: 0.0)"
+        ),
+    )
+    temporal_uncertain_unwrapping = traits.Bool(
+        argstr="--temporal-uncertain-unwrapping",
+        desc=(
+            "EXPERIMENTAL! "
+            "Uses spatial unwrapping on voxels that have high uncertainty values after "
+            "temporal unwrapping."
+        ),
+    )
+
+
+class _ROMEOOutputSpec(TraitedSpec):
+    """Output specification for ApplyAffine."""
+
+    out_file = File(exists=True, desc="output file")
+    quality_file = File(desc="Quality file. Only created if write_quality is True.")
+
+
+class ROMEO(CommandLine):
+    """Run ROMEO unwrapping."""
+
+    input_spec = _ROMEOInputSpec
+    output_spec = _ROMEOOutputSpec
+    _cmd = "romeo"
+
+    def _list_outputs(self):
+        outputs = self.output_spec().get()
+        outputs["out_file"] = os.path.abspath(self.inputs.out_file)
+        if self.inputs.write_quality:
+            outputs["quality_file"] = os.path.abspath("quality.nii")
+
+        return outputs
+
+
+class _MEDICB0InputSpec(TraitedSpec):
+    magnitude = traits.List(
+        File(exists=True),
+        mandatory=True,
+        desc="Echo-wise magnitude time series",
+    )
+    phase = traits.List(
+        File(exists=True),
+        mandatory=True,
+        desc="Echo-wise phase time series",
+    )
+    echo_times = traits.List(
+        traits.Float,
+        mandatory=True,
+        desc="the echo times of the EPI image",
+    )
+
+
+class _MEDICB0OutputSpec(TraitedSpec):
+    b0 = File(exists=True, desc="the B0 fieldmap time series")
+
+
+class MEDICB0(SimpleInterface):
+    """Run MEDIC B0 unwrapping."""
+
+    input_spec = _MEDICB0InputSpec
+    output_spec = _MEDICB0OutputSpec
+
+    def _run_interface(self, runtime):
+        import os
+
+        import nibabel as nb
+        import numpy as np
+        from nilearn import image
+
+        from sdcflows.utils.misc import weighted_regression
+
+        magnitude_files = self.inputs.magnitude
+        phase_files = self.inputs.phase
+        echo_times = np.array(self.inputs.echo_times)
+
+        assert len(magnitude_files) == len(phase_files) == len(echo_times)
+
+        temp_img = nb.load(magnitude_files[0])
+        n_volumes = temp_img.shape[3]
+        size = temp_img.shape[:3]
+        n_echoes = echo_times.size
+
+        out_b0 = np.zeros(temp_img.shape)
+        b0_file = os.path.abspath("b0.nii.gz")
+
+        # Split up and transpose the echo-wise data into volume-wise data
+        for i_vol in range(n_volumes):
+            magnitude_volume_imgs = []
+            phase_volume_imgs = []
+            for j_echo in range(n_echoes):
+                magnitude_volume_imgs.append(
+                    nb.load(magnitude_files[j_echo]).slicer[..., i_vol]
+                )
+                phase_volume_imgs.append(
+                    nb.load(phase_files[j_echo]).slicer[..., i_vol]
+                )
+
+            magnitude_volume_img = image.concat_imgs(magnitude_volume_imgs)
+            phase_volume_img = image.concat_imgs(phase_volume_imgs)
+
+            magnitude_volume_data = magnitude_volume_img.get_fdata()
+            phase_volume_data = phase_volume_img.get_fdata()
+
+            unwrapped_mat = phase_volume_data.reshape(-1, n_echoes).T
+            weights = magnitude_volume_data.reshape(-1, n_echoes).T
+            b0 = weighted_regression(
+                echo_times[:, np.newaxis],
+                unwrapped_mat,
+                weights,
+            )[0].T.reshape(*size)
+            b0 *= 1000 / (2 * np.pi)
+            out_b0[:, :, :, i_vol] = b0
+
+        b0_img = nb.Nifti1Image(out_b0, temp_img.affine, temp_img.header)
+        b0_img.to_filename(b0_file)
+        self._results["b0"] = b0_file
+
+        return runtime