Module medi_region_growing 

MEDI’s region-growing phase unwrapping: unwrapPhase / mexUnwrap.c (Tian Liu), a magnitude-guided, quality-ordered flood fill.

Liu T, Wisnieff C, Lou M, Chen W, Spincemaille P, Wang Y. Nonlinear formulation of the magnetic field to source relationship for robust quantitative susceptibility mapping. Magn Reson Med 2013;69(2):467–476, https://doi.org/10.1002/mrm.24272. The MEDI toolbox’s pipeline script unwraps Fit_ppm_complex’s field with it (README.m:80-95):

iMag = sqrt(sum(abs(iField).^2,4));
[iFreq_raw, N_std] = Fit_ppm_complex(iField);
iFreq = unwrapPhase(iMag, iFreq_raw, matrix_size);

Reference code

The original is the MEDI toolbox 2024.11.26’s functions/unwrapPhase.m and the C it calls, functions/mexUnwrap.c (“written by Tian Liu”, mex interface by Dong Zhou, last update 7.11.13). This module is a line-by-line port of that C, and reproduces it bit for bit: tests/medi_unwrap_cornell_parity.rs checks every voxel against the verified QSM-CI fit-ppm-complex-cornell container (ghcr.io/astewartau/qsm-ci/matlab-medi@sha256:7efa4bd6…, mexUnwrap compiled from that mexUnwrap.c with R2023b’s mex and identical to licensed MATLAB; QSMxT/QSM-CI#300), on phantoms built to exercise heap ties, the noise level, a region cut off from the seed and even and odd matrix sizes. From the echoes, the field-mapping stage with B0EstimationMethod::ComplexFit and this unwrapper gives that container’s total field to 8.9e-16 ppm (the units’ rounding). The fixture is sensitive to the heap: making heap_pop take the left child on a tie moves 541 of its 2970 voxels, and making heap_add sift up on a tie 1017. Every decision below is the C’s:

• Precision: single. mexUnwrap copies the magnitude and the phase from MATLAB’s doubles into float arrays (mexFunction), and every comparison and update is done on those floats. The output is the float phase written back as double, so every output value, including the voxels it never touches, is the input rounded to f32. This is kept (not double precision, unlike other ports here) because the growth order is decided by comparing f32 magnitudes: two magnitudes that differ in double but round to the same float are a tie in the original, and the tie decides which neighbour a voxel is unwrapped against.
• Quality map: the magnitude, nothing else. The heap is ordered by the magnitude value alone (heap_add / heap_pop compare .mag). MEDI passes iMag, the root-sum-of-squares of the echo magnitudes; the caller passes whatever magnitude it wants used (see UnwrappingAlgorithm::MediRegionGrowing for what the pipeline passes).
• Seed: the middle of the field of view, voxel (nx/2, ny/2, nz/2) in integer division (mexFunction: dim[2]/2*dim[1]*dim[0] + dim[1]/2*dim[0] + dim[0]/2, “choose the middle of the FOV as the seeding point”). Not the brightest voxel: the brightest voxel only sets the noise level. The seed is unwrapped (left as it is) whatever its magnitude, even below the noise level; if none of its neighbours is above the noise level, nothing else is unwrapped.
• Noise level: 0.005 × the largest magnitude, in single precision (noise_ratio = 5e-3, a float; noise_level = MRimage[max_ind].mag * noise_ratio). The largest magnitude is found with a strict > from voxel 0, so the first of equal maxima is taken (it only matters through its value). Only voxels with magnitude strictly above the noise level are ever planned; the rest keep their wrapped phase. MEDI_NOISE_RATIO. MEDI’s unwrapPhase has no way to change it, and neither does this function.
• No mask. unwrapPhase(iMag, iFreq_raw, matrix_size) takes none: the growth runs over the whole array, limited only by the noise level, so the phase outside the brain can decide the path into it. Voxels that the growth never reaches (below the noise level, or above it but cut off from the seed by voxels below it) are returned as they came in (rounded to f32).
• Neighbours: the 6-neighbourhood, in the order x−1, x+1, y−1, y+1, z−1, z+1 (get_neighbors), with x fastest in memory (MATLAB’s column-major order; QSM.rs’s layout). Neighbours outside the array are skipped.
• Growth. A binary max-heap of voxel indices keyed on magnitude (heap_create, heap_add, heap_pop). Each step takes the top voxel, unwraps it against the voxel it was planned from (its pre; the seed has none), pops it, and plans each neighbour that is unplanned and above the noise level, in neighbour order: its pre becomes the popped voxel, it is pushed, and it is marked planned. So a voxel is unwrapped against the first unwrapped neighbour that found it, not the best one, and pre is fixed when the voxel is planned, not when it is popped.
• Heap order and ties, exactly. heap_add sifts up while the parent’s magnitude is strictly less than the child’s (equal magnitudes stay below their parent). heap_pop moves the last element to the root, clears the vacated slot, and sifts down while the parent is strictly less than the chosen child, choosing the left child when its magnitude is strictly greater than the right’s and the right child otherwise (so the right one on a tie), and the left one when the right is past the end. Ties therefore resolve by heap position, not by voxel index or insertion time; reproducing them needs this exact heap, so no standard priority queue is used.
• The unwrap step. phase = phase − round((phase − phase[pre]) / (2π)) · 2π with the C’s types: the difference of the two floats is a float; it is divided by the double 2*PI (PI = 3.1415926535897931, the same double as std::f64::consts::PI); the quotient is passed to the file’s own static float round(float), i.e. converted to float, then floor(num + 0.5) in double (round half up, not half away from zero), returned as float; that float times the int 2 is a float, times PI a double, subtracted from the float phase in double, and the result stored as float.

The C also prints "%d voxels are unwrapped"; that is not reproduced.

Pipeline rules

complex-fit.unwrap-region-growing: with the complex fit (B0EstimationMethod::ComplexFit, MEDI’s Fit_ppm_complex), the unwrapper MEDI’s pipeline pairs it with is this one (README.m:80-95). It is a Kind::Reference rule, the authors’ setting, not a recommendation: nobody has shown it unwraps better than ROMEO, so a host shows it when the complex fit is chosen but never preselects it, and ROMEO stays the default unwrapper (decided 2026-11-04). It was left out when the port landed (2026-11-03) because QSMxT then preselected Apply on every reference value, which would have switched the unwrapper from ROMEO on one keypress. It is triggered by the complex fit rather than by MEDI, because the pairing is the field map’s (medi.field-complex-fit leads from MEDI to the fit).

Use

mexUnwrap.c assumes the input is wrapped; in MEDI it is Fit_ppm_complex’s field per echo spacing. The pipeline offers it as UnwrappingAlgorithm::MediRegionGrowing, on every field-mapping path; with B0EstimationMethod::ComplexFit it gives MEDI’s field map (QSM-CI’s fit-ppm-complex-cornell pair).

Structs

Heap 🔒

mexUnwrap.c’s binary max-heap of voxel indices, keyed on the magnitude. Kept as the C writes it (module docs, Heap order and ties): its tie behaviour decides the growth order.

Constants

MEDI_NOISE_RATIO

mexUnwrap.c’s noise_ratio: voxels with magnitude at or below this fraction of the largest magnitude are never unwrapped. A float in the C, as here.

PI 🔒

mexUnwrap.c’s #define PI 3.1415926535897931: that literal is this same double.

Functions

neighbours 🔒

get_neighbors: x−1, x+1, y−1, y+1, z−1, z+1, None past the array’s edge (the C’s −1).

round_c 🔒

The file’s static float round(float num) { return floor(num+0.5); }: num + 0.5 and the floor in double, the result converted to float.

unwrap_against 🔒

phase - round((phase - reference)/(2*PI))*2*PI with mexUnwrap.c’s types (module docs, The unwrap step).

unwrap_medi_region_growing

Unwrap phase (radians) by MEDI’s region growing, ordered by magnitude, as MEDI’s unwrapPhase(magnitude, phase, matrix_size) does. See the module docs for every detail of the original that is reproduced.