ENH: Add optional warpkit MEDIC support for multi-echo BOLD distortion correction#3645
ENH: Add optional warpkit MEDIC support for multi-echo BOLD distortion correction#3645DVSneuro wants to merge 16 commits into
Conversation
|
Thanks for opening this pull request! It looks like this is your first time contributing to fMRIPrep. 😄 |
=== Do not change lines below ===
{
"chain": [],
"cmd": "bash -c '! pixi lock --check || git checkout .'",
"exit": 0,
"extra_inputs": [],
"inputs": [
"pixi.lock",
"pyproject.toml"
],
"outputs": [
"pixi.lock"
],
"pwd": "."
}
^^^ Do not change lines above ^^^
|
This is extremely interesting and coincides with some work by @vanandrew in sdcflows (note to self: link PR). I've been meaning to look into how it would integrate and haven't found the time. Have you run datasets with this workflow? Do the results look good to you? Do you have comparisons with other fieldmap types? I will try to look at this in some detail soon. |
|
I've only ran a couple of subjects from our dataset (https://openneuro.org/datasets/ds005123/versions/1.1.3), and here's one example visualization for a single time point, which looks pretty good to me. 
I haven't compared it directly with standard fieldmaps, though we have those in the dataset above. We were drawn to warpkit/medic since some of our participants move a lot (older adults) and we had a few cases where the participant needed to get out of the scanner for a break but we forgot to get a new fieldmap when they got back in. So, the idea of making the fieldmaps for each frame using the multi-echo phase data was very appealing to us! But when we first started using medic/warpkit a year or so ago, I couldn't figure out how to get it integrated into fmriprep, so we'd generate a fieldmap for the first time point of each run and use that in the normal fmriprep workflow and the results have been good. I'm assuming the framewise implementation should be even better, especially for high-motion participants -- but hopefully I didn't break anything or do anything too silly with this PR! |
|
I got you @effigies: nipreps/sdcflows#541 :) I haven't reviewed @DVSneuro PR yet, but my vague plan was to have MEDIC be called from On the fmriprep side, the biggest changes would need to be in how the field map is applied: in MEDIC, you get a field map for every frame of your scan, so the one step resampling would need to be adjusted to account for this. Given that @DVSneuro has some results, I assume your PR already makes a change like this? Awesome job getting this working btw :). |
|
Thanks @vanandrew ! I had plenty of help from Codex. And yes, this PR does update the fMRIPrep side so a 4D MEDIC fieldmap series is applied volume-by-volume during one-step resampling. The current implementation calls MEDIC directly from fMRIPrep, but I’d be open to trying to restructure that piece if the preferred long-term architecture is to generate the MEDIC estimates in sdcflows and then consume them from fMRIPrep. I'm not seeing any other multi-echo datasets on OpenNeuro that have the phase data, so the testing here has been somewhat limited. |
|
@DVSneuro, I suggest checking https://me-ica.github.io/open-multi-echo-data/. I try to keep it up to date, and you can filter by if there's complex-valued data or not. There are nine on OpenNeuro that I'm aware of:
I think the goal should be to have MEDIC within SDCFlows so that dynamic distortion correction is treated like static distortion correction, though it would be cool if fMRIPrep could use One change that might be necessary on fMRIPrep's side is (optionally?) swapping the order of HMC and SDC when MEDIC is used, since I think dynamic distortion correction should improve motion correction. |
|
For datasets, I've mostly been using this one for my own tests: https://openneuro.org/datasets/ds006926/versions/1.0.0 (Thanks @tsalo for the list of mag/phase datasets). Getting into some of the more technical details:
Side question: I noticed that both nipreps project currently support python version 3.10+, but I've made |
|
Thanks, @tsalo and @vanandrew ! I'm happy to try and help with the sdcflows approach. In the meantime, I ran additional validation with an Apptainer image built from this branch using Successful runs so far:
One remaining test case is |
`warpkit.api.medic` exposes three outputs: `fieldmap` (undistorted-space, the consumer-facing output), `fieldmap_native` (a distorted-space debug output), and `displacement_map`. The Stage-3 wiring used `fieldmap_native`, which makes fmriprep's pull-resampler look up the field at the wrong physical location at each target voxel (it's off by the SDC displacement itself). `resample_vol` documents fmap_hz as "sampled to the target space, in Hz" (interfaces/resampling.py:272-273) — i.e., the undistorted/output grid. That's `fieldmap`, not `fieldmap_native`. Verified empirically by running `resample_series_async` against both fields on two datasets: * warpkit's 3.5mm fixture: RMS Δ ≈ 358, max|Δ| ≈ 17.5k vs mean signal ~7050 * ds007637 sub-04 at 2mm: RMS Δ ≈ 211, max|Δ| ≈ 21k vs mean signal ~4175 Relative error grows with resolution, as expected — a fixed physical displacement spans more voxels in a finer grid.
`func_fit_reports_wf.inputnode.fieldmap` flows into FieldmapReportlet (workflows/bold/outputs.py:373-401), which expects a 3D static field. The Stage-3 warpkit branch was passing the 4D HMC-aligned per-frame field directly, so the reportlet either crashes on the 4D NIfTI or silently renders frame 0 only. Pipe the temporal mean (already computed by `warpkit_mean` on line 526) instead. The 4D series remains the outputnode.fieldmap that the bold native workflow consumes for per-frame correction.
`mem_gb['resampled']` is the memory footprint of one resampled BOLD volume (`filesize * 4` in `estimate_bold_mem_usage`). MEDIC actually holds all N echoes of phase + magnitude in memory as float64, plus warpkit's internal state (unwrapped phases, masks, per-frame fields). For typical multi-echo data this underestimates real usage by an order of magnitude. ds007637 sub-04 is 5 echoes × ~320 MB compressed — close to 3 GB resident before warpkit's working set. The nipype scheduler will happily oversubscribe and OOM. Use `2 * n_echoes * filesize` instead, which roughly doubles the input footprint to cover outputs and warpkit's internal arrays.
The Stage-3 warpkit branch cleared `fmapref_buffer.sbref_files` and gated `raw_sbref_wf` with `not warpkit_enabled`. SBRef is acquired separately from the BOLD run and is orthogonal to MEDIC — dropping it silently is a regression for sites that acquire it. If there's a concrete incompatibility between SBRef and the MEDIC field-application path, it should be a stated limitation with a clear error, not an invisible behavior change. Removing the gate so SBRef flows through the same way it does on every other distortion-correction path.
`bold_series` is sorted by EchoTime in `workflows/base.py:250` and `phase_files` is independently sorted by EchoTime in `collect_bold_part_files`. They *should* line up, but a missing or zero-default EchoTime sidecar on one side and not the other silently transposes the pair set fed to warpkit. Add an explicit cross-check so that misalignment fails fast with a readable message instead of producing garbage fieldmaps.
|
Thanks, @vanandrew! Just tying up this loose end: your suggestion about the I pushed a small follow-up commit exposing warpkit’s
This resolved the previous warpkit phase-unwrapping failure. The log shows the expected trimming message, both MEDIC nodes completed, and fMRIPrep finished successfully. ds006072_v1_0_8_sub-P1R_run-20260511-132325_fmriprep-warpkit.log I also saw the draft changes in DVSneuro#1, including the move toward routing MEDIC through SDCFlows. Thanks for doing that, and I'm happy to test whenever it's ready! I understand that move may make |
FIX: correctness fixes for MEDIC pipeline
|
Validation update after testing @vanandrew 's correctness-fix branch from DVSneuro#1. I rebuilt the container and reran the public validation datasets I’ve been using. MEDIC workflows ran and completed for:
For ds005085/sub-10017, which includes a mix of single-echo and multi-echo runs, the single-echo runs followed the standard PHASEDIFF fieldmap path, while the eligible 4-echo runs used MEDIC/warpkit. I previously tested with ds005085/sub-10015, which didn't have any multi-echo data (an oversight on my part), but that test case didn't fail, which may be undesirable with the I also separately reran ds006072/sub-P1R/ses-1 with Otherwise, across these runs, I did not see crashes, tracebacks, fieldmap-length mismatches, or OOM failures in the logs. One small reporting/QC question I noticed while reviewing the HTML: the MEDIC fieldmap reportlet uses the sdcflows FieldmapReportlet hover/toggle behavior, and the Thanks again for helping with this! @effigies -- looking forward to your comments and feedback when you get a chance! Best wishes, |
|
Hey @DVSneuro @vanandrew @tsalo I wonder if we could catch up on a call this week. I'm pretty slammed at the moment, so I'm not sure when I can review this thoroughly and don't want to leave you hanging. A call where you could generally walk us through it all might be helpful. |
|
@effigies that would be great. I can make time for a call. |
Codecov Report❌ Patch coverage is Additional details and impacted files@@ Coverage Diff @@
## master #3645 +/- ##
==========================================
+ Coverage 74.41% 74.84% +0.42%
==========================================
Files 62 64 +2
Lines 4987 5235 +248
Branches 637 676 +39
==========================================
+ Hits 3711 3918 +207
- Misses 1142 1172 +30
- Partials 134 145 +11 ☔ View full report in Codecov by Sentry. 🚀 New features to boost your workflow:
|
|
@effigies No worries. I'm also pretty busy this week, but can do next week or later if we want to do a call. Where I think things stand: MEDIC needs to land in SDCflows first (nipreps/sdcflows#541), or at least be stable enough to build off of. I wouldn't merge this PR as-is, it predates the SDCflows work and would lock in a parallel MEDIC path we're about to rework. Once nipreps/sdcflows#541 is in, we can start the fMRIPrep side. @DVSneuro has a good start here, and pieces of it carry over (the part-phase BIDS filtering, the framewise-4D ResampleSeries path, boilerplate/citations). The main thing to revisit is that with SDCflows owning MEDIC as a first-class estimator, we can simplify by routing it through the normal fieldmap pipeline rather than as a multi-echo-specific branch. |
|
Thanks, @effigies. I’d be happy to chat through it, though no pressure on timing from my end. This isn’t holding anything up for us, and I’m also pretty swamped with a couple of grants right now, so the first week of June or later would work well for me! I agree with @vanandrew’s read here. If the preferred long-term path is for MEDIC to land in SDCFlows first via nipreps/sdcflows#541, then I’m very happy to treat this PR as a working prototype/validation branch rather than something that should merge as-is. When the timing makes sense, I’d be happy to walk through what this branch currently does, the datasets I tested, and the pieces that might carry over to the fMRIPrep side, especially the |
|
@vanandrew @DVSneuro Sounds good. I've created a when2meet for next week and filled in my availability: https://www.when2meet.com/?36851069-sVt5P If we can't make that, I will next have some availability after OHBM. cc @mgxd @tsalo @oesteban (Oscar, I assume you're not available, but I think you'd be interested if you are.) |
|
Meeting invitation sent |
|
Sorry, I missed this. When's the meeting, finally?
…On Sat, May 30, 2026, 15:51 Chris Markiewicz ***@***.***> wrote:
*effigies* left a comment (nipreps/fmriprep#3645)
<#3645 (comment)>
Meeting invitation sent
—
Reply to this email directly, view it on GitHub
<#3645?email_source=notifications&email_token=AAESDRTHX3RVW6YY6TZMNHL45LRN5A5CNFSNUABFM5UWIORPF5TWS5BNNB2WEL2JONZXKZKDN5WW2ZLOOQXTINJYGMYDCMJXGIZKM4TFMFZW63VHNVSW45DJN5XKKZLWMVXHJLDGN5XXIZLSL5RWY2LDNM#issuecomment-4583011722>,
or unsubscribe
<https://github.com/notifications/unsubscribe-auth/AAESDRWXKQSECETMOMJIGTT45LRN5AVCNFSM6AAAAACYVLXPPKVHI2DSMVQWIX3LMV43OSLTON2WKQ3PNVWWK3TUHM2DKOBTGAYTCNZSGI>
.
Triage notifications, keep track of coding agent tasks and review pull
requests on the go with GitHub Mobile for iOS
<https://github.com/notifications/mobile/ios/AAESDRXR7HLHIGX3ZXJN5ML45LRN5A5CNFSNUABFM5UWIORPF5TWS5BNNB2WEL2JONZXKZKDN5WW2ZLOOQXTINJYGMYDCMJXGIZKM4TFMFZW63VHNVSW45DJN5XKKZLWMVXHJKTGN5XXIZLSL5UW64Y>
and Android
<https://github.com/notifications/mobile/android/AAESDRXXYZCURC5TU77R4ND45LRN5A5CNFSNUABFM5UWIORPF5TWS5BNNB2WEL2JONZXKZKDN5WW2ZLOOQXTINJYGMYDCMJXGIZKM4TFMFZW63VHNVSW45DJN5XKKZLWMVXHJLTGN5XXIZLSL5QW4ZDSN5UWI>.
Download it today!
You are receiving this because you were mentioned.Message ID:
***@***.***>
|
|
3pm ET on Monday. There is also a Thursday slot everyone can make. |
| jacobian, | ||
| hmc_xfms[0] if hmc_xfms else None, | ||
| fmap_hz, | ||
| fmap_hz[..., 0] if fmap_is_series else fmap_hz, |
There was a problem hiding this comment.
I don't think you need the inline if/else since line 382 raises an error if fmap_is_series and data.ndim == 3. So this should only be fmap_hz.
6 participants
Changes proposed in this pull request
I think this is a little related to #3336 and maybe some other work from @tsalo , but this PR addresses a narrower use case than that issue: optional MEDIC-based susceptibility distortion correction for compatible multi-echo BOLD runs with
part-phasedata, usingwarpkitfollowing the logic from Andrew Van's ( @vanandrew ) preprint on bioRxiv: https://www.biorxiv.org/content/10.1101/2023.11.28.568744v1When
--me-use-warpkitis enabled, fMRIPrep will:part-phasedataMEDIC (warpkit)in the functional summaryThis path is opt-in and does not change default behavior for existing workflows.
Main changes in this PR:
--me-use-warpkitwarpkit.api.medicCurrent scope / limitations:
part-phasecompanionsValidation performed:
warpkitIn the HTML report, warpkit-enabled runs showed:
Susceptibility distortion correction: MEDIC (warpkit)This is my first contribution to fMRIPrep, so I tried to keep the scope focused and follow existing workflow patterns closely.
I used OpenAI Codex to help draft and iterate on portions of the code and tests, and I manually reviewed and tested the resulting changes as carefully as I could, but I wouldn't say I'm an expert on this (maybe a little better than vibe coding, but not by much).
Documentation that should be reviewed
This PR does not add standalone user documentation pages, but it does update user-facing/report-facing text in a few places:
--me-use-warpkit