qsm_core/pipeline/
separation.rs

1//! Susceptibility source-separation stage.
2//!
3//! Config-driven dispatcher over the χ-separation methods (χ+ / χ−). Unlike the
4//! main phase→field→BFR→inversion pipeline, separation is a downstream stage that
5//! needs extra inputs (a conventional QSM plus a relaxometry map / multi-echo
6//! magnitude), so it is invoked explicitly with a [`SeparationInputs`] bundle
7//! rather than folded into the top-level runner.
8//!
9//! Units follow the library convention: local field & QSM in ppm, R2'/R2* in Hz,
10//! echo times in seconds. Output is `(chi_pos ≥ 0, chi_neg ≤ 0, chi_total)` in ppm.
11//!
12//! # Reference pipelines
13//!
14//! Each method's module docs have a "Reference pipeline" section: the masking, unwrapping,
15//! background removal and QSM method its authors fed it, and where this stage's inputs differ
16//! from theirs. The main differences, all recorded there with sources:
17//! - `Decompose` takes a QSM per echo ([`SeparationInputs::qsm_per_echo`]), as its authors
18//!   do (a STAR-QSM per echo); producing those is the caller's.
19//! - The analytic methods' D_r is calibrated at 3 T only. `ChiSepMedi` and `ChiSepIlsqr` default
20//!   to the χ-separation toolbox's `137·(B0/2.9)` (from the authors' code). `WaveSep` keeps
21//!   137 Hz/ppm at every field, while `R2starQsm` and `HcChisep` scale theirs by B0/3;
22//!   those choices do not come from the authors.
23//! - `ChiSepNet` runs the authors' route by default ([`ChiSepNetRoute::Toolbox`](crate::separation::ChiSepNetRoute)):
24//!   its own QSMnetME QSM (so [`SeparationInputs::qsm`] is not read), R2′ from a measured map or
25//!   R2PRIMEnet on R2*, in the 3 T, 1 mm training frame.
26
27use super::config::*;
28use crate::separation::{
29    ChiSepIlsqrInputs, ChiSepIlsqrParams, ChiSepMediInputs, ChiSepMediParams, DecomposeParams, HcChisepParams, R2starQsmParams,
30};
31
32/// Inputs for [`run_separation`]. Provide whatever the chosen algorithm needs;
33/// the dispatcher errors if a required input is missing.
34#[derive(Clone, Copy)]
35pub struct SeparationInputs<'a> {
36    /// Local (tissue) field in ppm — for `ChiSepIlsqr` / `ChiSepMedi`.
37    pub local_field_ppm: &'a [f64],
38    /// Conventional QSM χ_total in ppm — the QSM-init for `ChiSepIlsqr` / `ChiSepMedi` and the
39    /// input for `R2starQsm` / `WaveSep` / `Decompose` / `HcChisep`.
40    pub qsm: &'a [f64],
41    /// QSM of each echo in ppm, voxel-major `(n_voxels, n_echoes)` (`Decompose`, required).
42    pub qsm_per_echo: Option<&'a [f64]>,
43    /// Binary brain mask (`n_voxels`, 1 = inside).
44    pub mask: &'a [u8],
45    /// R2' map in Hz (`ChiSepIlsqr` / `ChiSepMedi` / `WaveSep` / `HcChisep`).
46    pub r2prime: Option<&'a [f64]>,
47    /// R2* map in Hz (`R2starQsm`); if absent, fit from `magnitude_multi`.
48    pub r2star: Option<&'a [f64]>,
49    /// Root-sum-of-squares magnitude (`ChiSepMedi`; the toolbox's `chi_sep_iLSQR` takes it but
50    /// never uses it, so `ChiSepIlsqr` does not).
51    pub magnitude_rss: Option<&'a [f64]>,
52    /// Noise standard deviation of the field, `N_std` (`ChiSepMedi`; see
53    /// [`ChiSepMediInputs::n_std`]).
54    pub n_std: Option<&'a [f64]>,
55    /// CSF mask, nonzero inside (`ChiSepMedi`; see [`ChiSepMediInputs::mask_csf`]).
56    pub mask_csf: Option<&'a [u8]>,
57    /// Multi-echo magnitude, voxel-major `(n_voxels, n_echoes)`
58    /// (`R2starQsm` / `Decompose` / `HcChisep`).
59    pub magnitude_multi: Option<&'a [f64]>,
60    /// Multi-echo spin-echo magnitude, voxel-major `(n_voxels, n_se)` (`HcChisep`, optional).
61    pub se_magnitude_multi: Option<&'a [f64]>,
62}
63
64/// Result of [`run_separation`]: paramagnetic / diamagnetic / total maps in ppm.
65#[derive(Clone, Debug)]
66pub struct SeparationResult {
67    /// χ+ ≥ 0 (paramagnetic).
68    pub chi_pos: Vec<f64>,
69    /// χ− ≤ 0 (diamagnetic, signed).
70    pub chi_neg: Vec<f64>,
71    /// χ_total = χ+ + χ−.
72    pub chi_total: Vec<f64>,
73}
74
75/// Run the configured χ-separation algorithm.
76///
77/// # Arguments
78/// * `inputs` — See [`SeparationInputs`]; provide what the algorithm requires.
79/// * `metadata` — Scan metadata (grid, B0 direction, field strength, echo times).
80/// * `config` — See [`SeparationConfig`].
81/// * `progress` — Progress callback `(current_iter, max_iter)`.
82pub fn run_separation(
83    inputs: SeparationInputs,
84    metadata: &ScanMetadata,
85    config: &SeparationConfig,
86    progress: &mut dyn FnMut(usize, usize),
87) -> Result<SeparationResult, PipelineError> {
88    metadata.require_contiguous_slices()?;
89    let grid = metadata.grid();
90    let n = grid.n_total();
91    let bdir = metadata.b0_direction;
92    let cf = crate::constants::GAMMA_BAR_HZ_PER_T * metadata.field_strength;
93    let b0 = metadata.field_strength;
94
95    if inputs.qsm.len() != n || inputs.mask.len() != n {
96        return Err(PipelineError::DimensionMismatch { expected: n, got: inputs.qsm.len() });
97    }
98    check_inputs(&inputs, metadata, config.algorithm, n)?;
99    let mask = inputs.mask;
100
101    let (chi_pos, chi_neg, chi_total) = match config.algorithm {
102        SeparationAlgorithm::ChiSepIlsqr => {
103            let ilsqr_inputs = ChiSepIlsqrInputs {
104                local_field: inputs.local_field_ppm,
105                r2prime: need(inputs.r2prime, "r2prime")?,
106                qsm: inputs.qsm,
107                mask,
108            };
109            let params = ChiSepIlsqrParams { cf, ..config.chi_sep_ilsqr.clone() };
110            crate::separation::chi_sep_ilsqr(&ilsqr_inputs, &grid, bdir, &params, progress)?
111        }
112        SeparationAlgorithm::ChiSepMedi => {
113            let medi_inputs = ChiSepMediInputs {
114                magnitude: need(inputs.magnitude_rss, "magnitude_rss")?,
115                local_field: inputs.local_field_ppm,
116                r2prime: need(inputs.r2prime, "r2prime")?,
117                n_std: need(inputs.n_std, "n_std")?,
118                qsm: inputs.qsm,
119                mask,
120                mask_csf: need(inputs.mask_csf, "mask_csf")?,
121            };
122            let params = ChiSepMediParams { cf, ..config.chi_sep_medi.clone() };
123            crate::separation::chi_sep_medi(&medi_inputs, &grid, bdir, &params, progress)?
124        }
125        SeparationAlgorithm::R2starQsm => {
126            let params = R2starQsmParams { b0, ..config.r2star_qsm.clone() };
127            if let Some(r2star) = inputs.r2star {
128                crate::separation::r2star_qsm(inputs.qsm, r2star, mask, &params)
129            } else {
130                let magnitude = need(inputs.magnitude_multi, "r2star or magnitude_multi")?;
131                crate::separation::r2star_qsm_from_magnitude(
132                    magnitude, &metadata.echo_times, inputs.qsm, mask, &params,
133                )
134            }
135        }
136        SeparationAlgorithm::WaveSep => {
137            let r2prime = need(inputs.r2prime, "r2prime")?;
138            crate::separation::wavesep(
139                inputs.qsm, r2prime, mask, &grid, &config.wavesep, progress,
140            )?
141        }
142        SeparationAlgorithm::Decompose => {
143            let magnitude = need(inputs.magnitude_multi, "magnitude_multi")?;
144            let qsm_per_echo = need(inputs.qsm_per_echo, "qsm_per_echo")?;
145            let params = DecomposeParams { b0, ..config.decompose.clone() };
146            crate::separation::decompose(
147                qsm_per_echo, magnitude, &metadata.echo_times, mask, &params, progress,
148            )?
149        }
150        SeparationAlgorithm::HcChisep => {
151            let r2prime = need(inputs.r2prime, "r2prime")?;
152            let magnitude = need(inputs.magnitude_multi, "magnitude_multi")?;
153            let params = HcChisepParams { b0, ..config.hc_chisep.clone() };
154            crate::separation::hc_chisep(
155                inputs.qsm, r2prime, magnitude, &metadata.echo_times,
156                inputs.se_magnitude_multi, mask, &grid, &params, progress,
157            )
158        }
159        SeparationAlgorithm::SusepNet => {
160            let r2prime = need(inputs.r2prime, "r2prime")?;
161            config.susep_net.validate().map_err(PipelineError::InvalidConfig)?;
162            run_susep_net(
163                inputs.local_field_ppm, inputs.qsm, r2prime, mask, &grid, &config.susep_net,
164                progress,
165            )?
166        }
167        SeparationAlgorithm::ChiSepNet => {
168            config.chi_sepnet.validate().map_err(PipelineError::InvalidConfig)?;
169            match config.chi_sepnet.route {
170                crate::separation::ChiSepNetRoute::Native => {
171                    let r2prime = need(inputs.r2prime, "r2prime")?;
172                    run_chi_sepnet(
173                        inputs.local_field_ppm, inputs.qsm, r2prime, mask, &grid, &config.chi_sepnet,
174                    )?
175                }
176                // The authors' route makes its own QSM (QSMnetME) and, from R2*, its own R2′
177                // (R2PRIMEnet), as the toolbox's script does: a measured R2′ when there is one,
178                // else R2* (`Chisep_script.m`:177-185).
179                crate::separation::ChiSepNetRoute::Toolbox => run_chi_sepnet_toolbox(
180                    inputs.local_field_ppm, inputs.r2prime, inputs.r2star, mask, &grid, bdir, b0,
181                    &config.chi_sepnet,
182                )?,
183            }
184        }
185    };
186
187    Ok(SeparationResult { chi_pos, chi_neg, chi_total })
188}
189
190/// Source the SUSEP-Net weights and run inference (requires the `onnx` feature).
191#[cfg(feature = "onnx")]
192fn run_susep_net(
193    local_field_ppm: &[f64],
194    qsm: &[f64],
195    r2prime: &[f64],
196    mask: &[u8],
197    grid: &crate::Grid,
198    params: &crate::separation::SusepNetParams,
199    progress: &mut dyn FnMut(usize, usize),
200) -> Result<crate::separation::ChiSepMaps, PipelineError> {
201    let spec = crate::models::find_model("susep-net")
202        .ok_or_else(|| PipelineError::InvalidConfig("susep-net not in model registry".into()))?;
203    let bytes = crate::models::primary_weight_bytes(spec).map_err(PipelineError::InvalidConfig)?;
204    crate::separation::susep_net(
205        local_field_ppm, qsm, r2prime, mask, grid, &bytes,
206        &crate::separation::SusepNetNorm::default(),
207        params,
208        progress,
209    )
210    .map_err(|e| PipelineError::AlgorithmError(e.to_string()))
211}
212
213#[cfg(not(feature = "onnx"))]
214fn run_susep_net(
215    _local_field_ppm: &[f64],
216    _qsm: &[f64],
217    _r2prime: &[f64],
218    _mask: &[u8],
219    _grid: &crate::Grid,
220    _params: &crate::separation::SusepNetParams,
221    _progress: &mut dyn FnMut(usize, usize),
222) -> Result<crate::separation::ChiSepMaps, PipelineError> {
223    Err(PipelineError::InvalidConfig(
224        "SUSEP-Net requires building qsm-core with the 'onnx' feature".into(),
225    ))
226}
227
228/// Source the χ-sepnet weights and run inference (requires the `onnx` feature).
229#[cfg(feature = "onnx")]
230fn run_chi_sepnet(
231    local_field_ppm: &[f64],
232    qsm: &[f64],
233    r2prime: &[f64],
234    mask: &[u8],
235    grid: &crate::Grid,
236    params: &crate::separation::ChiSepNetParams,
237) -> Result<crate::separation::ChiSepMaps, PipelineError> {
238    let bytes = crate::models::primary_weight("chi-sepnet").map_err(PipelineError::InvalidConfig)?;
239    crate::separation::chisepnet(
240        local_field_ppm, qsm, r2prime, mask, grid, &bytes,
241        &crate::separation::ChiSepNetNorm::default(),
242        params,
243    )
244    .map_err(|e| PipelineError::AlgorithmError(e.to_string()))
245}
246
247#[cfg(not(feature = "onnx"))]
248fn run_chi_sepnet(
249    _local_field_ppm: &[f64],
250    _qsm: &[f64],
251    _r2prime: &[f64],
252    _mask: &[u8],
253    _grid: &crate::Grid,
254    _params: &crate::separation::ChiSepNetParams,
255) -> Result<crate::separation::ChiSepMaps, PipelineError> {
256    Err(PipelineError::InvalidConfig(
257        "χ-sepnet requires building qsm-core with the 'onnx' feature".into(),
258    ))
259}
260
261/// Source the weights and run χ-sepnet by the authors' route ([`crate::separation::chisepnet_toolbox`]).
262#[cfg(feature = "onnx")]
263#[allow(clippy::too_many_arguments)]
264fn run_chi_sepnet_toolbox(
265    local_field_ppm: &[f64],
266    r2prime: Option<&[f64]>,
267    r2star: Option<&[f64]>,
268    mask: &[u8],
269    grid: &crate::Grid,
270    bdir: (f64, f64, f64),
271    b0_tesla: f64,
272    params: &crate::separation::ChiSepNetParams,
273) -> Result<crate::separation::ChiSepMaps, PipelineError> {
274    use crate::separation::{ChiSepNetR2, ChiSepNetToolboxModels};
275    let r2 = match (r2prime, r2star) {
276        (Some(v), _) => ChiSepNetR2::R2Prime(v),
277        (None, Some(v)) => ChiSepNetR2::R2Star(v),
278        (None, None) => {
279            return Err(PipelineError::InvalidInput(
280                "χ-sepnet needs a measured R2' (χ-sepnet-R2') or an R2* map (χ-sepnet-R2*)".into(),
281            ))
282        }
283    };
284    let weight = |id: &str| crate::models::primary_weight(id).map_err(PipelineError::InvalidConfig);
285    let chisepnet = weight("chi-sepnet")?;
286    let qsmnetme = weight("qsmnetme")?;
287    let r2primenet = match r2 {
288        ChiSepNetR2::R2Star(_) => Some(weight("r2primenet")?),
289        ChiSepNetR2::R2Prime(_) => None,
290    };
291    let models = ChiSepNetToolboxModels { chisepnet: &chisepnet, qsmnetme: &qsmnetme, r2primenet: r2primenet.as_deref() };
292    crate::separation::chisepnet_toolbox(
293        local_field_ppm, r2, mask, grid, bdir, b0_tesla, &models,
294        &crate::separation::ChiSepNetNorm::default(), params,
295    )
296    .map_err(|e| PipelineError::AlgorithmError(e.to_string()))
297}
298
299#[cfg(not(feature = "onnx"))]
300#[allow(clippy::too_many_arguments)]
301fn run_chi_sepnet_toolbox(
302    _local_field_ppm: &[f64],
303    _r2prime: Option<&[f64]>,
304    _r2star: Option<&[f64]>,
305    _mask: &[u8],
306    _grid: &crate::Grid,
307    _bdir: (f64, f64, f64),
308    _b0_tesla: f64,
309    _params: &crate::separation::ChiSepNetParams,
310) -> Result<crate::separation::ChiSepMaps, PipelineError> {
311    Err(PipelineError::InvalidConfig(
312        "χ-sepnet requires building qsm-core with the 'onnx' feature".into(),
313    ))
314}
315
316/// Check what the separation algorithms would otherwise `assert!` on, so bad input is an error
317/// rather than a panic: map sizes, the multi-echo magnitude's echo count, Decompose's minimum
318/// of three echoes, and an empty mask.
319fn check_inputs(
320    inputs: &SeparationInputs,
321    metadata: &ScanMetadata,
322    algorithm: SeparationAlgorithm,
323    n: usize,
324) -> Result<(), PipelineError> {
325    let maps = [
326        ("local_field_ppm", Some(inputs.local_field_ppm)),
327        ("r2prime", inputs.r2prime),
328        ("r2star", inputs.r2star),
329        ("magnitude_rss", inputs.magnitude_rss),
330    ];
331    let needs_local_field = matches!(
332        algorithm,
333        SeparationAlgorithm::ChiSepIlsqr | SeparationAlgorithm::ChiSepMedi
334            | SeparationAlgorithm::SusepNet | SeparationAlgorithm::ChiSepNet
335    );
336    for (name, map) in maps {
337        if name == "local_field_ppm" && !needs_local_field {
338            continue;
339        }
340        if let Some(m) = map {
341            crate::error::check_len(name, m.len(), n)?;
342        }
343    }
344    let n_echoes = metadata.echo_times.len();
345    if let Some(m) = inputs.magnitude_multi {
346        crate::error::check_len(
347            &format!("magnitude_multi ({n_echoes} echoes x {n} voxels)"), m.len(), n * n_echoes,
348        )?;
349    }
350    if let Some(m) = inputs.qsm_per_echo {
351        crate::error::check_len(
352            &format!("qsm_per_echo ({n_echoes} echoes x {n} voxels)"), m.len(), n * n_echoes,
353        )?;
354    }
355    if algorithm == SeparationAlgorithm::Decompose && n_echoes < 3 {
356        return Err(PipelineError::InvalidInput(format!(
357            "Decompose needs at least 3 echoes, got {n_echoes}"
358        )));
359    }
360    if !inputs.mask.iter().any(|&m| m != 0) {
361        return Err(PipelineError::InvalidInput("the mask is empty".into()));
362    }
363    Ok(())
364}
365
366/// Require an optional input, or return an `InvalidInput` error naming it.
367fn need<'a, T>(opt: Option<&'a [T]>, what: &str) -> Result<&'a [T], PipelineError> {
368    opt.ok_or_else(|| PipelineError::InvalidInput(format!("{what} required for this algorithm")))
369}
370
371#[cfg(test)]
372mod tests {
373    use super::*;
374
375    fn meta(nx: usize, ny: usize, nz: usize) -> ScanMetadata {
376        ScanMetadata {
377            dims: (nx, ny, nz),
378            voxel_size: (1.0, 1.0, 1.0),
379            echo_times: vec![0.004, 0.012, 0.020, 0.028],
380            field_strength: 3.0,
381            b0_direction: (0.0, 0.0, 1.0),
382            slice_geometry: None,
383            affine: None,
384        }
385    }
386
387    fn base_inputs<'a>(qsm: &'a [f64], mask: &'a [u8]) -> SeparationInputs<'a> {
388        SeparationInputs {
389            local_field_ppm: qsm,
390            qsm,
391            qsm_per_echo: None,
392            mask,
393            r2prime: None,
394            r2star: None,
395            magnitude_rss: None,
396            n_std: None,
397            mask_csf: None,
398            magnitude_multi: None,
399            se_magnitude_multi: None,
400        }
401    }
402
403    #[test]
404    fn r2star_qsm_dispatch_closed_form() {
405        let (nx, ny, nz) = (6, 6, 6);
406        let n = nx * ny * nz;
407        let qsm = vec![0.02; n];
408        let r2star = vec![10.0; n];
409        let mask = vec![1u8; n];
410        let m = meta(nx, ny, nz);
411        let mut inputs = base_inputs(&qsm, &mask);
412        inputs.r2star = Some(&r2star);
413        let cfg = SeparationConfig { algorithm: SeparationAlgorithm::R2starQsm, ..Default::default() };
414        let r = run_separation(inputs, &m, &cfg, &mut |_, _| {}).unwrap();
415        assert_eq!(r.chi_pos.len(), n);
416        for i in 0..n {
417            assert!(r.chi_pos[i] >= 0.0 && r.chi_neg[i] <= 0.0);
418            assert!((r.chi_total[i] - (r.chi_pos[i] + r.chi_neg[i])).abs() < 1e-9);
419        }
420    }
421
422    #[test]
423    fn wavesep_dispatch_runs() {
424        let (nx, ny, nz) = (8, 8, 8);
425        let n = nx * ny * nz;
426        let qsm = vec![0.01; n];
427        let r2prime = vec![2.0; n];
428        let mask = vec![1u8; n];
429        let m = meta(nx, ny, nz);
430        let mut inputs = base_inputs(&qsm, &mask);
431        inputs.r2prime = Some(&r2prime);
432        let cfg = SeparationConfig { algorithm: SeparationAlgorithm::WaveSep, ..Default::default() };
433        let r = run_separation(inputs, &m, &cfg, &mut |_, _| {}).unwrap();
434        assert_eq!(r.chi_pos.len(), n);
435    }
436
437    #[test]
438    fn chi_sep_ilsqr_and_medi_dispatch() {
439        let (nx, ny, nz) = (6, 6, 6);
440        let n = nx * ny * nz;
441        let field = vec![0.01; n];
442        let r2prime = vec![2.0; n];
443        let mag = vec![1.0; n];
444        let mask = vec![1u8; n];
445        let n_std = vec![0.01; n];
446        let no_csf = vec![0u8; n];
447        let m = meta(nx, ny, nz);
448        for alg in [SeparationAlgorithm::ChiSepIlsqr, SeparationAlgorithm::ChiSepMedi] {
449            let mut inputs = base_inputs(&field, &mask);
450            inputs.r2prime = Some(&r2prime);
451            inputs.magnitude_rss = Some(&mag);
452            inputs.n_std = Some(&n_std);
453            inputs.mask_csf = Some(&no_csf);
454            let cfg = SeparationConfig { algorithm: alg, ..Default::default() };
455            let r = run_separation(inputs, &m, &cfg, &mut |_, _| {}).unwrap();
456            assert_eq!(r.chi_pos.len(), n);
457        }
458        // chi-sep-ilsqr takes no magnitude, noise map or CSF mask (the toolbox's chi_sep_iLSQR
459        // never uses them)
460        let mut inputs = base_inputs(&field, &mask);
461        inputs.r2prime = Some(&r2prime);
462        let cfg = SeparationConfig { algorithm: SeparationAlgorithm::ChiSepIlsqr, ..Default::default() };
463        assert!(run_separation(inputs, &m, &cfg, &mut |_, _| {}).is_ok());
464        // chi-sep-medi takes the noise map and the CSF mask the toolbox takes, and says so
465        for missing in ["n_std", "mask_csf"] {
466            let mut inputs = base_inputs(&field, &mask);
467            inputs.r2prime = Some(&r2prime);
468            inputs.magnitude_rss = Some(&mag);
469            if missing != "n_std" {
470                inputs.n_std = Some(&n_std);
471            }
472            if missing != "mask_csf" {
473                inputs.mask_csf = Some(&no_csf);
474            }
475            let cfg = SeparationConfig { algorithm: SeparationAlgorithm::ChiSepMedi, ..Default::default() };
476            let e = run_separation(inputs, &m, &cfg, &mut |_, _| {}).unwrap_err();
477            assert!(e.to_string().contains(&format!("{missing} required")), "{e}");
478        }
479    }
480
481    #[test]
482    fn decompose_and_hc_dispatch() {
483        let (nx, ny, nz) = (6, 6, 6);
484        let n = nx * ny * nz;
485        let ne = 4;
486        let qsm = vec![0.01; n];
487        let r2prime = vec![2.0; n];
488        let mag_multi = vec![1.0; n * ne];
489        let mask = vec![1u8; n];
490        let m = meta(nx, ny, nz);
491
492        let qsm_per_echo = vec![0.01; n * ne];
493        let mut d_inputs = base_inputs(&qsm, &mask);
494        d_inputs.magnitude_multi = Some(&mag_multi);
495        let d_cfg = SeparationConfig { algorithm: SeparationAlgorithm::Decompose, ..Default::default() };
496        // The per-echo QSM is required.
497        assert!(run_separation(d_inputs, &m, &d_cfg, &mut |_, _| {}).is_err());
498        d_inputs.qsm_per_echo = Some(&qsm_per_echo);
499        let rd = run_separation(d_inputs, &m, &d_cfg, &mut |_, _| {}).unwrap();
500        assert_eq!(rd.chi_pos.len(), n);
501        // A per-echo QSM of the wrong length is an error, not a panic.
502        d_inputs.qsm_per_echo = Some(&qsm_per_echo[..n]);
503        assert!(run_separation(d_inputs, &m, &d_cfg, &mut |_, _| {}).is_err());
504
505        let mut h_inputs = base_inputs(&qsm, &mask);
506        h_inputs.r2prime = Some(&r2prime);
507        h_inputs.magnitude_multi = Some(&mag_multi);
508        let h_cfg = SeparationConfig { algorithm: SeparationAlgorithm::HcChisep, ..Default::default() };
509        let rh = run_separation(h_inputs, &m, &h_cfg, &mut |_, _| {}).unwrap();
510        assert_eq!(rh.chi_pos.len(), n);
511    }
512
513    #[test]
514    fn r2star_from_magnitude_path() {
515        let (nx, ny, nz) = (6, 6, 6);
516        let n = nx * ny * nz;
517        let ne = 4;
518        let qsm = vec![0.02; n];
519        let mag_multi = vec![1.0; n * ne];
520        let mask = vec![1u8; n];
521        let m = meta(nx, ny, nz);
522        // No r2star provided → fit from magnitude_multi.
523        let mut inputs = base_inputs(&qsm, &mask);
524        inputs.magnitude_multi = Some(&mag_multi);
525        let cfg = SeparationConfig { algorithm: SeparationAlgorithm::R2starQsm, ..Default::default() };
526        let r = run_separation(inputs, &m, &cfg, &mut |_, _| {}).unwrap();
527        assert_eq!(r.chi_pos.len(), n);
528    }
529
530    #[test]
531    fn missing_input_errors() {
532        let (nx, ny, nz) = (6, 6, 6);
533        let n = nx * ny * nz;
534        let qsm = vec![0.02; n];
535        let mask = vec![1u8; n];
536        let m = meta(nx, ny, nz);
537        // R2starQsm with neither r2star nor magnitude_multi → InvalidInput.
538        let inputs = base_inputs(&qsm, &mask);
539        let cfg = SeparationConfig { algorithm: SeparationAlgorithm::R2starQsm, ..Default::default() };
540        let err = run_separation(inputs, &m, &cfg, &mut |_, _| {}).unwrap_err();
541        assert!(matches!(err, PipelineError::InvalidInput(_)), "got {err:?}");
542    }
543
544    /// The deep-learning methods take their patch from the config, and a patch the network
545    /// cannot take is a config error, refused before any weights are fetched (so this needs
546    /// none). The defaults are the networks' own: the authors' χ-sepnet patch and SUSEP-Net's
547    /// whole volume.
548    #[test]
549    fn dl_patches_from_config() {
550        use crate::separation::{chisepnet, ChiSepNetParams, SusepNetParams};
551        let d = SeparationConfig::default();
552        assert_eq!(d.chi_sepnet, ChiSepNetParams { patch: chisepnet::AUTHORS_PATCH, ..Default::default() });
553        assert_eq!(d.chi_sepnet.route, crate::separation::ChiSepNetRoute::Toolbox);
554        assert_eq!(d.susep_net, SusepNetParams { patch: None });
555
556        let (nx, ny, nz) = (6, 6, 6);
557        let n = nx * ny * nz;
558        let v = vec![0.01; n];
559        let mask = vec![1u8; n];
560        let m = meta(nx, ny, nz);
561        let mut inputs = base_inputs(&v, &mask);
562        inputs.r2prime = Some(&v);
563        let bad = [
564            SeparationConfig {
565                algorithm: SeparationAlgorithm::ChiSepNet,
566                chi_sepnet: ChiSepNetParams { patch: (128, 120, 64), ..Default::default() },
567                ..Default::default()
568            },
569            SeparationConfig {
570                algorithm: SeparationAlgorithm::SusepNet,
571                susep_net: SusepNetParams { patch: Some((128, 128, 4)) },
572                ..Default::default()
573            },
574        ];
575        for cfg in bad {
576            let err = run_separation(inputs, &m, &cfg, &mut |_, _| {}).unwrap_err();
577            match err {
578                PipelineError::InvalidConfig(msg) => assert!(msg.contains("patch"), "{msg}"),
579                e => panic!("{:?}: expected InvalidConfig, got {e:?}", cfg.algorithm),
580            }
581        }
582    }
583
584    /// Inputs the separation algorithms would `assert!` on are refused by the stage instead
585    /// (QSMxT #273): too few echoes for Decompose, a multi-echo magnitude with the wrong echo
586    /// count, a map of the wrong size, an empty mask.
587    #[test]
588    fn bad_inputs_are_errors_not_panics() {
589        let (nx, ny, nz) = (6, 6, 6);
590        let n = nx * ny * nz;
591        let qsm = vec![0.01; n];
592        let mask = vec![1u8; n];
593        let mut m = meta(nx, ny, nz);
594        let run = |inputs: SeparationInputs, m: &ScanMetadata, alg| {
595            let cfg = SeparationConfig { algorithm: alg, ..Default::default() };
596            run_separation(inputs, m, &cfg, &mut |_, _| {}).map(|_| ())
597        };
598        let invalid = |r: Result<(), PipelineError>| matches!(r, Err(PipelineError::InvalidInput(_)));
599
600        // Decompose with two echoes
601        m.echo_times = vec![0.004, 0.012];
602        let mag2 = vec![1.0; n * 2];
603        let mut inputs = base_inputs(&qsm, &mask);
604        inputs.magnitude_multi = Some(&mag2);
605        assert!(invalid(run(inputs, &m, SeparationAlgorithm::Decompose)));
606
607        // multi-echo magnitude with 3 echoes for 4 echo times
608        let m4 = meta(nx, ny, nz);
609        let mag3 = vec![1.0; n * 3];
610        let mut inputs = base_inputs(&qsm, &mask);
611        inputs.magnitude_multi = Some(&mag3);
612        assert!(invalid(run(inputs, &m4, SeparationAlgorithm::R2starQsm)));
613
614        // R2' on another grid
615        let short = vec![1.0; n - 1];
616        let mut inputs = base_inputs(&qsm, &mask);
617        inputs.r2prime = Some(&short);
618        assert!(invalid(run(inputs, &m4, SeparationAlgorithm::WaveSep)));
619
620        // empty mask
621        let empty = vec![0u8; n];
622        let r2prime = vec![2.0; n];
623        let mut inputs = base_inputs(&qsm, &empty);
624        inputs.r2prime = Some(&r2prime);
625        inputs.magnitude_rss = Some(&qsm);
626        assert!(invalid(run(inputs, &m4, SeparationAlgorithm::ChiSepIlsqr)));
627    }
628}
qsm_core/pipeline/separation.rs

qsm_core/pipeline/
separation.rs
