use crate::definitions::Image;
use crate::map::{ChannelMap, WithChannel};
use image::{GenericImageView, GrayImage, Luma, Pixel, Primitive, Rgb, Rgba};
use std::ops::AddAssign;
pub fn integral_image<P, T>(image: &Image<P>) -> Image<ChannelMap<P, T>>
where
P: Pixel<Subpixel = u8> + WithChannel<T> + 'static,
T: From<u8> + Primitive + AddAssign + 'static,
{
integral_image_impl(image, false)
}
pub fn integral_squared_image<P, T>(image: &Image<P>) -> Image<ChannelMap<P, T>>
where
P: Pixel<Subpixel = u8> + WithChannel<T> + 'static,
T: From<u8> + Primitive + AddAssign + 'static,
{
integral_image_impl(image, true)
}
fn integral_image_impl<P, T>(image: &Image<P>, square: bool) -> Image<ChannelMap<P, T>>
where
P: Pixel<Subpixel = u8> + WithChannel<T> + 'static,
T: From<u8> + Primitive + AddAssign + 'static,
{
let (in_width, in_height) = image.dimensions();
let out_width = in_width + 1;
let out_height = in_height + 1;
let mut out = Image::<ChannelMap<P, T>>::new(out_width, out_height);
if in_width == 0 || in_height == 0 {
return out;
}
for y in 0..in_height {
let mut sum = vec![T::zero(); P::CHANNEL_COUNT as usize];
for x in 0..in_width {
let input = unsafe { image.unsafe_get_pixel(x, y) };
for (s, c) in sum.iter_mut().zip(input.channels()) {
let pix: T = (*c).into();
*s += if square { pix * pix } else { pix };
}
let above = unsafe { out.unsafe_get_pixel(x + 1, y) };
let current = out.get_pixel_mut(x + 1, y + 1);
for c in 0..P::CHANNEL_COUNT {
current.channels_mut()[c as usize] = above.channels()[c as usize] + sum[c as usize];
}
}
}
out
}
pub trait ArrayData {
type DataType;
fn data(&self) -> Self::DataType;
fn add(lhs: Self::DataType, other: Self::DataType) -> Self::DataType;
fn sub(lhs: Self::DataType, other: Self::DataType) -> Self::DataType;
}
impl<T: Primitive + 'static> ArrayData for Luma<T> {
type DataType = [T; 1];
fn data(&self) -> Self::DataType {
[self.channels()[0]]
}
fn add(lhs: Self::DataType, rhs: Self::DataType) -> Self::DataType {
[lhs[0] + rhs[0]]
}
fn sub(lhs: Self::DataType, rhs: Self::DataType) -> Self::DataType {
[lhs[0] - rhs[0]]
}
}
impl<T: Primitive + 'static> ArrayData for Rgb<T> {
type DataType = [T; 3];
fn data(&self) -> Self::DataType {
[self.channels()[0], self.channels()[1], self.channels()[2]]
}
fn add(lhs: Self::DataType, rhs: Self::DataType) -> Self::DataType {
[lhs[0] + rhs[0], lhs[1] + rhs[1], lhs[2] + rhs[2]]
}
fn sub(lhs: Self::DataType, rhs: Self::DataType) -> Self::DataType {
[lhs[0] - rhs[0], lhs[1] - rhs[1], lhs[2] - rhs[2]]
}
}
impl<T: Primitive + 'static> ArrayData for Rgba<T> {
type DataType = [T; 4];
fn data(&self) -> Self::DataType {
[
self.channels()[0],
self.channels()[1],
self.channels()[2],
self.channels()[3],
]
}
fn add(lhs: Self::DataType, rhs: Self::DataType) -> Self::DataType {
[
lhs[0] + rhs[0],
lhs[1] + rhs[1],
lhs[2] + rhs[2],
lhs[3] + rhs[3],
]
}
fn sub(lhs: Self::DataType, rhs: Self::DataType) -> Self::DataType {
[
lhs[0] - rhs[0],
lhs[1] - rhs[1],
lhs[2] - rhs[2],
lhs[3] - rhs[3],
]
}
}
pub fn sum_image_pixels<P>(
integral_image: &Image<P>,
left: u32,
top: u32,
right: u32,
bottom: u32,
) -> P::DataType
where
P: Pixel + ArrayData + Copy + 'static,
{
let (a, b, c, d) = (
integral_image.get_pixel(right + 1, bottom + 1).data(),
integral_image.get_pixel(left, top).data(),
integral_image.get_pixel(right + 1, top).data(),
integral_image.get_pixel(left, bottom + 1).data(),
);
P::sub(P::sub(P::add(a, b), c), d)
}
pub fn variance(
integral_image: &Image<Luma<u32>>,
integral_squared_image: &Image<Luma<u32>>,
left: u32,
top: u32,
right: u32,
bottom: u32,
) -> f64 {
let n = (right - left + 1) as f64 * (bottom - top + 1) as f64;
let sum_sq = sum_image_pixels(integral_squared_image, left, top, right, bottom)[0];
let sum = sum_image_pixels(integral_image, left, top, right, bottom)[0];
(sum_sq as f64 - (sum as f64).powi(2) / n) / n
}
pub fn row_running_sum(image: &GrayImage, row: u32, buffer: &mut [u32], padding: u32) {
let (width, height) = image.dimensions();
let (width, padding) = (width as usize, padding as usize);
assert!(
buffer.len() >= width + 2 * padding,
"Buffer length {} is less than {} + 2 * {}",
buffer.len(),
width,
padding
);
assert!(row < height, "row out of bounds: {} >= {}", row, height);
assert!(width > 0, "image is empty");
let row_data = &(**image)[width * row as usize..][..width];
let first = row_data[0] as u32;
let last = row_data[width - 1] as u32;
let mut sum = 0;
for b in &mut buffer[..padding] {
sum += first;
*b = sum;
}
for (b, p) in buffer[padding..].iter_mut().zip(row_data) {
sum += *p as u32;
*b = sum;
}
for b in &mut buffer[padding + width..] {
sum += last;
*b = sum;
}
}
pub fn column_running_sum(image: &GrayImage, column: u32, buffer: &mut [u32], padding: u32) {
let (width, height) = image.dimensions();
assert!(
buffer.len() >= height as usize + 2 * padding as usize,
"Buffer length {} is less than {} + 2 * {}",
buffer.len(),
height,
padding
);
assert!(
column < width,
"column out of bounds: {} >= {}",
column,
width
);
assert!(
height > 0,
"image is empty"
);
let first = image.get_pixel(column, 0)[0] as u32;
let last = image.get_pixel(column, height - 1)[0] as u32;
let mut sum = 0;
for b in &mut buffer[..padding as usize] {
sum += first;
*b = sum;
}
unsafe {
for y in 0..height {
sum += image.unsafe_get_pixel(column, y)[0] as u32;
*buffer.get_unchecked_mut(y as usize + padding as usize) = sum;
}
}
for b in &mut buffer[padding as usize + height as usize..] {
sum += last;
*b = sum;
}
}
#[cfg(test)]
mod tests {
use super::*;
use crate::definitions::Image;
use crate::property_testing::GrayTestImage;
use crate::utils::{gray_bench_image, pixel_diff_summary, rgb_bench_image};
use ::test;
use image::{GenericImage, ImageBuffer, Luma};
use quickcheck::{quickcheck, TestResult};
#[test]
fn test_integral_image_gray() {
let image = gray_image!(
1, 2, 3;
4, 5, 6);
let expected = gray_image!(type: u32,
0, 0, 0, 0;
0, 1, 3, 6;
0, 5, 12, 21);
assert_pixels_eq!(integral_image::<_, u32>(&image), expected);
}
#[test]
fn test_integral_image_rgb() {
let image = rgb_image!(
[1, 11, 21], [2, 12, 22], [3, 13, 23];
[4, 14, 24], [5, 15, 25], [6, 16, 26]);
let expected = rgb_image!(type: u32,
[0, 0, 0], [0, 0, 0], [ 0, 0, 0], [ 0, 0, 0];
[0, 0, 0], [1, 11, 21], [ 3, 23, 43], [ 6, 36, 66];
[0, 0, 0], [5, 25, 45], [12, 52, 92], [21, 81, 141]);
assert_pixels_eq!(integral_image::<_, u32>(&image), expected);
}
#[test]
fn test_sum_image_pixels() {
let image = gray_image!(
1, 2;
3, 4);
let integral = integral_image::<_, u32>(&image);
assert_eq!(sum_image_pixels(&integral, 0, 0, 0, 0)[0], 1);
assert_eq!(sum_image_pixels(&integral, 0, 0, 1, 0)[0], 3);
assert_eq!(sum_image_pixels(&integral, 0, 0, 0, 1)[0], 4);
assert_eq!(sum_image_pixels(&integral, 0, 0, 1, 1)[0], 10);
assert_eq!(sum_image_pixels(&integral, 1, 0, 1, 0)[0], 2);
assert_eq!(sum_image_pixels(&integral, 1, 0, 1, 1)[0], 6);
assert_eq!(sum_image_pixels(&integral, 0, 1, 0, 1)[0], 3);
assert_eq!(sum_image_pixels(&integral, 0, 1, 1, 1)[0], 7);
assert_eq!(sum_image_pixels(&integral, 1, 1, 1, 1)[0], 4);
}
#[test]
fn test_sum_image_pixels_rgb() {
let image = rgb_image!(
[1, 2, 3], [ 4, 5, 6];
[7, 8, 9], [10, 11, 12]);
let integral = integral_image::<_, u32>(&image);
assert_eq!(sum_image_pixels(&integral, 0, 0, 0, 0), [1, 2, 3]);
assert_eq!(sum_image_pixels(&integral, 0, 0, 1, 0), [5, 7, 9]);
assert_eq!(sum_image_pixels(&integral, 0, 0, 0, 1), [8, 10, 12]);
assert_eq!(sum_image_pixels(&integral, 0, 0, 1, 1), [22, 26, 30]);
assert_eq!(sum_image_pixels(&integral, 1, 0, 1, 0), [4, 5, 6]);
assert_eq!(sum_image_pixels(&integral, 1, 0, 1, 1), [14, 16, 18]);
assert_eq!(sum_image_pixels(&integral, 0, 1, 0, 1), [7, 8, 9]);
assert_eq!(sum_image_pixels(&integral, 0, 1, 1, 1), [17, 19, 21]);
assert_eq!(sum_image_pixels(&integral, 1, 1, 1, 1), [10, 11, 12]);
}
#[bench]
fn bench_integral_image_gray(b: &mut test::Bencher) {
let image = gray_bench_image(500, 500);
b.iter(|| {
let integral = integral_image::<_, u32>(&image);
test::black_box(integral);
});
}
#[bench]
fn bench_integral_image_rgb(b: &mut test::Bencher) {
let image = rgb_bench_image(500, 500);
b.iter(|| {
let integral = integral_image::<_, u32>(&image);
test::black_box(integral);
});
}
fn integral_image_ref<I>(image: &I) -> Image<Luma<u32>>
where
I: GenericImage<Pixel = Luma<u8>>,
{
let (in_width, in_height) = image.dimensions();
let (out_width, out_height) = (in_width + 1, in_height + 1);
let mut out = ImageBuffer::from_pixel(out_width, out_height, Luma([0u32]));
for y in 1..out_height {
for x in 0..out_width {
let mut sum = 0u32;
for iy in 0..y {
for ix in 0..x {
sum += image.get_pixel(ix, iy)[0] as u32;
}
}
out.put_pixel(x, y, Luma([sum]));
}
}
out
}
#[test]
fn test_integral_image_matches_reference_implementation() {
fn prop(image: GrayTestImage) -> TestResult {
let expected = integral_image_ref(&image.0);
let actual = integral_image(&image.0);
match pixel_diff_summary(&actual, &expected) {
None => TestResult::passed(),
Some(err) => TestResult::error(err),
}
}
quickcheck(prop as fn(GrayTestImage) -> TestResult);
}
#[bench]
fn bench_row_running_sum(b: &mut test::Bencher) {
let image = gray_bench_image(1000, 1);
let mut buffer = [0; 1010];
b.iter(|| {
row_running_sum(&image, 0, &mut buffer, 5);
});
}
#[bench]
fn bench_column_running_sum(b: &mut test::Bencher) {
let image = gray_bench_image(100, 1000);
let mut buffer = [0; 1010];
b.iter(|| {
column_running_sum(&image, 0, &mut buffer, 5);
});
}
}