Skip to main content

mtb_align/
threshold.rs

1use crate::{Bitmap, Gray};
2
3/// Distinct sample values an 8-bit image can take.
4const LEVELS: usize = 256;
5
6/// Which population split the threshold bitmap is cut at. Ward uses the median,
7/// falling back to the 17th or 83rd for exposures too dark or light for it.
8#[derive(Clone, Copy, Debug, PartialEq, Eq)]
9pub struct Percentile(u8);
10
11impl Percentile {
12    pub const MEDIAN: Self = Self(50);
13
14    pub const SEVENTEENTH: Self = Self(17);
15    pub const EIGHTY_THIRD: Self = Self(83);
16
17    /// Panics unless `percent <= 100`.
18    pub fn new(percent: u8) -> Self {
19        assert!(percent <= 100, "{percent} is not a percentile");
20        Self(percent)
21    }
22
23    pub fn percent(self) -> u8 {
24        self.0
25    }
26}
27
28impl Default for Percentile {
29    fn default() -> Self {
30        Self::MEDIAN
31    }
32}
33
34/// The sample value the bitmaps are cut at. A `u16` because it can land one past
35/// the range: an image of nothing but 255s gives 256, as it does in OpenCV.
36///
37/// A result of 0 or 256 means every sample fell on one side of the cut, so the
38/// threshold bitmap is empty and the exposure carries nothing to align on. More
39/// than half a frame clipped does this. Callers wanting to detect it can compare
40/// against those two values before aligning.
41pub fn threshold(gray: &Gray, percentile: Percentile) -> u16 {
42    let histogram = histogram(gray.as_slice());
43
44    // OpenCV's scan: the sum is tested before each bin, so the answer is one past
45    // the bin that tipped it.
46    let target = gray.as_slice().len() as u64 * percentile.percent() as u64 / 100;
47    let mut sum = 0;
48    let mut level = 0;
49    while sum < target && level < LEVELS {
50        sum += histogram[level];
51        level += 1;
52    }
53
54    level as u16
55}
56
57/// Counts samples at each of the 256 levels. Counting is commutative, so
58/// threads tally a slice each rather than locking one shared histogram.
59fn histogram(samples: &[u8]) -> [u64; LEVELS] {
60    let tally = |samples: &[u8]| {
61        let mut bins = [0u64; LEVELS];
62        for &sample in samples {
63            bins[sample as usize] += 1;
64        }
65        bins
66    };
67
68    #[cfg(feature = "rayon")]
69    {
70        use rayon::prelude::*;
71
72        // Big enough to amortise a per-slice histogram, small enough not to
73        // leave a long tail.
74        const CHUNK: usize = 1 << 16;
75
76        if samples.len() > CHUNK {
77            return samples.par_chunks(CHUNK).map(tally).reduce(
78                || [0; LEVELS],
79                |mut carried, bins| {
80                    for (total, count) in carried.iter_mut().zip(bins) {
81                        *total += count;
82                    }
83                    carried
84                },
85            );
86        }
87    }
88
89    tally(samples)
90}
91
92/// One exposure reduced to the two bitmaps the search compares. They travel
93/// together: the exclusion bitmap only means anything against its own threshold.
94#[derive(Clone, PartialEq, Eq)]
95pub struct Bitmaps {
96    /// Set where a sample is above the threshold.
97    pub threshold: Bitmap,
98    /// *Clear* within the noise tolerance, where a sample's side is untrustworthy.
99    pub exclusion: Bitmap,
100}
101
102/// Reduces one exposure to its threshold and exclusion bitmaps.
103pub fn compute_bitmaps(gray: &Gray, percentile: Percentile, tolerance: u8) -> Bitmaps {
104    let cut = threshold(gray, percentile);
105    let tolerance = tolerance as u16;
106
107    let (width, height) = (gray.width(), gray.height());
108    let sample = |x, y| gray.sample(x, y) as u16;
109
110    Bitmaps {
111        threshold: Bitmap::packed(width, height, |x, y| sample(x, y) > cut),
112        exclusion: Bitmap::packed(width, height, |x, y| sample(x, y).abs_diff(cut) > tolerance),
113    }
114}
115
116#[cfg(test)]
117mod tests {
118    use super::*;
119
120    /// One pixel at every level, so the threshold has one place it can land.
121    fn ramp() -> Gray {
122        Gray::from_vec((0..=255).collect(), 16, 16)
123    }
124
125    #[test]
126    fn the_median_splits_the_population() {
127        assert_eq!(threshold(&ramp(), Percentile::MEDIAN), 128);
128    }
129
130    #[test]
131    fn a_lower_percentile_cuts_lower() {
132        assert_eq!(threshold(&ramp(), Percentile::SEVENTEENTH), 43);
133        assert_eq!(threshold(&ramp(), Percentile::EIGHTY_THIRD), 212);
134    }
135
136    #[test]
137    fn a_uniform_image_pushes_the_threshold_past_every_sample() {
138        let white = Gray::from_vec(vec![255; 64], 8, 8);
139
140        assert_eq!(threshold(&white, Percentile::MEDIAN), 256);
141    }
142
143    #[test]
144    fn the_threshold_bitmap_marks_only_what_is_above() {
145        let tb = compute_bitmaps(&ramp(), Percentile::MEDIAN, 4).threshold;
146
147        assert_eq!(tb.count_ones(), 127);
148        assert!(!tb.get(128 % 16, 128 / 16), "the threshold itself is below");
149        assert!(tb.get(129 % 16, 129 / 16));
150    }
151
152    /// Ward zeroes it within the tolerance: here 124 to 132, nine values.
153    #[test]
154    fn the_exclusion_bitmap_clears_the_band_around_the_threshold() {
155        let eb = compute_bitmaps(&ramp(), Percentile::MEDIAN, 4).exclusion;
156
157        assert_eq!(eb.count_ones(), 256 - 9);
158        for value in 124..=132usize {
159            assert!(
160                !eb.get(value % 16, value / 16),
161                "{value} should be excluded"
162            );
163        }
164        assert!(eb.get(123 % 16, 123 / 16));
165        assert!(eb.get(133 % 16, 133 / 16));
166    }
167
168    #[test]
169    fn a_zero_tolerance_still_excludes_the_threshold_value() {
170        let eb = compute_bitmaps(&ramp(), Percentile::MEDIAN, 0).exclusion;
171
172        assert_eq!(eb.count_ones(), 255);
173        assert!(!eb.get(128 % 16, 128 / 16));
174    }
175
176    /// The threshold is 256, so every sample is within one step of it.
177    #[test]
178    fn a_uniform_image_excludes_itself_entirely() {
179        let white = Gray::from_vec(vec![255; 64], 8, 8);
180
181        let Bitmaps {
182            threshold: tb,
183            exclusion: eb,
184        } = compute_bitmaps(&white, Percentile::MEDIAN, 4);
185
186        assert_eq!(tb.count_ones(), 0);
187        assert_eq!(eb.count_ones(), 0);
188    }
189
190    #[test]
191    fn the_bitmaps_match_the_image_dimensions() {
192        let Bitmaps {
193            threshold: tb,
194            exclusion: eb,
195        } = compute_bitmaps(&Gray::from_vec(vec![7; 12], 4, 3), Percentile::MEDIAN, 4);
196
197        assert_eq!((tb.width(), tb.height()), (4, 3));
198        assert_eq!((eb.width(), eb.height()), (4, 3));
199    }
200
201    #[test]
202    fn a_percentile_reports_what_it_was_built_with() {
203        assert_eq!(Percentile::new(17), Percentile::SEVENTEENTH);
204        assert_eq!(Percentile::MEDIAN.percent(), 50);
205        assert_eq!(Percentile::default(), Percentile::MEDIAN);
206    }
207
208    #[test]
209    #[should_panic(expected = "101 is not a percentile")]
210    fn a_percentile_past_a_hundred_is_rejected() {
211        Percentile::new(101);
212    }
213}