qsm_core/bgremove/
sharp.rs

1//! SHARP background field removal
2//!
3//! Sophisticated Harmonic Artifact Reduction for Phase data.
4//! Uses the spherical mean value property of harmonic functions to
5//! separate local from background fields.
6//!
7//! Reference:
8//! Schweser, F., Deistung, A., Lehr, B.W., Reichenbach, J.R. (2011).
9//! "Quantitative imaging of intrinsic magnetic tissue properties using MRI signal phase."
10//! NeuroImage, 54(4):2789-2807. https://doi.org/10.1016/j.neuroimage.2010.10.070
11//!
12//! Reference implementation: https://github.com/kamesy/QSM.jl
13
14use num_complex::Complex64;
15use crate::Grid;
16use crate::fft::{fft3d, ifft3d, fft_real_kernel};
17use crate::kernels::smv::{smv_kernel, erode_mask_smv};
18
19/// SHARP algorithm parameters
20#[cfg_attr(feature = "introspection", derive(serde::Serialize))]
21#[derive(Clone, Debug)]
22pub struct SharpParams {
23    /// Deconvolution threshold
24    pub threshold: f64,
25    /// SMV kernel radius in mm (default: 6.0)
26    pub radius: f64,
27}
28
29impl Default for SharpParams {
30    fn default() -> Self {
31        Self {
32            threshold: 0.05,
33            radius: 6.0,
34        }
35    }
36}
37
38/// SHARP background field removal
39///
40/// Uses spherical mean value (SMV) filtering to remove background field.
41/// The local field is obtained by deconvolving the SMV-filtered field.
42///
43/// # Arguments
44/// * `field` - Unwrapped total field (nx * ny * nz)
45/// * `mask` - Binary mask (nx * ny * nz), 1 = inside ROI
46/// * `grid` - Volume dimensions and voxel sizes
47/// * `params` - SHARP parameters (threshold, radius in mm)
48///
49/// # Returns
50/// (local_field, eroded_mask)
51pub fn sharp(
52    field: &[f64],
53    mask: &[u8],
54    grid: &Grid,
55    params: &SharpParams,
56) -> (Vec<f64>, Vec<u8>) {
57    sharp_core(field, mask, grid, params.threshold, params.radius)
58}
59
60/// SHARP with an explicit absolute kernel radius (mm).
61///
62/// Internal entry point shared with V-SHARP; the public [`sharp`] wrapper
63/// derives the radius from [`SharpParams`].
64pub(crate) fn sharp_core(
65    field: &[f64],
66    mask: &[u8],
67    grid: &Grid,
68    threshold: f64,
69    radius: f64,
70) -> (Vec<f64>, Vec<u8>) {
71    let (nx, ny, nz) = grid.dims;
72    let n_total = nx * ny * nz;
73
74    // Generate SMV kernel and FFT it
75    let s_kernel = smv_kernel(grid, radius);
76    let s_fft = fft_real_kernel(&s_kernel, nx, ny, nz);
77
78    // Erode mask via SMV convolution
79    let eroded_mask = erode_mask_smv(mask, &s_fft, grid, 1.0 - 1e-7_f64.sqrt());
80
81    // Apply SHARP:
82    // 1. Multiply field by (1-S) in k-space (high-pass filter)
83    // 2. Apply eroded mask
84    // 3. Divide by (1-S) with threshold (deconvolution)
85    // 4. Apply eroded mask
86
87    // FFT of field
88    let mut field_complex: Vec<Complex64> = field.iter()
89        .map(|&x| Complex64::new(x, 0.0))
90        .collect();
91    fft3d(&mut field_complex, nx, ny, nz);
92
93    // High-pass filter: multiply by (1-S)
94    for i in 0..n_total {
95        field_complex[i] *= 1.0 - s_fft[i];
96    }
97
98    // IFFT
99    ifft3d(&mut field_complex, nx, ny, nz);
100
101    // Apply eroded mask
102    for i in 0..n_total {
103        if eroded_mask[i] == 0 {
104            field_complex[i] = Complex64::new(0.0, 0.0);
105        }
106    }
107
108    // FFT again for deconvolution
109    fft3d(&mut field_complex, nx, ny, nz);
110
111    // Deconvolution: divide by (1-S) with threshold
112    for i in 0..n_total {
113        let one_minus_s = 1.0 - s_fft[i];
114        if one_minus_s.abs() < threshold {
115            field_complex[i] = Complex64::new(0.0, 0.0);
116        } else {
117            field_complex[i] /= one_minus_s;
118        }
119    }
120
121    // Final IFFT
122    ifft3d(&mut field_complex, nx, ny, nz);
123
124    // Apply eroded mask and extract real part
125    let local_field: Vec<f64> = field_complex.iter()
126        .enumerate()
127        .map(|(i, c)| if eroded_mask[i] == 1 { c.re } else { 0.0 })
128        .collect();
129
130    (local_field, eroded_mask)
131}
132
133#[cfg(test)]
134mod tests {
135    use super::*;
136
137    #[test]
138    fn test_sharp_zero_field() {
139        // Zero field should give zero local field
140        let n = 16;
141        let field = vec![0.0; n * n * n];
142        let mask = vec![1u8; n * n * n];
143        let grid = Grid::new(n, n, n, 1.0, 1.0, 1.0);
144
145        // Use small radius for small test array
146        let (local, _) = sharp(&field, &mask, &grid, &SharpParams { threshold: 0.05, radius: 2.0 });
147
148        for &val in local.iter() {
149            assert!(val.abs() < 1e-8, "Zero field should give zero local field, got {}", val);
150        }
151    }
152
153    #[test]
154    fn test_sharp_finite() {
155        // Result should be finite (no NaN or Inf)
156        let n = 16;
157        let field: Vec<f64> = (0..n*n*n).map(|i| (i as f64) * 0.01).collect();
158        let mask = vec![1u8; n * n * n];
159        let grid = Grid::new(n, n, n, 1.0, 1.0, 1.0);
160
161        let (local, eroded) = sharp(&field, &mask, &grid, &SharpParams { threshold: 0.05, radius: 2.0 });
162
163        for (i, &val) in local.iter().enumerate() {
164            assert!(val.is_finite(), "Local field should be finite at index {}", i);
165        }
166
167        // Eroded mask should have at least some voxels
168        let eroded_count: usize = eroded.iter().map(|&m| m as usize).sum();
169        assert!(eroded_count > 0, "Eroded mask should have some voxels");
170    }
171}
qsm_core/bgremove/sharp.rs

qsm_core/bgremove/
sharp.rs
