Add Rayon implementation filter_clamped_parallel for Kernel filter function (formerly filter3x3, ...)

examples/gradients.rs

@@ -5,8 +5,11 @@
 //! `cargo run --example gradients ./examples/empire-state-building.jpg ./tmp`
 
 use image::{open, GrayImage};
+#[cfg(not(feature = "rayon"))]
+use imageproc::filter::filter_clamped;
+#[cfg(feature = "rayon")]
+use imageproc::filter::filter_clamped_parallel as filter_clamped;
 use imageproc::{
-    filter::filter_clamped,
     kernel::{self, Kernel},
     map::map_subpixels,
 };

src/filter/mod.rs

@@ -17,6 +17,9 @@ use num::Num;
 use std::cmp::{max, min};
 use std::f32;
 
+#[cfg(feature = "rayon")]
+use rayon::prelude::*;
+
 /// Denoise 8-bit grayscale image using bilateral filtering.
 ///
 /// # Arguments
@@ -248,6 +251,68 @@ where
     out
 }
 
+#[cfg(feature = "rayon")]
+#[cfg_attr(docsrs, doc(cfg(feature = "rayon")))]
+/// Returns 2d correlation of an image. Intermediate calculations are performed
+/// at type K, and the results converted to pixel Q via f. Pads by continuity.
+/// This version uses rayon to parallelize the computation.
+pub fn filter_parallel<P, K, F, Q>(image: &Image<P>, kernel: Kernel<K>, f: F) -> Image<Q>
+where
+    P: Pixel + Sync,
+    <P as Pixel>::Subpixel: Into<K> + Send + Sync,
+    Q: Pixel + Send + Sync,
+    <Q as Pixel>::Subpixel: Send,
+    K: num::Num + Copy + Send + Sync,
+    F: Fn(&mut Q::Subpixel, K) + Send + Sync,
+{
+    let (width, height) = image.dimensions();
+    let num_channels = P::CHANNEL_COUNT as usize;
+    let zero = K::zero();
+    let (k_width, k_height) = (kernel.width, kernel.height);
+    let (width, height) = (width as i64, height as i64);
+
+    let out = Image::<Q>::new(width as u32, height as u32);
+    if width == 0 || height == 0 {
+        return out;
+    }
+
+    let mut out_data = out.into_raw();
+
+    out_data
+        .par_chunks_mut(width as usize * num_channels)
+        .enumerate()
+        .for_each(|(y, out_row)| {
+            let mut acc = vec![zero; num_channels];
+
+            for x in 0..width {
+                for k_y in 0..k_height {
+                    let y_p = min(
+                        height - 1,
+                        max(0, y as i64 + k_y as i64 - k_height as i64 / 2),
+                    ) as u32;
+                    for k_x in 0..k_width {
+                        let x_p =
+                            min(width - 1, max(0, x + k_x as i64 - k_width as i64 / 2)) as u32;
+                        debug_assert!(image.in_bounds(x_p, y_p));
+                        let px = unsafe { &image.unsafe_get_pixel(x_p, y_p) };
+                        let coefficient = *kernel.at(k_x, k_y);
+                        accumulate(&mut acc, px, coefficient);
+                    }
+                }
+
+                let pixel_slice = &mut out_row[(x as usize * num_channels)..][..num_channels];
+
+                for (a, c) in acc.iter_mut().zip(pixel_slice.iter_mut()) {
+                    f(c, *a);
+                    *a = zero;
+                }
+            }
+        });
+
+    Image::<Q>::from_raw(width as u32, height as u32, out_data)
+        .expect("failed to create output from raw data")
+}
+
 #[inline]
 fn gaussian(x: f32, r: f32) -> f32 {
     ((2.0 * f32::consts::PI).sqrt() * r).recip() * (-x.powi(2) / (2.0 * r.powi(2))).exp()
@@ -321,6 +386,9 @@ where
 
 /// Returns 2d correlation of an image with a row-major kernel. Intermediate calculations are
 /// performed at type K, and the results clamped to subpixel type S. Pads by continuity.
+///
+/// A parallelized version of this function exists with [`filter_clamped_parallel`] when
+/// the crate `rayon` feature is enabled.
 pub fn filter_clamped<P, K, S>(image: &Image<P>, kernel: Kernel<K>) -> Image<ChannelMap<P, S>>
 where
     P::Subpixel: Into<K>,
@@ -331,6 +399,27 @@ where
     filter(image, kernel, |channel, acc| *channel = S::clamp(acc))
 }
 
+/// Returns 2d correlation of an image with a 3x3 row-major kernel. Intermediate calculations are
+/// performed at type K, and the results clamped to subpixel type S. Pads by continuity.
+/// This version uses rayon to parallelize the computation.
+#[must_use = "the function does not modify the original image"]
+#[cfg(feature = "rayon")]
+#[cfg_attr(docsrs, doc(cfg(feature = "rayon")))]
+pub fn filter_clamped_parallel<P, K, S>(
+    image: &Image<P>,
+    kernel: Kernel<K>,
+) -> Image<ChannelMap<P, S>>
+where
+    P: Sync,
+    P::Subpixel: Into<K> + Send + Sync,
+    <P as WithChannel<S>>::Pixel: Send + Sync,
+    S: Clamp<K> + Primitive + Send,
+    P: WithChannel<S>,
+    K: Num + Copy + Send + Sync,
+{
+    filter_parallel(image, kernel, |channel, acc| *channel = S::clamp(acc))
+}
+
 /// Returns horizontal correlations between an image and a 1d kernel.
 /// Pads by continuity. Intermediate calculations are performed at
 /// type K.
