qsm_core/separation/
r2star_qsm.rs

1//! R2*-QSM: susceptibility source separation from gradient-echo data alone.
2//!
3//! Dimov et al. separate paramagnetic (χ+, iron) and diamagnetic (χ−, myelin)
4//! susceptibility from GRE data alone — R2* from the multi-echo magnitude plus a
5//! conventional QSM (χ_total) from the phase — with **no** separate R2/R2'
6//! measurement. Per voxel the model is
7//!
8//! ```text
9//!   R2*      = 𝓇·(|χ+| + |χ−|)      (a single relaxometric constant for both sources)
10//!   χ_total  = χ+ + χ−               (χ+ ≥ 0, χ− ≤ 0)
11//! ```
12//!
13//! which inverts in closed form to
14//!
15//! ```text
16//!   χ+   = (χ_total + R2*/𝓇) / 2
17//!   |χ−| = (R2*/𝓇 − χ_total) / 2
18//! ```
19//!
20//! with the physical constraints `χ+ ≥ 0`, `|χ−| ≥ 0` enforced by clipping.
21//!
22//! **Relaxometric constant.** Dimov calibrated `𝓇 = 274 Hz/ppm` at 3 T (a single
23//! value for both sources). R2* susceptibility-induced decay scales linearly with
24//! B0 while χ does not, so at the acquisition field `𝓇 = 274·(B0/3)` Hz/ppm.
25//!
26//! This is the paper's voxel-level solve; the full method adds a field
27//! data-consistency term (auto-satisfied when χ_total is supplied) and an
28//! edge-masked L1 regularisation, both omitted here.
29//!
30//! Reference:
31//! Dimov, A.V., et al. (2022). "Magnetic susceptibility source separation solely
32//! from gradient echo data: histological validation." Tomography, 8(3):1544-1551
33//! (and J. Neuroimaging 2022, https://doi.org/10.1111/jon.13014). R2*
34//! susceptibility model: Yablonskiy & Haacke, Magn. Reson. Med. 1994.
35
36/// Parameters for [`r2star_qsm`] / [`r2star_qsm_from_magnitude`].
37#[cfg_attr(feature = "introspection", derive(serde::Serialize))]
38#[derive(Clone, Debug)]
39pub struct R2starQsmParams {
40    /// Main field strength in Tesla (the constant is scaled by `B0/3`).
41    pub b0: f64,
42    /// Relaxometric constant at 3 T in Hz/ppm (Dimov et al.: 274).
43    pub r_const_3t: f64,
44}
45
46impl Default for R2starQsmParams {
47    fn default() -> Self {
48        Self { b0: 3.0, r_const_3t: 274.0 }
49    }
50}
51
52/// Closed-form R2*-QSM separation from a QSM and an R2* map.
53///
54/// # Arguments
55/// * `chi_total` — Conventional QSM χ_total in **ppm** (`n` voxels).
56/// * `r2star` — R2* map in **Hz** (`n` voxels), e.g. from [`crate::r2star::r2star_arlo`].
57/// * `mask` — Binary brain mask (`n` voxels, 1 = inside).
58/// * `params` — See [`R2starQsmParams`].
59///
60/// # Returns
61/// `(chi_pos, chi_neg, chi_total)` in ppm, restricted to `mask` — matching the
62/// [`chi_sep_ilsqr`](super::chi_sep_ilsqr)/[`chi_sep_medi`](super::chi_sep_medi)
63/// convention: `chi_pos` ≥ 0 (paramagnetic), `chi_neg` ≤ 0 (diamagnetic, signed),
64/// and `chi_total = chi_pos + chi_neg`.
65pub fn r2star_qsm(
66    chi_total: &[f64],
67    r2star: &[f64],
68    mask: &[u8],
69    params: &R2starQsmParams,
70) -> (Vec<f64>, Vec<f64>, Vec<f64>) {
71    let n = chi_total.len();
72    assert_eq!(r2star.len(), n, "r2star length must match chi_total");
73    assert_eq!(mask.len(), n, "mask length must match chi_total");
74
75    // Field-scaled relaxometric constant (Hz/ppm).
76    let r = params.r_const_3t * (params.b0 / 3.0);
77
78    let mut chi_pos = vec![0.0_f64; n];
79    let mut chi_neg = vec![0.0_f64; n];
80    let mut chi_out = vec![0.0_f64; n];
81    for i in 0..n {
82        if mask[i] == 0 {
83            continue;
84        }
85        let chi = chi_total[i];
86        let s = r2star[i] / r; // R2*/𝓇 = |χ+| + |χ−|
87        let pos = ((chi + s) / 2.0).max(0.0); // χ+
88        let dia_mag = ((s - chi) / 2.0).max(0.0); // |χ−|
89        chi_pos[i] = pos;
90        chi_neg[i] = -dia_mag; // signed χ− ≤ 0
91        chi_out[i] = pos - dia_mag; // χ_total = χ+ + χ−
92    }
93    (chi_pos, chi_neg, chi_out)
94}
95
96/// R2*-QSM separation fitting R2* from multi-echo magnitude, then separating.
97///
98/// Fits R2* by magnitude-weighted log-linear least squares (weight ∝ magnitude²),
99/// matching the QSM-CI reference implementation, then calls [`r2star_qsm`].
100///
101/// # Arguments
102/// * `magnitude` — Multi-echo magnitude, flattened as `(n_voxels, n_echoes)` in
103///   row-major order (echo fastest per voxel) — the same layout as
104///   [`crate::r2star::r2star_arlo`].
105/// * `echo_times` — Echo times in **seconds** (`n_echoes`).
106/// * `chi_total` — Conventional QSM χ_total in **ppm** (`n_voxels`).
107/// * `mask` — Binary brain mask (`n_voxels`).
108/// * `params` — See [`R2starQsmParams`].
109///
110/// # Returns
111/// `(chi_pos, chi_neg, chi_total)` as in [`r2star_qsm`].
112pub fn r2star_qsm_from_magnitude(
113    magnitude: &[f64],
114    echo_times: &[f64],
115    chi_total: &[f64],
116    mask: &[u8],
117    params: &R2starQsmParams,
118) -> (Vec<f64>, Vec<f64>, Vec<f64>) {
119    let n = chi_total.len();
120    let ne = echo_times.len();
121    assert_eq!(magnitude.len(), n * ne, "magnitude must be n_voxels * n_echoes");
122    let r2star = fit_r2star_weighted_loglinear(magnitude, echo_times, mask, n);
123    r2star_qsm(chi_total, &r2star, mask, params)
124}
125
126/// R2* (Hz) by magnitude-weighted log-linear least squares over echoes.
127///
128/// Weight ∝ magnitude² (SNR) so noisy late echoes don't dominate — the QSM-CI
129/// reference's robust stand-in for ARLO. `R2* = −slope of log(magnitude) vs TE`,
130/// clamped to ≥ 0.
131fn fit_r2star_weighted_loglinear(
132    magnitude: &[f64],
133    echo_times: &[f64],
134    mask: &[u8],
135    n_voxels: usize,
136) -> Vec<f64> {
137    let ne = echo_times.len();
138    let mut r2s = vec![0.0_f64; n_voxels];
139    for v in 0..n_voxels {
140        if mask[v] == 0 {
141            continue;
142        }
143        // Weighted means over echoes.
144        let mut sw = 0.0;
145        let mut swt = 0.0;
146        let mut swl = 0.0;
147        for e in 0..ne {
148            let m = magnitude[v * ne + e].max(1e-9);
149            let w = m * m;
150            sw += w;
151            swt += w * echo_times[e];
152            swl += w * m.ln();
153        }
154        if sw <= 0.0 {
155            continue;
156        }
157        let tbar = swt / sw;
158        let lbar = swl / sw;
159        // Weighted covariance / variance.
160        let mut cov = 0.0;
161        let mut var = 0.0;
162        for e in 0..ne {
163            let m = magnitude[v * ne + e].max(1e-9);
164            let w = m * m;
165            let dt = echo_times[e] - tbar;
166            cov += w * dt * (m.ln() - lbar);
167            var += w * dt * dt;
168        }
169        if var > 0.0 {
170            r2s[v] = (-cov / var).max(0.0);
171        }
172    }
173    r2s
174}
175
176#[cfg(test)]
177mod tests {
178    use super::*;
179
180    #[test]
181    fn closed_form_recovers_sources() {
182        // Build χ_total and R2* from known χ+ / |χ−| via the forward model, then
183        // check the closed form inverts them exactly.
184        let params = R2starQsmParams { b0: 3.0, r_const_3t: 274.0 };
185        let r = params.r_const_3t; // B0 = 3 -> r = 274
186        let chi_pos = vec![0.10, 0.05, 0.00, 0.20];
187        let chi_neg_mag = vec![0.02, 0.10, 0.15, 0.00]; // |χ−|
188        let n = chi_pos.len();
189        let chi_total: Vec<f64> = (0..n).map(|i| chi_pos[i] - chi_neg_mag[i]).collect();
190        let r2star: Vec<f64> = (0..n).map(|i| r * (chi_pos[i] + chi_neg_mag[i])).collect();
191        let mask = vec![1u8; n];
192
193        let (para, neg, total) = r2star_qsm(&chi_total, &r2star, &mask, &params);
194        for i in 0..n {
195            assert!((para[i] - chi_pos[i]).abs() < 1e-9, "χ+ voxel {i}");
196            assert!((-neg[i] - chi_neg_mag[i]).abs() < 1e-9, "|χ−| voxel {i}");
197            assert!((total[i] - chi_total[i]).abs() < 1e-9, "χ_total voxel {i}");
198            assert!(neg[i] <= 0.0, "χ− must be ≤ 0 at voxel {i}");
199        }
200    }
201
202    #[test]
203    fn field_scaling_and_masking() {
204        // At B0 = 6, r doubles, so the same R2* yields half the source magnitude.
205        let p3 = R2starQsmParams { b0: 3.0, r_const_3t: 274.0 };
206        let p6 = R2starQsmParams { b0: 6.0, r_const_3t: 274.0 };
207        let chi_total = vec![0.0, 0.0];
208        let r2star = vec![274.0, 274.0];
209        let mask = vec![1u8, 0u8];
210
211        let (para3, neg3, _t3) = r2star_qsm(&chi_total, &r2star, &mask, &p3);
212        let (para6, _neg6, _t6) = r2star_qsm(&chi_total, &r2star, &mask, &p6);
213        // voxel 0: χ_total=0, R2*/r = 1 at 3T -> χ+ = |χ−| = 0.5 (χ− = −0.5)
214        assert!((para3[0] - 0.5).abs() < 1e-9);
215        assert!((neg3[0] + 0.5).abs() < 1e-9);
216        // at 6T r = 548 -> R2*/r = 0.5 -> χ+ = 0.25
217        assert!((para6[0] - 0.25).abs() < 1e-9);
218        // masked-out voxel stays zero
219        assert_eq!(para3[1], 0.0);
220        assert_eq!(neg3[1], 0.0);
221    }
222
223    #[test]
224    fn r2star_fit_from_magnitude() {
225        // Mono-exponential decay with known R2*, single voxel.
226        let r2star_true = 40.0_f64;
227        let s0 = 1000.0_f64;
228        let te = vec![0.004, 0.012, 0.020, 0.028];
229        let mag: Vec<f64> = te.iter().map(|&t| s0 * (-r2star_true * t).exp()).collect();
230        let mask = vec![1u8];
231        let r2s = fit_r2star_weighted_loglinear(&mag, &te, &mask, 1);
232        assert!((r2s[0] - r2star_true).abs() / r2star_true < 1e-6, "R2* fit {}", r2s[0]);
233    }
234}
qsm_core/separation/r2star_qsm.rs

qsm_core/separation/
r2star_qsm.rs
