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
7use crate::constants::GAMMA_BAR_HZ_PER_T;
8
9/// Scale phase data to [-pi, pi] range in-place.
10///
11/// If data is already approximately in [-pi, pi] (within 10% tolerance),
12/// leaves it unchanged. Otherwise linearly maps from [min, max] to [-pi, pi].
13/// NaN/Inf values are replaced with 0.
14pub fn scale_phase_to_pi(data: &mut [f64]) {
15    if data.is_empty() {
16        return;
17    }
18
19    let mut min_val = f64::INFINITY;
20    let mut max_val = f64::NEG_INFINITY;
21    for &v in data.iter() {
22        if v.is_finite() {
23            if v < min_val {
24                min_val = v;
25            }
26            if v > max_val {
27                max_val = v;
28            }
29        }
30    }
31
32    // Replace non-finite values with 0
33    for v in data.iter_mut() {
34        if !v.is_finite() {
35            *v = 0.0;
36        }
37    }
38
39    let range = max_val - min_val;
40    if range < 1e-10 {
41        return;
42    }
43
44    // Check if already approximately in [-pi, pi]
45    let tol = 0.1 * PI;
46    if (min_val + PI).abs() < tol && (max_val - PI).abs() < tol {
47        return;
48    }
49
50    // Linearly rescale to [-pi, pi]
51    let scale = 2.0 * PI / range;
52    for v in data.iter_mut() {
53        *v = (*v - min_val) * scale - PI;
54    }
55}
56
57/// Convert field from Hz to ppm.
58pub fn hz_to_ppm(field_hz: &[f64], field_strength: f64) -> Vec<f64> {
59    let gamma_hz = GAMMA_BAR_HZ_PER_T;
60    let scale = 1e6 / (gamma_hz * field_strength);
61    field_hz.iter().map(|&v| v * scale).collect()
62}
63
64/// Convert field from rad/s to ppm.
65pub fn rads_to_ppm(field_rads: &[f64], field_strength: f64) -> Vec<f64> {
66    let gamma_hz = GAMMA_BAR_HZ_PER_T;
67    let scale = 1e6 / (2.0 * PI * gamma_hz * field_strength);
68    field_rads.iter().map(|&v| v * scale).collect()
69}
70
71/// Convert field from ppm to Hz: the inverse of [`hz_to_ppm`].
72pub fn ppm_to_hz(field_ppm: &[f64], field_strength: f64) -> Vec<f64> {
73    let scale = GAMMA_BAR_HZ_PER_T * field_strength * 1e-6;
74    field_ppm.iter().map(|&v| v * scale).collect()
75}
76
77/// Convert field from ppm to rad/s: the inverse of [`rads_to_ppm`].
78pub fn ppm_to_rads(field_ppm: &[f64], field_strength: f64) -> Vec<f64> {
79    let scale = 2.0 * PI * GAMMA_BAR_HZ_PER_T * field_strength * 1e-6;
80    field_ppm.iter().map(|&v| v * scale).collect()
81}
82
83/// Root-sum-of-squares combination of multiple magnitude images.
84///
85/// Every echo must hold the same number of voxels. The combined volume inherits
86/// that size and is then indexed by the run's grid downstream (bias correction,
87/// masking), so an echo of a different matrix size would otherwise either be
88/// silently truncated or make the result the wrong size for every later stage.
89///
90/// # Errors
91/// If the echoes do not all have the same length. No echoes gives an empty volume.
92pub fn rss_combine(magnitudes: &[&[f64]]) -> Result<Vec<f64>, crate::error::InputError> {
93    if magnitudes.is_empty() {
94        return Ok(Vec::new());
95    }
96    let n = magnitudes[0].len();
97    for (e, mag) in magnitudes.iter().enumerate() {
98        if mag.len() != n {
99            return Err(crate::error::InputError(format!(
100                "rss_combine: echo {e} has {} voxels but echo 0 has {n} \
101                 (all echoes must come from the same acquisition)",
102                mag.len()
103            )));
104        }
105    }
106    let mut combined = vec![0.0f64; n];
107    for mag in magnitudes {
108        for (i, &v) in mag.iter().enumerate() {
109            combined[i] += v * v;
110        }
111    }
112    for v in &mut combined {
113        *v = v.sqrt();
114    }
115    Ok(combined)
116}
117
118#[cfg(test)]
119mod tests {
120    use super::*;
121
122    #[test]
123    fn test_scale_phase_already_in_range() {
124        let mut data = vec![-PI, 0.0, PI];
125        let original = data.clone();
126        scale_phase_to_pi(&mut data);
127        for (a, b) in data.iter().zip(original.iter()) {
128            assert!((a - b).abs() < 1e-10);
129        }
130    }
131
132    #[test]
133    fn test_scale_phase_rescales_0_to_4096() {
134        let mut data = vec![0.0, 2048.0, 4096.0];
135        scale_phase_to_pi(&mut data);
136        assert!((data[0] - (-PI)).abs() < 1e-10);
137        assert!((data[2] - PI).abs() < 1e-10);
138        assert!(data[1].abs() < 1e-10);
139    }
140
141    #[test]
142    fn test_scale_phase_nan_replaced() {
143        let mut data = vec![0.0, f64::NAN, 4096.0];
144        scale_phase_to_pi(&mut data);
145        assert!(data[1].is_finite());
146    }
147
148    #[test]
149    fn test_scale_phase_empty() {
150        let mut data: Vec<f64> = vec![];
151        scale_phase_to_pi(&mut data);
152        assert!(data.is_empty());
153    }
154
155    /// A constant volume has no range to map: it is left as it is.
156    #[test]
157    fn test_scale_phase_constant_is_left_alone() {
158        let mut data = vec![5.0, 5.0, 5.0];
159        scale_phase_to_pi(&mut data);
160        assert!(data.iter().all(|v| (v - 5.0).abs() < 1e-10), "{data:?}");
161    }
162
163    /// Nothing finite: every value becomes 0, and there is then no range to map.
164    #[test]
165    fn test_scale_phase_all_nan_is_zero() {
166        let mut data = vec![f64::NAN, f64::INFINITY, f64::NAN];
167        scale_phase_to_pi(&mut data);
168        assert_eq!(data, vec![0.0, 0.0, 0.0]);
169    }
170
171    #[test]
172    fn test_hz_to_ppm_3t() {
173        let gamma_hz = GAMMA_BAR_HZ_PER_T;
174        let field = vec![gamma_hz * 3.0];
175        let ppm = hz_to_ppm(&field, 3.0);
176        assert!((ppm[0] - 1e6).abs() < 1.0);
177    }
178
179    #[test]
180    fn test_rads_to_ppm_3t() {
181        let gamma_hz = GAMMA_BAR_HZ_PER_T;
182        let rads = vec![2.0 * PI * gamma_hz * 3.0];
183        let ppm = rads_to_ppm(&rads, 3.0);
184        assert!((ppm[0] - 1e6).abs() < 1.0);
185    }
186
187    #[test]
188    fn ppm_hz_rads_round_trip() {
189        let ppm = vec![-0.25, 0.0, 1.5];
190        for (a, b) in hz_to_ppm(&ppm_to_hz(&ppm, 7.0), 7.0).iter().zip(&ppm) {
191            assert!((a - b).abs() < 1e-12);
192        }
193        for (a, b) in rads_to_ppm(&ppm_to_rads(&ppm, 3.0), 3.0).iter().zip(&ppm) {
194            assert!((a - b).abs() < 1e-12);
195        }
196        // 1 ppm at 3 T is 127.728 Hz.
197        assert!((ppm_to_hz(&[1.0], 3.0)[0] - 127.728).abs() < 1e-9);
198    }
199
200    #[test]
201    fn test_rss_combine() {
202        let a = vec![3.0, 0.0];
203        let b = vec![4.0, 1.0];
204        let result = rss_combine(&[&a, &b]).unwrap();
205        assert!((result[0] - 5.0).abs() < 1e-10);
206        assert!((result[1] - 1.0).abs() < 1e-10);
207    }
208
209    #[test]
210    fn test_rss_combine_empty() {
211        let result = rss_combine(&[]).unwrap();
212        assert!(result.is_empty());
213    }
214
215    // Issue #70: the combined volume inherits echo 0's size and is then indexed
216    // by the run's grid, so a mismatched echo either panicked on an index or made
217    // every later stage read past the end of the result. Now it is an error.
218    #[test]
219    fn rss_combine_refuses_a_longer_echo() {
220        let e0 = vec![1.0; 512];
221        let e1 = vec![1.0; 576];
222        let err = rss_combine(&[&e0, &e1]).unwrap_err();
223        assert!(err.to_string().contains("rss_combine: echo 1 has 576 voxels but echo 0 has 512"), "{err}");
224    }
225
226    #[test]
227    fn rss_combine_refuses_a_shorter_echo() {
228        // The short case never panicked - it silently summed only the part of
229        // the stack that echo 1 covered.
230        let e0 = vec![1.0; 512];
231        let e1 = vec![1.0; 216];
232        let err = rss_combine(&[&e0, &e1]).unwrap_err();
233        assert!(err.to_string().contains("rss_combine: echo 1 has 216 voxels but echo 0 has 512"), "{err}");
234    }
235}
qsm_core/pipeline/phase_utils.rs

qsm_core/pipeline/
phase_utils.rs
