✨ Detect and correct global 2π branch errors in MCPC-3D-S

[vanandrew]

Replaces a cascade of five hardcoded field-magnitude thresholds with a two-stage rule, and fixes a real 40 Hz field-map artifact found on OpenNeuro data.

The bug

Spatial unwrapping recovers phase differences between voxels, never the absolute turn count, so the dual-echo field is determined only modulo 1/dTE. ROMEO's correct_global resolves that by preferring the smallest field, which is equivalent to requiring |f| < 1/(2·dTE).

It re-derives that choice on every frame, from an unweighted median of rounded wrap counts over a dilated mask. When a subject's field approaches the half-wrap the ballot sits on a knife edge, and ordinary frame-to-frame variation tips it. Because the turn count is an integer, a tipped frame is displaced by a full wrap rather than slightly.

Two subjects in ds006131 (TEs 14.2/38.93/63.66 ms, wrap 40.44 Hz, half-wrap 20.22 Hz) do this:

                     correct_global alone              with the selector
sub-24630   -24.69..+17.26 Hz, 41.44 Hz max step   +14.13..+19.17 Hz, 5.04 Hz max step
            158 of 243 frames displaced             0 frames displaced
sub-20828   -23.74..+17.78 Hz, 40.98 Hz max step   +15.24..+18.63 Hz, 1.54 Hz max step
             19 of 243 frames displaced             0 frames displaced

The damage is temporal. A uniformly wrong turn count would be a constant offset in a quantity already referenced to an arbitrary demodulation frequency, and undetectable. A count that varies between frames of one acquisition corrupts the temporal structure the acquisition exists to measure.

The rule

1. Consistency. A turn-count error of M shifts the phase offset by wrap(2π·k·M), k = TE0/dTE, landing as the same additive constant on every echo. Fitting a line through the echoes and reading its intercept measures it: a self-consistent reconstruction extrapolates through the origin. The rejection scale is analytic, |wrap(2π·k)|, so nothing is tuned to data. Fitting and non-fitting candidates separate by 6–9 orders of magnitude.
2. Prior, only on a tie. Several candidates can be exactly through-origin; that alias is real and no phase statistic resolves it. Among survivors, take the smallest weighted-median field.

This is detectable at all only because the offset and the field are computed in separate stages. The offset is formed first; the field comes from a second unwrap that applies its own global correction with no knowledge of which count the offset assumed. On a tipped frame the two disagree, and that mismatch is what the intercept measures.

Stage 2 applies the same prior as correct_global, on a better-conditioned statistic: the field itself, weighted, over an eroded mask, rather than rounded counts over a dilated one. On the failing frames both candidates score ~1e-7, so stage 1 detects the flip but cannot rank them; the prior settles it.

The weighting

The prior weights by mag1², the second echo. The field is read from phase1 − phase0, whose noise is dominated by the weaker echo; weighting on mag0 gives full weight to fast-T2* voxels (healthy echo 0, collapsed echo 1), which are exactly the air–tissue interfaces the prior should discount.

This does not move the decision boundary, which is fixed at 1/(2·dTE), and barely moves the margin. What it buys is stability, and that is what predicts a flip:

	total errors	worst-case margin	between-subject spread
mag0²	0	0.884 Hz	2.044 Hz
mag1²	0	1.052 Hz	0.667 Hz

Also documented in the docstrings: the estimator is a median rather than a mean or M-estimator because the field distribution is broad and right-skewed with the half-wrap boundary inside its body, so anything pulled toward the mean reads closer to the boundary without being more accurate. A Huber estimator cut the worst-case margin to 0.045 Hz; a half-sample mode was tried and reverted for discontinuity.

Changes

• Add _evaluate_branch, _branch_intercept_step, _weighted_median, _select_branch; remove _select_branch_cascade and all four heuristic constants.
• Return 0 when TE0/dTE is an integer, where the count cannot move the offset.
• Correct two annotations: echo times are numpy scalars at every call site, not arrays; rescale_phase's min/max are floats.
• 16 tests, including a regression fixture built from the real failing frames (tests/data/branch_flip/, 9.7 MB, excluded from the sdist). Cropped to the brain bbox; deliberately not decimated, since resampling tips the knife-edge ballot and would invert what the test checks.

Breaking

--wrap-limit is removed from wk-medic and wk-unwrap-phase, and wrap_limit= from medic() and unwrap_and_compute_field_maps(). It existed to disable the heuristics this PR deletes.

Validation

100 runs, 24,105 frames, 10 OpenNeuro ME-EPI datasets, six TE0/dTE ratios (0.51–0.78) and echo spacings 16.2–27.5 ms. Zero errors under either weighting.

71 of those runs were processed after the weighting was chosen and held out from the comparison. mag1²'s spread advantage was 2.3× on the data that selected it and 2.7× on the held-out set, so it replicated rather than regressing.

The prior changes an answer on 2 of 77 prior-live runs, both above |f|/half-wrap = 0.81; below 0.59 it was never load-bearing. The worst frame anywhere reached 0.948, leaving a 2.10 Hz decision margin on a 40.44 Hz wrap.

CI

Pure Python — no pybind11, ITK, or wheel-matrix impact. Cost is ~+14% per frame for three candidate reconstructions.

[codecov]

Codecov Report

❌ Patch coverage is 95.65217% with 2 lines in your changes missing coverage. Please review.
✅ Project coverage is 96.25%. Comparing base (764aa57) to head (d066cf7).
⚠️ Report is 1 commits behind head on main.

Files with missing lines	Patch %	Lines
warpkit/unwrap.py	95.65%	2 Missing ⚠️

Additional details and impacted files

@@            Coverage Diff             @@
##             main      #32      +/-   ##
==========================================
+ Coverage   95.01%   96.25%   +1.23%     
==========================================
  Files          18       18              
  Lines        1265     1280      +15     
==========================================
+ Hits         1202     1232      +30     
+ Misses         63       48      -15     

☔ 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.

[Copilot]

Pull request overview

This PR updates MEDIC’s MCPC-3D-S global (2\pi) branch selection by replacing the prior field-magnitude heuristic cascade with an intercept-based consistency test (no tuned constants), and removes the now-obsolete wrap_limit API/CLI surface. It also adjusts a couple of type annotations and adds targeted tests to validate the new decision rule and an end-to-end injected-wrap recovery case.

Changes:

• Implement intercept-based branch evaluation/selection (_evaluate_branch, _branch_intercept_step, _select_branch) and remove the prior heuristic cascade and its tuning constants.
• Remove --wrap-limit / wrap_limit= from CLIs and library entry points that previously exposed heuristic disabling.
• Add tests for the analytic intercept step, decision-rule branches, and end-to-end recovery on bundled data.

Reviewed changes

Copilot reviewed 6 out of 6 changed files in this pull request and generated 3 comments.

Show a summary per file

File	Description
warpkit/utilities.py	Adjusts rescale_phase parameter annotations (int → float).
warpkit/unwrap.py	Adds intercept-based branch selection and removes heuristic-based wrap limiting.
warpkit/scripts/unwrap_phase.py	Removes --wrap-limit plumbing from wk-unwrap-phase.
warpkit/scripts/medic.py	Removes --wrap-limit plumbing from wk-medic.
warpkit/distortion.py	Removes wrap_limit from the public medic() wrapper call into unwrapping.
tests/test_unwrap.py	Adds unit tests for the new branch selector and an end-to-end injected-wrap recovery test.

---

💡 Add Copilot custom instructions for smarter, more guided reviews. Learn how to get started.