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

qsm_core/bgremove/
smv.rs
