imageproc/
stats.rs

1//! Statistical properties of images.
2
3use crate::definitions::Image;
4use crate::math::cast;
5use conv::ValueInto;
6use image::{GenericImageView, GrayImage, Pixel, Primitive};
7use num::Bounded;
8
9/// A set of per-channel histograms from an image with 8 bits per channel.
10pub struct ChannelHistogram {
11    /// Per-channel histograms.
12    pub channels: Vec<[u32; 256]>,
13}
14
15/// Returns a vector of per-channel histograms.
16pub fn histogram<P>(image: &Image<P>) -> ChannelHistogram
17where
18    P: Pixel<Subpixel = u8>,
19{
20    let mut hist = vec![[0u32; 256]; P::CHANNEL_COUNT as usize];
21
22    for pix in image.pixels() {
23        for (i, c) in pix.channels().iter().enumerate() {
24            hist[i][*c as usize] += 1;
25        }
26    }
27
28    ChannelHistogram { channels: hist }
29}
30
31/// A set of per-channel cumulative histograms from an image with 8 bits per channel.
32pub struct CumulativeChannelHistogram {
33    /// Per-channel cumulative histograms.
34    pub channels: Vec<[u32; 256]>,
35}
36
37/// Returns per-channel cumulative histograms.
38pub fn cumulative_histogram<P>(image: &Image<P>) -> CumulativeChannelHistogram
39where
40    P: Pixel<Subpixel = u8>,
41{
42    let mut hist = histogram(image);
43
44    for c in 0..hist.channels.len() {
45        for i in 1..hist.channels[c].len() {
46            hist.channels[c][i] += hist.channels[c][i - 1];
47        }
48    }
49
50    CumulativeChannelHistogram {
51        channels: hist.channels,
52    }
53}
54
55/// Returns the `p`th percentile of the pixel intensities in an image.
56///
57/// We define the `p`th percentile intensity to be the least `x` such
58/// that at least `p`% of image pixels have intensity less than or
59/// equal to `x`.
60///
61/// # Panics
62/// If `p > 100`.
63///
64/// # Examples
65/// ```
66/// # extern crate image;
67/// # #[macro_use]
68/// # extern crate imageproc;
69/// # fn main() {
70/// use imageproc::stats::percentile;
71///
72/// let image = gray_image!(
73///     1, 2, 3, 4, 5;
74///     6, 7, 8, 9, 10);
75///
76/// // The 0th percentile is always 0
77/// assert_eq!(percentile(&image, 0), 0);
78///
79/// // Exactly 10% of pixels have intensity <= 1.
80/// assert_eq!(percentile(&image, 10), 1);
81///
82/// // Fewer than 15% of pixels have intensity <=1, so the 15th percentile is 2.
83/// assert_eq!(percentile(&image, 15), 2);
84///
85/// // All pixels have intensity <= 10.
86/// assert_eq!(percentile(&image, 100), 10);
87/// # }
88/// ```
89pub fn percentile(image: &GrayImage, p: u8) -> u8 {
90    assert!(p <= 100, "requested percentile must be <= 100");
91
92    let cum_hist = cumulative_histogram(image).channels[0];
93    let total = cum_hist[255] as u64;
94
95    for i in 0..256 {
96        if 100 * cum_hist[i] as u64 / total >= p as u64 {
97            return i as u8;
98        }
99    }
100
101    unreachable!();
102}
103
104/// Returns the square root of the mean of the squares of differences
105/// between all subpixels in left and right. All channels are considered
106/// equally. If you do not want this (e.g. if using RGBA) then change
107/// image formats first.
108pub fn root_mean_squared_error<I, J, P>(left: &I, right: &J) -> f64
109where
110    I: GenericImageView<Pixel = P>,
111    J: GenericImageView<Pixel = P>,
112    P: Pixel,
113    P::Subpixel: ValueInto<f64>,
114{
115    mean_squared_error(left, right).sqrt()
116}
117
118/// Returns the peak signal to noise ratio for a clean image and its noisy
119/// aproximation. All channels are considered equally. If you do not want this
120/// (e.g. if using RGBA) then change image formats first.
121/// See also [peak signal-to-noise ratio (wikipedia)](https://en.wikipedia.org/wiki/Peak_signal-to-noise_ratio).
122pub fn peak_signal_to_noise_ratio<I, J, P>(original: &I, noisy: &J) -> f64
123where
124    I: GenericImageView<Pixel = P>,
125    J: GenericImageView<Pixel = P>,
126    P: Pixel,
127    P::Subpixel: ValueInto<f64> + Primitive,
128{
129    let max: f64 = cast(<P::Subpixel as Bounded>::max_value());
130    let mse = mean_squared_error(original, noisy);
131    20f64 * max.log(10f64) - 10f64 * mse.log(10f64)
132}
133
134fn mean_squared_error<I, J, P>(left: &I, right: &J) -> f64
135where
136    I: GenericImageView<Pixel = P>,
137    J: GenericImageView<Pixel = P>,
138    P: Pixel,
139    P::Subpixel: ValueInto<f64>,
140{
141    assert_dimensions_match!(left, right);
142    let mut sum_squared_diffs = 0f64;
143    for (p, q) in left.pixels().zip(right.pixels()) {
144        for (c, d) in p.2.channels().iter().zip(q.2.channels().iter()) {
145            let fc: f64 = cast(*c);
146            let fd: f64 = cast(*d);
147            let diff = fc - fd;
148            sum_squared_diffs += diff * diff;
149        }
150    }
151    let count = (left.width() * left.height() * P::CHANNEL_COUNT as u32) as f64;
152    sum_squared_diffs / count
153}
154
155#[cfg(test)]
156mod tests {
157    use super::*;
158    use image::{GrayImage, Luma, Rgb, RgbImage};
159    use test::Bencher;
160
161    #[test]
162    fn test_cumulative_histogram() {
163        let image = gray_image!(1u8, 2u8, 3u8, 2u8, 1u8);
164        let hist = cumulative_histogram(&image).channels[0];
165
166        assert_eq!(hist[0..4], [0, 2, 4, 5]);
167        assert!(hist.iter().skip(4).all(|x| *x == 5));
168    }
169
170    #[test]
171    fn test_cumulative_histogram_rgb() {
172        let image = rgb_image!(
173            [1u8, 10u8, 1u8],
174            [2u8, 20u8, 2u8],
175            [3u8, 30u8, 3u8],
176            [2u8, 20u8, 2u8],
177            [1u8, 10u8, 1u8]
178        );
179
180        let hist = cumulative_histogram(&image);
181        let r = hist.channels[0];
182        let g = hist.channels[1];
183        let b = hist.channels[2];
184
185        assert_eq!(r[0..4], [0, 2, 4, 5]);
186        assert!(r.iter().skip(4).all(|x| *x == 5));
187
188        assert_eq!(g[0..10], [0; 10]);
189        assert_eq!(g[10..20], [2; 10]);
190        assert_eq!(g[20..30], [4; 10]);
191        assert_eq!(g[30], 5);
192        assert!(b.iter().skip(30).all(|x| *x == 5));
193
194        assert_eq!(b[0..4], [0, 2, 4, 5]);
195        assert!(b.iter().skip(4).all(|x| *x == 5));
196    }
197
198    #[test]
199    fn test_histogram() {
200        let image = gray_image!(1u8, 2u8, 3u8, 2u8, 1u8);
201        let hist = histogram(&image).channels[0];
202
203        assert_eq!(hist[0..4], [0, 2, 2, 1]);
204    }
205
206    #[test]
207    fn test_histogram_rgb() {
208        let image = rgb_image!(
209            [1u8, 10u8, 1u8],
210            [2u8, 20u8, 2u8],
211            [3u8, 30u8, 3u8],
212            [2u8, 20u8, 2u8],
213            [1u8, 10u8, 1u8]
214        );
215
216        let hist = histogram(&image);
217        let r = hist.channels[0];
218        let g = hist.channels[1];
219        let b = hist.channels[2];
220
221        assert_eq!(r[0..4], [0, 2, 2, 1]);
222        assert!(r.iter().skip(4).all(|x| *x == 0));
223
224        assert_eq!(g[0..10], [0; 10]);
225        assert_eq!(g[10], 2);
226        assert_eq!(g[11..20], [0; 9]);
227        assert_eq!(g[20], 2);
228        assert_eq!(g[21..30], [0; 9]);
229        assert_eq!(g[30], 1);
230        assert!(b.iter().skip(30).all(|x| *x == 0));
231
232        assert_eq!(b[0..4], [0, 2, 2, 1]);
233        assert!(b.iter().skip(4).all(|x| *x == 0));
234    }
235
236    #[test]
237    fn test_root_mean_squared_error_grayscale() {
238        let left = gray_image!(
239            1, 2, 3;
240            4, 5, 6);
241
242        let right = gray_image!(
243            8, 4, 7;
244            6, 9, 1);
245
246        let rms = root_mean_squared_error(&left, &right);
247        let expected = (114f64 / 6f64).sqrt();
248        assert_eq!(rms, expected);
249    }
250
251    #[test]
252    fn test_root_mean_squared_error_rgb() {
253        let left = rgb_image!([1, 2, 3], [4, 5, 6]);
254        let right = rgb_image!([8, 4, 7], [6, 9, 1]);
255        let rms = root_mean_squared_error(&left, &right);
256        let expected = (114f64 / 6f64).sqrt();
257        assert_eq!(rms, expected);
258    }
259
260    #[test]
261    #[should_panic]
262    fn test_root_mean_squares_rejects_mismatched_dimensions() {
263        let left = gray_image!(1, 2);
264        let right = gray_image!(8; 4);
265        let _ = root_mean_squared_error(&left, &right);
266    }
267
268    fn left_image_rgb(width: u32, height: u32) -> RgbImage {
269        RgbImage::from_fn(width, height, |x, y| Rgb([x as u8, y as u8, (x + y) as u8]))
270    }
271
272    fn right_image_rgb(width: u32, height: u32) -> RgbImage {
273        RgbImage::from_fn(width, height, |x, y| Rgb([(x + y) as u8, x as u8, y as u8]))
274    }
275
276    #[bench]
277    fn bench_root_mean_squared_error_rgb(b: &mut Bencher) {
278        let left = left_image_rgb(50, 50);
279        let right = right_image_rgb(50, 50);
280
281        b.iter(|| {
282            let error = root_mean_squared_error(&left, &right);
283            test::black_box(error);
284        });
285    }
286
287    fn left_image_gray(width: u32, height: u32) -> GrayImage {
288        GrayImage::from_fn(width, height, |x, _| Luma([x as u8]))
289    }
290
291    fn right_image_gray(width: u32, height: u32) -> GrayImage {
292        GrayImage::from_fn(width, height, |_, y| Luma([y as u8]))
293    }
294
295    #[bench]
296    fn bench_root_mean_squared_error_gray(b: &mut Bencher) {
297        let left = left_image_gray(50, 50);
298        let right = right_image_gray(50, 50);
299
300        b.iter(|| {
301            let error = root_mean_squared_error(&left, &right);
302            test::black_box(error);
303        });
304    }
305}
imageproc/stats.rs

imageproc/
stats.rs
