diff --git a/Dockerfile b/Dockerfile
index f8fe67c2fa..73ec97108c 100644
--- a/Dockerfile
+++ b/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
diff --git a/pyproject.toml b/pyproject.toml
index fbdaae9e94..1e61e8820e 100644
--- a/pyproject.toml
+++ b/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",
 ]
diff --git a/sdcflows/fieldmaps.py b/sdcflows/fieldmaps.py
index b724a26a8b..228c1aa989 100644
--- a/sdcflows/fieldmaps.py
+++ b/sdcflows/fieldmaps.py
@@ -49,6 +49,7 @@ class EstimatorType(Enum):
     PHASEDIFF = auto()
     MAPPED = auto()
     ANAT = auto()
+    MEDIC = auto()
 
 
 MODALITIES = {
diff --git a/sdcflows/interfaces/fmap.py b/sdcflows/interfaces/fmap.py
index 6327460faf..462056790b 100644
--- a/sdcflows/interfaces/fmap.py
+++ b/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
diff --git a/sdcflows/interfaces/utils.py b/sdcflows/interfaces/utils.py
index 6696c786d7..3678ae24f2 100644
--- a/sdcflows/interfaces/utils.py
+++ b/sdcflows/interfaces/utils.py
@@ -21,26 +21,35 @@
 #     https://www.nipreps.org/community/licensing/
 #
 """Utilities."""
+import os
 from itertools import product
+
+import nibabel as nb
+from nilearn import image
+from nipype import logging
 from nipype.interfaces.base import (
     BaseInterfaceInputSpec,
+    DynamicTraitedSpec,
     TraitedSpec,
     File,
     traits,
     SimpleInterface,
     InputMultiObject,
     OutputMultiObject,
+    Undefined,
     isdefined,
 )
 from nipype.interfaces.ants.segmentation import (
     DenoiseImage as _DenoiseImageBase,
     DenoiseImageInputSpec as _DenoiseImageInputSpecBase,
 )
+from nipype.interfaces.io import add_traits
 from nipype.interfaces.mixins import CopyHeaderInterface as _CopyHeaderInterface
 
 from sdcflows.utils.tools import reorient_pedir
 
 OBLIQUE_THRESHOLD_DEG = 0.5
+LOGGER = logging.getLogger("nipype.interface")
 
 
 class _FlattenInputSpec(BaseInterfaceInputSpec):
@@ -552,3 +561,337 @@ def _ensure_positive_cosines(in_file: str, newpath: str = None):
     reoriented, axcodes = ensure_positive_cosines(nb.load(in_file))
     reoriented.to_filename(out_file)
     return out_file, "".join(axcodes)
+
+
+class _TransposeImagesInputSpec(DynamicTraitedSpec):
+    pass
+
+
+class TransposeImages(SimpleInterface):
+    """Convert input filenames to BIDS URIs, based on links in the dataset.
+
+    This interface can combine multiple lists of inputs.
+    """
+
+    input_spec = _TransposeImagesInputSpec
+    output_spec = DynamicTraitedSpec
+
+    def __init__(self, numinputs=0, numoutputs=0, **inputs):
+        super().__init__(**inputs)
+        self._numinputs = numinputs
+        self._numoutputs = numoutputs
+        if numinputs >= 1:
+            input_names = [f"in{i + 1}" for i in range(numinputs)]
+        else:
+            input_names = []
+        add_traits(self.inputs, input_names)
+
+    def _outputs(self):
+        # Mostly copied from nipype.interfaces.dcmstack.LookupMeta.
+        outputs = super()._outputs()
+
+        undefined_traits = {}
+        if self._numoutputs >= 1:
+            output_names = [f"out{i + 1}" for i in range(self._numoutputs)]
+        else:
+            output_names = []
+
+        for output_name in output_names:
+            outputs.add_trait(output_name, traits.Any)
+            undefined_traits[output_name] = Undefined
+
+        outputs.trait_set(trait_change_notify=False, **undefined_traits)
+
+        return outputs
+
+    def _run_interface(self, runtime):
+        inputs = [getattr(self.inputs, f"in{i + 1}") for i in range(self._numinputs)]
+
+        for i_vol in range(self._numoutputs):
+            vol_imgs = []
+            for file_ in inputs:
+                vol_img = image.index_img(file_, i_vol)
+                vol_imgs.append(vol_img)
+            concat_vol_img = image.concat_imgs(vol_imgs)
+
+            vol_file = os.path.abspath(f"vol{i_vol + 1}.nii.gz")
+            concat_vol_img.to_filename(vol_file)
+            self._results[f"out{i_vol + 1}"] = vol_file
+
+        return runtime
+
+    def _list_outputs(self):
+        outputs = self._outputs().get()
+        outputs.update(self._results)
+        return outputs
+
+
+class _SVDFilterInputSpec(TraitedSpec):
+    fieldmap = File(exists=True, mandatory=True, desc="Fieldmap image")
+    mask = File(exists=True, desc="Mask image")
+    border_filter = traits.List(
+        traits.Float,
+        desc="Border filter parameters",
+        minlen=2,
+        maxlen=2,
+    )
+    svd_filter = traits.Int(desc="SVD filter parameter")
+
+
+class _SVDFilterOutputSpec(TraitedSpec):
+    fieldmap = File(exists=True, desc="Filtered fieldmap image")
+
+
+class SVDFilter(SimpleInterface):
+
+    input_spec = _SVDFilterInputSpec
+    output_spec = _SVDFilterOutputSpec
+
+    def _run_interface(self, runtime):
+        import numpy as np
+        from nipype.utils.filemanip import fname_presuffix
+        from scipy.ndimage import gaussian_filter
+
+        border_filt = self.inputs.border_filter
+        svd_filt = self.inputs.svd_filter
+
+        fieldmap_img = nb.load(self.inputs.fieldmap)
+        mask_img = nb.load(self.inputs.mask)
+
+        n_frames = fieldmap_img.shape[3]
+        voxel_size = fieldmap_img.header.get_zooms()[0]  # assuming isotropic voxels I guess
+
+        mask_data = mask_img.get_fdata().astype(np.uint8)
+        fieldmap_data = fieldmap_img.get_fdata()
+
+        if mask_data.max() == 2 and n_frames >= np.max(border_filt):
+            LOGGER.info("Performing spatial/temporal filtering of border voxels...")
+            smoothed_field_maps = np.zeros(fieldmap_data.shape, dtype=np.float32)
+            # smooth by 4 mm kernel
+            sigma = (4 / voxel_size) / 2.355
+            for i_vol in range(n_frames):
+                smoothed_field_maps[..., i_vol] = gaussian_filter(
+                    fieldmap_data[..., i_vol],
+                    sigma=sigma,
+                )
+
+            # compute the union of all the masks
+            union_mask = np.sum(mask_data, axis=-1) > 0
+
+            # do temporal filtering of border voxels with SVD
+            U, S, VT = np.linalg.svd(smoothed_field_maps[union_mask], full_matrices=False)
+
+            # first pass of SVD filtering
+            recon = np.dot(U[:, : border_filt[0]] * S[: border_filt[0]], VT[: border_filt[0], :])
+            recon_img = np.zeros(fieldmap_data.shape, dtype=np.float32)
+            recon_img[union_mask] = recon
+
+            # set the border voxels in the field map to the recon values
+            for i_vol in range(n_frames):
+                fieldmap_data[mask_data[..., i_vol] == 1, i_vol] = recon_img[
+                    mask_data[..., i_vol] == 1,
+                    i_vol,
+                ]
+
+            # do second SVD filtering pass
+            U, S, VT = np.linalg.svd(fieldmap_data[union_mask], full_matrices=False)
+
+            # second pass of SVD filtering
+            recon = np.dot(U[:, : border_filt[1]] * S[: border_filt[1]], VT[: border_filt[1], :])
+            recon_img = np.zeros(fieldmap_data.shape, dtype=np.float32)
+            recon_img[union_mask] = recon
+
+            # set the border voxels in the field map to the recon values
+            for i_vol in range(n_frames):
+                fieldmap_data[mask_data[..., i_vol] == 1, i_vol] = recon_img[
+                    mask_data[..., i_vol] == 1,
+                    i_vol,
+                ]
+
+        # use svd filter to denoise the field maps
+        if n_frames >= svd_filt:
+            LOGGER.info("Denoising field maps with SVD...")
+            LOGGER.info(f"Keeping {svd_filt} components...")
+
+            # compute the union of all the masks
+            union_mask = np.sum(mask_data, axis=-1) > 0
+
+            # compute SVD
+            U, S, VT = np.linalg.svd(fieldmap_data[union_mask], full_matrices=False)
+
+            # only keep the first n_components components
+            recon = np.dot(U[:, :svd_filt] * S[:svd_filt], VT[:svd_filt, :])
+            recon_img = np.zeros(fieldmap_data.shape, dtype=np.float32)
+            recon_img[union_mask] = recon
+
+            # set the voxel values in the mask to the recon values
+            for i_vol in range(n_frames):
+                fieldmap_data[mask_data[..., i_vol] > 0, i_vol] = recon_img[
+                    mask_data[..., i_vol] > 0,
+                    i_vol,
+                ]
+
+        out_file = fname_presuffix(self.inputs.fieldmap, suffix="_filtered", newpath=runtime.cwd)
+        nb.Nifti1Image(
+            fieldmap_data,
+            fieldmap_img.affine,
+            fieldmap_img.header,
+        ).to_filename(out_file)
+        self._results["fieldmap"] = out_file
+
+        return runtime
+
+
+class _EnforceTemporalConsistencyInputSpec(TraitedSpec):
+    phase_unwrapped = traits.List(
+        File(exists=True),
+        mandatory=True,
+        desc="Unwrapped phase data",
+    )
+    echo_times = traits.List(
+        traits.Float,
+        mandatory=True,
+        desc="Echo times",
+    )
+    magnitude = traits.List(
+        File(exists=True),
+        mandatory=True,
+        desc="Magnitude images",
+    )
+    mask = File(exists=True, mandatory=True, desc="Brain mask")
+    threshold = traits.Float(
+        0.98,
+        usedefault=True,
+        desc="Threshold for correlation similarity",
+    )
+
+
+class _EnforceTemporalConsistencyOutputSpec(TraitedSpec):
+    phase_unwrapped = traits.List(
+        File(exists=True),
+        desc="Unwrapped phase data after temporal consistency is enforced.",
+    )
+
+
+class EnforceTemporalConsistency(SimpleInterface):
+    """Ensure phase unwrapping solutions are temporally consistent.
+
+    This uses correlation as a similarity metric between frames to enforce temporal consistency.
+
+    This is derived from ``warpkit.unwrap.check_temporal_consistency_corr``.
+
+    XXX: Small change from warpkit:
+    I used ``weights_mat[:j_echo, i_vol, :]`` instead of ``weights_mat[:j_echo]``
+    Otherwise, ``coefficients`` is voxels x time instead of just voxels
+    Not sure what the source of the problem is yet.
+    """
+
+    input_spec = _EnforceTemporalConsistencyInputSpec
+    output_spec = _EnforceTemporalConsistencyOutputSpec
+
+    def _run_interface(self, runtime):
+        import numpy as np
+        from nipype.utils.filemanip import fname_presuffix
+
+        from sdcflows.utils.misc import corr2_coeff, create_brain_mask, weighted_regression
+
+        echo_times = np.array(self.inputs.echo_times)
+        n_echoes = echo_times.size
+        mag_shortest = self.inputs.magnitude[0]
+
+        # generate brain mask (with 1 voxel erosion)
+        brain_mask_file = create_brain_mask(mag_shortest, -1)
+        brain_mask = nb.load(brain_mask_file).get_fdata().astype(bool)
+
+        # Concate the unwrapped phase data into a 5D array
+        phase_imgs = [nb.load(phase) for phase in self.inputs.phase_unwrapped]
+        unwrapped_echo_1 = phase_imgs[0].get_fdata()
+        phase_unwrapped = np.stack([img.get_fdata() for img in phase_imgs], axis=3)
+        magnitude_imgs = [nb.load(mag) for mag in self.inputs.magnitude]
+        n_volumes = phase_imgs[0].shape[3]
+        mask_img = nb.load(self.inputs.mask)
+        mask_data = mask_img.get_fdata().astype(bool)
+
+        for i_vol in range(n_volumes):
+            LOGGER.info(f"Computing temporal consistency check for frame: {i_vol}")
+
+            # get the current frame phase
+            current_frame_data = unwrapped_echo_1[brain_mask, i_vol][:, np.newaxis]
+
+            # get the correlation between the current frame and all other frames
+            corr = corr2_coeff(current_frame_data, unwrapped_echo_1[brain_mask, :]).ravel()
+
+            # threhold the RD
+            tmask = corr > self.inputs.threshold
+
+            # get indices of mask
+            indices = np.where(tmask)[0]
+
+            # get mask for frame
+            volume_mask = mask_data[..., i_vol] > 0
+
+            # for each frame compute the mean value along the time axis
+            # (masked by indices and mask)
+            mean_voxels = np.mean(unwrapped_echo_1[volume_mask][:, indices], axis=-1)
+
+            # for this frame figure out the integer multiple that minimizes the value to the
+            # mean voxel
+            int_map = np.round(
+                (mean_voxels - unwrapped_echo_1[volume_mask, i_vol]) / (2 * np.pi)
+            ).astype(int)
+
+            # correct the first echo's data using the integer map
+            phase_unwrapped[volume_mask, 0, i_vol] += 2 * np.pi * int_map
+
+            # format weight matrix
+            weights_mat = np.stack([m.dataobj for m in magnitude_imgs], axis=-1)[volume_mask].T
+
+            # form design matrix
+            X = echo_times[:, np.newaxis]
+
+            # fit subsequent echos to the weighted linear regression from the first echo
+            for j_echo in range(1, n_echoes):
+                # form response matrix
+                Y = phase_unwrapped[volume_mask, :j_echo, i_vol].T
+
+                # fit model to data
+                # XXX: Small change from warpkit:
+                # weights_mat[:j_echo, i_vol, :] instead of weights_mat[:j_echo]
+                # Otherwise, coefficients is voxels x time instead of just voxels
+                # Not sure what the issue is.
+                coefficients, _ = weighted_regression(
+                    X[:j_echo],
+                    Y,
+                    weights_mat[:j_echo, i_vol, :],
+                )
+
+                # get the predicted values for this echo
+                Y_pred = coefficients * echo_times[j_echo]
+
+                # compute the difference and get the integer multiple map
+                int_map = np.round(
+                    (Y_pred - phase_unwrapped[volume_mask, j_echo, i_vol]) / (2 * np.pi)
+                ).astype(int)
+
+                # correct the data using the integer map
+                phase_unwrapped[volume_mask, j_echo, i_vol] += 2 * np.pi * int_map
+
+        out_files = []
+        for i_echo in range(n_echoes):
+            phase_unwrapped_echo = phase_unwrapped[:, :, :, i_echo, :]
+            phase_unwrapped_echo_img = nb.Nifti1Image(
+                phase_unwrapped_echo,
+                phase_imgs[0].affine,
+                phase_imgs[0].header,
+            )
+            out_file = fname_presuffix(
+                self.inputs.phase_unwrapped[i_echo],
+                suffix="_temporallyconsistent",
+                newpath=runtime.cwd,
+            )
+            phase_unwrapped_echo_img.to_filename(out_file)
+            out_files.append(out_file)
+
+        self._results["phase_unwrapped"] = out_files
+
+        return runtime
diff --git a/sdcflows/tests/data/dsD/dataset_description.json b/sdcflows/tests/data/dsD/dataset_description.json
new file mode 100644
index 0000000000..5882169df5
--- /dev/null
+++ b/sdcflows/tests/data/dsD/dataset_description.json
@@ -0,0 +1,11 @@
+{
+    "Name": "Test Dataset D for MEDIC, only empty files",
+    "BIDSVersion": "",
+    "License": "CC0",
+    "Authors": ["Salo T."],
+    "Acknowledgements": "",
+    "HowToAcknowledge": "",
+    "Funding": "",
+    "ReferencesAndLinks": [""],
+    "DatasetDOI": ""
+}
diff --git a/sdcflows/tests/data/dsD/sub-01/func/sub-01_task-rest_echo-1_part-mag_bold.nii.gz b/sdcflows/tests/data/dsD/sub-01/func/sub-01_task-rest_echo-1_part-mag_bold.nii.gz
new file mode 100644
index 0000000000..e69de29bb2
diff --git a/sdcflows/tests/data/dsD/sub-01/func/sub-01_task-rest_echo-1_part-phase_bold.json b/sdcflows/tests/data/dsD/sub-01/func/sub-01_task-rest_echo-1_part-phase_bold.json
new file mode 100644
index 0000000000..fbe4bbcf8f
--- /dev/null
+++ b/sdcflows/tests/data/dsD/sub-01/func/sub-01_task-rest_echo-1_part-phase_bold.json
@@ -0,0 +1,94 @@
+{
+    "B0FieldIdentifier": "medic",
+    "B0FieldSource": "medic",
+    "EchoTime": 0.0142,
+    "FlipAngle": 68,
+    "ImageType": [
+        "ORIGINAL",
+        "PRIMARY",
+        "P",
+        "MB",
+        "TE1",
+        "ND",
+        "MOSAIC",
+        "PHASE"
+    ],
+    "PhaseEncodingDirection": "j",
+    "RepetitionTime": 1.761,
+    "SliceTiming": [
+        0,
+        0.725,
+        1.45,
+        0.435,
+        1.16,
+        0.145,
+        0.87,
+        1.595,
+        0.58,
+        1.305,
+        0.29,
+        1.015,
+        0,
+        0.725,
+        1.45,
+        0.435,
+        1.16,
+        0.145,
+        0.87,
+        1.595,
+        0.58,
+        1.305,
+        0.29,
+        1.015,
+        0,
+        0.725,
+        1.45,
+        0.435,
+        1.16,
+        0.145,
+        0.87,
+        1.595,
+        0.58,
+        1.305,
+        0.29,
+        1.015,
+        0,
+        0.725,
+        1.45,
+        0.435,
+        1.16,
+        0.145,
+        0.87,
+        1.595,
+        0.58,
+        1.305,
+        0.29,
+        1.015,
+        0,
+        0.725,
+        1.45,
+        0.435,
+        1.16,
+        0.145,
+        0.87,
+        1.595,
+        0.58,
+        1.305,
+        0.29,
+        1.015,
+        0,
+        0.725,
+        1.45,
+        0.435,
+        1.16,
+        0.145,
+        0.87,
+        1.595,
+        0.58,
+        1.305,
+        0.29,
+        1.015
+    ],
+    "SpacingBetweenSlices": 2,
+    "Units": "arbitrary"
+}
\ No newline at end of file
diff --git a/sdcflows/tests/data/dsD/sub-01/func/sub-01_task-rest_echo-1_part-phase_bold.nii.gz b/sdcflows/tests/data/dsD/sub-01/func/sub-01_task-rest_echo-1_part-phase_bold.nii.gz
new file mode 100644
index 0000000000..e69de29bb2
diff --git a/sdcflows/tests/data/dsD/sub-01/func/sub-01_task-rest_echo-2_part-mag_bold.json b/sdcflows/tests/data/dsD/sub-01/func/sub-01_task-rest_echo-2_part-mag_bold.json
new file mode 100644
index 0000000000..3475655916
--- /dev/null
+++ b/sdcflows/tests/data/dsD/sub-01/func/sub-01_task-rest_echo-2_part-mag_bold.json
@@ -0,0 +1,93 @@
+{
+    "B0FieldIdentifier": "medic",
+    "B0FieldSource": "medic",
+    "EchoTime": 0.03893,
+    "FlipAngle": 68,
+    "ImageType": [
+        "ORIGINAL",
+        "PRIMARY",
+        "M",
+        "MB",
+        "TE2",
+        "ND",
+        "NORM",
+        "MOSAIC"
+    ],
+    "PhaseEncodingDirection": "j",
+    "RepetitionTime": 1.761,
+    "SliceTiming": [
+        0,
+        0.725,
+        1.45,
+        0.435,
+        1.16,
+        0.145,
+        0.87,
+        1.595,
+        0.58,
+        1.305,
+        0.29,
+        1.015,
+        0,
+        0.725,
+        1.45,
+        0.435,
+        1.16,
+        0.145,
+        0.87,
+        1.595,
+        0.58,
+        1.305,
+        0.29,
+        1.015,
+        0,
+        0.725,
+        1.45,
+        0.435,
+        1.16,
+        0.145,
+        0.87,
+        1.595,
+        0.58,
+        1.305,
+        0.29,
+        1.015,
+        0,
+        0.725,
+        1.45,
+        0.435,
+        1.16,
+        0.145,
+        0.87,
+        1.595,
+        0.58,
+        1.305,
+        0.29,
+        1.015,
+        0,
+        0.725,
+        1.45,
+        0.435,
+        1.16,
+        0.145,
+        0.87,
+        1.595,
+        0.58,
+        1.305,
+        0.29,
+        1.015,
+        0,
+        0.725,
+        1.45,
+        0.435,
+        1.16,
+        0.145,
+        0.87,
+        1.595,
+        0.58,
+        1.305,
+        0.29,
+        1.015
+    ],
+    "SpacingBetweenSlices": 2
+}
\ No newline at end of file
diff --git a/sdcflows/tests/data/dsD/sub-01/func/sub-01_task-rest_echo-2_part-mag_bold.nii.gz b/sdcflows/tests/data/dsD/sub-01/func/sub-01_task-rest_echo-2_part-mag_bold.nii.gz
new file mode 100644
index 0000000000..e69de29bb2
diff --git a/sdcflows/tests/data/dsD/sub-01/func/sub-01_task-rest_echo-2_part-phase_bold.json b/sdcflows/tests/data/dsD/sub-01/func/sub-01_task-rest_echo-2_part-phase_bold.json
new file mode 100644
index 0000000000..1691b161de
--- /dev/null
+++ b/sdcflows/tests/data/dsD/sub-01/func/sub-01_task-rest_echo-2_part-phase_bold.json
@@ -0,0 +1,94 @@
+{
+    "B0FieldIdentifier": "medic",
+    "B0FieldSource": "medic",
+    "EchoTime": 0.03893,
+    "FlipAngle": 68,
+    "ImageType": [
+        "ORIGINAL",
+        "PRIMARY",
+        "P",
+        "MB",
+        "TE2",
+        "ND",
+        "MOSAIC",
+        "PHASE"
+    ],
+    "PhaseEncodingDirection": "j",
+    "RepetitionTime": 1.761,
+    "SliceTiming": [
+        0,
+        0.725,
+        1.45,
+        0.435,
+        1.16,
+        0.145,
+        0.87,
+        1.595,
+        0.58,
+        1.305,
+        0.29,
+        1.015,
+        0,
+        0.725,
+        1.45,
+        0.435,
+        1.16,
+        0.145,
+        0.87,
+        1.595,
+        0.58,
+        1.305,
+        0.29,
+        1.015,
+        0,
+        0.725,
+        1.45,
+        0.435,
+        1.16,
+        0.145,
+        0.87,
+        1.595,
+        0.58,
+        1.305,
+        0.29,
+        1.015,
+        0,
+        0.725,
+        1.45,
+        0.435,
+        1.16,
+        0.145,
+        0.87,
+        1.595,
+        0.58,
+        1.305,
+        0.29,
+        1.015,
+        0,
+        0.725,
+        1.45,
+        0.435,
+        1.16,
+        0.145,
+        0.87,
+        1.595,
+        0.58,
+        1.305,
+        0.29,
+        1.015,
+        0,
+        0.725,
+        1.45,
+        0.435,
+        1.16,
+        0.145,
+        0.87,
+        1.595,
+        0.58,
+        1.305,
+        0.29,
+        1.015
+    ],
+    "SpacingBetweenSlices": 2,
+    "Units": "arbitrary"
+}
\ No newline at end of file
diff --git a/sdcflows/tests/data/dsD/sub-01/func/sub-01_task-rest_echo-2_part-phase_bold.nii.gz b/sdcflows/tests/data/dsD/sub-01/func/sub-01_task-rest_echo-2_part-phase_bold.nii.gz
new file mode 100644
index 0000000000..e69de29bb2
diff --git a/sdcflows/tests/data/dsD/sub-01/func/sub-01_task-rest_echo-3_part-mag_bold.json b/sdcflows/tests/data/dsD/sub-01/func/sub-01_task-rest_echo-3_part-mag_bold.json
new file mode 100644
index 0000000000..773894cc81
--- /dev/null
+++ b/sdcflows/tests/data/dsD/sub-01/func/sub-01_task-rest_echo-3_part-mag_bold.json
@@ -0,0 +1,93 @@
+{
+    "B0FieldIdentifier": "medic",
+    "B0FieldSource": "medic",
+    "EchoTime": 0.06366,
+    "FlipAngle": 68,
+    "ImageType": [
+        "ORIGINAL",
+        "PRIMARY",
+        "M",
+        "MB",
+        "TE3",
+        "ND",
+        "NORM",
+        "MOSAIC"
+    ],
+    "PhaseEncodingDirection": "j",
+    "RepetitionTime": 1.761,
+    "SliceTiming": [
+        0,
+        0.725,
+        1.45,
+        0.435,
+        1.16,
+        0.145,
+        0.87,
+        1.595,
+        0.58,
+        1.305,
+        0.29,
+        1.015,
+        0,
+        0.725,
+        1.45,
+        0.435,
+        1.16,
+        0.145,
+        0.87,
+        1.595,
+        0.58,
+        1.305,
+        0.29,
+        1.015,
+        0,
+        0.725,
+        1.45,
+        0.435,
+        1.16,
+        0.145,
+        0.87,
+        1.595,
+        0.58,
+        1.305,
+        0.29,
+        1.015,
+        0,
+        0.725,
+        1.45,
+        0.435,
+        1.16,
+        0.145,
+        0.87,
+        1.595,
+        0.58,
+        1.305,
+        0.29,
+        1.015,
+        0,
+        0.725,
+        1.45,
+        0.435,
+        1.16,
+        0.145,
+        0.87,
+        1.595,
+        0.58,
+        1.305,
+        0.29,
+        1.015,
+        0,
+        0.725,
+        1.45,
+        0.435,
+        1.16,
+        0.145,
+        0.87,
+        1.595,
+        0.58,
+        1.305,
+        0.29,
+        1.015
+    ],
+    "SpacingBetweenSlices": 2
+}
\ No newline at end of file
diff --git a/sdcflows/tests/data/dsD/sub-01/func/sub-01_task-rest_echo-3_part-mag_bold.nii.gz b/sdcflows/tests/data/dsD/sub-01/func/sub-01_task-rest_echo-3_part-mag_bold.nii.gz
new file mode 100644
index 0000000000..e69de29bb2
diff --git a/sdcflows/tests/data/dsD/sub-01/func/sub-01_task-rest_echo-3_part-phase_bold.json b/sdcflows/tests/data/dsD/sub-01/func/sub-01_task-rest_echo-3_part-phase_bold.json
new file mode 100644
index 0000000000..523a7f28cd
--- /dev/null
+++ b/sdcflows/tests/data/dsD/sub-01/func/sub-01_task-rest_echo-3_part-phase_bold.json
@@ -0,0 +1,94 @@
+{
+    "B0FieldIdentifier": "medic",
+    "B0FieldSource": "medic",
+    "EchoTime": 0.06366,
+    "FlipAngle": 68,
+    "ImageType": [
+        "ORIGINAL",
+        "PRIMARY",
+        "P",
+        "MB",
+        "TE3",
+        "ND",
+        "MOSAIC",
+        "PHASE"
+    ],
+    "PhaseEncodingDirection": "j",
+    "RepetitionTime": 1.761,
+    "SliceTiming": [
+        0,
+        0.725,
+        1.45,
+        0.435,
+        1.16,
+        0.145,
+        0.87,
+        1.595,
+        0.58,
+        1.305,
+        0.29,
+        1.015,
+        0,
+        0.725,
+        1.45,
+        0.435,
+        1.16,
+        0.145,
+        0.87,
+        1.595,
+        0.58,
+        1.305,
+        0.29,
+        1.015,
+        0,
+        0.725,
+        1.45,
+        0.435,
+        1.16,
+        0.145,
+        0.87,
+        1.595,
+        0.58,
+        1.305,
+        0.29,
+        1.015,
+        0,
+        0.725,
+        1.45,
+        0.435,
+        1.16,
+        0.145,
+        0.87,
+        1.595,
+        0.58,
+        1.305,
+        0.29,
+        1.015,
+        0,
+        0.725,
+        1.45,
+        0.435,
+        1.16,
+        0.145,
+        0.87,
+        1.595,
+        0.58,
+        1.305,
+        0.29,
+        1.015,
+        0,
+        0.725,
+        1.45,
+        0.435,
+        1.16,
+        0.145,
+        0.87,
+        1.595,
+        0.58,
+        1.305,
+        0.29,
+        1.015
+    ],
+    "SpacingBetweenSlices": 2,
+    "Units": "arbitrary"
+}
\ No newline at end of file
diff --git a/sdcflows/tests/data/dsD/sub-01/func/sub-01_task-rest_echo-3_part-phase_bold.nii.gz b/sdcflows/tests/data/dsD/sub-01/func/sub-01_task-rest_echo-3_part-phase_bold.nii.gz
new file mode 100644
index 0000000000..e69de29bb2
diff --git a/sdcflows/tests/data/dsD/sub-01/func/sub-01_task-rest_echo-4_part-mag_bold.json b/sdcflows/tests/data/dsD/sub-01/func/sub-01_task-rest_echo-4_part-mag_bold.json
new file mode 100644
index 0000000000..3e022762f8
--- /dev/null
+++ b/sdcflows/tests/data/dsD/sub-01/func/sub-01_task-rest_echo-4_part-mag_bold.json
@@ -0,0 +1,93 @@
+{
+    "B0FieldIdentifier": "medic",
+    "B0FieldSource": "medic",
+    "EchoTime": 0.08839,
+    "FlipAngle": 68,
+    "ImageType": [
+        "ORIGINAL",
+        "PRIMARY",
+        "M",
+        "MB",
+        "TE4",
+        "ND",
+        "NORM",
+        "MOSAIC"
+    ],
+    "PhaseEncodingDirection": "j",
+    "RepetitionTime": 1.761,
+    "SliceTiming": [
+        0,
+        0.725,
+        1.45,
+        0.435,
+        1.16,
+        0.145,
+        0.87,
+        1.595,
+        0.58,
+        1.305,
+        0.29,
+        1.015,
+        0,
+        0.725,
+        1.45,
+        0.435,
+        1.16,
+        0.145,
+        0.87,
+        1.595,
+        0.58,
+        1.305,
+        0.29,
+        1.015,
+        0,
+        0.725,
+        1.45,
+        0.435,
+        1.16,
+        0.145,
+        0.87,
+        1.595,
+        0.58,
+        1.305,
+        0.29,
+        1.015,
+        0,
+        0.725,
+        1.45,
+        0.435,
+        1.16,
+        0.145,
+        0.87,
+        1.595,
+        0.58,
+        1.305,
+        0.29,
+        1.015,
+        0,
+        0.725,
+        1.45,
+        0.435,
+        1.16,
+        0.145,
+        0.87,
+        1.595,
+        0.58,
+        1.305,
+        0.29,
+        1.015,
+        0,
+        0.725,
+        1.45,
+        0.435,
+        1.16,
+        0.145,
+        0.87,
+        1.595,
+        0.58,
+        1.305,
+        0.29,
+        1.015
+    ],
+    "SpacingBetweenSlices": 2
+}
\ No newline at end of file
diff --git a/sdcflows/tests/data/dsD/sub-01/func/sub-01_task-rest_echo-4_part-mag_bold.nii.gz b/sdcflows/tests/data/dsD/sub-01/func/sub-01_task-rest_echo-4_part-mag_bold.nii.gz
new file mode 100644
index 0000000000..e69de29bb2
diff --git a/sdcflows/tests/data/dsD/sub-01/func/sub-01_task-rest_echo-4_part-phase_bold.json b/sdcflows/tests/data/dsD/sub-01/func/sub-01_task-rest_echo-4_part-phase_bold.json
new file mode 100644
index 0000000000..bdbaee1286
--- /dev/null
+++ b/sdcflows/tests/data/dsD/sub-01/func/sub-01_task-rest_echo-4_part-phase_bold.json
@@ -0,0 +1,94 @@
+{
+    "B0FieldIdentifier": "medic",
+    "B0FieldSource": "medic",
+    "EchoTime": 0.08839,
+    "FlipAngle": 68,
+    "ImageType": [
+        "ORIGINAL",
+        "PRIMARY",
+        "P",
+        "MB",
+        "TE4",
+        "ND",
+        "MOSAIC",
+        "PHASE"
+    ],
+    "PhaseEncodingDirection": "j",
+    "RepetitionTime": 1.761,
+    "SliceTiming": [
+        0,
+        0.725,
+        1.45,
+        0.435,
+        1.16,
+        0.145,
+        0.87,
+        1.595,
+        0.58,
+        1.305,
+        0.29,
+        1.015,
+        0,
+        0.725,
+        1.45,
+        0.435,
+        1.16,
+        0.145,
+        0.87,
+        1.595,
+        0.58,
+        1.305,
+        0.29,
+        1.015,
+        0,
+        0.725,
+        1.45,
+        0.435,
+        1.16,
+        0.145,
+        0.87,
+        1.595,
+        0.58,
+        1.305,
+        0.29,
+        1.015,
+        0,
+        0.725,
+        1.45,
+        0.435,
+        1.16,
+        0.145,
+        0.87,
+        1.595,
+        0.58,
+        1.305,
+        0.29,
+        1.015,
+        0,
+        0.725,
+        1.45,
+        0.435,
+        1.16,
+        0.145,
+        0.87,
+        1.595,
+        0.58,
+        1.305,
+        0.29,
+        1.015,
+        0,
+        0.725,
+        1.45,
+        0.435,
+        1.16,
+        0.145,
+        0.87,
+        1.595,
+        0.58,
+        1.305,
+        0.29,
+        1.015
+    ],
+    "SpacingBetweenSlices": 2,
+    "Units": "arbitrary"
+}
\ No newline at end of file
diff --git a/sdcflows/tests/data/dsD/sub-01/func/sub-01_task-rest_echo-4_part-phase_bold.nii.gz b/sdcflows/tests/data/dsD/sub-01/func/sub-01_task-rest_echo-4_part-phase_bold.nii.gz
new file mode 100644
index 0000000000..e69de29bb2
diff --git a/sdcflows/tests/data/dsD/sub-01/func/sub-01_task-rest_echo-5_part-mag_bold.json b/sdcflows/tests/data/dsD/sub-01/func/sub-01_task-rest_echo-5_part-mag_bold.json
new file mode 100644
index 0000000000..54a616bbf5
--- /dev/null
+++ b/sdcflows/tests/data/dsD/sub-01/func/sub-01_task-rest_echo-5_part-mag_bold.json
@@ -0,0 +1,93 @@
+{
+    "B0FieldIdentifier": "medic",
+    "B0FieldSource": "medic",
+    "EchoTime": 0.11312,
+    "FlipAngle": 68,
+    "ImageType": [
+        "ORIGINAL",
+        "PRIMARY",
+        "M",
+        "MB",
+        "TE5",
+        "ND",
+        "NORM",
+        "MOSAIC"
+    ],
+    "PhaseEncodingDirection": "j",
+    "RepetitionTime": 1.761,
+    "SliceTiming": [
+        0,
+        0.725,
+        1.45,
+        0.435,
+        1.16,
+        0.145,
+        0.87,
+        1.595,
+        0.58,
+        1.305,
+        0.29,
+        1.015,
+        0,
+        0.725,
+        1.45,
+        0.435,
+        1.16,
+        0.145,
+        0.87,
+        1.595,
+        0.58,
+        1.305,
+        0.29,
+        1.015,
+        0,
+        0.725,
+        1.45,
+        0.435,
+        1.16,
+        0.145,
+        0.87,
+        1.595,
+        0.58,
+        1.305,
+        0.29,
+        1.015,
+        0,
+        0.725,
+        1.45,
+        0.435,
+        1.16,
+        0.145,
+        0.87,
+        1.595,
+        0.58,
+        1.305,
+        0.29,
+        1.015,
+        0,
+        0.725,
+        1.45,
+        0.435,
+        1.16,
+        0.145,
+        0.87,
+        1.595,
+        0.58,
+        1.305,
+        0.29,
+        1.015,
+        0,
+        0.725,
+        1.45,
+        0.435,
+        1.16,
+        0.145,
+        0.87,
+        1.595,
+        0.58,
+        1.305,
+        0.29,
+        1.015
+    ],
+    "SpacingBetweenSlices": 2
+}
\ No newline at end of file
diff --git a/sdcflows/tests/data/dsD/sub-01/func/sub-01_task-rest_echo-5_part-mag_bold.nii.gz b/sdcflows/tests/data/dsD/sub-01/func/sub-01_task-rest_echo-5_part-mag_bold.nii.gz
new file mode 100644
index 0000000000..e69de29bb2
diff --git a/sdcflows/tests/data/dsD/sub-01/func/sub-01_task-rest_echo-5_part-phase_bold.json b/sdcflows/tests/data/dsD/sub-01/func/sub-01_task-rest_echo-5_part-phase_bold.json
new file mode 100644
index 0000000000..7714ed273a
--- /dev/null
+++ b/sdcflows/tests/data/dsD/sub-01/func/sub-01_task-rest_echo-5_part-phase_bold.json
@@ -0,0 +1,94 @@
+{
+    "B0FieldIdentifier": "medic",
+    "B0FieldSource": "medic",
+    "EchoTime": 0.11312,
+    "FlipAngle": 68,
+    "ImageType": [
+        "ORIGINAL",
+        "PRIMARY",
+        "P",
+        "MB",
+        "TE5",
+        "ND",
+        "MOSAIC",
+        "PHASE"
+    ],
+    "PhaseEncodingDirection": "j",
+    "RepetitionTime": 1.761,
+    "SliceTiming": [
+        0,
+        0.725,
+        1.45,
+        0.435,
+        1.16,
+        0.145,
+        0.87,
+        1.595,
+        0.58,
+        1.305,
+        0.29,
+        1.015,
+        0,
+        0.725,
+        1.45,
+        0.435,
+        1.16,
+        0.145,
+        0.87,
+        1.595,
+        0.58,
+        1.305,
+        0.29,
+        1.015,
+        0,
+        0.725,
+        1.45,
+        0.435,
+        1.16,
+        0.145,
+        0.87,
+        1.595,
+        0.58,
+        1.305,
+        0.29,
+        1.015,
+        0,
+        0.725,
+        1.45,
+        0.435,
+        1.16,
+        0.145,
+        0.87,
+        1.595,
+        0.58,
+        1.305,
+        0.29,
+        1.015,
+        0,
+        0.725,
+        1.45,
+        0.435,
+        1.16,
+        0.145,
+        0.87,
+        1.595,
+        0.58,
+        1.305,
+        0.29,
+        1.015,
+        0,
+        0.725,
+        1.45,
+        0.435,
+        1.16,
+        0.145,
+        0.87,
+        1.595,
+        0.58,
+        1.305,
+        0.29,
+        1.015
+    ],
+    "SpacingBetweenSlices": 2,
+    "Units": "arbitrary"
+}
\ No newline at end of file
diff --git a/sdcflows/tests/data/dsD/sub-01/func/sub-01_task-rest_echo-5_part-phase_bold.nii.gz b/sdcflows/tests/data/dsD/sub-01/func/sub-01_task-rest_echo-5_part-phase_bold.nii.gz
new file mode 100644
index 0000000000..e69de29bb2
diff --git a/sdcflows/utils/misc.py b/sdcflows/utils/misc.py
index b87f37bac4..3cedd4548d 100644
--- a/sdcflows/utils/misc.py
+++ b/sdcflows/utils/misc.py
@@ -22,6 +22,11 @@
 #
 """Basic miscellaneous utilities."""
 import logging
+from scipy.stats import mode
+from typing import Tuple, cast
+
+import numpy as np
+import numpy.typing as npt
 
 
 def front(inlist):
@@ -82,3 +87,696 @@ def create_logger(name: str, level: int = 40) -> logging.Logger:
     handler.setFormatter(formatter)
     logger.addHandler(handler)
     return logger
+
+
+def get_largest_connected_component(mask_data: npt.NDArray[np.bool_]) -> npt.NDArray[np.bool_]:
+    """Get the largest connected component of a mask
+
+    Parameters
+    ----------
+    mask_data : :obj:`numpy.ndarray`
+        Mask to get the largest connected component of
+
+    Returns
+    -------
+    :obj:`numpy.ndarray`
+        Mask with only the largest connected component
+    """
+    from skimage.measure import label, regionprops
+
+    # get the largest connected component
+    labelled_mask = label(mask_data)
+    props = regionprops(labelled_mask)
+    sorted_props = sorted(props, key=lambda x: x.area, reverse=True)
+    mask_data = labelled_mask == sorted_props[0].label
+
+    # return the mask
+    return mask_data
+
+
+def medic_automask(mag_file, voxel_quality, echo_times, automask_dilation, out_file="mask.nii.gz"):
+    import os
+    from typing import cast
+
+    import nibabel as nb
+    import numpy as np
+    import numpy.typing as npt
+    from scipy.ndimage import (
+        binary_dilation,
+        binary_erosion,
+        binary_fill_holes,
+        generate_binary_structure,
+    )
+    from skimage.filters import threshold_otsu
+
+    from sdcflows.utils.misc import create_brain_mask, get_largest_connected_component
+
+    out_file = os.path.abspath(out_file)
+
+    mag_img = nb.load(mag_file)
+
+    vq = nb.load(voxel_quality).get_fdata()
+    vq[np.isnan(vq)] = 0  # I (TS) am seeing NaNs in this file, which breaks things.
+    vq_mask = vq > threshold_otsu(vq)
+    strel = generate_binary_structure(3, 2)
+    vq_mask = binary_fill_holes(vq_mask, strel).astype(bool)
+    # get largest connected component
+    vq_mask = get_largest_connected_component(vq_mask)
+
+    # combine masks
+    brain_mask_file = create_brain_mask(mag_file)
+    brain_mask = nb.load(brain_mask_file).get_fdata().astype(bool)
+    combined_mask = brain_mask | vq_mask
+    combined_mask = get_largest_connected_component(combined_mask)
+
+    # erode then dilate
+    combined_mask = cast(npt.NDArray[np.bool_], binary_erosion(combined_mask, strel, iterations=2))
+    combined_mask = get_largest_connected_component(combined_mask)
+    combined_mask = cast(
+        npt.NDArray[np.bool_],
+        binary_dilation(combined_mask, strel, iterations=2),
+    )
+
+    # get a dilated verision of the mask
+    combined_mask_dilated = cast(
+        npt.NDArray[np.bool_],
+        binary_dilation(combined_mask, strel, iterations=automask_dilation),
+    )
+
+    # get sum of masks (we can select dilated vs original version by indexing)
+    mask_data_select = combined_mask.astype(np.int8) + combined_mask_dilated.astype(np.int8)
+
+    # let mask_data be the dilated version
+    mask_data = mask_data_select > 0
+
+    mask_img = nb.Nifti1Image(mask_data, mag_img.affine, mag_img.header)
+    mask_img.to_filename(out_file)
+
+    return out_file, None
+
+
+def calculate_diffs(mag_file, phase_file):
+    """Calculate magnitude and phase differences between first two echoes.
+
+    Parameters
+    ----------
+    mag_file : str
+        Magnitude image file, with shape (x, y, z, n_echoes).
+    phase_file : str
+        Phase image file, with shape (x, y, z, n_echoes).
+
+    Returns
+    -------
+    mag_diff_file : str
+        Magnitude difference between echoes 2 and 1, with shape (x, y, z).
+    phase_diff_file : str
+        Phase difference between echoes 2 and 1, with shape (x, y, z).
+
+    Notes
+    -----
+    This calculates a signal difference map between the first two echoes using both the
+    magnitude and phase data,
+    then separates the magnitude and phase differences out.
+    """
+    import os
+
+    import nibabel as nb
+    import numpy as np
+
+    mag_img = nb.load(mag_file)
+    phase_img = nb.load(phase_file)
+    mag_data = mag_img.get_fdata()
+    phase_data = phase_img.get_fdata()
+    signal_diff = (
+        mag_data[..., 0]
+        * mag_data[..., 1]
+        * np.exp(1j * (phase_data[..., 1] - phase_data[..., 0]))
+    )
+    mag_diff = np.abs(signal_diff)
+    phase_diff = np.angle(signal_diff)
+    mag_diff_img = nb.Nifti1Image(mag_diff, mag_img.affine, mag_img.header)
+    phase_diff_img = nb.Nifti1Image(phase_diff, phase_img.affine, phase_img.header)
+    mag_diff_file = os.path.abspath("mag_diff.nii.gz")
+    phase_diff_file = os.path.abspath("phase_diff.nii.gz")
+    mag_diff_img.to_filename(mag_diff_file)
+    phase_diff_img.to_filename(phase_diff_file)
+    return mag_diff_file, phase_diff_file
+
+
+def compute_offset(echo_ind: int, W: npt.NDArray, X: npt.NDArray, Y: npt.NDArray) -> int:
+    """Method for computing the global mode offset for echoes > 1.
+
+    TODO: Rename this and calculate_offset.
+
+    Parameters
+    ----------
+    echo_ind : int
+        Echo index
+    W : :obj:`numpy.ndarray`
+        Weights
+    X : :obj:`numpy.ndarray`
+        TEs in 2d matrix
+    Y : :obj:`numpy.ndarray`
+        Masked unwrapped data weighted by magnitude
+
+    Returns
+    -------
+    best_offset : int
+    """
+    # fit the model to the up to previous echo
+    coefficients, _ = weighted_regression(X[:echo_ind], Y[:echo_ind], W[:echo_ind])
+
+    # compute the predicted phase for the current echo
+    Y_pred = X[echo_ind] * coefficients
+
+    # compute the difference between the predicted phase and the unwrapped phase
+    Y_diff = Y_pred - Y[echo_ind]
+
+    # compute closest multiple of 2pi to the difference
+    int_map = np.round(Y_diff / (2 * np.pi)).astype(int)
+
+    # compute the most often occuring multiple
+    best_offset = mode(int_map, axis=0, keepdims=False).mode
+    best_offset = cast(int, best_offset)
+
+    return best_offset
+
+
+def weighted_regression(
+    X: npt.NDArray,
+    Y: npt.NDArray,
+    W: npt.NDArray,
+) -> Tuple[npt.NDArray, npt.NDArray]:
+    """Single parameter weighted regression.
+
+    This function does a single parameter weighted regression using elementwise matrix operations.
+
+    Parameters
+    ----------
+    X : :obj:`numpy.ndarray`
+        Design matrix, rows are each echo, columns are voxels.
+        If column is size 1, this array will be broadcasted across all voxels.
+    Y : :obj:`numpy.ndarray`
+        Ordinate matrix, rows are each echo, columns are voxels.
+        (This is meant for phase values)
+    W : :obj:`numpy.ndarray`
+        Weight matrix (usually the magnitude image) echos x voxels
+
+    Returns
+    -------
+    :obj:`numpy.ndarray`
+        model weights
+    :obj:`numpy.ndarray`
+        residuals
+    """
+    # compute weighted X and Y
+    WY = W * Y
+    WX = W * X
+
+    # now compute the weighted squared magnitude of the X
+    sq_WX = np.sum(WX**2, axis=0, keepdims=True)
+
+    # get the inverse X
+    inv_WX = np.zeros(WX.shape)
+    np.divide(WX, sq_WX, out=inv_WX, where=(sq_WX != 0))
+
+    # compute the weights
+    model_weights = (inv_WX * WY).sum(axis=0)
+
+    # now compute residuals (on unweighted data)
+    residuals = np.sum((Y - model_weights * X) ** 2, axis=0)
+
+    # return model weights and residuals
+    return model_weights, residuals
+
+
+def calculate_diffs2(magnitude, phase):
+    """Calculate the magnitude and phase differences between two complex-valued images.
+
+    Parameters
+    ----------
+    magnitude : :obj:`str`
+        The path to the magnitude image (concatenated across echoes).
+    phase : :obj:`str`
+        The path to the phase image (concatenated across echoes).
+
+    Returns
+    -------
+    mag_diff_file : :obj:`str`
+        The path to the magnitude difference image.
+    phase_diff_file : :obj:`str`
+        The path to the phase difference image.
+    """
+    import os
+
+    import nibabel as nb
+    import numpy as np
+
+    mag_diff_file = os.path.abspath("magnitude_diff.nii.gz")
+    phase_diff_file = os.path.abspath("phase_diff.nii.gz")
+
+    magnitude_img = nb.load(magnitude)
+    phase_img = nb.load(phase)
+    magnitude_data = magnitude_img.get_fdata()
+    phase_data = phase_img.get_fdata()
+
+    signal_diff = (
+        magnitude_data[..., 0]
+        * magnitude_data[..., 1]
+        * np.exp(1j * (phase_data[..., 1] - phase_data[..., 0]))
+    )
+    mag_diff = np.abs(signal_diff)
+    phase_diff = np.angle(signal_diff)
+    mag_diff_img = nb.Nifti1Image(mag_diff, magnitude_img.affine, magnitude_img.header)
+    phase_diff_img = nb.Nifti1Image(phase_diff, phase_img.affine, phase_img.header)
+    mag_diff_img.to_filename(mag_diff_file)
+    phase_diff_img.to_filename(phase_diff_file)
+
+    return mag_diff_file, phase_diff_file
+
+
+def create_brain_mask(magnitude, extra_dilation=0):
+    """Create a quick brain mask for a single frame.
+
+    Parameters
+    ----------
+    magnitude : str
+        Path to magnitude file, where last dimension is the number of echoes.
+    extra_dilation : int
+        Number of extra dilations (or erosions if negative) to perform, by default 0
+
+    Returns
+    -------
+    mask_file : str
+        Mask of voxels to use for unwrapping
+    """
+    import os
+
+    import nibabel as nb
+    import numpy as np
+    from scipy.ndimage import (
+        binary_dilation,
+        binary_fill_holes,
+        binary_erosion,
+        generate_binary_structure,
+    )
+    from skimage.filters import threshold_otsu
+
+    from sdcflows.utils.misc import get_largest_connected_component
+
+    mag_img = nb.load(magnitude)
+    mag_data = mag_img.get_fdata()
+    mag_data = mag_data[..., 0]
+
+    mask_file = os.path.abspath("mask.nii.gz")
+
+    # create structuring element
+    strel = generate_binary_structure(3, 2)
+
+    # get the otsu threshold
+    threshold = threshold_otsu(mag_data)
+    mask_data = mag_data > threshold
+    if mask_data.ndim != strel.ndim:
+        raise ValueError(f"{mask_data.shape} ({mask_data.ndim}) != {strel.shape} ({strel.ndim})")
+    mask_data = binary_fill_holes(mask_data, strel).astype(np.float32)
+
+    # erode mask
+    mask_data = binary_erosion(
+        mask_data,
+        structure=strel,
+        iterations=2,
+        border_value=1,
+    ).astype(bool)
+
+    # get largest connected component
+    mask_data = get_largest_connected_component(mask_data)
+
+    # dilate the mask
+    mask_data = binary_dilation(mask_data, structure=strel, iterations=2)
+
+    # extra dilation to get areas on the edge of the brain
+    if extra_dilation > 0:
+        mask_data = binary_dilation(mask_data, structure=strel, iterations=extra_dilation)
+    # if negative, erode instead
+    if extra_dilation < 0:
+        mask_data = binary_erosion(mask_data, structure=strel, iterations=abs(extra_dilation))
+
+    # since we can't have a completely empty mask, set all zeros to ones
+    # if the mask is all empty
+    if np.all(np.isclose(mask_data, 0)):
+        mask_data = np.ones(mask_data.shape)
+
+    mask_data = mask_data.astype(np.bool_).astype(np.uint8)
+    mask_img = nb.Nifti1Image(mask_data, mag_img.affine, mag_img.header)
+    mask_img.to_filename(mask_file)
+    return mask_file
+
+
+def calculate_offset(phase, unwrapped_diff, echo_times):
+    """Calculate the phase offset between two echoes.
+
+    TODO: Rename this and compute_offset.
+
+    Parameters
+    ----------
+    phase : :obj:`str`
+        The path to the phase image (concatenated across echoes).
+    unwrapped_diff : :obj:`str`
+        The path to the unwrapped phase difference image from the first two echoes.
+    echo_times : :obj:`list` of :obj:`float`
+        The echo times.
+
+    Returns
+    -------
+    offset : :obj:`str`
+        The path to the phase offset image.
+    """
+    import os
+
+    import nibabel as nb
+    import numpy as np
+
+    offset_file = os.path.abspath("offset.nii.gz")
+
+    phase_img = nb.load(phase)
+    unwrapped_diff_img = nb.load(unwrapped_diff)
+    phase_data = phase_img.get_fdata()
+    unwrapped_diff_data = unwrapped_diff_img.get_fdata()
+
+    proposed_offset = np.angle(
+        np.exp(
+            1j
+            * (
+                phase_data[..., 0]
+                - ((echo_times[0] * unwrapped_diff_data) / (echo_times[1] - echo_times[0]))
+            )
+        )
+    )
+    proposed_offset_img = nb.Nifti1Image(proposed_offset, phase_img.affine, phase_img.header)
+    proposed_offset_img.to_filename(offset_file)
+    return offset_file
+
+
+def subtract_offset(phase, offset):
+    """Subtract the offset from the phase image.
+
+    Parameters
+    ----------
+    phase : :obj:`str`
+        The path to the phase image.
+    offset : :obj:`str`
+        The path to the offset image.
+
+    Returns
+    -------
+    updated_phase : :obj:`str`
+        The path to the updated phase image.
+    """
+    import os
+
+    import nibabel as nb
+    import numpy as np
+
+    updated_phase_file = os.path.abspath("updated_phase.nii.gz")
+
+    phase_img = nb.load(phase)
+    offset_img = nb.load(offset)
+    phase_data = phase_img.get_fdata()
+    offset_data = offset_img.get_fdata()
+    updated_phase = phase_data - offset_data[..., np.newaxis]
+    updated_phase_img = nb.Nifti1Image(updated_phase, phase_img.affine, phase_img.header)
+    updated_phase_img.to_filename(updated_phase_file)
+    return updated_phase_file
+
+
+def calculate_dual_echo_fieldmap(unwrapped_phase, echo_times):
+    """Calculate the fieldmap from the unwrapped phase difference of the first two echoes.
+
+    Parameters
+    ----------
+    unwrapped_phase : :obj:`str`
+        The path to the unwrapped phase difference image.
+    echo_times : :obj:`list` of :obj:`float`
+        The echo times.
+
+    Returns
+    -------
+    fieldmap : :obj:`str`
+        The path to the fieldmap image.
+    """
+    import os
+
+    import nibabel as nb
+    import numpy as np
+
+    fieldmap_file = os.path.abspath("fieldmap.nii.gz")
+
+    unwrapped_phase_img = nb.load(unwrapped_phase)
+    unwrapped_phase_data = unwrapped_phase_img.get_fdata()
+
+    phase_diff = unwrapped_phase_data[..., 1] - unwrapped_phase_data[..., 0]
+    fieldmap = (1000 / (2 * np.pi)) * phase_diff / (echo_times[1] - echo_times[0])
+    fieldmap_img = nb.Nifti1Image(
+        fieldmap,
+        unwrapped_phase_img.affine,
+        unwrapped_phase_img.header,
+    )
+    fieldmap_img.to_filename(fieldmap_file)
+
+    return fieldmap_file
+
+
+def modify_unwrapped_diff(phase, unwrapped_diff, echo_times):
+    """Add 2*pi to unwrapped difference data and recalculate the offset."""
+    import os
+
+    import nibabel as nb
+    import numpy as np
+
+    updated_phase_file = os.path.abspath("updated_phase.nii.gz")
+
+    phase_img = nb.load(phase)
+    unwrapped_diff_img = nb.load(unwrapped_diff)
+    phase_data = phase_img.get_fdata()
+    unwrapped_diff_data = unwrapped_diff_img.get_fdata()
+
+    new_proposed_offset = np.angle(
+        np.exp(
+            1j
+            * (
+                phase_data[..., 0]
+                - (
+                    (echo_times[0] * (unwrapped_diff_data + 2 * np.pi))
+                    / (echo_times[1] - echo_times[0])
+                )
+            )
+        )
+    )
+    new_proposed_phase = phase_data - new_proposed_offset[..., np.newaxis]
+    new_proposed_phase_img = nb.Nifti1Image(
+        new_proposed_phase,
+        phase_img.affine,
+        phase_img.header,
+    )
+    new_proposed_phase_img.to_filename(updated_phase_file)
+
+    return updated_phase_file
+
+
+def select_fieldmap(
+    original_fieldmap,
+    original_unwrapped_phase,
+    original_offset,
+    modified_fieldmap,
+    modified_unwrapped_phase,
+    unwrapped_diff,
+    voxel_mask,
+    echo_times,
+    wrap_limit,
+):
+    import os
+
+    import nibabel as nb
+    import numpy as np
+
+    FMAP_PROPORTION_HEURISTIC = 0.25
+    FMAP_AMBIGUIOUS_HEURISTIC = 0.5
+
+    new_unwrapped_diff_file = os.path.abspath("new_unwrapped_diff.nii.gz")
+
+    orig_fmap_img = nb.load(original_fieldmap)
+    orig_unwrapped_phase_img = nb.load(original_unwrapped_phase)
+    orig_offset_img = nb.load(original_offset)
+    mod_fmap_img = nb.load(modified_fieldmap)
+    mod_unwrapped_phase_img = nb.load(modified_unwrapped_phase)
+    unwrapped_diff_img = nb.load(unwrapped_diff)
+    voxel_mask_img = nb.load(voxel_mask)
+
+    orig_fmap_data = orig_fmap_img.get_fdata()
+    orig_unwrapped_phase_data = orig_unwrapped_phase_img.get_fdata()
+    orig_offset_data = orig_offset_img.get_fdata()
+    mod_fmap_data = mod_fmap_img.get_fdata()
+    mod_unwrapped_phase_data = mod_unwrapped_phase_img.get_fdata()
+    unwrapped_diff_data = unwrapped_diff_img.get_fdata()
+    voxel_mask_data = voxel_mask_img.get_fdata().astype(bool)
+
+    masked_orig_fmap = orig_fmap_data[voxel_mask_data]
+    if masked_orig_fmap.mean() < -10:
+        # check if the proposed fieldmap is below -10
+        unwrapped_diff_data += 2 * np.pi
+
+    elif masked_orig_fmap.mean() < 0 and not wrap_limit:
+        # check if the proposed fieldmap is between -10 and 0
+        # look at proportion of voxels that are positive
+        voxel_prop = np.count_nonzero(masked_orig_fmap > 0) / masked_orig_fmap.shape[0]
+
+        # if the proportion of positive voxels is less than 0.25, then add 2pi
+        if voxel_prop < FMAP_PROPORTION_HEURISTIC:
+            unwrapped_diff_data += 2 * np.pi
+        elif voxel_prop < FMAP_AMBIGUIOUS_HEURISTIC:
+            # compute mean of phase offset
+            mean_phase_offset = orig_offset_data[voxel_mask_data].mean()
+            if mean_phase_offset < -1:
+                phase_fits = np.concatenate(
+                    (
+                        np.zeros((*orig_unwrapped_phase_data.shape[:-1], 1)),
+                        orig_unwrapped_phase_data,
+                    ),
+                    axis=-1,
+                )
+                _, residuals_1, _, _, _ = np.polyfit(
+                    echo_times,
+                    phase_fits[voxel_mask_data, :].T,
+                    1,
+                    full=True,
+                )
+
+                masked_mod_fmap = mod_fmap_data[voxel_mask_data]
+                # fit linear model to the proposed phase
+                new_phase_fits = np.concatenate(
+                    (
+                        np.zeros((*mod_unwrapped_phase_data.shape[:-1], 1)),
+                        mod_unwrapped_phase_data,
+                    ),
+                    axis=-1,
+                )
+                _, residuals_2, _, _, _ = np.polyfit(
+                    echo_times,
+                    new_phase_fits[voxel_mask_data, :].T,
+                    1,
+                    full=True,
+                )
+
+                if (
+                    np.isclose(residuals_1.mean(), residuals_2.mean(), atol=1e-3, rtol=1e-3)
+                    and masked_mod_fmap.mean() > 0
+                ):
+                    unwrapped_diff_data += 2 * np.pi
+                else:
+                    unwrapped_diff_data -= 2 * np.pi
+
+    new_unwrapped_diff_img = nb.Nifti1Image(
+        unwrapped_diff_data,
+        unwrapped_diff_img.affine,
+        unwrapped_diff_img.header,
+    )
+    new_unwrapped_diff_img.to_filename(new_unwrapped_diff_file)
+    return new_unwrapped_diff_file
+
+
+def global_mode_correction(magnitude, unwrapped, mask, echo_times):
+    """Compute the global mode offset for the first echo and apply it to the subsequent echoes.
+
+    Parameters
+    ----------
+    magnitude : :obj:`str`
+    unwrapped : :obj:`str`
+    mask : :obj:`str`
+    echo_times : :obj:`list` of :obj:`float`
+
+    Returns
+    -------
+    new_unwrapped_file : :obj:`str`
+
+    Notes
+    -----
+    This computes the global mode offset for the first echo then tries to find the offset
+    that minimizes the residuals for each subsequent echo
+    """
+    import os
+
+    import nibabel as nb
+    import numpy as np
+
+    from sdcflows.utils.misc import compute_offset, create_brain_mask
+
+    new_unwrapped_file = os.path.abspath("unwrapped.nii.gz")
+
+    mag_img = nb.load(magnitude)
+    unwrapped_img = nb.load(unwrapped)
+    mask_img = nb.load(mask)
+    mag_data = mag_img.get_fdata()
+    unwrapped_data = unwrapped_img.get_fdata()
+    mask_data = mask_img.get_fdata().astype(bool)
+    echo_times = np.array(echo_times)
+
+    brain_mask = create_brain_mask(magnitude)
+    brain_mask_img = nb.load(brain_mask)
+    brain_mask_data = brain_mask_img.get_fdata().astype(bool)
+
+    # for each of these matrices TEs are on rows, voxels are columns
+    # get design matrix
+    X = echo_times[:, np.newaxis]
+
+    # get magnitude weight matrix
+    W = mag_data[brain_mask_data, :].T
+
+    # loop over each index past 1st echo
+    for i in range(1, echo_times.shape[0]):
+        # get matrix with the masked unwrapped data (weighted by magnitude)
+        Y = unwrapped_data[brain_mask_data, :].T
+
+        # Compute offset through linear regression method
+        best_offset = compute_offset(i, W, X, Y)
+
+        # apply the offset
+        unwrapped_data[..., i] += 2 * np.pi * best_offset
+
+    # set anything outside of mask_data to 0
+    unwrapped_data[~mask_data] = 0
+
+    new_unwrapped_img = nb.Nifti1Image(unwrapped_data, mag_img.affine, mag_img.header)
+    new_unwrapped_img.to_filename(new_unwrapped_file)
+
+    return new_unwrapped_file
+
+
+def corr2_coeff(A, B):
+    """Efficiently calculates correlation coefficient between the columns of two 2D arrays
+
+    Parameters
+    ----------
+    A : :obj:`numpy.ndarray`
+        1st array to correlate
+    B : :obj:`numpy.ndarray`
+        2nd array to correlate
+
+    Returns
+    -------
+    :obj:`numpy.ndarray`
+        array of correlation coefficients
+    """
+    # Transpose A and B
+    A = A.T
+    B = B.T
+
+    # Rowwise mean of input arrays & subtract from input arrays themeselves
+    A_mA = A - A.mean(axis=1, keepdims=True)
+    B_mB = B - B.mean(axis=1, keepdims=True)
+
+    # Sum of squares across rows
+    ssA = (A_mA**2).sum(axis=1, keepdims=True)
+    ssB = (B_mB**2).sum(axis=1, keepdims=True).T
+
+    # Finally get corr coeff
+    return np.dot(A_mA, B_mB.T) / np.sqrt(np.dot(ssA, ssB))
diff --git a/sdcflows/utils/phasemanip.py b/sdcflows/utils/phasemanip.py
index 0cc77d75cd..d5538011aa 100644
--- a/sdcflows/utils/phasemanip.py
+++ b/sdcflows/utils/phasemanip.py
@@ -46,6 +46,30 @@ 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)."""
+    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()))
+    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
diff --git a/sdcflows/workflows/fit/medic.py b/sdcflows/workflows/fit/medic.py
new file mode 100644
index 0000000000..6b89a83c4f
--- /dev/null
+++ b/sdcflows/workflows/fit/medic.py
@@ -0,0 +1,729 @@
+# emacs: -*- mode: python; py-indent-offset: 4; indent-tabs-mode: nil -*-
+# vi: set ft=python sts=4 ts=4 sw=4 et:
+#
+# Copyright 2021 The NiPreps Developers <nipreps@gmail.com>
+#
+# Licensed under the Apache License, Version 2.0 (the "License");
+# you may not use this file except in compliance with the License.
+# You may obtain a copy of the License at
+#
+#     http://www.apache.org/licenses/LICENSE-2.0
+#
+# Unless required by applicable law or agreed to in writing, software
+# distributed under the License is distributed on an "AS IS" BASIS,
+# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+# See the License for the specific language governing permissions and
+# limitations under the License.
+#
+# We support and encourage derived works from this project, please read
+# about our expectations at
+#
+#     https://www.nipreps.org/community/licensing/
+#
+"""Processing of dynamic field maps from complex-valued multi-echo BOLD data."""
+
+from nipype.pipeline import engine as pe
+from nipype.interfaces import utility as niu
+from niworkflows.engine.workflows import LiterateWorkflow as Workflow
+from niworkflows.interfaces.nibabel import MergeSeries
+
+from sdcflows.interfaces.fmap import MEDICB0, PhaseMap2rads2, ROMEO
+from sdcflows.interfaces.utils import EnforceTemporalConsistency, SVDFilter, TransposeImages
+from sdcflows.utils.misc import (
+    calculate_diffs2,
+    calculate_dual_echo_fieldmap,
+    calculate_offset,
+    create_brain_mask,
+    global_mode_correction,
+    medic_automask,
+    modify_unwrapped_diff,
+    select_fieldmap,
+    subtract_offset,
+)
+
+INPUT_FIELDS = ("magnitude", "phase")
+
+
+def init_medic_wf(
+    n_volumes,
+    echo_times,
+    automask,
+    automask_dilation,
+    border_filter=[1, 5],
+    svd_filter=10,
+    name="medic_wf",
+):
+    """Create the MEDIC dynamic field estimation workflow.
+
+    Workflow Graph
+        .. workflow ::
+            :graph2use: orig
+            :simple_form: yes
+
+            from sdcflows.workflows.fit.medic import init_medic_wf
+
+            wf = init_medic_wf(
+                n_volumes=4,
+                echo_times=[0.015, 0.030, 0.045, 0.06],
+                automask=True,
+                automask_dilation=3,
+                border_filter=[1, 5],
+                svd_filter=10,
+            )   # doctest: +SKIP
+
+    Parameters
+    ----------
+    n_volumes : :obj:`int`
+        The number of volumes in the multi-echo data.
+    echo_times : :obj:`list` of :obj:`float`
+        The echo times of the multi-echo data, in milliseconds.
+    automask : :obj:`bool`
+        Whether to automatically generate a mask for the fieldmap.
+    automask_dilation : :obj:`int`
+        The number of voxels by which to dilate the automatically generated mask.
+    border_filter : :obj:`list` of :obj:`int`
+        The border filter to apply to the fieldmap.
+    svd_filter : :obj:`int`
+        The number of singular values to retain in the SVD filter.
+    name : :obj:`str`
+        Name for this workflow
+
+    Inputs
+    ------
+    magnitude : :obj:`list` of :obj:`str`
+        A list of echo-wise magnitude EPI files that will be fed into MEDIC.
+    phase : :obj:`list` of :obj:`str`
+        A list of echo-wise phase EPI files that will be fed into MEDIC.
+
+    Outputs
+    -------
+    fieldmap : :obj:`str`
+        The path of the estimated fieldmap time series file.
+    phase_unwrapped : :obj:`list` of :obj:`str`
+        Path to the unwrapped phase time series files.
+    method: :obj:`str`
+        Short description of the estimation method that was run.
+
+    Notes
+    -----
+    This is a translation of the MEDIC algorithm, as implemented in ``vandandrew/warpkit``
+    (specifically the function ``unwrap_and_compute_field_maps``), into a Nipype workflow.
+
+    """
+    workflow = Workflow(name=name)
+
+    workflow.__desc__ = """\
+A dynamic fieldmap was estimated from multi-echo EPI data using the MEDIC algorithm (@medic).
+"""
+
+    inputnode = pe.Node(niu.IdentityInterface(fields=INPUT_FIELDS), name="inputnode")
+    outputnode = pe.Node(
+        niu.IdentityInterface(fields=["fieldmap", "phase_unwrapped", "method"]),
+        name="outputnode",
+    )
+    outputnode.inputs.method = "MEDIC"
+
+    n_echoes = len(echo_times)
+
+    # Convert phase to radians (-pi to pi, not 0 to 2pi)
+    phase2rad = pe.MapNode(
+        PhaseMap2rads2(),
+        iterfield=["in_file"],
+        name="phase2rad",
+    )
+    workflow.connect([(inputnode, phase2rad, [("phase", "in_file")])])
+
+    # Transpose the magnitude data into volume-wise images concatenated across echoes
+    transpose_magnitude = pe.Node(
+        TransposeImages(
+            numinputs=n_echoes,
+            numoutputs=n_volumes,
+        ),
+        name="transpose_magnitude",
+    )
+
+    # Transpose the phase data into volume-wise images concatenated across echoes
+    transpose_phase = pe.Node(
+        TransposeImages(
+            numinputs=n_echoes,
+            numoutputs=n_volumes,
+        ),
+        name="transpose_phase",
+    )
+
+    for i_echo in range(n_echoes):
+        # Split magnitude data into single-frame, multi-echo 4D files
+        select_mag_echo = pe.Node(
+            niu.Select(index=i_echo),
+            name=f"select_mag_echo_{i_echo:02d}",
+        )
+        workflow.connect([
+            (inputnode, select_mag_echo, [("magnitude", "inlist")]),
+            (select_mag_echo, transpose_magnitude, [("out", f"in{i_echo + 1}")]),
+        ])  # fmt:skip
+
+        select_phase_echo = pe.Node(
+            niu.Select(index=i_echo),
+            name=f"select_phase_echo_{i_echo:02d}",
+        )
+        workflow.connect([
+            (phase2rad, select_phase_echo, [("out_file", "inlist")]),
+            (select_phase_echo, transpose_phase, [("out", f"in{i_echo + 1}")]),
+        ])  # fmt:skip
+
+    # Prepare to recombine into echo-wise time series
+    transpose_phase_unwrapped = pe.Node(
+        TransposeImages(
+            numinputs=n_volumes,
+            numoutputs=n_echoes,
+        ),
+        name="transpose_phase_unwrapped",
+    )
+    aggregate_masks = pe.Node(
+        niu.Merge(n_volumes),
+        name="aggregate_masks",
+    )
+
+    for i_volume in range(n_volumes):
+        process_volume_wf = init_process_volume_wf(
+            echo_times=echo_times,
+            automask=automask,
+            automask_dilation=automask_dilation,
+            name=f"process_volume_{i_volume:02d}_wf",
+        )
+
+        workflow.connect([
+            (transpose_magnitude, process_volume_wf, [
+                (f"out{i_volume + 1}", "inputnode.magnitude"),
+            ]),
+            (transpose_phase, process_volume_wf, [(f"out{i_volume + 1}", "inputnode.phase")]),
+            (process_volume_wf, transpose_phase_unwrapped, [
+                ("outputnode.phase_unwrapped", f"in{i_volume + 1}"),
+            ]),
+            (process_volume_wf, aggregate_masks, [("outputnode.mask", f"in{i_volume + 1}")]),
+        ])  # fmt:skip
+
+    concatenate_masks = pe.Node(
+        MergeSeries(allow_4D=False),
+        name="concatenate_masks",
+    )
+    workflow.connect([(aggregate_masks, concatenate_masks, [("out", "in_files")])])
+
+    # Compute field maps
+    merge_phase_unwrapped = pe.Node(
+        niu.Merge(n_echoes),
+        name="merge_phase_unwrapped",
+    )
+    for i_echo in range(n_echoes):
+        workflow.connect([
+            (transpose_phase_unwrapped, merge_phase_unwrapped, [
+                (f"out{i_echo + 1}", f"in{i_echo + 1}"),
+            ]),
+        ])  # fmt:skip
+
+    # Enforce temporal consistency of phase unwrapping
+    enforce_consistency = pe.Node(
+        EnforceTemporalConsistency(echo_times=echo_times, threshold=0.98),
+        name="enforce_consistency",
+    )
+    workflow.connect([
+        (inputnode, enforce_consistency, [("magnitude", "magnitude")]),
+        (concatenate_masks, enforce_consistency, [("out_file", "mask")]),
+        (merge_phase_unwrapped, enforce_consistency, [("out", "phase_unwrapped")]),
+        (enforce_consistency, outputnode, [("phase_unwrapped", "phase_unwrapped")]),
+    ])  # fmt:skip
+
+    compute_fieldmap = pe.Node(
+        MEDICB0(echo_times=echo_times),
+        name="compute_fieldmap",
+    )
+    workflow.connect([
+        (inputnode, compute_fieldmap, [("magnitude", "magnitude")]),
+        (enforce_consistency, compute_fieldmap, [("phase_unwrapped", "phase")]),
+    ])  # fmt:skip
+
+    # Apply SVD filter to field maps
+    apply_svd_filter = pe.Node(
+        SVDFilter(
+            svd_filter=svd_filter,
+            border_filter=list(border_filter),
+        ),
+        name="apply_svd_filter",
+    )
+    workflow.connect([
+        (compute_fieldmap, apply_svd_filter, [("b0", "fieldmap")]),
+        (concatenate_masks, apply_svd_filter, [("out_file", "mask")]),
+        (apply_svd_filter, outputnode, [("fieldmap", "fieldmap")]),
+    ])  # fmt:skip
+
+    return workflow
+
+
+def init_process_volume_wf(
+    echo_times,
+    automask=True,
+    automask_dilation=3,
+    name="process_volume_wf",
+):
+    """Process a single volume of multi-echo data according to the MEDIC method.
+
+    Workflow Graph
+        .. workflow ::
+            :graph2use: orig
+            :simple_form: yes
+
+            from sdcflows.workflows.fit.medic import init_process_volume_wf
+
+            wf = init_process_volume_wf(
+                echo_times=[0.015, 0.030, 0.045, 0.06],
+                automask=True,
+                automask_dilation=3,
+            )   # doctest: +SKIP
+
+    Parameters
+    ----------
+    echo_times : :obj:`list` of :obj:`float`
+        The echo times of the multi-echo data.
+    automask : :obj:`bool`
+        Whether to automatically generate a mask for the fieldmap.
+    automask_dilation : :obj:`int`
+        The number of voxels by which to dilate the automatically generated mask.
+    name : :obj:`str`
+        The name of the workflow.
+
+    Inputs
+    ------
+    magnitude : :obj:`str`
+        One volume of magnitude EPI data, concatenated across echoes.
+    phase : :obj:`str`
+        One volume of phase EPI data, concatenated across echoes.
+        Must already be scaled from -pi to pi.
+    mask : :obj:`str`
+        The brain mask that will be used to constrain the fieldmap estimation.
+        If ``automask`` is True, this mask will be modified.
+        Otherwise, it will be returned in the outputnode unmodified.
+
+    Outputs
+    -------
+    phase_unwrapped : :obj:`str`
+        Unwrapped phase in radians.
+    mask : :obj:`str`
+        Path to a brain mask that can be used to constrain the fieldmap estimation.
+    """
+    workflow = Workflow(name=name)
+
+    inputnode = pe.Node(
+        niu.IdentityInterface(fields=["magnitude", "phase", "mask"]),
+        name="inputnode",
+    )
+    outputnode = pe.Node(
+        niu.IdentityInterface(fields=["phase_unwrapped", "mask"]),
+        name="outputnode",
+    )
+
+    mask_buffer = pe.Node(
+        niu.IdentityInterface(fields=["mask"]),
+        name="mask_buffer",
+    )
+    if automask:
+        # the theory goes like this, the magnitude/otsu base mask can be too aggressive
+        # occasionally and the voxel quality mask can get extra voxels that are not brain,
+        # but is noisy so we combine the two masks to get a better mask
+
+        # Use ROMEO's voxel-quality command
+        # XXX: In warpkit Andrew creates a "mag" image of all ones.
+        # With the current version (at least on some test data),
+        # there are NaNs in the voxel quality map.
+        voxqual = pe.Node(
+            ROMEO(write_quality=True, echo_times=echo_times, mask="robustmask"),
+            name="voxqual",
+        )
+        workflow.connect([
+            (inputnode, voxqual, [
+                ("magnitude", "mag_file"),
+                ("phase", "phase_file"),
+            ]),
+        ])  # fmt:skip
+
+        # Then use skimage's otsu thresholding to get a mask and do a bunch of other stuff
+        automask_medic = pe.Node(
+            niu.Function(
+                input_names=["mag_file", "voxel_quality", "echo_times", "automask_dilation"],
+                output_names=["mask", "masksum_file"],
+                function=medic_automask,
+            ),
+            name="automask_medic",
+        )
+        automask_medic.inputs.echo_times = echo_times
+        automask_medic.inputs.automask_dilation = automask_dilation
+        workflow.connect([
+            (inputnode, automask_medic, [("magnitude", "mag_file")]),
+            (voxqual, automask_medic, [("quality_file", "voxel_quality")]),
+            (automask_medic, mask_buffer, [("mask", "mask")]),
+        ])  # fmt:skip
+    else:
+        workflow.connect([(inputnode, mask_buffer, [("mask", "mask")])])
+
+    workflow.connect([(mask_buffer, outputnode, [("mask", "mask")])])
+
+    # Do MCPC-3D-S algo to compute phase offset
+    mcpc_3d_s_wf = init_mcpc_3d_s_wf(wrap_limit=False, name="mcpc_3d_s_wf")
+    mcpc_3d_s_wf.inputs.inputnode.echo_times = echo_times
+    workflow.connect([
+        (inputnode, mcpc_3d_s_wf, [
+            ("magnitude", "inputnode.magnitude"),
+            ("phase", "inputnode.phase"),
+        ]),
+        (mask_buffer, mcpc_3d_s_wf, [("mask", "inputnode.mask")]),
+    ])  # fmt:skip
+
+    # remove offset from phase data
+    remove_offset = pe.Node(
+        niu.Function(
+            input_names=["phase", "offset"],
+            output_names=["phase_modified"],
+            function=subtract_offset,
+        ),
+        name="remove_offset",
+    )
+    workflow.connect([
+        (inputnode, remove_offset, [("phase", "phase")]),
+        (mcpc_3d_s_wf, remove_offset, [("outputnode.offset", "offset")]),
+    ])  # fmt:skip
+
+    # Unwrap the modified phase data with ROMEO
+    unwrap_phase = pe.Node(
+        ROMEO(
+            echo_times=echo_times,
+            weights="romeo",
+            correct_global=True,
+            max_seeds=1,
+            merge_regions=False,
+            correct_regions=False,
+        ),
+        name="unwrap_phase",
+    )
+    workflow.connect([
+        (inputnode, unwrap_phase, [("magnitude", "mag_file")]),
+        (mask_buffer, unwrap_phase, [("mask", "mask")]),
+        (remove_offset, unwrap_phase, [("phase_modified", "phase_file")]),
+    ])  # fmt:skip
+
+    # Global mode correction
+    global_mode_corr = pe.Node(
+        niu.Function(
+            input_names=["magnitude", "unwrapped", "mask", "echo_times"],
+            output_names=["unwrapped"],
+            function=global_mode_correction,
+        ),
+        name="global_mode_corr",
+    )
+    global_mode_corr.inputs.echo_times = echo_times
+    workflow.connect([
+        (inputnode, global_mode_corr, [("magnitude", "magnitude")]),
+        (unwrap_phase, global_mode_corr, [("out_file", "unwrapped")]),
+        (mask_buffer, global_mode_corr, [("mask", "mask")]),
+        (global_mode_corr, outputnode, [("unwrapped", "phase_unwrapped")]),
+    ])  # fmt:skip
+
+    return workflow
+
+
+def init_mcpc_3d_s_wf(wrap_limit, name):
+    """Estimate and remove phase offset with MCPC-3D-S algorithm.
+
+    Parameters
+    ----------
+    wrap_limit : bool
+        If True, this turns off some heuristics for phase unwrapping.
+    name : str
+        The name of the workflow.
+
+    Inputs
+    ------
+    magnitude : str
+        The path to the magnitude image. A single volume, concatenated across echoes.
+    phase : str
+        The path to the phase image. A single volume, concatenated across echoes.
+    echo_times : list of float
+        The echo times of the multi-echo data.
+    mask : str
+        The path to the brain mask mask.
+
+    Outputs
+    -------
+    offset : str
+        The path to the estimated phase offset.
+    unwrapped_diff : str
+        The path to the unwrapped phase difference image.
+    """
+    workflow = Workflow(name=name)
+
+    inputnode = pe.Node(
+        niu.IdentityInterface(
+            fields=[
+                "magnitude",
+                "phase",
+                "echo_times",
+                "mask",
+            ],
+        ),
+        name="inputnode",
+    )
+    outputnode = pe.Node(
+        niu.IdentityInterface(fields=["offset", "unwrapped_diff"]),
+        name="outputnode",
+    )
+
+    # Calculate magnitude and phase differences from first two echoes
+    calc_diffs = pe.Node(
+        niu.Function(
+            input_names=["magnitude", "phase"],
+            output_names=["mag_diff_file", "phase_diff_file"],
+            function=calculate_diffs2,
+        ),
+        name="calc_diffs",
+    )
+    workflow.connect([
+        (inputnode, calc_diffs, [
+            ("magnitude", "magnitude"),
+            ("phase", "phase"),
+        ]),
+    ])  # fmt:skip
+
+    # Unwrap difference images
+    unwrap_diffs = pe.Node(
+        ROMEO(
+            weights="romeo",
+            correct_global=True,
+        ),
+        name="unwrap_diffs",
+    )
+    workflow.connect([
+        (inputnode, unwrap_diffs, [("mask", "mask")]),
+        (calc_diffs, unwrap_diffs, [
+            ("mag_diff_file", "mag_file"),
+            ("phase_diff_file", "phase_file"),
+        ]),
+    ])  # fmt:skip
+
+    # Calculate voxel mask
+    create_mask = pe.Node(
+        niu.Function(
+            input_names=["magnitude", "extra_dilation"],
+            output_names=["mask"],
+            function=create_brain_mask,
+        ),
+        name="create_mask",
+    )
+    create_mask.inputs.extra_dilation = -2  # hardcoded in warpkit
+    workflow.connect([(inputnode, create_mask, [("magnitude", "magnitude")])])
+
+    # Calculate initial offset estimate
+    calc_offset = pe.Node(
+        niu.Function(
+            input_names=["phase", "unwrapped_diff", "echo_times"],
+            output_names=["offset"],
+            function=calculate_offset,
+        ),
+        name="calc_offset",
+    )
+    workflow.connect([
+        (inputnode, calc_offset, [
+            ("phase", "phase"),
+            ("echo_times", "echo_times"),
+        ]),
+        (unwrap_diffs, calc_offset, [("out_file", "unwrapped_diff")]),
+    ])  # fmt:skip
+
+    # Get the new phase
+    calc_proposed_phase = pe.Node(
+        niu.Function(
+            input_names=["phase", "offset"],
+            output_names=["proposed_phase"],
+            function=subtract_offset,
+        ),
+        name="calc_proposed_phase",
+    )
+    workflow.connect([
+        (inputnode, calc_proposed_phase, [("phase", "phase")]),
+        (calc_offset, calc_proposed_phase, [("offset", "offset")]),
+    ])  # fmt:skip
+
+    # Compute the dual-echo field map
+    dual_echo_wf = init_dual_echo_wf(name="dual_echo_wf")
+    workflow.connect([
+        (inputnode, dual_echo_wf, [
+            ("mask", "inputnode.mask"),
+            ("echo_times", "inputnode.echo_times"),
+            ("magnitude", "inputnode.magnitude"),
+        ]),
+        (calc_proposed_phase, dual_echo_wf, [("proposed_phase", "inputnode.phase")]),
+    ])  # fmt:skip
+
+    # Calculate a modified field map with 2pi added to the unwrapped difference image
+    add_2pi = pe.Node(
+        niu.Function(
+            input_names=["phase", "unwrapped_diff", "echo_times"],
+            output_names=["phase"],
+            function=modify_unwrapped_diff,
+        ),
+        name="add_2pi",
+    )
+    workflow.connect([
+        (inputnode, add_2pi, [
+            ("phase", "phase"),
+            ("echo_times", "echo_times"),
+        ]),
+        (unwrap_diffs, add_2pi, [("out_file", "unwrapped_diff")]),
+    ])  # fmt:skip
+
+    modified_dual_echo_wf = init_dual_echo_wf(name="modified_dual_echo_wf")
+    workflow.connect([
+        (inputnode, modified_dual_echo_wf, [
+            ("mask", "inputnode.mask"),
+            ("echo_times", "inputnode.echo_times"),
+            ("magnitude", "inputnode.magnitude"),
+        ]),
+        (add_2pi, modified_dual_echo_wf, [("phase", "inputnode.phase")]),
+    ])  # fmt:skip
+
+    # Select the fieldmap
+    select_unwrapped_diff = pe.Node(
+        niu.Function(
+            input_names=[
+                "original_fieldmap",
+                "original_unwrapped_phase",
+                "original_offset",
+                "modified_fieldmap",
+                "modified_unwrapped_phase",
+                "unwrapped_diff",
+                "voxel_mask",
+                "echo_times",
+                "wrap_limit",
+            ],
+            output_names=["new_unwrapped_diff"],
+            function=select_fieldmap,
+        ),
+        name="select_unwrapped_diff",
+    )
+    select_unwrapped_diff.inputs.wrap_limit = wrap_limit
+    workflow.connect([
+        (inputnode, select_unwrapped_diff, [("echo_times", "echo_times")]),
+        (unwrap_diffs, select_unwrapped_diff, [("out_file", "unwrapped_diff")]),
+        (create_mask, select_unwrapped_diff, [("mask", "voxel_mask")]),
+        (calc_offset, select_unwrapped_diff, [("offset", "original_offset")]),
+        (dual_echo_wf, select_unwrapped_diff, [
+            ("outputnode.fieldmap", "original_fieldmap"),
+            ("outputnode.unwrapped_phase", "original_unwrapped_phase"),
+        ]),
+        (modified_dual_echo_wf, select_unwrapped_diff, [
+            ("outputnode.fieldmap", "modified_fieldmap"),
+            ("outputnode.unwrapped_phase", "modified_unwrapped_phase"),
+        ]),
+        (select_unwrapped_diff, outputnode, [("new_unwrapped_diff", "unwrapped_diff")]),
+    ])  # fmt:skip
+
+    # Compute the updated phase offset
+    calc_updated_offset = pe.Node(
+        niu.Function(
+            input_names=["phase", "unwrapped_diff", "echo_times"],
+            output_names=["offset"],
+            function=calculate_offset,
+        ),
+        name="calc_updated_offset",
+    )
+    workflow.connect([
+        (inputnode, calc_updated_offset, [
+            ("phase", "phase"),
+            ("echo_times", "echo_times"),
+        ]),
+        (select_unwrapped_diff, calc_updated_offset, [("new_unwrapped_diff", "unwrapped_diff")]),
+        (calc_updated_offset, outputnode, [("offset", "offset")]),
+    ])  # fmt:skip
+
+    return workflow
+
+
+def init_dual_echo_wf(name="dual_echo_wf"):
+    """Estimate a field map from the first two echoes of multi-echo data.
+
+    Parameters
+    ----------
+    name : str
+        The name of the workflow.
+
+    Inputs
+    ------
+    magnitude : str
+        The path to the magnitude image.
+    phase : str
+        The path to the phase image.
+    echo_times : list of float
+        The echo times of the multi-echo data.
+    mask : str
+        The path to the brain mask mask.
+
+    Outputs
+    -------
+    fieldmap : str
+        The path to the estimated fieldmap.
+    unwrapped_phase : str
+        The path to the unwrapped phase image.
+    """
+    workflow = Workflow(name=name)
+
+    inputnode = pe.Node(
+        niu.IdentityInterface(
+            fields=[
+                "magnitude",
+                "phase",
+                "echo_times",
+                "mask",
+            ],
+        ),
+        name="inputnode",
+    )
+
+    outputnode = pe.Node(
+        niu.IdentityInterface(fields=["fieldmap", "unwrapped_phase"]),
+        name="outputnode",
+    )
+
+    # Unwrap the phase with ROMEO
+    unwrap_phase = pe.Node(
+        ROMEO(
+            weights="romeo",
+            correct_global=True,
+            max_seeds=1,
+            merge_regions=False,
+            correct_regions=False,
+        ),
+        name="unwrap_phase",
+    )
+    workflow.connect([
+        (inputnode, unwrap_phase, [
+            ("magnitude", "mag_file"),
+            ("phase", "phase_file"),
+            ("mask", "mask"),
+            ("echo_times", "echo_times"),
+        ]),
+        (unwrap_phase, outputnode, [("out_file", "unwrapped_phase")]),
+    ])  # fmt:skip
+
+    # Calculate the fieldmap
+    calc_fieldmap = pe.Node(
+        niu.Function(
+            input_names=["unwrapped_phase", "echo_times"],
+            output_names=["fieldmap"],
+            function=calculate_dual_echo_fieldmap,
+        ),
+        name="calc_fieldmap",
+    )
+    workflow.connect([
+        (unwrap_phase, calc_fieldmap, [("out_file", "unwrapped_phase")]),
+        (inputnode, calc_fieldmap, [("echo_times", "echo_times")]),
+        (calc_fieldmap, outputnode, [("fieldmap", "fieldmap")]),
+    ])  # fmt:skip
+
+    return workflow