Merge pull request #2302 from torfmaster/feature/fast-blur

Fast Blur

Cargo.toml

@@ -105,3 +105,13 @@ harness = false
 [[bench]]
 name = "copy_from"
 harness = false
+
+[[bench]]
+path = "benches/fast_blur.rs"
+name = "fast_blur"
+harness = false
+
+[[bench]]
+path = "benches/blur.rs"
+name = "blur"
+harness = false

benches/blur.rs

@@ -0,0 +1,13 @@
+use criterion::{criterion_group, criterion_main, Criterion};
+use image::{imageops::blur, ImageBuffer, Rgb};
+
+pub fn bench_fast_blur(c: &mut Criterion) {
+    let src = ImageBuffer::from_pixel(1024, 768, Rgb([255u8, 0, 0]));
+
+    c.bench_function("blur", |b| {
+        b.iter(|| blur(&src, 50.0));
+    });
+}
+
+criterion_group!(benches, bench_fast_blur);
+criterion_main!(benches);

benches/fast_blur.rs

@@ -0,0 +1,13 @@
+use criterion::{criterion_group, criterion_main, Criterion};
+use image::{imageops::fast_blur, ImageBuffer, Rgb};
+
+pub fn bench_fast_blur(c: &mut Criterion) {
+    let src = ImageBuffer::from_pixel(1024, 768, Rgb([255u8, 0, 0]));
+
+    c.bench_function("fast_blur", |b| {
+        b.iter(|| fast_blur(&src, 50.0));
+    });
+}
+
+criterion_group!(benches, bench_fast_blur);
+criterion_main!(benches);

examples/fast_blur/.gitignore

@@ -0,0 +1,2 @@
+mandril_color_blurred.tif
+

examples/fast_blur/main.rs

@@ -0,0 +1,14 @@
+use image::ImageReader;
+
+fn main() {
+    let path = concat!(
+        env!("CARGO_MANIFEST_DIR"),
+        "/tests/images/tiff/testsuite/mandrill.tiff"
+    );
+    let img = ImageReader::open(path).unwrap().decode().unwrap();
+
+    let img2 = img.blur(10.0);
+
+    img2.save("examples/fast_blur/mandril_color_blurred.tif")
+        .unwrap();
+}

src/dynimage.rs

