Module prelude 

PRELUDE phase unwrapping — best-pair-first region merging.

After M. Jenkinson, “A Fast, Automated, N-Dimensional Phase Unwrapping Algorithm”, FMRIB Technical Report TR01MJ1 (2001), https://www.fmrib.ox.ac.uk/datasets/techrep/tr01mj1/tr01mj1.pdf, condensed as Jenkinson M., Magn Reson Med 2003;49(1):193–197, https://doi.org/10.1002/mrm.10354.

PRELUDE is the reference unwrapper that the QSM and field-mapping literature benchmarks against, which is the reason it is here: it makes this crate’s unwrapper comparisons readable against the published ones.

The literature generally reports it as slower and less accurate than unwrap_romeo — Dymerska et al., Magn Reson Med 2021;85(4):2294–2308, https://doi.org/10.1002/mrm.28563. That comparison ran FSL’s implementation. The accuracy half carries over, since this module gives FSL’s output (below); the speed half does not. FSL’s prelude at its defaults took 31 min for the KKI travelling head’s five echoes (256×256×128, 1.52M in-mask voxels); this module takes 4.6 s for the same output (release build, parallel, 8 cores; 0.8–1.1 s per echo). On the bids/ phantom (164×205×205, 1.31M in-mask voxels) it takes 0.40 s for the first echo and 1.45 s for the last at the defaults. With full-3D labelling at 6 partitions it took, per echo:

	first echo	last echo
PRELUDE (full 3D, 6 partitions)	0.19 s	0.34 s
ROMEO	0.34 s	0.29 s
best path	0.32 s	0.33 s
Laplacian	0.46 s	0.55 s

and the local-field correlations after an identical V-SHARP were 0.615 (PRELUDE), 0.608 (Laplacian), 0.602 (best path), 0.592 (ROMEO) — a spread narrower than any of them is from ground truth, which is the agreement the comparison is for rather than a ranking.

Reference pipeline

The original is FSL’s prelude (Jenkinson’s own implementation, fsl-fugue), and the defaults (PreludeParams::default) are its defaults: 8 phase partitions (its -n; the report’s worked example uses 6) and initial regions labelled slice by slice, then merged in 3D, on any volume over 64 voxels along either in-plane axis, and over the whole volume on smaller ones (PreludeLabelling::Auto). --labelslices labels slice by slice whatever the size (PreludeLabelling::Slices), and -f over the whole volume (PreludeLabelling::Volume). On the in-vivo and simulated volumes measured its output has a median in-mask wrap count of 0, which PreludeParams::correct_global (on by default) reproduces (see “The global constant” below). FSL’s own masking from the magnitude (when no mask is given) is not reproduced; a mask is required here. Nothing in the algorithm depends on voxel size, field strength or orientation. The QSM pipeline offers it as UnwrappingAlgorithm::Prelude (each echo on its own, then whole turns aligned between echoes); it is also UnwrapMethod::Prelude for the slice-wise unwrapper and MCPC-3D-S’s HIP unwrapping.

