qsm_core/pipeline/
per_echo.rs

1//! Per-echo QSM: a susceptibility map reconstructed from each echo on its own.
2//!
3//! The usual pipeline combines the echoes into one field ([`super::run_field_mapping`]) and
4//! inverts that. A per-echo reconstruction ([`QsmPipelineConfig::qsm_per_echo`](super::config::QsmPipelineConfig::qsm_per_echo)) instead runs
5//! every stage once per echo:
6//!
7//! 1. **Field**, [`run_field_mapping_per_echo`]: echo `e`'s phase unwrapped and divided by its own
8//!    echo time. The echo combination, phase-offset removal and bipolar correction combine or
9//!    compare echoes, so they do not run.
10//! 2. **Background removal, inversion and referencing**: the host calls
11//!    [`super::run_bg_removal`], [`super::run_dipole_inversion`] and [`super::apply_reference`]
12//!    on that echo's field with [`ScanMetadata::echo`](super::config::ScanMetadata::echo)`(e)`, so a method that converts at an
13//!    echo time (STAR-QSM's `te`, MEDI's, …) uses that echo's.
14//! 3. **Outputs**: the per-echo χ maps, and their voxel-wise mean across echoes
15//!    ([`mean_across_echoes`]) as the single QSM. [`voxel_major`] lays the per-echo maps out as
16//!    [`SeparationInputs::qsm_per_echo`](super::SeparationInputs::qsm_per_echo) takes them.
17//!
18//! # Why, and the reference
19//!
20//! DECOMPOSE's signal model synthesises each echo's phase from that echo's χ(t): "Note that QSM
21//! of each echo, χ(t), is used instead of the average QSM across all echoes" (Chen et al. 2021,
22//! NeuroImage 242:118477, §2.2), and "it is necessary to use TE-dependent QSM as the input"
23//! (§4.1). Their pipeline (§3.6, and `DECOMPOSE_analysis.m`) is the one above: Laplacian
24//! unwrapping of each echo's raw phase, V-SHARP, and STAR-QSM "to compute the susceptibility
25//! maps for each echo". The pipeline rule `decompose.qsm-per-echo` recommends this mode when
26//! DECOMPOSE is chosen.
27//!
28//! # The mean across echoes
29//!
30//! [`mean_across_echoes`] is the plain voxel-wise arithmetic mean with equal weights: the
31//! "average QSM across all echoes" Chen et al. contrast χ(t) with (§2.2), which names no
32//! weighting. Referencing commutes with it: the mean of maps each referenced to their mean over
33//! one mask is referenced to its mean over that mask.
34
35pub use super::field_mapping::run_field_mapping_per_echo;
36
37/// The voxel-wise mean of `maps` (one map per echo, each `n_voxels` long), with equal weights.
38/// See the [module docs](self#the-mean-across-echoes). Empty when `maps` is; panics if the maps'
39/// lengths differ.
40pub fn mean_across_echoes<M: AsRef<[f64]>>(maps: &[M]) -> Vec<f64> {
41    let Some(first) = maps.first() else { return Vec::new() };
42    let n = first.as_ref().len();
43    let mut out = vec![0.0; n];
44    for m in maps {
45        let m = m.as_ref();
46        assert_eq!(m.len(), n, "per-echo maps differ in length");
47        for (o, &v) in out.iter_mut().zip(m) {
48            *o += v;
49        }
50    }
51    let k = maps.len() as f64;
52    out.iter_mut().for_each(|o| *o /= k);
53    out
54}
55
56/// `maps` (one per echo, each `n_voxels` long) interleaved voxel-major, `(n_voxels, n_echoes)`
57/// with the echo fastest, as [`SeparationInputs::qsm_per_echo`](super::SeparationInputs::qsm_per_echo)
58/// and the multi-echo magnitude take them. Panics if the maps' lengths differ.
59pub fn voxel_major<M: AsRef<[f64]>>(maps: &[M]) -> Vec<f64> {
60    let ne = maps.len();
61    let n = maps.first().map_or(0, |m| m.as_ref().len());
62    let mut out = vec![0.0; n * ne];
63    for (e, m) in maps.iter().enumerate() {
64        let m = m.as_ref();
65        assert_eq!(m.len(), n, "per-echo maps differ in length");
66        for (v, &x) in m.iter().enumerate() {
67            out[v * ne + e] = x;
68        }
69    }
70    out
71}
72
73#[cfg(test)]
74mod tests {
75    use super::*;
76    use crate::pipeline::config::ScanMetadata;
77    use crate::pipeline::config::{FieldMappingConfig, UnwrappingAlgorithm};
78    use crate::pipeline::run_field_mapping;
79
80    #[test]
81    fn mean_and_layout() {
82        let maps = [vec![1.0, 2.0, 3.0], vec![3.0, 4.0, 5.0]];
83        assert_eq!(mean_across_echoes(&maps), vec![2.0, 3.0, 4.0]);
84        assert_eq!(voxel_major(&maps), vec![1.0, 3.0, 2.0, 4.0, 3.0, 5.0]);
85        assert!(mean_across_echoes::<Vec<f64>>(&[]).is_empty());
86    }
87
88    #[test]
89    fn echo_metadata_keeps_one_echo() {
90        let m = ScanMetadata {
91            dims: (2, 2, 2), voxel_size: (1.0, 1.0, 1.0), echo_times: vec![0.004, 0.008, 0.012],
92            field_strength: 3.0, b0_direction: (0.0, 0.0, 1.0), slice_geometry: None, affine: None,
93        };
94        let e = m.echo(1);
95        assert_eq!(e.echo_times, vec![0.008]);
96        assert_eq!((e.dims, e.field_strength), (m.dims, m.field_strength));
97    }
98
99    /// Each echo's field is that echo run alone through [`run_field_mapping`], whatever the echo
100    /// combination is set to (here a robust fit, which one echo could not run).
101    #[test]
102    fn each_echo_is_a_single_echo_field() {
103        let (nx, ny, nz) = (8, 8, 8);
104        let n = nx * ny * nz;
105        let tes = [0.005, 0.010, 0.015];
106        let metadata = ScanMetadata {
107            dims: (nx, ny, nz), voxel_size: (1.0, 1.0, 1.0), echo_times: tes.to_vec(),
108            field_strength: 3.0, b0_direction: (0.0, 0.0, 1.0), slice_geometry: None, affine: None,
109        };
110        let freq: Vec<f64> = (0..n).map(|i| 20.0 * ((i % nx) as f64 - 3.5)).collect(); // rad/s
111        let phases: Vec<Vec<f64>> = tes.iter().map(|&t| freq.iter().map(|f| f * t).collect()).collect();
112        let refs: Vec<&[f64]> = phases.iter().map(|p| p.as_slice()).collect();
113        let mask = vec![1u8; n];
114        for alg in [UnwrappingAlgorithm::Romeo, UnwrappingAlgorithm::Laplacian, UnwrappingAlgorithm::Prelude] {
115            let single = FieldMappingConfig { unwrapping_algorithm: alg, ..Default::default() };
116            let config = FieldMappingConfig {
117                b0_estimation: crate::pipeline::config::B0EstimationMethod::RobustFit(Default::default()),
118                phase_offset_removal: true,
119                ..single.clone()
120            };
121            let per_echo = run_field_mapping_per_echo(&refs, None, &mask, &metadata, &config, &mut |_, _| {}).unwrap();
122            assert_eq!(per_echo.len(), tes.len());
123            for (e, r) in per_echo.iter().enumerate() {
124                let alone = run_field_mapping(&refs[e..=e], None, &mask, &metadata.echo(e), &single, &mut |_, _| {}).unwrap();
125                assert_eq!(r.b0_field_ppm, alone.b0_field_ppm, "{alg:?} echo {e}");
126            }
127        }
128    }
129}
qsm_core/pipeline/per_echo.rs

qsm_core/pipeline/
per_echo.rs