@@ -560,6 +649,22 @@ pub fn laplacian_filter(image: &GrayImage) -> Image<Luma<i16>> {
     filter_clamped(image, kernel::FOUR_LAPLACIAN_3X3)
 }
 
+/// Calculates the Laplacian of an image.
+/// This version uses rayon to parallelize the computation.
+///
+/// The Laplacian is computed by filtering the image using the following 3x3 kernel:
+/// ```notrust
+/// 0, 1, 0,
+/// 1, -4, 1,
+/// 0, 1, 0
+/// ```
+#[must_use = "the function does not modify the original image"]
+#[cfg(feature = "rayon")]
+#[cfg_attr(docsrs, doc(cfg(feature = "rayon")))]
+pub fn laplacian_filter_parallel(image: &GrayImage) -> Image<Luma<i16>> {
+    filter_clamped_parallel(image, kernel::FOUR_LAPLACIAN_3X3)
+}
+
 #[cfg(test)]
 mod tests {
     use super::*;
@@ -710,7 +815,7 @@ mod tests {
             #[test]
             fn $test_name() {
                 // I think the interesting edge cases here are determined entirely
-                // by the relative sizes of the kernel and the image side length, so
+                // by the relative sizes of the kernel and the image side length, soz
                 // I'm just enumerating over small values instead of generating random
                 // examples via quickcheck.
                 for height in 0..5 {
@@ -823,6 +928,31 @@ mod tests {
         assert_pixels_eq!(filtered, expected);
     }
 
+    #[test]
+    #[cfg(feature = "rayon")]
+    fn test_filter_clamped_parallel_with_results_outside_input_channel_range() {
+        #[rustfmt::skip]
+            let kernel: Kernel<i32> = Kernel::new(&[
+            -1, 0, 1,
+            -2, 0, 2,
+            -1, 0, 1
+        ], 3, 3);
+
+        let image = gray_image!(
+            3, 2, 1;
+            6, 5, 4;
+            9, 8, 7);
+
+        let expected = gray_image!(type: i16,
+            -4, -8, -4;
+            -4, -8, -4;
+            -4, -8, -4
+        );
+
+        let filtered = filter_clamped_parallel(&image, kernel);
+        assert_pixels_eq!(filtered, expected);
+    }
+
     #[test]
     #[should_panic]
     fn test_kernel_must_be_nonempty() {
@@ -1072,6 +1202,23 @@ mod benches {
         });
     }
 
+    #[bench]
+    fn bench_filter_clamped_parallel_i32_filter(b: &mut Bencher) {
+        let image = gray_bench_image(500, 500);
+        #[rustfmt::skip]
+        let kernel: Kernel<i32> = Kernel::new(&[
+            -1, 0, 1,
+            -2, 0, 2,
+            -1, 0, 1
+        ], 3, 3);
+
+        b.iter(|| {
+            let filtered: ImageBuffer<Luma<i16>, Vec<i16>> =
+                filter_clamped_parallel::<_, _, i16>(&image, kernel);
+            black_box(filtered);
+        });
+    }
+
     /// Baseline implementation of Gaussian blur is that provided by image::imageops.
     /// We can also use this to validate correctness of any implementations we add here.
     fn gaussian_baseline_rgb<I>(image: &I, stdev: f32) -> Image<Rgb<u8>>

src/filter/sharpen.rs

@@ -1,3 +1,5 @@
+#[cfg(feature = "rayon")]
+use super::filter_clamped_parallel;
 use super::{filter_clamped, gaussian_blur_f32};
 use crate::{
     definitions::{Clamp, Image},
@@ -13,6 +15,16 @@ pub fn sharpen3x3(image: &GrayImage) -> GrayImage {
     filter_clamped(image, identity_minus_laplacian)
 }
 
+/// Sharpens a grayscale image by applying a 3x3 approximation to the Laplacian.
+/// This version uses rayon to parallelize the computation.
+#[must_use = "the function does not modify the original image"]
+#[cfg(feature = "rayon")]
+#[cfg_attr(docsrs, doc(cfg(feature = "rayon")))]
+pub fn sharpen3x3_parallel(image: &GrayImage) -> GrayImage {
+    let identity_minus_laplacian = Kernel::new(&[0, -1, 0, -1, 5, -1, 0, -1, 0], 3, 3);
+    filter_clamped_parallel(image, identity_minus_laplacian)
+}
+
 /// Sharpens a grayscale image using a Gaussian as a low-pass filter.
 ///
 /// * `sigma` is the standard deviation of the Gaussian filter used.