Function bipolar_correction 

pub fn bipolar_correction<P: AsMut<[f64]> + AsRef<[f64]>>(
    phases: &mut [P],
    mags: &[impl AsRef<[f64]>],
    tes: &[f64],
    mask: &[u8],
    sigma: [f64; 3],
    grid: &Grid,
)

Bipolar gradient correction for multi-echo phase data.

Removes linear phase artefact caused by bipolar readout gradients. Requires at least 3 echoes.

Reference: Eckstein PhD thesis (2021), Section 3.1.3 https://doi.org/10.34726/hss.2021.43447

MriResearchTools’ bipolar_correction! (src/mcpc3ds.jl:97-129, 771eadf), step for step:

1. artefact: with k = (TE₁ + TE₃) / TE₂, fG = φ₁ + φ₃ − k·φ₂, where φ₂ is first unwrapped by romeo(φ₂; mag = mag₂) when k is within 0.01 of an integer. Over the whole volume: no mask, and fG is not rewrapped.
2. gaussiansmooth3d_phase(fG, sigma; mask): MriResearchTools’ masked box-filter smoothing of exp(i·fG) (nan_box_smooth_3d_phase), NaN where it cannot extrapolate.
3. romeo!(fG; mag = mag₁, correctglobal = true): ROMEO on the whole volume, again with no mask (edges touching a NaN are invalid there and here), then the global n·2π correction over every finite voxel.
4. remove_artefact!: echo e loses t·fG (t = m/f for odd, (m + 1)/f for even echoes; m = TE₁/(TE₂ − TE₁), f = (2 − k)m − k) and is rewrapped (rem2pi).

Deviation: MriResearchTools writes NaN wherever fG is NaN (far from the mask); here such voxels keep their input phase. Everywhere fG is finite the result is the original’s.

Agreement (astewartau/QSM.rs#319): until 2026-10-30 this smoothed fG with a true Gaussian kernel (radius 3σ) instead of MriResearchTools’ 4-pass NaN-aware box filter, unwrapped φ₂ and fG inside the mask only, and corrected only inside the mask. On the 4-echo, 20×18×14, 1×1×1.5 mm phantom of tests/mritools_parity.rs (verified QSM-CI julia-qsm image, QSMxT/QSM-CI#299) the artefact estimate was up to 0.48 rad from the original’s (corrected phase 0.24 rad on odd, 0.48 rad on even echoes). Now it agrees to 4.4e-16 / 8.9e-16 rad (measured 2026-10-30). The smoothing was the whole gap there: with the new smoothing but ROMEO still masked, the phantom also agrees to 8.9e-16 rad. The unmasked ROMEO is ported because the original does it; it matters only where the growth inside the mask depends on paths outside it.

Arguments

• phases - Mutable phase data per echo (modified in-place)
• mags - Magnitude data per echo (echo 1’s and echo 2’s guide ROMEO; empty for none)
• tes - Echo times (any consistent unit; only ratios are used)
• mask - Binary mask: where the artefact is smoothed (MCPC-3D-S’s robustmask of √|HIP| in MriResearchTools’ mcpc3ds)
• sigma - Smoothing sigma in voxels for the artefact estimation
• grid - Volume grid (dimensions and voxel sizes)