[FEAT] MEDIC distortion correction via warpkit

[vanandrew]

Summary

Adds multi-echo dynamic distortion correction (MEDIC) to sdcflows, backed by warpkit. MEDIC estimates a per-volume B0 fieldmap directly from a multi-echo, mag+phase BOLD series, capturing breathing- and motion-driven field changes that a single static fieldmap can't.

Revives #435 / #438 with the simpler pure-Python warpkit (post vanandrew/warpkit#16, no Julia/C++ setup required). Closes #36.

What's new

Workflows

• init_medic_wf (sdcflows/workflows/fit/medic.py) — multi-echo phase + magnitude → 4D Hz fieldmap (one volume per timepoint, on the EPI grid) + brain-extracted reference + brain mask. Two-stage warpkit call (UnwrapPhase → ComputeFieldmap) so the per-frame masks stay accessible. Single fmap output is 4D for MEDIC, leaving the 3D-vs-4D dispatch to the apply consumer.
• init_dynamic_unwarp_wf (sdcflows/workflows/apply/dynamic.py) — per-volume apply path built on sdcflows.transform.apply_dynamic_unwarp, a per-frame extension of the same scipy/nitransforms-backed resampling that powers the static init_unwarp_wf. No warpkit needed at runtime. Includes Jacobian determinant intensity correction (jacobian=True default) sharing transform.fieldmap_jacobian with the static path.

Interfaces

• sdcflows/interfaces/warpkit.py — thin LibraryBaseInterface wrappers around the two MEDIC stages SDCFlows actually drives: UnwrapPhase (ROMEO) and ComputeFieldmap. Lazy import — sdcflows imports warpkit only when these interfaces actually run.

Detection / dispatch

