qsm_core/kernels/
laplacian.rs

1//! Laplacian kernel for QSM processing
2//!
3//! Discrete 7-point stencil Laplacian kernel used for Laplacian-based
4//! phase unwrapping and regularization.
5
6/// Generate Laplacian kernel in image space
7///
8/// Creates a 7-point stencil Laplacian centered at (0,0,0) for FFT compatibility.
9/// The stencil is: [-2(1/dx² + 1/dy² + 1/dz²)] at center,
10/// [1/dx²] at ±1 in x, [1/dy²] at ±1 in y, [1/dz²] at ±1 in z.
11///
12/// # Arguments
13/// * `nx`, `ny`, `nz` - Array dimensions
14/// * `vsx`, `vsy`, `vsz` - Voxel sizes in mm
15/// * `negative` - If true, return negative Laplacian
16///
17/// # Returns
18/// Flattened Laplacian kernel array of size nx*ny*nz in C order
19pub fn laplacian_kernel(
20    nx: usize, ny: usize, nz: usize,
21    vsx: f64, vsy: f64, vsz: f64,
22    negative: bool,
23) -> Vec<f64> {
24    let n_total = nx * ny * nz;
25    let mut l = vec![0.0; n_total];
26
27    // Stencil coefficients
28    let hx = 1.0 / (vsx * vsx);
29    let hy = 1.0 / (vsy * vsy);
30    let hz = 1.0 / (vsz * vsz);
31    let center = -2.0 * (hx + hy + hz);
32
33    let sign = if negative { -1.0 } else { 1.0 };
34
35    // Fortran order: idx(i, j, k) = i + j*nx + k*nx*ny
36
37    // Center (0, 0, 0)
38    l[0] = sign * center;
39
40    // x neighbors: (1, 0, 0) and (nx-1, 0, 0)
41    if nx > 1 {
42        l[1] = sign * hx;  // (1, 0, 0)
43        l[nx - 1] = sign * hx;  // (nx-1, 0, 0) wraps to -1
44    }
45
46    // y neighbors: (0, 1, 0) and (0, ny-1, 0)
47    if ny > 1 {
48        l[nx] = sign * hy;  // (0, 1, 0)
49        l[(ny - 1) * nx] = sign * hy;  // (0, ny-1, 0)
50    }
51
52    // z neighbors: (0, 0, 1) and (0, 0, nz-1)
53    if nz > 1 {
54        l[nx * ny] = sign * hz;  // (0, 0, 1)
55        l[(nz - 1) * nx * ny] = sign * hz;  // (0, 0, nz-1)
56    }
57
58    l
59}
60
61/// Generate negative Laplacian kernel (commonly used for gradient regularization)
62pub fn negative_laplacian_kernel(
63    nx: usize, ny: usize, nz: usize,
64    vsx: f64, vsy: f64, vsz: f64,
65) -> Vec<f64> {
66    laplacian_kernel(nx, ny, nz, vsx, vsy, vsz, true)
67}
68
69#[cfg(test)]
70mod tests {
71    use super::*;
72
73    #[test]
74    fn test_laplacian_kernel_sum() {
75        // Laplacian kernel should sum to 0 (for periodic BC)
76        let l = laplacian_kernel(8, 8, 8, 1.0, 1.0, 1.0, false);
77        let sum: f64 = l.iter().sum();
78        assert!(sum.abs() < 1e-10, "Laplacian should sum to 0, got {}", sum);
79    }
80
81    #[test]
82    fn test_laplacian_kernel_center() {
83        let l = laplacian_kernel(8, 8, 8, 1.0, 1.0, 1.0, false);
84        // With voxel size 1, center should be -6
85        assert!((l[0] - (-6.0)).abs() < 1e-10, "Center should be -6, got {}", l[0]);
86    }
87
88    #[test]
89    fn test_negative_laplacian() {
90        let l_pos = laplacian_kernel(8, 8, 8, 1.0, 1.0, 1.0, false);
91        let l_neg = laplacian_kernel(8, 8, 8, 1.0, 1.0, 1.0, true);
92
93        for (p, n) in l_pos.iter().zip(l_neg.iter()) {
94            assert!((p + n).abs() < 1e-10, "Negative should be opposite sign");
95        }
96    }
97}
qsm_core/kernels/laplacian.rs

qsm_core/kernels/
laplacian.rs
