sci_gaussfilt_rs/
lib.rs

1pub mod error;
2pub mod types;
3
4use crate::error::FilterError;
5use crate::types::BoundaryMode;
6use num_traits::{Float, FloatConst, NumAssign};
7
8// --- KERNEL GENERATION ---
9
10pub fn gaussian_kernel<T>(sigma: T, truncate: T) -> Vec<T>
11where
12    T: Float + FloatConst + NumAssign,
13{
14    // SciPy Logic: int(truncate * sigma + 0.5)
15    // We add 0.5 and floor to simulate "Round Half Up" behavior perfectly
16    let radius_f = (truncate * sigma + T::from(0.5).unwrap()).floor();
17    let radius = radius_f.to_isize().unwrap();
18    let size = (2 * radius + 1) as usize;
19
20    if sigma <= T::zero() {
21        let mut k = vec![T::zero(); size];
22        k[radius as usize] = T::one();
23        return k;
24    }
25
26    let mut kernel = Vec::with_capacity(size);
27    let two = T::from(2.0).unwrap();
28    let denom = two * sigma * sigma;
29
30    for i in -radius..=radius {
31        let x = T::from(i).unwrap();
32        kernel.push((-x * x / denom).exp());
33    }
34
35    // Normalize
36    let sum: T = kernel.iter().fold(T::zero(), |acc, &x| acc + x);
37    for v in &mut kernel {
38        *v /= sum;
39    }
40    kernel
41}
42
43// --- PADDING LOGIC ---
44
45fn pad_data<T>(data: &[T], radius: usize, mode: BoundaryMode<T>) -> Vec<T>
46where
47    T: Float + Copy,
48{
49    let len = data.len();
50    let new_len = len + 2 * radius;
51    let mut padded = Vec::with_capacity(new_len);
52
53    // 1. Left Padding
54    match mode {
55        BoundaryMode::Constant(val) => {
56            padded.extend(std::iter::repeat(val).take(radius));
57        }
58        BoundaryMode::Nearest => {
59            let val = data[0];
60            padded.extend(std::iter::repeat(val).take(radius));
61        }
62        BoundaryMode::Reflect => {
63            // SciPy-style reflect: d c b a | a b c ...
64            for i in (1..=radius).rev() {
65                 let src = reflect_index_helper(-(i as isize), len);
66                 padded.push(data[src]);
67            }
68        }
69    }
70
71    // 2. Middle (Copy Data)
72    padded.extend_from_slice(data);
73
74    // 3. Right Padding
75    match mode {
76        BoundaryMode::Constant(val) => {
77            padded.extend(std::iter::repeat(val).take(radius));
78        }
79        BoundaryMode::Nearest => {
80            let val = data[len - 1];
81            padded.extend(std::iter::repeat(val).take(radius));
82        }
83        BoundaryMode::Reflect => {
84             for i in 0..radius {
85                 let src = reflect_index_helper((len + i) as isize, len);
86                 padded.push(data[src]);
87             }
88        }
89    }
90
91    padded
92}
93
94#[inline(always)]
95fn reflect_index_helper(idx: isize, len: usize) -> usize {
96    if len == 0 { return 0; }
97    let n_i = len as isize;
98    let mut i = idx;
99    if i < 0 { i = -i - 1; }
100    let period = 2 * n_i;
101    i = ((i % period) + period) % period;
102    if i >= n_i { (period - 1 - i) as usize } else { i as usize }
103}
104
105// --- CONVOLUTION HELPER (HOT LOOP) ---
106
107/// Inner loop broken out for LLVM optimization.
108/// Marked inline(always) to ensure it dissolves into the caller for SIMD.
109#[inline(always)]
110fn convolve_sample<T>(window: &[T], kernel: &[T]) -> T
111where 
112    T: Float + NumAssign
113{
114    let mut acc = T::zero();
115    let len = kernel.len(); 
116    
117    // SAFETY: The caller ensures window and kernel are the same length.
118    // This removes bounds checks to allow aggressive SIMD auto-vectorization.
119    unsafe {
120        for j in 0..len {
121            acc += *window.get_unchecked(j) * *kernel.get_unchecked(j);
122        }
123    }
124    acc
125}
126
127// --- MAIN ENTRY POINT ---
128
129pub fn gaussian_filter1d<T>(
130    data: &[T],
131    sigma: T,
132    truncate: T,
133    mode: BoundaryMode<T>,
134) -> Result<Vec<T>, FilterError>
135where
136    T: Float + FloatConst + NumAssign + std::fmt::Debug,
137{
138    if data.is_empty() {
139        return Err(FilterError::EmptyInput);
140    }
141    if sigma < T::zero() {
142        return Err(FilterError::InvalidSigma(sigma.to_f64().unwrap()));
143    }
144    // Optimization: If sigma is effectively zero, return copy
145    if sigma.abs() < T::epsilon() {
146        return Ok(data.to_vec());
147    }
148
149    // 1. Prepare Kernel
150    let kernel = gaussian_kernel(sigma, truncate);
151    let radius = (kernel.len() - 1) / 2;
152
153    // 2. Pad Data
154    // Pre-calculating padding allows the inner loop to be branchless
155    let padded = pad_data(data, radius, mode);
156    
157    let mut out = Vec::with_capacity(data.len());
158    let k_len = kernel.len();
159
160    // 3. Convolution
161    // Iterate only over the valid data region
162    for i in 0..data.len() {
163        let window = &padded[i..i + k_len];
164        out.push(convolve_sample(window, &kernel));
165    }
166
167    Ok(out)
168}
sci_gaussfilt_rs/lib.rs

sci_gaussfilt_rs/
lib.rs
