Function multi_echo_robust_fit 

pub fn multi_echo_robust_fit(
    unwrapped_phases: &[impl AsRef<[f64]>],
    mags: &[impl AsRef<[f64]>],
    tes: &[f64],
    mask: &[u8],
    params: &RobustFitParams,
) -> RobustFitResult

Multi-echo fit that downweights echoes which disagree with the rest, and reports which.

Iteratively reweighted least squares around multi_echo_linear_fit’s model (phase = intercept + slope·TE, magnitude-weighted), with a robust loss applied to each echo’s residual. An echo that is corrupted — through-plane motion, a spin-history dropout, a shot artefact, anything that is not a straight line in TE — is pulled out of the fit at the voxels where it misbehaves, instead of dragging the slope at every voxel it touches.

Supported by Biondetti et al., MRM 2022, which found weighted and nonlinear-fit echo combination clearly outperforms unweighted averaging or plain linear fitting.

The scale is global, and it has to be

IRLS needs a scale to measure “large” against. The textbook choice is a per-observation-set MAD, which here means the MAD of one voxel’s n_echoes residuals — and that is unusable on this data. A four-echo fit with an intercept has two degrees of freedom, so a per-voxel MAD is the MAD of four numbers of which two are structurally near zero; it collapses toward zero at a large fraction of voxels, and dividing by it rejects healthy echoes essentially at random.

So the scale is estimated once per iteration over the whole mask, from residuals standardised by magnitude. Phase noise goes as 1/SNR, so residual × magnitude is roughly homoscedastic across the volume and a single scale is meaningful for all of it. That is also what makes the per-echo fractions in EchoQuality comparable between echoes: every echo is being measured against the same ruler.

What it does not detect

Most of what follows is measured at four echoes, which is this function’s own minimum. With an intercept the model has two parameters, so four is the fewest that leaves anything to detect with, and the tests sit exactly on that boundary. More echoes improve both halves markedly — on the same corruption, going from four to six tightens the clean-series score spread from 1.88 to 1.34 (more headroom under the flagging threshold) and takes the slope correction from 9% to 28% of the plain fit’s error. Read the limitations below as the four-echo case, not as the method.

It reports more reliably than it repairs. Detection and correction are separate capabilities with different sensitivities, and the gap is worth knowing about. On the phantom, a dropout over 12% of the volume is flagged (score ratio 9.8) while the field map it produces is neither better nor worse than the plain fit (0.0259 against 0.0246). The same is true of corruption severe enough to pull the per-voxel fit at the voxels it affects: the echo is named every time, but the weighting recovers progressively less of the field as the corruption grows. Treat a flag as “look at this echo”, not as “this has been handled”.

Corruption of the first or last echo, as readily as a middle one. The endpoints of a linear fit carry the most leverage, so the fit follows them and their residuals shrink. An echo in the middle of the train is the easiest to catch and the first and last the hardest.

A series in which an echo carries no signal. This wants n_params + 2 echoes carrying signal, not merely present, and it cannot tell the difference: a near-dead echo scores ~0 because the score is magnitude-weighted, and the remaining three points then put all their leverage on the middle echo, which gets flagged with nothing wrong with it. Drop a dead echo before believing this. a_dead_echo_causes_a_known_false_positive pins the behaviour and records the two fixes that do not work.

Anything that is still linear in TE. A per-echo scaling of the whole field, or a corruption that happens to look like a different slope, fits the model perfectly well and leaves no residual to find. What this detects is departure from a straight line, which is what motion, dropout and shot artefacts produce — not every way an echo can be wrong.

How many echoes this needs

Robustness costs degrees of freedom that a short echo train does not have. With estimate_offset the model has two parameters, so n_echoes must be at least 4 for a corrupted echo to be distinguishable from the fit; without it, at least 3. Below that there is nothing to detect with and this returns the ordinary magnitude-weighted fit unchanged, with every quality figure zeroed and nothing flagged — not an error, because a short series is not a mistake, but the caller should not read a clean report as evidence of a clean series. Check EchoQuality::len against the echo count if that distinction matters.

Even at 4 echoes this is working near its limit: one corrupted echo is 25% of the data, and the fit it is being measured against is itself pulled by that echo. Expect it to catch an echo that is badly wrong and to be unreliable about one that is marginally wrong.

Arguments

Identical to multi_echo_linear_fit except that the two fit options are carried in params: unwrapped phase and magnitude per echo, echo times in seconds, a mask, and RobustFitParams.