• EstimatorType.MEDIC added. FieldmapEstimation.__attrs_post_init__ detects MEDIC inputs (part-{phase,mag} on bold/epi/sbref sources), enforces matched echo cardinality (≥2 phase, equal mag count), and rejects partial or mixed part sets. get_workflow instantiates init_medic_wf with EchoTime-sorted lists (BIDS doesn't guarantee echo entity == numeric order).
• Wrangler seeds MEDIC detection from both part='phase' and part='mag' queries (some datasets carry IntendedFor only on one side); dedup walks the full sibling set, no reliance on pybids ordering.
• force_medic opt-in on find_estimators — auto-discover MEDIC estimators from complex multi-echo BOLD even when neither IntendedFor nor B0FieldIdentifier is set. Pairing is unambiguous because the part-mag/part-phase echoes of the same run are MEDIC sources by construction. Intended for public datasets that ship the required echoes without the metadata the default discovery path needs.

Plumbing

• init_fmap_preproc_wf skips fmap_coeff for MEDIC (the fieldmap is on the EPI grid by construction, no B-spline rep).
• init_fmap_derivatives_wf uses MergeSeries(allow_4D=True) so MEDIC's 4D fmap passes through the same ds_fieldmap sink as the 3D static fmaps.

Packaging

• New warpkit extra in pyproject.toml, explicitly excluded from [all] because warpkit ships under a non-commercial WUSTL license. Default pip install sdcflows stays Apache-clean.
• tox.ini: the veryslow env pulls the warpkit extra so MEDIC end-to-end tests only run there.
• CI: datalad-fetches sub-04/ses-2 of ds007637 and sub-a01 of ds006926 as MEDIC fixtures (cache key bumped to v3).

Validation

• Three-layer multi-echo guard: wrangler seed query (echo=Query.REQUIRED), FieldmapEstimation cardinality check (len(phase_files) < 2), and _unpack_metadata runtime guard inside init_medic_wf.
• Construction tests for both workflows run on every CI lane (no warpkit needed).
• test_apply_dynamic_unwarp_matches_static pins apply_dynamic_unwarp to the same Hz→VSM + scipy.ndimage convention as the static _sdc_unwarp path — catches drift in sign / pe_info handling.
• test_wrangler_filter / test_wrangler_URIs parametrized with a 3-session × 3-echo × {mag,phase} BIDS skeleton.
• End-to-end runs on ds006926 (3-echo, 64×64×40, 205 frames) and ds007637 (5-echo, 110×110×72, 237 frames) confirm Jacobian factor sits at ~1.0 in the brain mask (std 0.087 / 0.118) with ~1–2% of voxels in compression/expansion regions and total signal change <2% — physically sensible.

Known compromises

• Default wrangler MEDIC discovery still keys off IntendedFor / B0FieldIdentifier — same constraint as the existing single-PE EPI branch. Datasets missing both can opt into discovery via force_medic=True on find_estimators (see Detection / dispatch above).
• The single fmap output is 4D for MEDIC. Downstream tools that expect 3D field maps need to either dispatch on dimensionality or block MEDIC-based estimators until they do.

Test plan

• pytest sdcflows/utils/tests/test_wrangler.py — wrangler MEDIC detection paths (including force_medic)
• pytest sdcflows/workflows/fit/tests/test_medic.py — init_medic_wf construction + _unpack_metadata guards
• pytest sdcflows/workflows/apply/tests/test_dynamic.py — init_dynamic_unwarp_wf construction + jacobian flag + per-frame resampling vs. static path
• pytest -m veryslow with pip install sdcflows[warpkit] and ds006926 / ds007637 fixtures present — full MEDIC fit + dynamic apply

[codecov]

Codecov Report

❌ Patch coverage is 77.35043% with 53 lines in your changes missing coverage. Please review.
✅ Project coverage is 81.07%. Comparing base (16a6d75) to head (b799d17).
⚠️ Report is 1 commits behind head on main.

Files with missing lines	Patch %	Lines
sdcflows/fieldmaps.py	56.25%	14 Missing ⚠️
sdcflows/interfaces/warpkit.py	76.36%	13 Missing ⚠️
sdcflows/workflows/base.py	0.00%	9 Missing ⚠️
sdcflows/workflows/apply/dynamic.py	75.00%	7 Missing ⚠️
sdcflows/transform.py	86.95%	3 Missing and 3 partials ⚠️
sdcflows/utils/wrangler.py	85.71%	3 Missing and 1 partial ⚠️

Additional details and impacted files

@@            Coverage Diff             @@
##             main     #541      +/-   ##
==========================================
- Coverage   84.10%   81.07%   -3.04%     
==========================================
  Files          30       33       +3     
  Lines        2938     3138     +200     
  Branches      391      432      +41     
==========================================
+ Hits         2471     2544      +73     
- Misses        389      513     +124     
- Partials       78       81       +3     

☔ View full report in Codecov by Harness.
📢 Have feedback on the report? Share it here.

🚀 New features to boost your workflow:

• ❄️ Test Analytics: Detect flaky tests, report on failures, and find test suite problems.

[vanandrew]

I have no idea why test (3.13, pre, fast) is failing on the cache step.

It looks like an out-of-space error on the runner (but the other jobs with the same steps succeeded, so very confused):

[test (3.13, pre, fast)](https://github.com/nipreps/sdcflows/actions/runs/25840682073/job/75990934202)
System.IO.IOException: No space left on device : '/home/runner/actions-runner/cached/2.334.0/_diag/Worker_20260514-125354-utc.log'
   at System.IO.RandomAccess.WriteAtOffset(SafeFileHandle handle, ReadOnlySpan`1 buffer, Int64 fileOffset)
   at System.IO.StreamWriter.Flush(Boolean flushStream, Boolean flushEncoder)
   at System.Diagnostics.TextWriterTraceListener.Flush()
   at GitHub.Runner.Common.HostTraceListener.WriteHeader(String source, TraceEventType eventType, Int32 id)
   at System.Diagnostics.TraceSource.TraceEvent(TraceEventType eventType, Int32 id, String message)
   at GitHub.Runner.Worker.Worker.RunAsync(String pipeIn, String pipeOut)
   at GitHub.Runner.Worker.Program.MainAsync(IHostContext context, String[] args)
System.IO.IOException: No space left on device : '/home/runner/actions-runner/cached/2.334.0/_diag/Worker_20260514-125354-utc.log'
   at System.IO.RandomAccess.WriteAtOffset(SafeFileHandle handle, ReadOnlySpan`1 buffer, Int64 fileOffset)
   at System.IO.StreamWriter.Flush(Boolean flushStream, Boolean flushEncoder)
   at System.Diagnostics.TextWriterTraceListener.Flush()
   at GitHub.Runner.Common.HostTraceListener.WriteHeader(String source, TraceEventType eventType, Int32 id)
   at System.Diagnostics.TraceSource.TraceEvent(TraceEventType eventType, Int32 id, String message)
   at GitHub.Runner.Common.Tracing.Error(Exception exception)
   at GitHub.Runner.Worker.Program.MainAsync(IHostContext context, String[] args)
Unhandled exception. System.IO.IOException: No space left on device : '/home/runner/actions-runner/cached/2.334.0/_diag/Worker_20260514-125354-utc.log'
   at System.IO.RandomAccess.WriteAtOffset(SafeFileHandle handle, ReadOnlySpan`1 buffer, Int64 fileOffset)
   at System.IO.StreamWriter.Flush(Boolean flushStream, Boolean flushEncoder)
   at System.Diagnostics.TextWriterTraceListener.Flush()
   at System.Diagnostics.TraceSource.Flush()
   at GitHub.Runner.Common.Tracing.Dispose(Boolean disposing)
   at GitHub.Runner.Common.Tracing.Dispose()
   at GitHub.Runner.Common.TraceManager.Dispose(Boolean disposing)
   at GitHub.Runner.Common.TraceManager.Dispose()
   at GitHub.Runner.Common.HostContext.Dispose(Boolean disposing)
   at GitHub.Runner.Common.HostContext.Dispose()
   at GitHub.Runner.Worker.Program.Main(String[] args)

[tsalo]

@vanandrew I think someone must have rerun that failing job and now it's passing.

[tsalo]

Thanks for this! I made a quick pass through the code.

[vanandrew]

Pushed four commits addressing your feedback (noise_frames, ctor inputs, nitransforms switch + interface cleanup, unified fmap output).

[vanandrew]

I can't tell if these test failures are something wrong with my dataset additions or if it's just another one-off error:

update(error): . (dataset) [Fetch failed: CommandError(CommandError: 'git -c diff.ignoreSubmodules=none -c core.quotepath=false fetch --verbose --progress --no-recurse-submodules --
prune origin' failed with exitcode 128 under /home/runner/sdcflows-tests/brain-extraction-tests [err: 'fatal: unable to access 'https://gin.g-node.org/nipreps-data/brain-extraction-
tests/': The requested URL returned error: 403'])] [CommandError: 'git -c diff.ignoreSubmodules=none -c core.quotepath=false fetch --verbose --progress --no-recurse-submodules --
prune origin' failed with exitcode 128 under /home/runner/sdcflows-tests/brain-extraction-tests [err: 'fatal: unable to access 'https://gin.g-node.org/nipreps-data/brain-extraction-
tests/': The requested URL returned error: 403']]
Error: Process completed with exit code 1.

[vanandrew]

Looks like the CI keeps failing because it's running out of space (I added two datasets for tests so that might have pushed it over the edge).

I pushed 6f5cfff2 adding a jlumbroso/free-disk-space@main step at the top of the test job:

- name: Free disk space
  uses: jlumbroso/free-disk-space@main
  with:
    tool-cache: false
    android: true
    dotnet: true
    haskell: true
    large-packages: false
    swap-storage: false

This gets rid of some unneeded development tools on the github runner to free room for data.

[vanandrew]

@tsalo Alright. I made changes based on our convo today:

• find_estimators discovers a fieldmap estimator for every fieldmap declared in BIDS intent metadata: grouping inputs by B0FieldIdentifier (Step 1), or falling back to IntendedFor heuristics (Step 2) only when no B0FieldIdentifier exists. (no_medic suppresses MEDIC via either route). It adds a sort that always returns dynamic estimators (i.e. MEDIC) first.
• Simplified fmap_ref and fmap_mask to be the first echo magntude and MEDIC fmap masks respectively. These likely won't be used by downstream consumers.
• Removed force_medic and replaced with no_medic

[vanandrew]

@tsalo I did a rebase to the latest main, hopefully to quash the 3.13 test error that happened in the last CI run: https://github.com/nipreps/sdcflows/actions/runs/26792754492/job/79155816998#step:22:320

[effigies]

There's a significant conceptual problem here. I'm not sure if it had been done correctly and then the distinction was lost after the latest refactor. If so, apologies for not getting to this sooner and saving you the time.

There is a comment that says the caller is responsible for ensuring that the fieldmap has already been projected into the moving (source) space, but we project them into the fixed (target) space.

To explain: We are resampling into the target space, so we need a VSM to correspond with every coordinate in the target space, and to then apply this shift in the source space. For a fixed fieldmap, this value is the same, so we can save computation by interpolating the fieldmap exactly once into the target space. For a dynamic fieldmap, there are no savings to be had here (and we don't know how to map the fieldmap into the target space without using the fieldmap to unwarp itself) and we need to interpolate the values in the target space. So if we want to have a single resampling function, it needs to have the logic:

if not dynamic:
    vsm = fmap_hz * pe_info[1]
else:
    vsm = ndi.map_coordinates(fmap_hz, coordinates, ...) * pe_info[1]

coordinates[pe_info[0], ...] += vsm

resampled = ndi.map_coordinates(data, coordinates)

As it's a simple function, it might be cleaner to do two separate functions, one with a fieldmap in target space and one with it in source space. I don't know if that makes the wrapping code easier or harder.

---

I started making some other comments, but stopped when I got to this point. Happy to discuss further on this thread or on a call.

[vanandrew]

@effigies Thanks for the review! I'll take a closer look at the comments later this week.

As for your conceptual concern, I think that comment might be a typo. I'm pretty sure I'm doing everything in the target space, not the moving space. But I'll take a look again and let you know for sure.

[vanandrew]

There is a comment that says the caller is responsible for ensuring that the fieldmap has already been projected into the moving (source) space, but we project them into the fixed (target) space.

@effigies Can you point out the line numbers of this comment?

[effigies]

sdcflows/sdcflows/transform.py

Lines 594 to 598 in a68a1e4

	Shared core of :meth:`B0FieldTransform.apply`. The caller is responsible
	for producing ``fmap_hz`` already on the ``moving`` grid — B-spline
	reconstruction plus coregistration for static estimators, or a pre-gridded
	per-frame field for dynamic estimators — and for ensuring ``moving`` has
	positive cosines.

[vanandrew]

@effigies I took a closer look at this. It's more of a badly phrased docstring. I reworded it in 7b02d9a to make it clearer, but there are two separate concepts:

• Grid (lattice): fmap_hz is co-gridded with moving — same voxel lattice — so the per-voxel scaling and coordinate shift broadcast element-wise. The sampling coordinates come straight from moving's own grid (transform.py:642, nt.linear.Affine(reference=moving)), which is why the field has to share that lattice. This is a precondition for the pull-resample, not a statement about distortion space.
• Semantics: the field is valued in the undistorted (target) sense. In _sdc_unwarp the shift is applied exactly once, on the target lattice (transform.py:95-97):

  # voxcoords is the deformation field, i.e., the target position of each voxel
  vsm = fmap_hz * pe_info[1]
  coordinates[pe_info[0], ...] += vsm

  No source-space treatment, no inversion.

So we agree the apply path works in target space — the old docstring just described it as "moving/source," which is exactly the contradiction you flagged. MEDIC hands in a field that's already target-valued and on the EPI grid, so it skips fit() (transform.py:508); the static estimators reach the same state via B-spline reconstruction + coregistration inside fit() (transform.py:521).