image_compare/
histogram.rs

1use crate::prelude::*;
2
3const BINS: u8 = u8::MAX;
4
5fn correlation(first_hist: &Histogram, second_hist: &Histogram) -> Option<f64> {
6    let first_mean = first_hist.mean();
7    let second_mean = second_hist.mean();
8    let numerator = (0..=BINS)
9        .map(|i| {
10            (first_hist.get_bin_content(i) - first_mean)
11                * (second_hist.get_bin_content(i) - second_mean)
12        })
13        .sum::<f64>();
14    let denominator = (first_hist.variance() * second_hist.variance()).sqrt();
15
16    if denominator != 0. {
17        Some(numerator / denominator)
18    } else {
19        None
20    }
21}
22
23fn chi_square(first_hist: &Histogram, second_hist: &Histogram) -> Option<f64> {
24    let score = (0..=BINS)
25        .map(|i| {
26            let num = (first_hist.get_bin_content(i) - second_hist.get_bin_content(i)).powi(2);
27            let den = first_hist.get_bin_content(i);
28            if num == 0. {
29                0.
30            } else {
31                num / den
32            }
33        })
34        .sum::<f64>();
35
36    if score.is_nan() || score.is_infinite() {
37        None
38    } else {
39        Some(score)
40    }
41}
42
43fn intersection(first_hist: &Histogram, second_hist: &Histogram) -> f64 {
44    (0..=BINS)
45        .map(|i| {
46            first_hist
47                .get_bin_content(i)
48                .min(second_hist.get_bin_content(i))
49        })
50        .sum::<f64>()
51}
52
53fn hellinger(first_hist: &Histogram, second_hist: &Histogram) -> Option<f64> {
54    let bc = (0..=BINS)
55        .map(|i| (first_hist.get_bin_content(i) * second_hist.get_bin_content(i)).sqrt())
56        .sum::<f64>();
57    let normalization = (first_hist.integral() * second_hist.integral()).sqrt();
58    if normalization == 0. {
59        return None;
60    }
61    let bc_normalized = bc / normalization;
62    Some((1. - bc_normalized).sqrt())
63}
64
65/// The distance metric choices for histogram comparisons
66pub enum Metric {
67    /// <img src="https://render.githubusercontent.com/render/math?math=d(H_1,H_2) = \frac{\sum_I (H_1(I) - \bar{H_1}) (H_2(I) - \bar{H_2})}{\sqrt{\sum_I(H_1(I) - \bar{H_1})^2 \sum_I(H_2(I) - \bar{H_2})^2}}">
68    Correlation,
69    /// <img src="https://render.githubusercontent.com/render/math?math=d(H_1,H_2) = \sum _I \frac{\left(H_1(I)-H_2(I)\right)^2}{H_1(I)}">
70    /// First histogram may not have empty bins
71    ChiSquare,
72    /// <img src="https://render.githubusercontent.com/render/math?math=d(H_1,H_2) = \sum _I \min (H_1(I), H_2(I))">
73    Intersection,
74    /// <img src="https://render.githubusercontent.com/render/math?math=d(H_1,H_2) = \sqrt{1 - \frac{1}{\sqrt{\int{H_1} \int{H_2}}} \sum_I \sqrt{H_1(I) \cdot H_2(I)}}">
75    /// Both histograms need to be normalizable
76    Hellinger,
77}
78
79pub fn img_compare(
80    first: &GrayImage,
81    second: &GrayImage,
82    metric: Metric,
83) -> Result<f64, CompareError> {
84    let first_hist = Histogram::from_gray_image(first);
85    let second_hist = Histogram::from_gray_image(second);
86    let score = match metric {
87        Metric::Correlation => correlation(&first_hist, &second_hist).ok_or_else(|| {
88            CompareError::CalculationFailed(
89                "One or both histograms' variances were zero!".to_owned(),
90            )
91        })?,
92        Metric::ChiSquare => chi_square(&first_hist, &second_hist).ok_or_else(|| {
93            CompareError::CalculationFailed(
94                "First histogram contained empty bins - relative error calculation impossible!"
95                    .to_owned(),
96            )
97        })?,
98        Metric::Intersection => intersection(&first_hist, &second_hist),
99        Metric::Hellinger => hellinger(&first_hist, &second_hist).ok_or_else(|| {
100            CompareError::CalculationFailed(
101                "One or both histograms were not normalizable!".to_owned(),
102            )
103        })?,
104    };
105    Ok(score)
106}
107
108struct Histogram {
109    data: Vec<f64>,
110}
111
112impl Histogram {
113    pub fn get_bin_content(&self, bin: u8) -> f64 {
114        self.data[bin as usize]
115    }
116
117    pub fn from_gray_image(image: &GrayImage) -> Histogram {
118        let mut data = vec![0.; 256];
119        image.pixels().for_each(|p| data[p[0] as usize] += 1.);
120        Histogram { data }
121    }
122
123    pub fn mean(&self) -> f64 {
124        self.data.iter().sum::<f64>() / self.data.len() as f64
125    }
126
127    pub fn variance(&self) -> f64 {
128        let mean = self.mean();
129        self.data.iter().map(|v| (v - mean).powi(2)).sum()
130    }
131
132    pub fn integral(&self) -> f64 {
133        self.data.iter().sum()
134    }
135
136    #[cfg(test)]
137    pub fn from_vec(data: Vec<f64>) -> Option<Histogram> {
138        if data.len() != 256 {
139            None
140        } else {
141            Some(Histogram { data })
142        }
143    }
144}
145
146#[cfg(test)]
147mod tests {
148    use super::*;
149
150    #[test]
151    fn correlation_simple() {
152        let mut first_vec = vec![255.; 256];
153        first_vec[0] = 0.;
154        let mut second_vec = vec![0.; 256];
155        second_vec[0] = 255.;
156
157        let first = Histogram::from_vec(first_vec).unwrap();
158        let second = Histogram::from_vec(second_vec).unwrap();
159
160        assert_eq!(correlation(&first, &first).unwrap(), 1.);
161        assert_eq!(correlation(&first, &second).unwrap(), -1.);
162    }
163
164    #[test]
165    fn correlation_zero_variance() {
166        let first_vec = vec![0.; 256];
167
168        let first = Histogram::from_vec(first_vec).unwrap();
169
170        assert!(correlation(&first, &first).is_none());
171    }
172
173    #[test]
174    fn chi_square_simple() {
175        let first_vec = vec![10.; 256];
176        let second_vec = vec![0.; 256];
177
178        let first = Histogram::from_vec(first_vec).unwrap();
179        let second = Histogram::from_vec(second_vec).unwrap();
180
181        assert_eq!(chi_square(&first, &first).unwrap(), 0.);
182        assert_eq!(chi_square(&first, &second).unwrap(), 10. * 256.);
183    }
184
185    #[test]
186    fn chi_square_edge_cases() {
187        let first_vec = vec![0.; 256];
188        let second_vec = vec![10.; 256];
189
190        let first = Histogram::from_vec(first_vec).unwrap();
191        let second = Histogram::from_vec(second_vec).unwrap();
192
193        assert_eq!(chi_square(&first, &first).unwrap(), 0.);
194        assert!(chi_square(&first, &second).is_none());
195    }
196
197    #[test]
198    fn intersection_simple() {
199        let first_vec = vec![1.; 256];
200        let second_vec = vec![10.; 256];
201
202        let first = Histogram::from_vec(first_vec).unwrap();
203        let second = Histogram::from_vec(second_vec).unwrap();
204
205        assert_eq!(intersection(&first, &first), first.integral());
206        assert_eq!(intersection(&first, &second), first.integral());
207    }
208
209    #[test]
210    fn hellinger_tests() {
211        let first_vec = vec![1.; 256];
212        let second_vec = vec![10.; 256];
213        let mut third_vec = vec![0.; 256];
214        for i in 0..127 {
215            third_vec[i] = 100.;
216        }
217
218        let zeros = vec![0.; 256];
219        let zeros = Histogram::from_vec(zeros).unwrap();
220
221        let first = Histogram::from_vec(first_vec).unwrap();
222        let second = Histogram::from_vec(second_vec).unwrap();
223        let third = Histogram::from_vec(third_vec).unwrap();
224
225        assert_eq!(hellinger(&first, &first).unwrap(), 0.0);
226        assert_eq!(hellinger(&first, &second).unwrap(), 5.1619136559035694e-8);
227        assert_eq!(hellinger(&first, &third).unwrap(), 0.5437469730039513);
228        assert!(hellinger(&first, &zeros).is_none());
229    }
230}