qsm_core/pipeline/
referencing.rs

1//! QSM referencing
2//!
3//! Adjusts the susceptibility map by subtracting the mean value
4//! within the mask (mean referencing) or zeroing outside the mask.
5
6use super::config::QsmReference;
7
8/// Apply QSM referencing to a susceptibility map.
9///
10/// - `Mean`: subtract the mean of masked voxels, zero outside mask
11/// - `None`: zero outside mask only
12pub fn apply_reference(chi: &[f64], mask: &[u8], method: QsmReference) -> Vec<f64> {
13    let n = chi.len();
14    let mut result = vec![0.0; n];
15
16    match method {
17        QsmReference::Mean => {
18            let mut sum = 0.0;
19            let mut count = 0usize;
20            for i in 0..n {
21                if mask[i] > 0 {
22                    sum += chi[i];
23                    count += 1;
24                }
25            }
26            let mean = if count > 0 { sum / count as f64 } else { 0.0 };
27            for i in 0..n {
28                if mask[i] > 0 {
29                    result[i] = chi[i] - mean;
30                }
31            }
32        }
33        QsmReference::None => {
34            for i in 0..n {
35                if mask[i] > 0 {
36                    result[i] = chi[i];
37                }
38            }
39        }
40    }
41
42    result
43}
44
45#[cfg(test)]
46mod tests {
47    use super::*;
48
49    #[test]
50    fn test_reference_mean() {
51        let chi = vec![1.0, 2.0, 3.0, 4.0];
52        let mask = vec![1, 1, 1, 0];
53        let result = apply_reference(&chi, &mask, QsmReference::Mean);
54        // Mean of masked = (1+2+3)/3 = 2.0
55        assert!((result[0] - (-1.0)).abs() < 1e-10);
56        assert!((result[1] - 0.0).abs() < 1e-10);
57        assert!((result[2] - 1.0).abs() < 1e-10);
58        assert_eq!(result[3], 0.0); // masked out
59    }
60
61    #[test]
62    fn test_reference_none() {
63        let chi = vec![1.0, 2.0, 3.0];
64        let mask = vec![1, 0, 1];
65        let result = apply_reference(&chi, &mask, QsmReference::None);
66        assert_eq!(result[0], 1.0);
67        assert_eq!(result[1], 0.0);
68        assert_eq!(result[2], 3.0);
69    }
70
71    #[test]
72    fn test_reference_mean_empty_mask() {
73        let chi = vec![1.0, 2.0];
74        let mask = vec![0, 0];
75        let result = apply_reference(&chi, &mask, QsmReference::Mean);
76        assert_eq!(result, vec![0.0, 0.0]);
77    }
78}
qsm_core/pipeline/referencing.rs

qsm_core/pipeline/
referencing.rs