@@ -824,11 +824,21 @@ impl DynamicImage {
 
     /// Performs a Gaussian blur on this image.
     /// `sigma` is a measure of how much to blur by.
+    /// Use [DynamicImage::fast_blur()] for a faster but less
+    /// accurate version.
     #[must_use]
     pub fn blur(&self, sigma: f32) -> DynamicImage {
         dynamic_map!(*self, ref p => imageops::blur(p, sigma))
     }
 
+    /// Performs a fast blur on this image.
+    /// `sigma` is the standard deviation of the
+    /// (approximated) Gaussian
+    #[must_use]
+    pub fn fast_blur(&self, sigma: f32) -> DynamicImage {
+        dynamic_map!(*self, ref p => imageops::fast_blur(p, sigma))
+    }
+
     /// Performs an unsharpen mask on this image.
     /// `sigma` is the amount to blur the image by.
     /// `threshold` is a control of how much to sharpen.

src/imageops/fast_blur.rs

@@ -0,0 +1,152 @@
+use num_traits::clamp;
+
+use crate::{ImageBuffer, Pixel, Primitive};
+
+/// Approximation of Gaussian blur after
+/// Kovesi, P.:  Fast Almost-Gaussian Filtering The Australian Pattern
+/// Recognition Society Conference: DICTA 2010. December 2010. Sydney.
+pub fn fast_blur<P: Pixel>(
+    image_buffer: &ImageBuffer<P, Vec<P::Subpixel>>,
+    sigma: f32,
+) -> ImageBuffer<P, Vec<P::Subpixel>> {
+    let (width, height) = image_buffer.dimensions();
+
+    if width == 0 || height == 0 {
+        return image_buffer.clone();
+    }
+    let mut samples = image_buffer.as_flat_samples().samples.to_vec();
+    let num_passes = 3;
+
+    let boxes = boxes_for_gauss(sigma, num_passes);
+
+    for radius in boxes.iter().take(num_passes) {
+        let horizontally_blurred_transposed = horizontal_fast_blur_half::<P::Subpixel>(
+            &samples,
+            width as usize,
+            height as usize,
+            (*radius - 1) / 2,
+            P::CHANNEL_COUNT as usize,
+        );
+        samples = horizontal_fast_blur_half::<P::Subpixel>(
+            &horizontally_blurred_transposed,
+            height as usize,
+            width as usize,
+            (*radius - 1) / 2,
+            P::CHANNEL_COUNT as usize,
+        );
+    }
+    ImageBuffer::from_raw(width, height, samples).unwrap()
+}
+
+fn boxes_for_gauss(sigma: f32, n: usize) -> Vec<usize> {
+    let w_ideal = f32::sqrt((12.0 * sigma.powi(2) / (n as f32)) + 1.0);
+    let mut w_l = w_ideal.floor();
+    if w_l % 2.0 == 0.0 {
+        w_l -= 1.0
+    };
+    let w_u = w_l + 2.0;
+
+    let m_ideal = 0.25 * (n as f32) * (w_l + 3.0) - 3.0 * sigma.powi(2) * (w_l + 1.0).recip();
+
+    let m = f32::round(m_ideal) as usize;
+
+    (0..n)
+        .map(|i| if i < m { w_l as usize } else { w_u as usize })
+        .collect::<Vec<_>>()
+}
+
+fn channel_idx(channel: usize, idx: usize, channel_num: usize) -> usize {
+    channel_num * idx + channel
+}
+
+fn horizontal_fast_blur_half<P: Primitive>(
+    samples: &[P],
+    width: usize,
+    height: usize,
+    r: usize,
+    channel_num: usize,
+) -> Vec<P> {
+    let channel_size = width * height;
+
+    let mut out_samples = vec![P::from(0).unwrap(); channel_size * channel_num];
+    let mut vals = vec![0.0; channel_num];
+
+    let min_value = P::DEFAULT_MIN_VALUE.to_f32().unwrap();
+    let max_value = P::DEFAULT_MAX_VALUE.to_f32().unwrap();
+
+    for row in 0..height {
+        for (channel, value) in vals.iter_mut().enumerate().take(channel_num) {
+            *value = ((-(r as isize))..(r + 1) as isize)
+                .map(|x| {
+                    extended_f(
+                        samples,
+                        width,
+                        height,
+                        x,
+                        row as isize,
+                        channel,
+                        channel_num,
+                    )
+                    .to_f32()
+                    .unwrap_or(0.0)
+                })
+                .sum()
+        }
+
+        for column in 0..width {
+            for (channel, channel_val) in vals.iter_mut().enumerate() {
+                let val = *channel_val / (2.0 * r as f32 + 1.0);
+                let val = clamp(val, min_value, max_value);
+                let val = P::from(val).unwrap();
+
+                let destination_row = column;
+                let destination_column = row;
+                let destination_sample_index = channel_idx(
+                    channel,
+                    destination_column + destination_row * height,
+                    channel_num,
+                );
+                out_samples[destination_sample_index] = val;
+                *channel_val = *channel_val
+                    - extended_f(
+                        samples,
+                        width,
+                        height,
+                        column as isize - r as isize,
+                        row as isize,
+                        channel,
+                        channel_num,
+                    )
+                    .to_f32()
+                    .unwrap_or(0.0)
+                    + extended_f(
+                        samples,
+                        width,
+                        height,
+                        { column + r + 1 } as isize,
+                        row as isize,
+                        channel,
+                        channel_num,
+                    )
+                    .to_f32()
+                    .unwrap_or(0.0)
+            }
+        }
+    }
+
+    out_samples
+}
+
+fn extended_f<P: Primitive>(
+    samples: &[P],
+    width: usize,
+    height: usize,
+    x: isize,
+    y: isize,
+    channel: usize,
+    channel_num: usize,
+) -> P {
+    let x = clamp(x, 0, width as isize - 1) as usize;
+    let y = clamp(y, 0, height as isize - 1) as usize;
+    samples[channel_idx(channel, y * width + x, channel_num)]
+}

src/imageops/mod.rs

@@ -31,8 +31,11 @@ mod affine;
 // Public only because of Rust bug:
 // https://github.com/rust-lang/rust/issues/18241
 pub mod colorops;
+mod fast_blur;
 mod sample;
 
+pub use fast_blur::fast_blur;
+
 /// Return a mutable view into an image
 /// The coordinates set the position of the top left corner of the crop.
 pub fn crop<I: GenericImageView>(
@@ -353,9 +356,12 @@ where
 #[cfg(test)]
 mod tests {
 
-    use super::{overlay, overlay_bounds_ext};
+    use super::*;
     use crate::color::Rgb;
+    use crate::GrayAlphaImage;
+    use crate::GrayImage;
     use crate::ImageBuffer;
+    use crate::RgbImage;
     use crate::RgbaImage;
 
     #[test]
@@ -478,9 +484,86 @@ mod tests {
     }
 
     #[test]
-    /// Test blur doesn't panick when passed 0.0
+    /// Test blur doesn't panic when passed 0.0
     fn test_blur_zero() {
         let image = RgbaImage::new(50, 50);
-        let _ = super::blur(&image, 0.0);
+        let _ = blur(&image, 0.0);
+    }
+
+    #[test]
+    /// Test fast blur doesn't panic when passed 0.0
+    fn test_fast_blur_zero() {
+        let image = RgbaImage::new(50, 50);
+        let _ = fast_blur(&image, 0.0);
+    }
+
+    #[test]
+    /// Test fast blur doesn't panic when passed negative numbers
+    fn test_fast_blur_negative() {
+        let image = RgbaImage::new(50, 50);
+        let _ = fast_blur(&image, -1.0);
+    }
+
+    #[test]
+    /// Test fast blur doesn't panic when sigma produces boxes larger than the image
+    fn test_fast_large_sigma() {
+        let image = RgbaImage::new(1, 1);
+        let _ = fast_blur(&image, 50.0);
+    }
+
+    #[test]
+    /// Test blur doesn't panic when passed an empty image (any direction)
+    fn test_fast_blur_empty() {
+        let image = RgbaImage::new(0, 0);
+        let _ = fast_blur(&image, 1.0);
+        let image = RgbaImage::new(20, 0);
+        let _ = fast_blur(&image, 1.0);
+        let image = RgbaImage::new(0, 20);
+        let _ = fast_blur(&image, 1.0);
+    }
+
+    #[test]
+    /// Test fast blur works with 3 channels
+    fn test_fast_blur_3_channels() {
+        let image = RgbImage::new(50, 50);
+        let _ = fast_blur(&image, 1.0);
+    }
+
+    #[test]
+    /// Test fast blur works with 2 channels
+    fn test_fast_blur_2_channels() {
+        let image = GrayAlphaImage::new(50, 50);
+        let _ = fast_blur(&image, 1.0);
+    }
+
+    #[test]
+    /// Test fast blur works with 1 channel
+    fn test_fast_blur_1_channels() {
+        let image = GrayImage::new(50, 50);
+        let _ = fast_blur(&image, 1.0);
+    }
+
+    #[test]
+    #[cfg(feature = "tiff")]
+    fn fast_blur_approximates_gaussian_blur_well() {
+        let path = concat!(
+            env!("CARGO_MANIFEST_DIR"),
+            "/tests/images/tiff/testsuite/rgb-3c-16b.tiff"
+        );
+        let image = crate::open(path).unwrap();
+        let image_blurred_gauss = image.blur(50.0).to_rgb8();
+        let image_blurred_gauss_samples = image_blurred_gauss.as_flat_samples();
+        let image_blurred_gauss_bytes = image_blurred_gauss_samples.as_slice();
+        let image_blurred_fast = image.fast_blur(50.0).to_rgb8();
+        let image_blurred_fast_samples = image_blurred_fast.as_flat_samples();
+        let image_blurred_fast_bytes = image_blurred_fast_samples.as_slice();
+
+        let error = image_blurred_gauss_bytes
+            .iter()
+            .zip(image_blurred_fast_bytes.iter())
+            .map(|(a, b)| ((*a as f32 - *b as f32) / (*a as f32)))
+            .sum::<f32>()
+            / (image_blurred_gauss_bytes.len() as f32);
+        assert!(error < 0.05);
     }
 }

src/imageops/sample.rs

@@ -953,6 +953,8 @@ where
 
 /// Performs a Gaussian blur on the supplied image.
 /// ```sigma``` is a measure of how much to blur by.
+/// Use [crate::imageops::fast_blur()] for a faster but less
+/// accurate version.
 pub fn blur<I: GenericImageView>(
     image: &I,
     sigma: f32,