Function enforce_inter_echo_consistency 

pub fn enforce_inter_echo_consistency(
    unwrapped: &mut [Vec<f64>],
    wrapped: &[&[f64]],
    mask: &[u8],
)

Re-seat every echo on the first, so the echoes of a voxel differ by their measured phase evolution and nothing else.

Slice-wise unwrapping leaves echo e as phi + 2*pi*k, and the field is recoverable from a fit across echoes only if k does not depend on e. This makes it so. It does not, and cannot, say what k is — the leftover wrap stays in the intercept, along with the physical slice offset, and the fit discards both.

The relative wrap is pinned by the wrapped input, without unwrapping anything and without knowing the slice offset. Writing p for the wrapped phase and u for the unwrapped,

u_e - u_{e-1} - wrap(p_e - p_{e-1}) = 2*pi * (an integer)

because wrap(p_e - p_{e-1}) is the echo-to-echo evolution, wrapped. Rounding that to the nearest multiple of 2π and subtracting it leaves u_e - u_{e-1} equal to the measured step exactly. Formulating it on the wrapped difference rather than against a TE-scaled reference is what lets it survive an echo-independent offset, which the TE-scaled form of correct_multi_echo_wraps would read as field.

This is per voxel, deliberately

An earlier version of this did it per connected component of each slice, taking a median, on the reasoning that slice-wise unwrapping offsets a whole component at a time. That reasoning is wrong in practice: a region grower that cannot route through the slice axis has to cross every in-plane fringe head-on, and where one is ambiguous it leaves a part of the component offset by 2π. Measured on the qsm-forward 2D phantom, the wrap count varies within a single slice of a single echo (2-3 distinct values), so a single median per component has nothing it can do. Correcting per voxel repairs those plateaus as well, and took the fitted field from r = -0.44 to r = 0.98 against ground truth during development. The per-slice variant was not kept, so that figure records the decision rather than a test.

This is the half that fixes a linear-fit B0 map

Public, and deliberately not tied to slice-wise unwrapping: it repairs any unwrapper’s output, including none at all. On the phantom’s interleaved-offset session it takes 3D ROMEO from 0.20 to 0.983 and the raw wrapped phase from -0.01 to 0.983. What it does not do is make the phase spatially continuous — it propagates the first echo’s wrap state faithfully into the rest, which on that session leaves more in-plane 2π jumps than not unwrapping at all. See the module docs for the division of labour.

What it assumes

That the true phase evolution between consecutive echoes stays inside ±π, which is the usual design constraint on echo spacing and the same assumption correct_multi_echo_wraps makes. A voxel whose evolution exceeds it is not recoverable from this data by any means that only sees the wrapped steps, and will carry a slope error of 2*pi / delta_TE.

Arguments

• unwrapped - Unwrapped phase per echo, modified in place
• wrapped - The wrapped input phase per echo
• mask - Binary mask (nx * ny * nz), 1 = inside ROI; voxels outside are left alone