qsm_core/separation/mod.rs
1//! Susceptibility source separation
2//!
3//! This module provides algorithms for separating total magnetic susceptibility
4//! into paramagnetic (χ+, primarily iron) and diamagnetic (χ-, primarily myelin)
5//! components using local field maps and R2' relaxation data.
6//!
7//! # Methods
8//! - `chi_sep_ilsqr`: the SNU χ-separation toolbox's `chi_sep_iLSQR` (LSQR, streaking removal
9//! and NNLS from a QSM)
10//! - `chi_sep_medi`: the SNU χ-separation toolbox's `chi_sep_MEDI` (edge-masked L1 Gauss-Newton
11//! from a QSM, with a CSF term)
12//! - `r2star_qsm`: closed-form separation from a QSM + R2* (Dimov 2022, no R2' needed)
13//! - `wavesep`: wavelet-L1 proximal-gradient separation from a QSM + R2' (Fang 2023)
14//! - `decompose`: signal-domain 3-compartment per-voxel fit from a QSM + multi-echo
15//! magnitude (Chen 2021)
16//! - `hc_chisep`: hollow-cylinder χ-separation with signal-derived fiber orientation
17//! from a QSM + R2' + multi-echo magnitude (Wharton & Bowtell model)
18//!
19//! # Reference
20//! Shin, H., et al. (2021). "χ-separation: Magnetic susceptibility source separation
21//! toward iron and myelin mapping in the brain." NeuroImage, 240:118371.
22
23/// χ-separation's relaxometric constant D_r, in Hz/ppm, as published by Shin et al. 2021.
24///
25/// The forward model is `R2' = D_r,pos·|χ+| + D_r,neg·|χ−|`, and both published derivations give
26/// the two sources the *same* relaxivity — neither paper fits them separately. 137 Hz/ppm comes
27/// from linear regression of R2' against **single-orientation** QSM in iron-rich deep grey matter
28/// (Shin et al. 2021, NeuroImage 240:118371), which is the reconstruction every analytic method
29/// in this module is given, so it is their default.
30///
31/// [`DR_KIM_2025_COSMOS`] is the same quantity re-measured against COSMOS. Which is correct for a
32/// given dataset is an open question in the field; what is not open is that a value measured
33/// against one kind of QSM belongs with that kind of QSM. Changing one method's default without
34/// the others makes two variants of the same method incomparable on the same data — which is the
35/// state QSM.rs#130 was filed about — so these are shared constants rather than four literals.
36pub const DR_SHIN_2021: f64 = 137.0;
37
38/// χ-separation's relaxometric constant D_r re-measured against COSMOS, in Hz/ppm.
39///
40/// Kim et al. 2025 (Human Brain Mapping 46(4):e70136) repeated Shin's regression using
41/// COSMOS-reconstructed χ maps instead of single-orientation QSM and obtained 114 Hz/ppm,
42/// attributing the difference to the QSM algorithm behind the training data. The SNU-LIST
43/// toolbox's own script carries the same note against this value.
44///
45/// It is the right constant for a COSMOS-referenced χ — which is what the χ-sepnet family of
46/// networks was trained on, and why
47/// [`ChiSepNetNorm`](chisepnet::ChiSepNetNorm) and
48/// [`R2PrimeNetNorm`](crate::relaxometry::R2PrimeNetNorm) bake 114 in as a fixed training
49/// normalisation rather than a tunable. Those are deliberately *not* defined in terms of this
50/// constant: a network's normalisation cannot move when the field's best estimate of the physical
51/// relaxivity does, and tying them together would invite exactly that edit.
52pub const DR_KIM_2025_COSMOS: f64 = 114.0;
53
54/// The three maps every separation method returns, in the crate's standard order:
55/// `(chi_pos, chi_neg, chi_total)` in ppm, where `chi_pos >= 0`, `chi_neg <= 0`, and
56/// `chi_total` is their sum. Named because the bare triple is otherwise indistinguishable
57/// from any other three volumes at a call site.
58pub type ChiSepMaps = (Vec<f64>, Vec<f64>, Vec<f64>);
59
60/// One voxel of a network's de-normalised source magnitudes, as `(χ+, χ−, χ_total)` in the
61/// crate's order: χ+ and |χ−| are set to 0 where they are negative, χ− is returned signed (≤ 0),
62/// and χ_total is the sum of the two after that.
63///
64/// A source magnitude below zero has no physical meaning; it comes from the de-normalisation
65/// (`z·std + mean`) of a network output with no output activation. χ-sepnet's own `test.py` sets
66/// those voxels to 0 ("zero-truncation"); for SUSEP-Net, whose `test_demo.py` does not, it is a
67/// deliberate deviation.
68#[cfg_attr(not(feature = "onnx"), allow(dead_code))]
69pub(crate) fn sources_from_magnitudes(pos: f64, neg_mag: f64) -> (f64, f64, f64) {
70 let pos = pos.max(0.0);
71 let neg = -neg_mag.max(0.0);
72 (pos, neg, pos + neg)
73}
74
75pub mod chi_sep_ilsqr;
76pub mod chi_sep_medi;
77pub mod decompose;
78pub mod hc_chisep;
79pub mod r2star_qsm;
80pub mod wavesep;
81// The modules are always built, so configs can carry the networks' parameters; the inference
82// functions themselves need `onnx`.
83pub mod susep_net;
84pub mod chisepnet;
85
86pub use chi_sep_ilsqr::{chi_sep_ilsqr, chisep_toolbox_ilsqr, ChiSepIlsqrInputs, ChiSepIlsqrParams, NNLS_ITER_MAX, NNLS_TOL};
87pub use chi_sep_medi::{chi_sep_medi, ChiSepMediInputs, ChiSepMediParams, CHI_SEP_MEDI_MIN_ITER, toolbox_default_dr, TOOLBOX_DR_REFERENCE_TESLA, TOOLBOX_GYRO_HZ_PER_T};
88pub use decompose::{decompose, DecomposeChiBound, DecomposeParams, DecomposeStageSolver, DecomposeStageSolvers};
89pub use hc_chisep::{hc_chisep, HcChisepParams};
90pub use r2star_qsm::{r2star_qsm, r2star_qsm_from_magnitude, R2starQsmParams};
91pub use wavesep::{wavesep, WaveSepParams};
92pub use susep_net::{SusepNetNorm, SusepNetParams};
93pub use chisepnet::{ChiSepNetNorm, ChiSepNetParams, ChiSepNetR2, ChiSepNetRoute, ChiSepNetToolboxModels};
94#[cfg(feature = "onnx")]
95pub use susep_net::susep_net;
96#[cfg(feature = "onnx")]
97pub use chisepnet::{chisepnet, chisepnet_toolbox};
98
99#[cfg(test)]
100mod tests {
101 use super::*;
102
103 /// The analytic methods must default to the same relaxivity.
104 ///
105 /// `chi_sep_medi` and `chi_sep_ilsqr` follow the χ-separation toolbox, which scales
106 /// 137 Hz/ppm with the field (`137·(B0/2.9)`); [`toolbox_default_dr_is_shin_2021_at_2_9_tesla`]
107 /// checks it.
108 ///
109 /// D_r is a physical constant of the forward model, not a tuning knob: when two variants of
110 /// χ-separation default to different values, the same data analysed with each is not
111 /// comparable, and nothing in the output says why. `chi_sep_medi` sat at 114/30 against the
112 /// others' 137/137 until QSM.rs#130 — undetected because no test asserted they agreed.
113 #[test]
114 fn analytic_methods_share_one_relaxivity() {
115 let wave = WaveSepParams::default();
116 let hc = HcChisepParams::default();
117
118 for (name, dr) in [
119 ("wavesep.dr_pos", wave.dr_pos),
120 ("wavesep.dr_neg", wave.dr_neg),
121 ("hc_chisep.dr_pos_3t", hc.dr_pos_3t),
122 ] {
123 assert_eq!(
124 dr, DR_SHIN_2021,
125 "{name} defaults to {dr} Hz/ppm, not the shared {DR_SHIN_2021}. If that is \
126 deliberate, say which paper it comes from here; if it is not, it makes this \
127 method's output incomparable with the others' on the same data."
128 );
129 }
130 }
131
132 /// The toolbox's default `D_r` is Shin et al. 2021's 137 Hz/ppm at 2.9 T and scales with
133 /// the field.
134 #[test]
135 fn toolbox_default_dr_is_shin_2021_at_2_9_tesla() {
136 let at = |b0: f64| toolbox_default_dr(TOOLBOX_GYRO_HZ_PER_T * b0);
137 assert!((at(2.9) - DR_SHIN_2021).abs() < 1e-12);
138 assert_eq!(at(3.0), 141.7241379310345, "the toolbox's value at CF = 127.74 MHz");
139 let (dp, dn) = ChiSepMediParams { cf: TOOLBOX_GYRO_HZ_PER_T * 7.0, ..Default::default() }.dr();
140 assert!((dp - 137.0 * 7.0 / 2.9).abs() < 1e-9 && dp == dn);
141 let p = ChiSepIlsqrParams { cf: TOOLBOX_GYRO_HZ_PER_T * 7.0, ..Default::default() };
142 assert_eq!(p.dr(), (dp, dn), "chi_sep_ilsqr shares chi_sep_medi's default");
143 }
144
145 /// Negative source magnitudes become 0; non-negative ones pass through unchanged, and
146 /// χ_total stays χ+ + χ−.
147 #[test]
148 fn network_sources_are_never_the_wrong_sign() {
149 assert_eq!(sources_from_magnitudes(0.03, 0.02), (0.03, -0.02, 0.03 - 0.02));
150 assert_eq!(sources_from_magnitudes(-0.01, 0.02), (0.0, -0.02, -0.02));
151 assert_eq!(sources_from_magnitudes(0.03, -0.01), (0.03, 0.0, 0.03));
152 let (p, n, t) = sources_from_magnitudes(-0.01, -0.01);
153 assert_eq!((p, t), (0.0, 0.0));
154 assert!(n == 0.0);
155 }
156
157 /// The two published values are distinct and neither drifted.
158 ///
159 /// They are not interchangeable: 137 was regressed against single-orientation QSM, 114
160 /// against COSMOS. A future edit that "unifies" them would silently re-reference every
161 /// analytic method's χ+/χ− to a reconstruction its input is not.
162 #[test]
163 fn published_relaxivities_are_distinct() {
164 assert_eq!(DR_SHIN_2021, 137.0, "Shin et al. 2021, NeuroImage 240:118371");
165 assert_eq!(DR_KIM_2025_COSMOS, 114.0, "Kim et al. 2025, Hum Brain Mapp 46(4):e70136");
166 assert_ne!(DR_SHIN_2021, DR_KIM_2025_COSMOS);
167 }
168
169 /// The χ-sepnet family's 114 is a baked training normalisation, not the physical relaxivity,
170 /// so it must keep its own literal rather than following [`DR_SHIN_2021`] when that moves.
171 #[cfg(feature = "onnx")]
172 #[test]
173 fn network_normalisation_is_not_the_analytic_default() {
174 assert_eq!(ChiSepNetNorm::default().dr, DR_KIM_2025_COSMOS);
175 assert_eq!(crate::relaxometry::R2PrimeNetNorm::default().dr, DR_KIM_2025_COSMOS);
176 }
177}