Function field_mapping_direct 

fn field_mapping_direct(
    phases: &[&[f64]],
    mag_slices: &[&[f64]],
    mask: &[u8],
    metadata: &ScanMetadata,
    config: &FieldMappingConfig,
    _n_voxels: usize,
    nx: usize,
    ny: usize,
    nz: usize,
    vsx: f64,
    vsy: f64,
    vsz: f64,
    progress: &mut dyn FnMut(usize, usize),
) -> Result<FieldMappingResult, PipelineError>

Path B: Multi-echo without phase offset removal (per-echo unwrap + B0 estimation).

Taken whenever offset removal is off: always with Laplacian unwrapping, and with ROMEO when phase_offset_removal is false. Each echo is unwrapped on its own (for Laplacian, with the phase zeroed outside mask first unless config.laplacian_unmasked, and config.laplacian_solver as the Poisson solver), and the echoes are then combined with config.b0_estimation, as on the offset-removal path:

• B0EstimationMethod::WeightedAvg (the default): each echo is unwrapped only up to a constant of its own, which a weighted mean of φ/TE would turn into a field offset (one that varies over the brain for magnitude-dependent weights). So that constant is aligned across echoes first, by the least that the unwrapper leaves undetermined:

  • Laplacian fixes each echo only up to an arbitrary constant, so the masked mean is subtracted from every echo. Echo e is TEₑ·ω + cₑ inside the mask; demeaned it is TEₑ·(ω − ω̄) for every echo, which loses only a global constant field ω̄ (removed by background removal and referencing anyway).
  • ROMEO (unwrap_romeo_multi_echo with romeo_params, as on the offset-removal path) fixes each echo up to a whole number of turns, so each echo’s median wrap count relative to the TE-scaled previous echo is removed (correct_multi_echo_wraps); this is a no-op when the echoes are already consistent (temporal unwrapping, or individual unwrapping with correct_global, which already ends with it). Kept for individual mode without global correction, where the authors’ calculateB0_unwrapped would average echoes a turn apart.

  Then B0 = Σ w (φ/TE) / Σ w with the variant’s weighting.

  This alignment is QSM.rs’s own (kept, Ashley 2026-10-28): each echo’s Laplacian unwrap is determined only up to a constant of its own, so the constants are removed before the echoes are combined, which loses only a global constant field. No reference combines Laplacian-unwrapped echoes this way. The MATLAB pipelines QSM.rs reproduces with other settings: UK Biobank (STI Suite MRPhaseUnwrap on the masked phase, then the T2*-weighted average; LaplacianSolver::Fft { pad: [64; 3] } with B0WeightType::AssumedDecay, matched to ≈1e-14 rad up to a global constant, see run_field_mapping), QSMART (Hongfu Sun’s lapunwrap per echo on the raw phase, then echofit; LaplacianSolver::Mirror with laplacian_unmasked and B0EstimationMethod::QsmartEchofit, which aligns nothing) and STI Suite’s GUI (MRPhaseUnwrap on the raw phase, echoes summed: laplacian_unmasked with B0WeightType::TEs). QSM.jl’s README unwraps the raw phase and inverts each echo’s field separately, so it has no echo combination to compare with. Measured on the QSM-CI in-vivo travelling-head scan (Prisma 2.89 T, 1 mm, 5 echoes) against QSM.jl’s unwrap_laplacian(…; solver = :dct) with each echo scaled by its TE and combined by QSM.jl’s multi_echo_average (W = magnitude), in the QSM-CI julia-qsm image (sha256:b3c048ac…): with the phase masked and each echo’s masked mean removed on the QSM.jl side too, the field maps agree to NRMSE 1.4e-9 % (float32 output rounding); with QSM.jl’s raw phase and no mean removal, NRMSE 46 % against this default and 28 % against laplacian_unmasked: true, the rest being the per-echo constants this alignment removes.

• B0EstimationMethod::LinearFit: the magnitude-weighted fit of phase against TE with the variant’s LinearFitParams. With ROMEO it runs on the per-echo unwrapped phases as they come out of the unwrapper, which is what this path computed unconditionally before it honoured b0_estimation (bit-identical). With Laplacian the per-echo constants are aligned first, as for the weighted mean: a constant cₑ that differs between echoes is not linear in TE, so the fit turns it into a slope error, and because the fit is magnitude-weighted that error varies over the brain rather than being a constant that background removal would take out. With an intercept and two echoes the fit is the echo difference (φ₂ − φ₁)/(TE₂ − TE₁), which discards the absolute phase, so it is only worth choosing when the echoes carry a phase offset that has not been removed.

• B0EstimationMethod::RobustFit: the per-echo constants are aligned exactly as for the weighted mean (an unaligned constant is a departure from linearity in TE that the robust loss would blame on the echo), then the robust fit runs on the aligned echoes. Keep estimate_offset on unless the echoes are known to be offset-free (e.g. MCPC-3D-S combined): the offset removal that would otherwise take the intercept out does not run here.

• B0EstimationMethod::QsmartEchofit: QSMART’s echofit.m on the echoes as unwrapped. With Laplacian nothing is aligned, as in QSMART (whose unwrap_phase → echofit removes no per-echo constant), and with laplacian_unmasked every echo is kept over the whole volume, because echofit.m fits every voxel and takes its R_0 threshold from all of them. With ROMEO, whole turns are aligned as for the weighted mean: a fit through the origin turns a turn between echoes into a slope error.