[FEAT] MEDIC distortion correction via warpkit

.github/workflows/build-test-publish.yml

@@ -55,6 +55,15 @@ jobs:
       TEST_DATA_HOME: /home/runner/sdcflows-tests
 
     steps:
+    - name: Free disk space
+      uses: jlumbroso/free-disk-space@54081f138730dfa15788a46383842cd2f914a1be # v1.3.1
+      with:
+        tool-cache: false
+        android: true
+        dotnet: true
+        haskell: true
+        large-packages: false
+        swap-storage: false
     - uses: actions/checkout@3d3c42e5aac5ba805825da76410c181273ba90b1 # v7.0.1
     - name: Git settings (pacify DataLad)
       run: |
@@ -72,7 +81,7 @@ jobs:
     - uses: actions/cache@55cc8345863c7cc4c66a329aec7e433d2d1c52a9 # v6.1.0
       with:
         path: ${{ env.TEST_DATA_HOME }}
-        key: data-cache-v3
+        key: data-cache-v4
     - name: Install test data
       run: |
         mkdir -p ${{ env.TEST_DATA_HOME }}
@@ -115,6 +124,16 @@ jobs:
         datalad update -r --merge -d hcph-pilot_fieldmaps/
         datalad get -r -J 2 -d hcph-pilot_fieldmaps/ hcph-pilot_fieldmaps/*
 
+        # ds006926 — MEDIC multi-echo mag+phase BOLD (sub-a01 only)
+        datalad install -r https://github.com/OpenNeuroDatasets/ds006926.git
+        datalad update -r --merge -d ds006926/
+        datalad get -r -J 2 -d ds006926/ ds006926/sub-a01/func/sub-a01_task-VisMot_acq-tr1800_*
+
+        # ds007637 — MEDIC multi-echo mag+phase BOLD (sub-04/ses-2 fracback only)
+        datalad install -r https://github.com/OpenNeuroDatasets/ds007637.git
+        datalad update -r --merge -d ds007637/
+        datalad get -r -J 2 -d ds007637/ ds007637/sub-04/ses-2/func/sub-04_ses-2_task-fracback_acq-MBME_echo-*_part-mag_bold.nii.gz ds007637/sub-04/ses-2/func/sub-04_ses-2_task-fracback_acq-MBME_echo-*_part-phase_bold.nii.gz
+
   test:
     if: "!contains(github.event.head_commit.message, '[skip ci]' && (github.event_name == 'push' || github.event.pull_request.head.repo.full_name != 'nipreps/sdcflows'))"
     runs-on: ubuntu-latest
@@ -148,6 +167,15 @@ jobs:
             marks: "veryslow"
 
     steps:
+    - name: Free disk space
+      uses: jlumbroso/free-disk-space@54081f138730dfa15788a46383842cd2f914a1be # v1.3.1
+      with:
+        tool-cache: false
+        android: true
+        dotnet: true
+        haskell: true
+        large-packages: false
+        swap-storage: false
     - uses: actions/checkout@3d3c42e5aac5ba805825da76410c181273ba90b1 # v7.0.1
     - name: Install dependencies
       run: |
@@ -207,7 +235,7 @@ jobs:
     - uses: actions/cache/restore@55cc8345863c7cc4c66a329aec7e433d2d1c52a9 # v6.1.0
       with:
         path: ${{ env.TEST_DATA_HOME }}
-        key: data-cache-v3
+        key: data-cache-v4
 
     - name: Set FreeSurfer variables
       run: |

.zenodo.json

@@ -135,6 +135,12 @@
       "affiliation": "Department of Psychology, Stanford University, CA, USA",
       "name": "Russell A. Poldrack",
       "type": "Researcher"
+    },
+    {
+      "orcid": "0000-0002-8787-0943",
+      "affiliation": "Department of Biomedical Engineering, Washington University in St. Louis, MO, USA",
+      "name": "Andrew Van",
+      "type": "Researcher"
     }
   ],
   "keywords": [

pyproject.toml

@@ -34,6 +34,7 @@ dependencies = [
     "scipy >= 1.15",
     "templateflow >= 23.1",
     "toml >= 0.10",
+    "warpkit >= 1.5.0",
 ]
 dynamic = ["version"]
 

sdcflows/cli/main.py

@@ -70,6 +70,7 @@ def main(argv=None):
                 layout=config.execution.layout,
                 subject=subject,
                 fmapless=config.workflow.fmapless,
+                no_medic=config.workflow.no_medic,
                 logger=config.loggers.cli,
             )
 

sdcflows/cli/parser.py

@@ -247,6 +247,14 @@ def _bids_filter(value):
         default=True,
         help='Allow fieldmap-less estimation',
     )
+    g_outputs.add_argument(
+        '--no-medic',
+        action='store_true',
+        dest='no_medic',
+        default=False,
+        help='Disable MEDIC discovery (by default MEDIC takes priority for '
+        'complex multi-echo BOLD)',
+    )
     g_outputs.add_argument(
         '--use-plugin',
         action='store',

sdcflows/config.py

@@ -498,6 +498,8 @@ class workflow(_Config):
     """Level of analysis."""
     fmapless = False
     """Allow fieldmap-less estimation"""
+    no_medic = False
+    """Disable MEDIC discovery (otherwise MEDIC takes priority for complex multi-echo BOLD)"""
     species = 'human'
     """Subject species to choose most appropriate template"""
     template_id = 'MNI152NLin2009cAsym'

sdcflows/conftest.py

@@ -41,10 +41,17 @@
 test_workdir = os.getenv('TEST_WORK_DIR')
 _sloppy_mode = os.getenv('TEST_PRODUCTION', 'off').lower() not in ('on', '1', 'true', 'yes', 'y')
 
+# MEDIC fixtures live in full OpenNeuro trees (tens of thousands of JSON
+# sidecars) but only a few files are actually fetched via ``datalad get``.
+# Indexing those trees with ``BIDSLayout(derivatives=True)`` at collection
+# time stalled CI past the 20-minute tox watchdog. The MEDIC tests reach
+# their files via the ``datadir`` fixture directly, not via ``layouts``.
+_SKIP_LAYOUTS = {'ds006926', 'ds007637'}
+
 layouts = {
     p.name: BIDSLayout(str(p), validate=False, derivatives=True)
     for p in Path(test_data_env).glob('*')
-    if p.is_dir()
+    if p.is_dir() and p.name not in _SKIP_LAYOUTS
 }
 
 data_dir = Path(__file__).parent / 'tests' / 'data'
@@ -128,3 +135,52 @@ def dsA_dir():
 @pytest.fixture
 def sloppy_mode():
     return _sloppy_mode
+
+
+# MEDIC end-to-end fixtures, shared by the fit (``test_medic``) and apply
+# (``test_dynamic``) test modules. A handful of timepoints is enough to
+# exercise the full per-volume path; the source datasets ship 200+ volumes ×
+# 5 echoes × mag+phase, which OOM-kills CI runners when xdist schedules these
+# in parallel.
+_MEDIC_DATASETS = [
+    pytest.param(
+        (
+            'ds007637',
+            'sub-04/ses-2/func/sub-04_ses-2_task-fracback_acq-MBME_echo-*_part-mag_bold.nii.gz',
+        ),
+        id='ds007637',
+    ),
+    pytest.param(
+        ('ds006926', 'sub-a01/func/sub-a01_task-VisMot_acq-tr1800_echo-*_part-mag_bold.nii.gz'),
+        id='ds006926',
+    ),
+]
+
+
+@pytest.fixture
+def medic_test_volumes():
+    return 3
+
+
+@pytest.fixture(params=_MEDIC_DATASETS)
+def medic_fixture(request):
+    """Yield ``(dataset, mag_glob_under_dataset)`` for each MEDIC fixture."""
+    return request.param
+
+
+@pytest.fixture
+def truncate_to_volumes():
+    """Return a helper that slices 4D NIfTIs down to ``volumes`` timepoints."""
+
+    def _truncate(in_files, volumes, dest):
+        out = []
+        for f in in_files:
+            img = nibabel.load(str(f))
+            if img.shape[-1] > volumes:
+                img = img.slicer[..., :volumes]
+            new = dest / f.name
+            img.to_filename(new)
+            out.append(new)
+        return out
+
+    return _truncate

sdcflows/fieldmaps.py

@@ -55,6 +55,7 @@ class EstimatorType(Enum):
     PHASEDIFF = auto()
     MAPPED = auto()
     ANAT = auto()
+    MEDIC = auto()
 
 
 MODALITIES = {
@@ -75,6 +76,12 @@ class EstimatorType(Enum):
     'T2w': EstimatorType.ANAT,
 }
 
+# Estimator types that emit a per-volume 4D fieldmap on the EPI grid (and
+# therefore do not produce B-spline coefficients). Add new dynamic methods
+# here so consumers — ``init_fmap_preproc_wf`` in particular — pick them up
+# without per-method branching.
+_DYNAMIC_METHODS = frozenset({EstimatorType.MEDIC})
+
 
 def _type_setter(obj, attribute, value):
     """Make sure the type of estimation is not changed."""
@@ -88,6 +95,7 @@ def _type_setter(obj, attribute, value):
         EstimatorType.PHASEDIFF,
         EstimatorType.MAPPED,
         EstimatorType.ANAT,
+        EstimatorType.MEDIC,
     ):
         raise ValueError(f'Invalid estimation method type {value}.')
 
@@ -338,6 +346,59 @@ def __attrs_post_init__(self):
         suffix_list = [f.suffix for f in self.sources]
         suffix_set = set(suffix_list)
 
+        # Fieldmap option 0: MEDIC — multi-echo phase + magnitude
+        # ``bold`` / ``epi`` sources tagged with the BIDS ``part-{phase,mag}``
+        # entity. PEPOLAR uses ``dir-`` instead, so the part entity is the
+        # cleanest way to disambiguate.
+        parts = {f.entities.get('part') for f in self.sources}
+        medic_parts = parts & {'phase', 'mag'}
+        if suffix_set <= {'bold', 'epi', 'sbref'} and medic_parts:
+            # Any sources is ``part``-tagged: this is a MEDIC-shaped input.
+            # Reject incomplete sets explicitly rather than letting them slip
+            # through to the PEPOLAR branch and produce a confusing failure.
+            if parts != {'phase', 'mag'}:
+                raise ValueError(
+                    'MEDIC requires every source to be tagged ``part-mag`` or '
+                    '``part-phase``, with both present; got '
+                    f'parts={sorted(str(p) for p in parts)!r}.'
+                )
+            phase_files = [f for f in self.sources if f.entities.get('part') == 'phase']
+            mag_files = [f for f in self.sources if f.entities.get('part') == 'mag']
+            if len(phase_files) < 2:
+                raise ValueError(
+                    f'MEDIC requires at least two echoes of phase data; got {len(phase_files)}.'
+                )
+            if len(phase_files) != len(mag_files):
+                raise ValueError(
+                    f'MEDIC requires matched magnitude/phase pairs per echo; '
+                    f'got {len(phase_files)} phase and {len(mag_files)} '
+                    'magnitude file(s).'
+                )
+            # Counts alone would let phase echoes 1/2 pair with magnitude
+            # echoes 1/3: require one file per part for each echo, with
+            # matching EchoTime, so ``get_workflow`` can align them.
+            phase_echoes = {f.entities.get('echo'): f for f in phase_files}
+            mag_echoes = {f.entities.get('echo'): f for f in mag_files}
+            if len(phase_echoes) != len(phase_files) or len(mag_echoes) != len(mag_files):
+                raise ValueError(
+                    'MEDIC requires exactly one phase and one magnitude file per ``echo`` entity.'
+                )
+            if phase_echoes.keys() != mag_echoes.keys():
+                raise ValueError(
+                    'MEDIC requires matched magnitude/phase pairs per echo; got '
+                    f'phase echoes {sorted(map(str, phase_echoes))} and magnitude '
+                    f'echoes {sorted(map(str, mag_echoes))}.'
+                )
+            for echo, phase_file in phase_echoes.items():
+                phase_te = phase_file.metadata.get('EchoTime')
+                mag_te = mag_echoes[echo].metadata.get('EchoTime')
+                if phase_te is None or mag_te is None or abs(phase_te - mag_te) > 1e-6:
+                    raise ValueError(
+                        f'MEDIC echo {echo} requires matching EchoTime on its phase '
+                        f'and magnitude files; got {phase_te} and {mag_te}.'
+                    )
+            self.method = EstimatorType.MEDIC
+
         # Fieldmap option 1: actual field-mapping sequences
         fmap_types = suffix_set.intersection(('fieldmap', 'phasediff', 'phase1', 'phase2'))
         if len(fmap_types) > 1 and fmap_types - {'phase1', 'phase2'}:
@@ -399,7 +460,7 @@ def __attrs_post_init__(self):
             > 1
         )
 
-        if _pepolar_estimation and not anat_types:
+        if self.method == EstimatorType.UNKNOWN and _pepolar_estimation and not anat_types:
             self.method = MODALITIES[pepolar_types.pop()]
             _pe = {f.metadata['PhaseEncodingDirection'] for f in self.sources}
             if len(_pe) == 1:
@@ -455,6 +516,11 @@ def __attrs_post_init__(self):
         # special characters are not allowed.
         self.sanitized_id = re.sub(r'[^a-zA-Z0-9]', '_', self.bids_id)
 
+    @property
+    def is_dynamic(self) -> bool:
+        """The estimator emits a per-volume 4D fieldmap and no B-spline coefficients."""
+        return self.method in _DYNAMIC_METHODS
+
     def paths(self):
         """Return a tuple of paths that are sorted."""
         return tuple(sorted(str(f.path) for f in self.sources))
@@ -502,6 +568,33 @@ def get_workflow(self, set_inputs=True, **kwargs):
             from .workflows.fit.syn import init_syn_sdc_wf
 
             self._wf = init_syn_sdc_wf(**kwargs)
+        elif self.method == EstimatorType.MEDIC:
+            from .workflows.fit.medic import init_medic_wf
+
+            for f in self.sources:
+                if not f.path.is_file():
+                    raise FileNotFoundError(
+                        f'File path <{f.path}> does not exist, '
+                        'is a broken link, or it is not a file'
+                    )
+
+            self._wf = init_medic_wf(**kwargs)
+
+            if set_inputs:
+                phase_files = [f for f in self.sources if f.entities.get('part') == 'phase']
+                mag_echoes = {
+                    f.entities.get('echo'): f
+                    for f in self.sources
+                    if f.entities.get('part') == 'mag'
+                }
+                # Order by EchoTime (BIDS does not guarantee echo entity ==
+                # numeric order) and pair each magnitude with its phase echo,
+                # as validated in ``__attrs_post_init__``.
+                phase_files = sorted(phase_files, key=lambda f: f.metadata['EchoTime'])
+                mag_files = [mag_echoes[f.entities.get('echo')] for f in phase_files]
+                self._wf.inputs.inputnode.phase = [str(f.path.absolute()) for f in phase_files]
+                self._wf.inputs.inputnode.magnitude = [str(f.path.absolute()) for f in mag_files]
+                self._wf.inputs.inputnode.metadata = [f.metadata for f in phase_files]
 
         return self._wf
 

sdcflows/interfaces/reportlets.py

@@ -76,18 +76,26 @@ def _generate_report(self):
 
         fmapnii = nb.squeeze_image(rotate_affine(load_img(self.inputs.fieldmap), rot=canonical_r))
 
+        frame = 0
         if fmapnii.dataobj.ndim == 4:
-            for tstep in nb.four_to_three(fmapnii):
+            for idx, tstep in enumerate(nb.four_to_three(fmapnii)):
                 if np.any(np.asanyarray(tstep.dataobj) != 0):
-                    fmapnii = tstep
+                    fmapnii, frame = tstep, idx
                     break
 
+        # Dynamic estimators (e.g., MEDIC) pair a 4D fieldmap with a per-frame
+        # reference and mask; show the frame that matches the fieldmap snapshot.
+        movnii = refnii = _pick_frame(movnii, frame)
+
         if isdefined(self.inputs.moving):
             movnii = rotate_affine(load_img(self.inputs.moving), rot=canonical_r)
+            movnii = _pick_frame(movnii, frame)
 
         contour_nii = mask_nii = None
         if isdefined(self.inputs.mask):
-            contour_nii = rotate_affine(load_img(self.inputs.mask), rot=canonical_r)
+            contour_nii = _pick_frame(
+                rotate_affine(load_img(self.inputs.mask), rot=canonical_r), frame
+            )
             maskdata = contour_nii.get_fdata() > 0
         else:
             mask_nii = threshold_img(refnii, 1e-3)
@@ -140,3 +148,11 @@ def _generate_report(self):
             ),
             out_file=self._out_report,
         )
+
+
+def _pick_frame(img, frame):
+    """Return ``frame`` of a 4D image (clamped to its length), or ``img`` if 3D."""
+    img = nb.squeeze_image(img)
+    if img.ndim < 4:
+        return img
+    return img.slicer[..., min(frame, img.shape[3] - 1)]