qsm_core/bgremove/
smv.rs

1//! Simple Spherical Mean Value (SMV) background field removal
2//!
3//! Basic SMV filtering: subtracts the spherical mean of the field.
4//! Simpler than SHARP (no deconvolution step).
5//!
6//! local_field = field - SMV(field)
7//!
8//! Reference:
9//! Schweser, F., Deistung, A., Lehr, B.W., Reichenbach, J.R. (2011).
10//! "Quantitative imaging of intrinsic magnetic tissue properties using MRI signal phase."
11//! NeuroImage, 54(4):2789-2807. https://doi.org/10.1016/j.neuroimage.2010.10.070
12//!
13//! Reference implementation: https://github.com/kamesy/QSM.jl
14
15use crate::Grid;
16use crate::fft::fft_real_kernel;
17use crate::kernels::smv::{smv_kernel, erode_mask_smv};
18
19/// Simple SMV background field removal
20///
21/// Computes: local_field = field - SMV(field)
22///
23/// # Arguments
24/// * `field` - Unwrapped total field (nx * ny * nz)
25/// * `mask` - Binary mask (nx * ny * nz), 1 = inside ROI
26/// * `grid` - Volume dimensions and voxel sizes
27/// * `radius` - SMV kernel radius in mm
28///
29/// # Returns
30/// (local_field, eroded_mask)
31pub fn smv(
32    field: &[f64],
33    mask: &[u8],
34    grid: &Grid,
35    radius: f64,
36) -> (Vec<f64>, Vec<u8>) {
37    let (nx, ny, nz) = grid.dims;
38
39    // Generate SMV kernel and FFT it
40    let s_kernel = smv_kernel(grid, radius);
41    let s_fft = fft_real_kernel(&s_kernel, nx, ny, nz);
42
43    // Erode mask via SMV convolution
44    let eroded_mask = erode_mask_smv(mask, &s_fft, grid, 1.0 - 1e-10);
45
46    // Compute SMV(field) = background field estimate
47    let smv_field = crate::fft::apply_real_kernel(field, &s_fft, nx, ny, nz);
48
49    // Local field = field - SMV(field), within eroded mask
50    let local_field: Vec<f64> = field.iter()
51        .zip(smv_field.iter())
52        .enumerate()
53        .map(|(i, (&f, &smv_f))| {
54            if eroded_mask[i] == 1 { f - smv_f } else { 0.0 }
55        })
56        .collect();
57
58    (local_field, eroded_mask)
59}
60
61#[cfg(test)]
62mod tests {
63    use super::*;
64
65    #[test]
66    fn test_smv_zero_field() {
67        let n = 16;
68        let field = vec![0.0; n * n * n];
69        let mask = vec![1u8; n * n * n];
70        let grid = Grid::new(n, n, n, 1.0, 1.0, 1.0);
71
72        let (local, _) = smv(&field, &mask, &grid, 2.0);
73
74        for &val in local.iter() {
75            assert!(val.abs() < 1e-8, "Zero field should give zero local field, got {}", val);
76        }
77    }
78
79    #[test]
80    fn test_smv_finite() {
81        let n = 16;
82        let field: Vec<f64> = (0..n*n*n).map(|i| (i as f64) * 0.01).collect();
83        let mask = vec![1u8; n * n * n];
84        let grid = Grid::new(n, n, n, 1.0, 1.0, 1.0);
85
86        let (local, eroded) = smv(&field, &mask, &grid, 2.0);
87
88        for (i, &val) in local.iter().enumerate() {
89            assert!(val.is_finite(), "Local field should be finite at index {}", i);
90        }
91
92        let eroded_count: usize = eroded.iter().map(|&m| m as usize).sum();
93        assert!(eroded_count > 0, "Eroded mask should have some voxels");
94    }
95}
qsm_core/bgremove/smv.rs

qsm_core/bgremove/
smv.rs
