Function phase_offset_removal 

pub fn phase_offset_removal(
    phases: &[impl AsRef<[f64]>],
    mags: &[impl AsRef<[f64]>],
    tes: &[f64],
    sigma: [f64; 3],
    echoes: [usize; 2],
    unwrap_method: UnwrapMethod,
    grid: &Grid,
) -> PhaseOffsetRemoval

Remove phase offset from multi-echo phase data using HIP (Hermitian Inner Product): MCPC-3D-S, as MriResearchTools.jl’s mcpc3ds does it for one channel (src/mcpc3ds.jl:38-57, 771eadf; v3.8.0).

1. hip = conj(echo₁)·echo₂ over the whole volume, weight = sqrt(|hip|).
2. mask = robustmask(weight) (robust_mask); mcpc3ds makes its own mask and takes none from the caller.
3. The HIP phase unwrapped inside that mask (ROMEO with the default weights and weight as the magnitude, romeo(angle.(hip); mag=weight, mask)), scaled by TE₁/ΔTE to the phase that evolved by TE₁; the offset is φ₁ minus that, at every voxel.
4. The offset smoothed inside the same mask with MriResearchTools’ gaussiansmooth3d_phase (nan_box_smooth_3d_phase), and removed from every echo at every voxel (combinewithPO).

Where the offset is undefined

The masked smoothing leaves the offset NaN outside its mask, beyond the margin it extrapolates into, and so the reference’s corrected phase is NaN there too. ROMEO.jl then never grows through those voxels: every edge that touches a NaN phase has a NaN weight, which rescale makes invalid (ROMEO.jl src/weights.jl), so they stay NaN and so does their B0 (calculateB0_unwrapped). This function returns exactly that: NaN phase and offset there, and PhaseOffsetRemoval::valid marking the voxels that are defined.

Until 2026-11 QSM.rs set the offset to 0 there instead and kept the raw phase. That was QSM.rs’s own deviation and it was wrong: a brain mask is not inside MCPC-3D-S’s mask (scalp, the periphery and, on Philips, the scanner-masked background lie outside it), and the raw phase there carries the full coil-combination offset. On GE (real/imaginary, 7 echoes) and Philips (8 echoes) 3 T data those few voxels (0.6 % and 1.7 % of the brain mask) bridged the ROMEO growth, gave each echo a path-dependent number of turns and left ring artefacts in the field and in χ; given the same NaN→raw semantics, ROMEO.jl reproduces those rings almost voxel for voxel, and with its own NaN semantics it does not (χ correlation with Laplacian-unwrapped χ 0.21 → 0.78 on the GE scan; Q0748 audit, workstream C, ROMEO.jl 1.6.0 / MriResearchTools 3.8.0, 2026-11-04). On Siemens adaptive-combine data, with small offsets, those voxels were harmless. No option keeps the old behaviour.

Arguments

• phases - Wrapped phase per echo (n_echoes arrays of nxnynz)
• mags - Magnitude per echo
• tes - Echo times (any consistent unit)
• sigma - Gaussian smoothing sigma in voxels [x, y, z]
• echoes - Which two echoes to use for HIP [e1, e2]
• unwrap_method - Method to unwrap the HIP phase (mcpc3ds uses ROMEO; the others are QSM.rs options)
• grid - Volume grid (dimensions and voxel sizes)