qsm_core/pipeline/
phase_utils.rs

1//! Phase and field conversion utilities
2//!
3//! Canonical implementations shared by all pipeline consumers.
4
5use std::f64::consts::PI;
6
7/// Scale phase data to [-pi, pi] range in-place.
8///
9/// If data is already approximately in [-pi, pi] (within 10% tolerance),
10/// leaves it unchanged. Otherwise linearly maps from [min, max] to [-pi, pi].
11/// NaN/Inf values are replaced with 0.
12pub fn scale_phase_to_pi(data: &mut [f64]) {
13    if data.is_empty() {
14        return;
15    }
16
17    let mut min_val = f64::INFINITY;
18    let mut max_val = f64::NEG_INFINITY;
19    for &v in data.iter() {
20        if v.is_finite() {
21            if v < min_val {
22                min_val = v;
23            }
24            if v > max_val {
25                max_val = v;
26            }
27        }
28    }
29
30    // Replace non-finite values with 0
31    for v in data.iter_mut() {
32        if !v.is_finite() {
33            *v = 0.0;
34        }
35    }
36
37    let range = max_val - min_val;
38    if range < 1e-10 {
39        return;
40    }
41
42    // Check if already approximately in [-pi, pi]
43    let tol = 0.1 * PI;
44    if (min_val + PI).abs() < tol && (max_val - PI).abs() < tol {
45        return;
46    }
47
48    // Linearly rescale to [-pi, pi]
49    let scale = 2.0 * PI / range;
50    for v in data.iter_mut() {
51        *v = (*v - min_val) * scale - PI;
52    }
53}
54
55/// Convert field from Hz to ppm.
56pub fn hz_to_ppm(field_hz: &[f64], field_strength: f64) -> Vec<f64> {
57    let gamma_hz = 42.576e6; // Hz/T (proton gyromagnetic ratio)
58    let scale = 1e6 / (gamma_hz * field_strength);
59    field_hz.iter().map(|&v| v * scale).collect()
60}
61
62/// Convert field from rad/s to ppm.
63pub fn rads_to_ppm(field_rads: &[f64], field_strength: f64) -> Vec<f64> {
64    let gamma_hz = 42.576e6;
65    let scale = 1e6 / (2.0 * PI * gamma_hz * field_strength);
66    field_rads.iter().map(|&v| v * scale).collect()
67}
68
69/// Root-sum-of-squares combination of multiple magnitude images.
70pub fn rss_combine(magnitudes: &[&[f64]]) -> Vec<f64> {
71    if magnitudes.is_empty() {
72        return Vec::new();
73    }
74    let n = magnitudes[0].len();
75    let mut combined = vec![0.0f64; n];
76    for mag in magnitudes {
77        for (i, &v) in mag.iter().enumerate() {
78            combined[i] += v * v;
79        }
80    }
81    for v in &mut combined {
82        *v = v.sqrt();
83    }
84    combined
85}
86
87// Re-export mask operations from canonical location
88pub use crate::utils::mask::{erode_mask, dilate_mask};
89
90#[cfg(test)]
91mod tests {
92    use super::*;
93
94    #[test]
95    fn test_scale_phase_already_in_range() {
96        let mut data = vec![-PI, 0.0, PI];
97        let original = data.clone();
98        scale_phase_to_pi(&mut data);
99        for (a, b) in data.iter().zip(original.iter()) {
100            assert!((a - b).abs() < 1e-10);
101        }
102    }
103
104    #[test]
105    fn test_scale_phase_rescales_0_to_4096() {
106        let mut data = vec![0.0, 2048.0, 4096.0];
107        scale_phase_to_pi(&mut data);
108        assert!((data[0] - (-PI)).abs() < 1e-10);
109        assert!((data[2] - PI).abs() < 1e-10);
110        assert!(data[1].abs() < 1e-10);
111    }
112
113    #[test]
114    fn test_scale_phase_nan_replaced() {
115        let mut data = vec![0.0, f64::NAN, 4096.0];
116        scale_phase_to_pi(&mut data);
117        assert!(data[1].is_finite());
118    }
119
120    #[test]
121    fn test_scale_phase_empty() {
122        let mut data: Vec<f64> = vec![];
123        scale_phase_to_pi(&mut data);
124        assert!(data.is_empty());
125    }
126
127    #[test]
128    fn test_hz_to_ppm_3t() {
129        let gamma_hz = 42.576e6;
130        let field = vec![gamma_hz * 3.0];
131        let ppm = hz_to_ppm(&field, 3.0);
132        assert!((ppm[0] - 1e6).abs() < 1.0);
133    }
134
135    #[test]
136    fn test_rads_to_ppm_3t() {
137        let gamma_hz = 42.576e6;
138        let rads = vec![2.0 * PI * gamma_hz * 3.0];
139        let ppm = rads_to_ppm(&rads, 3.0);
140        assert!((ppm[0] - 1e6).abs() < 1.0);
141    }
142
143    #[test]
144    fn test_rss_combine() {
145        let a = vec![3.0, 0.0];
146        let b = vec![4.0, 1.0];
147        let result = rss_combine(&[&a, &b]);
148        assert!((result[0] - 5.0).abs() < 1e-10);
149        assert!((result[1] - 1.0).abs() < 1e-10);
150    }
151
152    #[test]
153    fn test_rss_combine_empty() {
154        let result = rss_combine(&[]);
155        assert!(result.is_empty());
156    }
157
158    #[test]
159    fn test_erode_mask_cube() {
160        let mask = vec![1u8; 27]; // 3x3x3
161        let grid = crate::Grid::new(3, 3, 3, 1.0, 1.0, 1.0);
162        let result = erode_mask(&mask, &grid, 1);
163        let center = 1 + 3 + 9;
164        assert_eq!(result[center], 1);
165        let total: u32 = result.iter().map(|&v| v as u32).sum();
166        assert_eq!(total, 1);
167    }
168
169    #[test]
170    fn test_dilate_mask_single() {
171        let mut mask = vec![0u8; 27];
172        mask[1 + 3 + 9] = 1; // center
173        let grid = crate::Grid::new(3, 3, 3, 1.0, 1.0, 1.0);
174        let result = dilate_mask(&mask, &grid, 1);
175        let total: u32 = result.iter().map(|&v| v as u32).sum();
176        assert_eq!(total, 7); // center + 6 neighbors
177    }
178}
qsm_core/pipeline/phase_utils.rs

qsm_core/pipeline/
phase_utils.rs
