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.
8/// Matches MATLAB's graythresh: operates on all values including zeros,
9/// normalizes to [0,1] range, and returns threshold at bin edge.
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, matching MATLAB's graythresh)
79    min_val + optimal_threshold_bin as f64 * bin_width
80}
81
82#[cfg(test)]
83mod tests {
84    use super::*;
85
86    #[test]
87    fn test_otsu_threshold_bimodal() {
88        // Bimodal distribution with spread around two clusters
89        let mut data = Vec::new();
90        for i in 0..100 {
91            data.push(0.1 + 0.2 * (i as f64 / 100.0)); // 0.1 to 0.3
92        }
93        for i in 0..100 {
94            data.push(0.7 + 0.2 * (i as f64 / 100.0)); // 0.7 to 0.9
95        }
96
97        let threshold = otsu_threshold(&data, 256);
98        assert!(threshold > 0.2 && threshold < 0.8,
99            "Threshold {} should be between the two clusters", threshold);
100    }
101
102    #[test]
103    fn test_otsu_threshold_empty() {
104        assert_eq!(otsu_threshold(&[], 256), 0.0);
105    }
106
107    #[test]
108    fn test_otsu_threshold_constant() {
109        let data = vec![5.0; 100];
110        assert_eq!(otsu_threshold(&data, 256), 5.0);
111    }
112}
qsm_core/utils/threshold.rs

qsm_core/utils/
threshold.rs
