qsm_core/inversion/
tkd.rs

1//! Truncated k-space division (TKD) / TSVD for QSM
2//!
3//! TKD is the simplest dipole inversion method. It directly divides the
4//! field in k-space by the dipole kernel, with truncation to avoid
5//! division by small values near the magic angle.
6//!
7//! Reference:
8//! Shmueli, K., de Zwart, J.A., van Gelderen, P., Li, T.Q., Dodd, S.J., Duyn, J.H. (2009).
9//! "Magnetic susceptibility mapping of brain tissue in vivo using MRI phase data."
10//! Magnetic Resonance in Medicine, 62(6):1510-1522. https://doi.org/10.1002/mrm.22135
11//!
12//! Reference implementation: https://github.com/kamesy/QSM.jl
13
14use num_complex::Complex64;
15use crate::fft::{fft3d, ifft3d};
16use crate::kernels::dipole::dipole_kernel;
17use crate::utils::apply_mask_zero;
18use crate::Grid;
19
20/// TKD algorithm parameters
21#[cfg_attr(feature = "introspection", derive(serde::Serialize))]
22#[derive(Clone, Debug)]
23pub struct TkdParams {
24    /// Truncation threshold (typically 0.1-0.2)
25    pub threshold: f64,
26}
27
28impl Default for TkdParams {
29    fn default() -> Self {
30        Self { threshold: 0.15 }
31    }
32}
33
34/// Truncated k-space division (TKD) for dipole inversion
35///
36/// Computes susceptibility map from local field using direct k-space division
37/// with threshold-based truncation.
38///
39/// # Arguments
40/// * `local_field` - Local field values (unwrapped phase / TE / gamma / B0)
41/// * `mask` - Binary mask (1 = inside ROI, 0 = outside)
42/// * `grid` - Volume grid (dimensions and voxel sizes)
43/// * `bdir` - B0 field direction (bx, by, bz)
44/// * `params` - TKD parameters (truncation threshold)
45///
46/// # Returns
47/// Susceptibility map (same size as input)
48pub fn tkd(
49    local_field: &[f64],
50    mask: &[u8],
51    grid: &Grid,
52    bdir: (f64, f64, f64),
53    params: &TkdParams,
54) -> Vec<f64> {
55    let (nx, ny, nz) = grid.dims;
56    let n_total = grid.n_total();
57    let threshold = params.threshold;
58
59    // Generate dipole kernel
60    let d = dipole_kernel(grid, bdir);
61
62    // Compute inverse dipole kernel with truncation
63    // TKD: if |D| <= threshold, use sign(D)/threshold; else use 1/D
64    let inv_threshold = 1.0 / threshold;
65    let inv_d: Vec<f64> = d.iter().map(|&dval| {
66        if dval.abs() <= threshold {
67            // Truncate: use sign(D)/threshold
68            if dval >= 0.0 { inv_threshold } else { -inv_threshold }
69        } else {
70            // Normal inverse
71            1.0 / dval
72        }
73    }).collect();
74
75    // Convert local field to complex
76    let mut field_complex: Vec<Complex64> = local_field.iter()
77        .map(|&x| Complex64::new(x, 0.0))
78        .collect();
79
80    // FFT of local field
81    fft3d(&mut field_complex, nx, ny, nz);
82
83    // Multiply by inverse dipole kernel
84    for i in 0..n_total {
85        field_complex[i] *= inv_d[i];
86    }
87
88    // IFFT to get susceptibility
89    ifft3d(&mut field_complex, nx, ny, nz);
90
91    // Extract real part and apply mask
92    let mut chi: Vec<f64> = field_complex.iter()
93        .map(|c| c.re)
94        .collect();
95
96    apply_mask_zero(&mut chi, mask);
97
98    chi
99}
100
101/// Truncated singular value decomposition (TSVD) variant
102///
103/// Similar to TKD but zeros out values below threshold instead of truncating.
104/// This produces smoother results but may have more artifacts at the magic angle.
105pub fn tsvd(
106    local_field: &[f64],
107    mask: &[u8],
108    grid: &Grid,
109    bdir: (f64, f64, f64),
110    params: &TkdParams,
111) -> Vec<f64> {
112    let (nx, ny, nz) = grid.dims;
113    let n_total = grid.n_total();
114    let threshold = params.threshold;
115
116    // Generate dipole kernel
117    let d = dipole_kernel(grid, bdir);
118
119    // Compute inverse dipole kernel with TSVD truncation
120    // TSVD: if |D| <= threshold, use 0; else use 1/D
121    let inv_d: Vec<f64> = d.iter().map(|&dval| {
122        if dval.abs() <= threshold {
123            0.0  // Zero out small values
124        } else {
125            1.0 / dval
126        }
127    }).collect();
128
129    // Convert local field to complex
130    let mut field_complex: Vec<Complex64> = local_field.iter()
131        .map(|&x| Complex64::new(x, 0.0))
132        .collect();
133
134    // FFT of local field
135    fft3d(&mut field_complex, nx, ny, nz);
136
137    // Multiply by inverse dipole kernel
138    for i in 0..n_total {
139        field_complex[i] *= inv_d[i];
140    }
141
142    // IFFT to get susceptibility
143    ifft3d(&mut field_complex, nx, ny, nz);
144
145    // Extract real part and apply mask
146    let mut chi: Vec<f64> = field_complex.iter()
147        .map(|c| c.re)
148        .collect();
149
150    apply_mask_zero(&mut chi, mask);
151
152    chi
153}
154
155#[cfg(test)]
156mod tests {
157    use super::*;
158
159    #[test]
160    fn test_tkd_zero_field() {
161        // Zero field should give zero susceptibility
162        let n = 8;
163        let field = vec![0.0; n * n * n];
164        let mask = vec![1u8; n * n * n];
165        let grid = Grid::new(n, n, n, 1.0, 1.0, 1.0);
166
167        let chi = tkd(&field, &mask, &grid, (0.0, 0.0, 1.0), &TkdParams { threshold: 0.15 });
168
169        for val in chi.iter() {
170            assert!(val.abs() < 1e-10, "Zero field should give zero chi");
171        }
172    }
173
174    #[test]
175    fn test_tkd_mask() {
176        // Values outside mask should be zero
177        let n = 8;
178        let field = vec![1.0; n * n * n];
179        let mut mask = vec![1u8; n * n * n];
180        let grid = Grid::new(n, n, n, 1.0, 1.0, 1.0);
181
182        // Set some voxels outside mask
183        mask[0] = 0;
184        mask[1] = 0;
185
186        let chi = tkd(&field, &mask, &grid, (0.0, 0.0, 1.0), &TkdParams { threshold: 0.15 });
187
188        assert_eq!(chi[0], 0.0, "Outside mask should be 0");
189        assert_eq!(chi[1], 0.0, "Outside mask should be 0");
190    }
191
192    #[test]
193    fn test_tkd_finite() {
194        // Result should be finite (no NaN or Inf)
195        let n = 8;
196        let field: Vec<f64> = (0..n*n*n).map(|i| (i as f64) * 0.01).collect();
197        let mask = vec![1u8; n * n * n];
198        let grid = Grid::new(n, n, n, 1.0, 1.0, 1.0);
199
200        let chi = tkd(&field, &mask, &grid, (0.0, 0.0, 1.0), &TkdParams { threshold: 0.15 });
201
202        for (i, val) in chi.iter().enumerate() {
203            assert!(val.is_finite(), "Chi should be finite at index {}", i);
204        }
205    }
206
207    #[test]
208    fn test_tsvd_vs_tkd() {
209        // TSVD and TKD should give different results near magic angle
210        let n = 16;
211        let field: Vec<f64> = (0..n*n*n).map(|i| ((i as f64) * 0.1).sin()).collect();
212        let mask = vec![1u8; n * n * n];
213        let grid = Grid::new(n, n, n, 1.0, 1.0, 1.0);
214
215        let chi_tkd = tkd(&field, &mask, &grid, (0.0, 0.0, 1.0), &TkdParams { threshold: 0.15 });
216        let chi_tsvd = tsvd(&field, &mask, &grid, (0.0, 0.0, 1.0), &TkdParams { threshold: 0.15 });
217
218        // They should be different
219        let diff: f64 = chi_tkd.iter().zip(chi_tsvd.iter())
220            .map(|(a, b)| (a - b).abs())
221            .sum();
222
223        assert!(diff > 1e-10, "TKD and TSVD should give different results");
224    }
225}
qsm_core/inversion/tkd.rs

qsm_core/inversion/
tkd.rs
