qsm_core/inversion/
iqsm_echo_weighting.rs

1//! Echo weighting for the multi-echo iQSM, iQFM and iQSM+ combinations.
2//!
3//! Not behind the `onnx` feature, so callers can name it (e.g. in configuration) whatever features
4//! they build with; the networks themselves need `onnx`.
5
6/// How the multi-echo iQSM and iQFM entry points combine their per-echo reconstructions:
7/// `Σ w·χ / Σ w` over the echoes, voxel by voxel, with the weight `w` below. Without magnitudes
8/// both are `TE²` (uniform magnitude), as in both of the authors' implementations.
9///
10/// The authors' two implementations disagree (checked at `sunhongfu/iQSM` 81f5462 and
11/// `sunhongfu/iQSM_Plus` 3ba25f5):
12///
13/// - **MATLAB** (published with the paper): `sunhongfu/iQSM` `legacy/matlab/demo_multi_echo.m:152-160`
14///   and `sunhongfu/iQSM_Plus` `legacy/matlab/iQSM.m:260-267` / `iQFM.m:263-270` multiply each
15///   echo's output by its TE and fit that against TE with `echofit.m` (`sunhongfu/iQSM`
16///   `legacy/matlab/iQSM_fcns/echofit.m`, zero intercept): `Σ mag·TE·(TE·χ) / Σ mag·TE²`, i.e.
17///   weights `mag·TE²`.
18/// - **Python** (the current drivers): `sunhongfu/iQSM` `run.py:170` and `app.py:416` weight by
19///   `(mag·TE)²`, i.e. magnitude squared, although they print "Magnitude × TE²". QSM-CI's iQSM and
20///   iQFM originals (bb163ef) run this Python.
21///
22/// The default is [`MagnitudeSquaredTeSquared`](Self::MagnitudeSquaredTeSquared), `(mag·TE)²`
23/// (Ashley's decision, 2026-10-31). Each echo's χ (or local field) estimate is its phase divided
24/// by TE, and the phase noise has standard deviation ∝ 1/mag, so the estimate's noise variance is
25/// ∝ 1/(mag·TE)²: `(mag·TE)²` is the inverse-variance weight, which gives the combination of least
26/// variance. It is also the authors' current Python, and iQSM+'s weighting in both of its
27/// authors' implementations (`sunhongfu/iQSM_Plus` `legacy/matlab/iQSM_plus.m:273-275`,
28/// `run.py:199`; see `iqsm_plus_multi_echo`, which has no choice). The MATLAB's `mag·TE²`
29/// ([`MagnitudeTeSquared`](Self::MagnitudeTeSquared)) stays as an option, for parity with the
30/// implementation published with the paper.
31#[derive(Clone, Copy, Debug, Default, PartialEq, Eq)]
32pub enum IqsmEchoWeighting {
33    /// `w = mag·TE²`: the authors' MATLAB echo fit (an option, for parity with it).
34    MagnitudeTeSquared,
35    /// `w = (mag·TE)² = mag²·TE²`: inverse-variance; the authors' Python drivers (the default).
36    #[default]
37    MagnitudeSquaredTeSquared,
38}
39
40impl IqsmEchoWeighting {
41    /// The weight of one echo with magnitude `mag` and echo time `te`.
42    pub fn weight(self, mag: f64, te: f64) -> f64 {
43        match self {
44            Self::MagnitudeTeSquared => mag * te * te,
45            Self::MagnitudeSquaredTeSquared => (mag * te) * (mag * te),
46        }
47    }
48}
49
50/// Reconstruct each echo with `echo(phase, te)` and combine them with `weighting`'s weights.
51#[cfg_attr(not(feature = "onnx"), allow(dead_code))]
52pub(crate) fn combine_echoes<E>(
53    phases: &[&[f64]],
54    magnitudes: &[&[f64]],
55    tes: &[f64],
56    n: usize,
57    weighting: IqsmEchoWeighting,
58    mut echo: impl FnMut(&[f64], f64) -> Result<Vec<f64>, E>,
59) -> Result<Vec<f64>, E> {
60    let mut acc = vec![0.0f64; n];
61    let mut wsum = vec![0.0f64; n];
62    for (e, &phase) in phases.iter().enumerate() {
63        let chi = echo(phase, tes[e])?;
64        for i in 0..n {
65            let w = weighting.weight(magnitudes.get(e).map(|m| m[i]).unwrap_or(1.0), tes[e]);
66            acc[i] += w * chi[i];
67            wsum[i] += w;
68        }
69    }
70    for i in 0..n {
71        acc[i] = if wsum[i] > 0.0 { acc[i] / wsum[i] } else { 0.0 };
72    }
73    Ok(acc)
74}
75
76#[cfg(test)]
77mod tests {
78    use super::*;
79
80    /// The combination's exact weights on two voxels and two echoes, each echo's "reconstruction"
81    /// a constant per echo, so the result is the weighted mean of those constants.
82    #[test]
83    fn echo_weightings_are_the_authors() {
84        let tes = [0.01, 0.03];
85        let (m1, m2) = ([2.0, 1.0], [1.0, 0.0]);
86        let mags: [&[f64]; 2] = [&m1, &m2];
87        let phases: [&[f64]; 2] = [&[0.0; 2], &[0.0; 2]];
88        let chi = |te: f64| if te < 0.02 { 1.0 } else { 5.0 };
89        let run = |w, mags: &[&[f64]]| {
90            combine_echoes(&phases, mags, &tes, 2, w, |_, te| Ok::<_, ()>(vec![chi(te); 2])).unwrap()
91        };
92        let mean = |w1: f64, w2: f64| (w1 * 1.0 + w2 * 5.0) / (w1 + w2);
93        let close = |a: &[f64], b: &[f64]| a.iter().zip(b).all(|(x, y)| (x - y).abs() < 1e-12);
94
95        // MATLAB echo fit: mag·TE².
96        let matlab = run(IqsmEchoWeighting::MagnitudeTeSquared, &mags);
97        assert!(close(&matlab, &[mean(2.0 * 1e-4, 1.0 * 9e-4), 1.0]), "{matlab:?}");
98        // Python: (mag·TE)².
99        let python = run(IqsmEchoWeighting::MagnitudeSquaredTeSquared, &mags);
100        assert!(close(&python, &[mean(4.0 * 1e-4, 1.0 * 9e-4), 1.0]), "{python:?}");
101        // Without magnitudes both are TE².
102        for w in [IqsmEchoWeighting::MagnitudeTeSquared, IqsmEchoWeighting::MagnitudeSquaredTeSquared] {
103            let r = run(w, &[]);
104            assert!(close(&r, &[mean(1e-4, 9e-4); 2]), "{w:?}: {r:?}");
105        }
106        // All weights zero: zero.
107        let zero = run(IqsmEchoWeighting::MagnitudeTeSquared, &[&[0.0; 2], &[0.0; 2]]);
108        assert_eq!(zero, vec![0.0; 2]);
109        assert_eq!(IqsmEchoWeighting::default(), IqsmEchoWeighting::MagnitudeSquaredTeSquared);
110    }
111}
qsm_core/inversion/iqsm_echo_weighting.rs

qsm_core/inversion/
iqsm_echo_weighting.rs
