phasesmith_core/
background.rs

1//! Deterministic background-estimation kernels.
2//!
3//! Background estimation is preprocessing, not a refinable profile term. The
4//! Bruckner smoother intentionally preserves the scan and boundary semantics of
5//! xypattern's pinned Cython implementation.
6
7use std::error::Error;
8use std::fmt::{Display, Formatter};
9
10/// Validation failures for background-estimation kernels.
11#[derive(Clone, Debug, PartialEq, Eq)]
12pub enum BackgroundError {
13    /// At least one intensity sample is required for endpoint padding.
14    EmptyInput,
15    /// One input sample is NaN or infinite.
16    NonFiniteSample {
17        /// Zero-based location of the invalid intensity.
18        index: usize,
19    },
20    /// Requested padding cannot be represented by the platform index type.
21    SizeOverflow,
22    /// The padded work buffer cannot be allocated.
23    AllocationFailed,
24}
25
26impl Display for BackgroundError {
27    fn fmt(&self, formatter: &mut Formatter<'_>) -> std::fmt::Result {
28        match self {
29            Self::EmptyInput => formatter.write_str("background input must not be empty"),
30            Self::NonFiniteSample { index } => {
31                write!(formatter, "background input sample {index} must be finite")
32            }
33            Self::SizeOverflow => formatter.write_str("background smoothing size overflow"),
34            Self::AllocationFailed => {
35                formatter.write_str("background smoothing work buffer allocation failed")
36            }
37        }
38    }
39}
40
41impl Error for BackgroundError {}
42
43/// Reproduce the pinned xypattern Cython Smooth Bruckner estimator.
44///
45/// `smooth_points` is the half-window in samples, so each moving mean contains
46/// `2 * smooth_points + 1` values. Each iteration scans extended indices
47/// `smooth_points..y.len() - smooth_points - 2`; this intentionally leaves the
48/// same trailing samples unchanged as the compatibility implementation.
49///
50/// # Errors
51///
52/// Returns [`BackgroundError`] when `y` is empty, contains a non-finite sample,
53/// or the padded allocation length overflows `usize`.
54pub fn smooth_bruckner(
55    y: &[f64],
56    smooth_points: usize,
57    iterations: usize,
58) -> Result<Vec<f64>, BackgroundError> {
59    if y.is_empty() {
60        return Err(BackgroundError::EmptyInput);
61    }
62    if let Some((index, _)) = y.iter().enumerate().find(|(_, value)| !value.is_finite()) {
63        return Err(BackgroundError::NonFiniteSample { index });
64    }
65
66    let double_padding = smooth_points
67        .checked_mul(2)
68        .ok_or(BackgroundError::SizeOverflow)?;
69    let extended_len = y
70        .len()
71        .checked_add(double_padding)
72        .ok_or(BackgroundError::SizeOverflow)?;
73    let window_len = double_padding
74        .checked_add(1)
75        .ok_or(BackgroundError::SizeOverflow)?;
76    let window_len_u32 = u32::try_from(window_len).map_err(|_| BackgroundError::SizeOverflow)?;
77
78    let mut extended = Vec::new();
79    extended
80        .try_reserve_exact(extended_len)
81        .map_err(|_| BackgroundError::AllocationFailed)?;
82    extended.resize(extended_len, 0.0);
83    extended[..smooth_points].fill(y[0]);
84    extended[smooth_points..smooth_points + y.len()].copy_from_slice(y);
85    extended[smooth_points + y.len()..].fill(y[y.len() - 1]);
86
87    let scan_end = y.len().saturating_sub(smooth_points.saturating_add(2));
88    let window_scale = f64::from(window_len_u32);
89    for _ in 0..iterations {
90        let mut window_sum = 0.0;
91        for value in &extended[..window_len] {
92            window_sum += value;
93        }
94        let mut window_average = window_sum / window_scale;
95
96        for index in smooth_points..scan_end {
97            let outgoing = extended[index - smooth_points];
98            let incoming = extended[index + smooth_points + 1];
99            if extended[index] > window_average {
100                let old_value = extended[index];
101                extended[index] = window_average;
102                window_average +=
103                    ((window_average - old_value) + (incoming - outgoing)) / window_scale;
104            } else {
105                window_average += (incoming - outgoing) / window_scale;
106            }
107        }
108    }
109
110    Ok(extended[smooth_points..smooth_points + y.len()].to_vec())
111}
112
113#[cfg(test)]
114mod tests {
115    use super::*;
116
117    #[test]
118    fn zero_iterations_only_round_trips_the_padded_center() {
119        let y = [1.0, 5.0, 2.0, 9.0];
120        assert_eq!(smooth_bruckner(&y, 3, 0), Ok(y.to_vec()));
121    }
122
123    #[test]
124    fn zero_half_window_preserves_every_sample() {
125        let y = [1.0, 5.0, 2.0, 9.0];
126        assert_eq!(smooth_bruckner(&y, 0, 20), Ok(y.to_vec()));
127    }
128
129    #[test]
130    fn clipping_and_trailing_compatibility_range_are_explicit() {
131        let y = [0.0, 0.0, 9.0, 0.0, 0.0, 7.0, 8.0, 9.0, 10.0];
132        let actual = smooth_bruckner(&y, 1, 1).expect("valid smoother input");
133        assert_eq!(actual, vec![0.0, 0.0, 3.0, 0.0, 0.0, 7.0, 8.0, 9.0, 10.0]);
134    }
135
136    #[test]
137    fn window_larger_than_signal_is_a_safe_noop() {
138        let y = [1.0, 4.0, 2.0];
139        assert_eq!(smooth_bruckner(&y, 20, 5), Ok(y.to_vec()));
140    }
141
142    #[test]
143    fn invalid_samples_are_rejected_with_their_index() {
144        assert_eq!(smooth_bruckner(&[], 1, 1), Err(BackgroundError::EmptyInput));
145        assert_eq!(
146            smooth_bruckner(&[1.0, f64::NAN], 1, 1),
147            Err(BackgroundError::NonFiniteSample { index: 1 })
148        );
149    }
150}
phasesmith_core/background.rs

phasesmith_core/
background.rs