Why slice labelling is the default (QSM.rs#397, Ashley’s decision 2026-11-02): it is the original’s default on every in-vivo volume (FSL’s switch to full 3D below 65 voxels in-plane is kept too, as PreludeLabelling::Auto, so that small volumes also give FSL’s output), and on the KKI travelling head’s late echoes full-3D labelling lets regions straddle wraps that slice labelling does not. Measured as the share of in-mask voxels whose echo-to-echo step u_e − u_{e−1} is not wrap(φ_e − φ_{e−1}) plus the most common 2π multiple (a temporal-consistency check that needs no ground truth; lower is better), KKI echoes 2→5:

labelling, partitions	e2	e3	e4	e5
full 3D, 6 (prelude -f -n 6)	3.21 %	3.48 %	12.36 %	6.87 %
full 3D, 8 (prelude -f)	3.21 %	3.47 %	4.66 %	11.93 %
slices, 6 (prelude -n 6)	3.21 %	3.47 %	3.45 %	3.45 %
slices, 8 (defaults)	3.21 %	3.47 %	3.45 %	3.46 %

On the simulation (4 echoes at 7 T) the labelling makes no difference to it (5.5–5.7 % for all; its 8 ms echo spacing sets that floor). Full-3D labelling is kept for parity with prelude -f.

Against FSL’s prelude

FSL’s binary was run as a black box (QSM-CI’s algorithms/_fsl-prelude: fsl-fugue 2603.0 from FSL’s conda channel; its source was not read) and compared with this module inside the mask, each echo on its own, on QSM-CI’s simulation (164×205×205, 4 echoes, 1.31M in-mask voxels), the KKI travelling head (256×256×128, 5 echoes, 1.52M) and the Prisma travelling head’s first echo (184×256×128, 1.70M) (QSMxT/QSM-CI#302, QSM.rs#397; measured 2026-11-02). Counted is every in-mask voxel whose wrap count differs from FSL’s, beyond one offset per connected component of the mask:

here	FSL	sim	KKI	Prisma e1
defaults (auto = slices here, 8)	prelude	0 on all 4 echoes	0 on all 5 echoes	0
slices, 6	prelude -n 6	0 on all 4 echoes	0 on echoes 1–4; 7 voxels (0.0005 %) on echo 5	—
full 3D, 8	prelude -f	0 on all 4 echoes	0 on all 5 echoes	0
full 3D, 6	prelude -f -n 6	0 on all 4 echoes	0 on all 5 echoes	—

At the defaults the output is FSL’s exactly, global constant included (largest difference 0 after FSL’s single-precision output), and the region counts agree: FSL’s -v reports 13111, 23565, 40070, 59635 and 82032 initial regions on KKI’s echoes at its defaults (4739 to 40265 with -f), which is what 4-connected labelling of the third axis’s slices gives, and not the other two axes’. FSL writes 0 outside the mask. FSL’s default run took 31 min on KKI and, with the simulation’s four echoes in one 4-D run, did not finish; run one echo at a time it took 72 s, 492 s and 1964 s on echoes 2 to 4 (31618 to 118208 regions); this module takes 0.4 to 1.5 s per echo there.

When FSL’s default labels slices (PreludeLabelling::Auto): only when the volume is over 64 voxels along its first or second axis. 64×64×5, 64×64×64, 64×64×128, 40×36×200 and 20×20×300 volumes were labelled in full 3D (-v region count equal to -f‘s), and 65×10×20, 10×65×20, 65×65×20, 128×32×20 and 256×256×8 slice by slice (equal to --labelslices’), whatever the voxel size (0.4, 1, 3 mm) or the mask’s extent (a 128-wide volume masked to its first 38 columns is still labelled by slices).

The global constant. Each connected set of regions comes out of the merge with an arbitrary whole number of turns. FSL’s output on the three volumes above has a median in-mask wrap count of 0, which PreludeParams::correct_global gives; on the small noisy phantoms of tests/prelude_fsl_parity.rs FSL’s constant is up to two turns off that median (and is not the mean, the mode, the first voxel’s, the centre voxel’s or the largest region’s either). No rule was found for it, so the fixture compares wrap counts up to that one offset: every voxel of all eight runs (defaults, --labelslices, -f, -n 6, on a 72- and a 64-wide phantom) matches.

prelude -s (everything slice by slice, no 3D merge) differs from unwrap_slicewise with UnwrapMethod::Prelude: FSL writes 0 in small in-slice components (14k and 23k in-mask voxels on the simulation and KKI) and, beyond a per-component offset, 0.1–0.6 % of the remaining voxels take another wrap count. Slice-wise unwrapping is this crate’s own, so that is information only.

Provenance

This is a clean-room implementation written from the technical report alone. FSL’s prelude.cc was not read, fetched, or consulted. That is a licensing requirement, not a stylistic one: FSL is licensed for non-commercial use only and the licence reaches “the Software or any derivative of it”, so transliterating it would encumber this MIT crate. The report specifies the whole algorithm, down to the data-structure choices, so nothing had to be reverse-engineered. Equation numbers in the comments below refer to it. The slice labelling is the planar restriction of the report’s §4.1; that it is FSL’s default, with 8 partitions, and that it labels slices of the third axis, was established by running FSL’s binary as a black box (its -v output and region counts; see above), not from its code.

The algorithm

The unwrapped phase is the wrapped phase plus a whole number of turns, θ_j = φ_j + 2π·m_j (eq 1). Solving for a per-voxel m_j is hopeless, so the volume is first cut into regions that are each assumed wrap-free, and every region gets a single offset.

1. Initial regions (§2.4, §4.1). Split [-π, π) into num_phase_partitions equal bins. For each bin, mask the in-mask voxels whose wrapped phase falls in it and label the connected components: within each slice of the third axis, 4-connected (PreludeLabelling::Slices; what the default PreludeLabelling::Auto does on in-vivo volumes), or over the volume, 6-connected (PreludeLabelling::Volume). Every component across every bin (and slice) becomes a region with its own id. Two voxels in one region therefore differ in wrapped phase by less than one bin width, so a region straddling a wrap is unlikely — and the narrower the bin, the less likely still. Limiting a region to one slice also cuts the chance of two distant areas at the same phase being labelled one region (§4.1), which is how a full-3D region comes to straddle a wrap. The merge below is 3D either way.

2. Cost (eq 2–4). Across the interface between regions A and B the cost is Σ (φ_Aj − φ_Bk + 2π·M_AB)², summed over interfacing voxel pairs, with M_AB = M_A − M_B. Totalled over all interfaces that is the function to minimise. Differentiating (eq 5–6) gives the continuous minimum M_AB = −P_AB / (2π·N_AB) for N_AB interfacing pairs and P_AB their summed phase difference. But M_AB has to be an integer, which makes the exact problem integer programming with 2^N neighbouring solutions — hence the greedy scheme below.

3. Best-pair-first merging (§2.3). Treat each interface as the constraint M_AB = round(K_AB), K_AB = −P_AB/(2π·N_AB) (eq 7). Getting one of these wrong by ±1 costs ΔC_AB = 8π²·N_AB·(½ − |K_AB − L_AB|) (eq 8–9). Merge the pair that maximises it (eq 10) — in the report’s words, “choose the pair of regions where getting the offset ‘wrong’ is the most disastrous”. A wide interface whose mean phase step sits close to a whole number of turns is both the most confident call available and the most expensive one to get wrong, so it is settled first and the ambiguous interfaces inherit what it decided. Merging sums the statistics over the merged region’s remaining neighbours (eq 11–12). Repeat until no interfaces are left.

Deviations from the report

• The merge statistics are kept in the current frame. Equations 11–12 read N_QW = N_AW + N_BW, P_QW = P_AW + P_BW. The counts add directly; the summed differences do not, until both sides are expressed in one frame. Merging at L_AB only fixes the two sets relative to each other, so one of them has to move — and its P to every other region shifts by ±2π·L_AB·N before the sum. Which side moves is not free to choose here: the union-find keeps the larger set’s representative, so the absorbed set is the one that moves, and the sign follows from that rather than always being B’s. The report says to merge “using the ‘optimal’ offset L_AB” and leaves the bookkeeping implicit; merge_pair writes it out. Getting the direction wrong is silent — it biases every later merge by one turn, and only shows up on volumes where the two sets differ in size enough for the union-find to pick the other root.
• A lazily-invalidated binary heap replaces the O(C) scan. The report searches its map of interfaces linearly for the maximum ΔC, at most R times, which it costs as a small penalty next to the cheap insert and delete. That holds at its worked example’s scale (1920 regions); it does not at the hundreds of thousands a high-resolution 3D volume can produce, where O(R·C) dominates everything else. The merge order is identical — the heap returns the same maximum — only the way it is found differs.
• Masked voxels keep their original wrapped phase on output, matching unwrap_bestpath and unwrap_romeo.

Each connected set of regions is unwrapped independently and carries its own arbitrary multiple of 2π — nothing links regions across a gap in the mask. Set PreludeParams::correct_global for a defined choice of that constant.

Structs

Candidate 🔒

One entry in the merge queue: an interface, with the ΔC it had when queued.

Interface 🔒

Running statistics for one interface between two regions.

PreludeParams

Parameters for unwrap_prelude.

Enums

PreludeLabelling

How the initial regions are labelled (§2.4, §4.1): which voxels of one phase partition count as connected.

Constants

AUTO_SLICES_ABOVE

The in-plane size above which PreludeLabelling::Auto labels slice by slice: FSL prelude labels a 64×64×N volume in full 3D and a 65×N×M or N×65×M one slice by slice, whatever the slice count, voxel size or mask.

NO_REGION 🔒

Sentinel in the region map for a voxel that belongs to no region — outside the mask, or carrying a non-finite phase.

PI 🔒

TWO_PI 🔒

Functions

build_interfaces 🔒

Accumulate N_AB and P_AB for every interface, as a symmetric adjacency list.

label_regions 🔒

Partition the phase range, label the connected components of each partition, and hand every component across every partition a unique region id (§2.4).

merge_pair 🔒

Merge regions a and b at the offset their interface implies, fold their adjacency together, and queue the interfaces that changed.

merge_regions 🔒

Run the best-pair-first merge to exhaustion (§2.3).

phase_partition 🔒

Which phase partition a voxel falls in, or None if it has no usable phase.

unwrap_prelude

Unwrap a 3D phase volume with PRELUDE’s region-merging algorithm.