Function fit_ppm_complex 

pub fn fit_ppm_complex(
    magnitudes: &[&[f64]],
    phases: &[&[f64]],
    echo_times: &[f64],
    params: &FitPpmComplexParams,
) -> Result<FitPpmComplexResult, InputError>

Complex multi-echo field fit, a port of Fit_ppm_complex.m (MEDI toolbox 2024.11.26; Liu et al., MRM 2013;69:467; Kressler et al., IEEE TMI 2010;29:273; de Rochefort et al., MRM 2008;60:1003).

Fits Sₑ = |Sₑ| exp(i(p0 + e·p1)) over echo indices e = 0, 1, … to the complex echoes by magnitude-weighted Gauss-Newton, starting from a linear fit to the temporally unwrapped phase of the first three echoes, and returns p1, the phase per echo spacing, wrapped to [-π, π]. MEDI then unwraps it spatially (unwrapPhase, a region-growing unwrapper) and divides by delta_TE; run_field_mapping does the same with QSM.rs’s own unwrappers when B0EstimationMethod::ComplexFit is chosen.

As in the reference:

• the echoes must be evenly spaced (echo_times is only used to check that, with echo_spacing_problem, and to report the spacing; the model is in echo index);
• the stopping test and the iteration count are global over the volume, not per voxel;
• a single echo is fitted as cat(4, abs(M), M), i.e. with zero phase at TE = 0, so the field is the echo’s phase over TE₁;
• a zero magnitude has phase 0 (MATLAB’s angle(0)), and NaN updates are set to 0.

Not ported, because a QSM.rs caller does not need them: the coil combination of a 5th (channel) dimension (pass coil-combined echoes) and the splitting of inputs over 1e9 bytes into slabs (memory only; the result is the same, apart from the global stopping test being evaluated per slab there).

Sign: MATLAB fits conj(M), so its p1 and p0 are minus FitPpmComplexResult::field and FitPpmComplexResult::phase0. Here the conjugate is taken as the first step and undone as the last, so every intermediate is MATLAB’s and the outputs follow QSM.rs’s convention (the sign of the input phase, as the ROMEO and Laplacian field maps).

Reference pipeline

MEDI’s pipeline script (README.m:80-103 of the 2024.11.26 release) computes iMag = sqrt(sum(abs(iField).^2,4)), then [iFreq_raw, N_std] = Fit_ppm_complex(iField) on the complex (coil-combined or per-channel) echoes with the default opts, unwraps iFreq_raw with unwrapPhase(iMag, iFreq_raw, matrix_size) (region-growing), and hands the field and N_std (both in radians per delta_TE) to PDF and MEDI_L1. For uneven echo spacing the toolbox has Fit_ppm_complex_TE(iField, TE) and for bipolar readouts Fit_ppm_complex_bipolar (unwrapped at twice the phase and halved, README.m:98-102); neither is ported. Nothing in the fit depends on the field strength.

magnitudes[e] and phases[e] are echo e (phase in radians, wrapped or not).