qsm_core/utils/
threshold.rs

1//! Automatic thresholding algorithms
2//!
3//! Provides Otsu's method for computing optimal thresholds on bimodal histograms.
4
5/// Otsu's method for automatic threshold selection
6///
7/// Finds the threshold that maximizes inter-class variance over a histogram spanning the data's
8/// own range (all values, zeros included), and returns the threshold at a bin edge. For MATLAB's
9/// `graythresh` on double data, which does not rescale, see [`graythresh`].
10///
11/// # Arguments
12/// * `data` - Input data (e.g. flattened 3D image)
13/// * `num_bins` - Number of histogram bins (typically 256)
14///
15/// # Returns
16/// The optimal threshold value
17pub fn otsu_threshold(data: &[f64], num_bins: usize) -> f64 {
18    if data.is_empty() {
19        return 0.0;
20    }
21
22    let min_val = data.iter().fold(f64::MAX, |a, &b| a.min(b));
23    let max_val = data.iter().fold(f64::MIN, |a, &b| a.max(b));
24
25    if (max_val - min_val).abs() < 1e-10 {
26        return min_val;
27    }
28
29    // Build histogram over full range including zeros
30    let bin_width = (max_val - min_val) / num_bins as f64;
31    let mut histogram = vec![0usize; num_bins];
32
33    for &v in data {
34        let bin = ((v - min_val) / bin_width).floor() as usize;
35        let bin = bin.min(num_bins - 1);
36        histogram[bin] += 1;
37    }
38
39    let total_pixels = data.len() as f64;
40
41    // Compute cumulative sums
42    let mut sum_total = 0.0;
43    for (i, &count) in histogram.iter().enumerate() {
44        sum_total += i as f64 * count as f64;
45    }
46
47    let mut sum_background = 0.0;
48    let mut weight_background = 0.0;
49    let mut max_variance = 0.0;
50    let mut optimal_threshold_bin = 0;
51
52    for (t, &count) in histogram.iter().enumerate() {
53        weight_background += count as f64;
54        if weight_background == 0.0 {
55            continue;
56        }
57
58        let weight_foreground = total_pixels - weight_background;
59        if weight_foreground == 0.0 {
60            break;
61        }
62
63        sum_background += t as f64 * count as f64;
64
65        let mean_background = sum_background / weight_background;
66        let mean_foreground = (sum_total - sum_background) / weight_foreground;
67
68        // Inter-class variance
69        let variance = weight_background * weight_foreground
70            * (mean_background - mean_foreground).powi(2);
71
72        if variance > max_variance {
73            max_variance = variance;
74            optimal_threshold_bin = t;
75        }
76    }
77
78    // Convert bin to threshold value (bin edge)
79    min_val + optimal_threshold_bin as f64 * bin_width
80}
81
82/// MATLAB's `graythresh` on double data: Otsu's threshold, in [0, 1].
83///
84/// As `graythresh` does, the data are converted with `im2uint8` (clipped to [0, 1], scaled by
85/// 255 and rounded; NaN becomes 0), not rescaled to their own range, and Otsu's method runs on
86/// the 256-bin histogram of every value (`otsuthresh`: ties take the mean of the maximising bins;
87/// a non-finite maximum gives 0).
88pub fn graythresh(data: &[f64]) -> f64 {
89    let mut counts = [0f64; 256];
90    for &v in data {
91        let q = if v.is_nan() { 0.0 } else { (v.clamp(0.0, 1.0) * 255.0).round() };
92        counts[q as usize] += 1.0;
93    }
94    let total: f64 = counts.iter().sum();
95    if total == 0.0 {
96        return 0.0;
97    }
98    let (mut omega, mut mu) = ([0f64; 256], [0f64; 256]);
99    let (mut o, mut m) = (0.0, 0.0);
100    for i in 0..256 {
101        let p = counts[i] / total;
102        o += p;
103        m += p * (i + 1) as f64;
104        omega[i] = o;
105        mu[i] = m;
106    }
107    let mu_t = mu[255];
108    let sigma: Vec<f64> = (0..256)
109        .map(|i| (mu_t * omega[i] - mu[i]).powi(2) / (omega[i] * (1.0 - omega[i])))
110        .collect();
111    let max = sigma.iter().copied().filter(|v| !v.is_nan()).fold(f64::NEG_INFINITY, f64::max);
112    if !max.is_finite() {
113        return 0.0;
114    }
115    let idx: Vec<f64> = (0..256).filter(|&i| sigma[i] == max).map(|i| (i + 1) as f64).collect();
116    let idx = idx.iter().sum::<f64>() / idx.len() as f64;
117    (idx - 1.0) / 255.0
118}
119
120#[cfg(test)]
121mod tests {
122    use super::*;
123
124    #[test]
125    fn graythresh_matches_matlab() {
126        // MATLAB R2023b graythresh on the same vectors.
127        let a = [0.0, 0.0, 0.0, 0.1, 0.2, 0.8, 0.9, 0.95, 1.2, -0.1, f64::NAN, 0.5, 0.51];
128        assert!((graythresh(&a) - 0.349019607843137).abs() < 1e-14);
129        let b: Vec<f64> = std::iter::repeat_n(0.0, 1000)
130            .chain((1..=50).map(|i| 0.02 * i as f64))
131            .chain((1..=30).map(|i| 0.6 + 0.003 * i as f64))
132            .collect();
133        assert!((graythresh(&b) - 0.329411764705882).abs() < 1e-14);
134        assert_eq!(graythresh(&[0.3; 10]), 0.0);
135        // A tie across bins 1..254 takes their mean.
136        assert!((graythresh(&[0.004, 0.005, 0.006, 0.996]) - 0.5).abs() < 1e-14);
137    }
138
139    #[test]
140    fn test_otsu_threshold_bimodal() {
141        // Bimodal distribution with spread around two clusters
142        let mut data = Vec::new();
143        for i in 0..100 {
144            data.push(0.1 + 0.2 * (i as f64 / 100.0)); // 0.1 to 0.3
145        }
146        for i in 0..100 {
147            data.push(0.7 + 0.2 * (i as f64 / 100.0)); // 0.7 to 0.9
148        }
149
150        let threshold = otsu_threshold(&data, 256);
151        assert!(threshold > 0.2 && threshold < 0.8,
152            "Threshold {} should be between the two clusters", threshold);
153    }
154
155    #[test]
156    fn test_otsu_threshold_empty() {
157        assert_eq!(otsu_threshold(&[], 256), 0.0);
158    }
159
160    #[test]
161    fn test_otsu_threshold_constant() {
162        let data = vec![5.0; 100];
163        assert_eq!(otsu_threshold(&data, 256), 5.0);
164    }
165}
qsm_core/utils/threshold.rs

qsm_core/utils/
threshold.rs
