Add filter_map_parallel

src/filter/mod.rs

@@ -9,6 +9,7 @@ mod sharpen;
 pub use self::sharpen::*;
 
 use image::{GenericImage, GenericImageView, GrayImage, Luma, Pixel, Primitive};
+use itertools::Itertools;
 
 use crate::definitions::{Clamp, Image};
 use crate::integral_image::{column_running_sum, row_running_sum};
@@ -83,52 +84,111 @@ pub fn box_filter(image: &GrayImage, x_radius: u32, y_radius: u32) -> Image<Luma
     out
 }
 
+/// Calculates the new pixel value for a particular pixel and kernel.
+fn filter_pixel<P, K, F, Q>(x: u32, y: u32, kernel: Kernel<K>, f: F, image: &Image<P>) -> Q
+where
+    P: Pixel,
+    Q: Pixel,
+    F: Fn(K) -> Q::Subpixel,
+    K: num::Num + Copy + From<P::Subpixel>,
+{
+    let (width, height) = (image.width() as i64, image.height() as i64);
+    let (k_width, k_height) = (kernel.width as i64, kernel.height as i64);
+    let (x, y) = (i64::from(x), i64::from(y));
+
+    let weighted_pixels = (0..kernel.height as i64)
+        .cartesian_product(0..kernel.width as i64)
+        .map(|(k_y, k_x)| {
+            let kernel_weight = *kernel.at(k_x as u32, k_y as u32);
+
+            let window_y = (y + k_y - k_height / 2).clamp(0, height - 1);
+            let window_x = (x + k_x - k_width / 2).clamp(0, width - 1);
+
+            debug_assert!(image.in_bounds(window_x as u32, window_y as u32));
+
+            // Safety: we clamped `window_x` and `window_y` to be in bounds.
+            let window_pixel = unsafe { image.unsafe_get_pixel(window_x as u32, window_y as u32) };
+
+            //optimisation: remove allocation when `Pixel::map` allows mapping to a different
+            //type
+            #[allow(clippy::unnecessary_to_owned)]
+            let weighted_pixel = window_pixel
+                .channels()
+                .to_vec()
+                .into_iter()
+                .map(move |c| kernel_weight * K::from(c));
+
+            weighted_pixel
+        });
+
+    let final_channel_sum = weighted_pixels.fold(
+        //optimisation: do this without allocation when `Pixel` gains a method of constant initialization
+        vec![K::zero(); Q::CHANNEL_COUNT as usize],
+        |mut accumulator, weighted_pixel| {
+            for (i, weighted_subpixel) in weighted_pixel.enumerate() {
+                accumulator[i] = accumulator[i] + weighted_subpixel;
+            }
+
+            accumulator
+        },
+    );
+
+    *Q::from_slice(&final_channel_sum.into_iter().map(f).collect_vec())
+}
+
 /// 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.
-pub fn filter<P, K, F, Q>(image: &Image<P>, kernel: Kernel<K>, mut f: F) -> Image<Q>
+pub fn filter<P, K, F, Q>(image: &Image<P>, kernel: Kernel<K>, f: F) -> Image<Q>
 where
     P: Pixel,
-    <P as Pixel>::Subpixel: Into<K>,
     Q: Pixel,
-    F: FnMut(&mut Q::Subpixel, K),
-    K: num::Num + Copy,
+    F: Fn(K) -> Q::Subpixel,
+    K: num::Num + Copy + From<P::Subpixel>,
 {
-    let (width, height) = image.dimensions();
-    let mut out = Image::<Q>::new(width, height);
-    let num_channels = P::CHANNEL_COUNT as usize;
-    let zero = K::zero();
-    let mut acc = vec![zero; num_channels];
+    //TODO refactor this to use the `crate::maps` functions once that lands
 
-    let (k_width, k_height) = (kernel.width, kernel.height);
+    let (width, height) = image.dimensions();
 
-    let (width, height, k_width, k_height) =
-        (width as i64, height as i64, k_width as i64, k_height as i64);
+    let mut out = Image::<Q>::new(width, height);
 
     for y in 0..height {
         for x in 0..width {
-            for k_y in 0..k_height {
-                let y_p = min(height - 1, max(0, y + k_y - k_height / 2)) as u32;
-                for k_x in 0..k_width {
-                    let x_p = min(width - 1, max(0, x + k_x - k_width / 2)) as u32;
-                    accumulate(
-                        &mut acc,
-                        unsafe { &image.unsafe_get_pixel(x_p, y_p) },
-                        *kernel.at(k_x as u32, k_y as u32),
-                    );
-                }
-            }
-            let out_channels = out.get_pixel_mut(x as u32, y as u32).channels_mut();
-            for (a, c) in acc.iter_mut().zip(out_channels.iter_mut()) {
-                f(c, *a);
-                *a = zero;
-            }
+            out.put_pixel(x, y, filter_pixel(x, y, kernel, &f, image));
         }
     }
 
     out
 }
+#[cfg(feature = "rayon")]
+#[doc = generate_parallel_doc_comment!("filter")]
+pub fn filter_parallel<P, K, F, Q>(image: &Image<P>, kernel: Kernel<K>, f: F) -> Image<Q>
+where
+    P: Pixel + Sync,
+    Q: Pixel + Send + Sync,
+    P::Subpixel: Sync,
+    Q::Subpixel: Send + Sync,
+    F: Fn(K) -> Q::Subpixel + Send + Sync,
+    K: Num + Copy + From<P::Subpixel> + Sync,
+{
+    //TODO refactor this to use the `crate::maps` functions once that lands
+
+    use rayon::iter::IndexedParallelIterator;
+    use rayon::iter::ParallelIterator;
+
+    let (width, height) = image.dimensions();
+
+    let mut out: Image<Q> = Image::new(width, height);
+
+    image
+        .par_enumerate_pixels()
+        .zip_eq(out.par_pixels_mut())
+        .for_each(move |((x, y, _), output_pixel)| {
+            *output_pixel = filter_pixel(x, y, kernel, &f, image);
+        });
+
+    out
+}
 
-/// The gaussian probability density function with a mean of zero.
 #[inline]
 fn gaussian_pdf(x: f32, sigma: f32) -> f32 {
     // Equation derived from <https://en.wikipedia.org/wiki/Gaussian_function>
@@ -203,14 +263,33 @@ 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>,
     S: Clamp<K> + Primitive,
     P: WithChannel<S>,
-    K: Num + Copy,
+    K: Num + Copy + From<<P as image::Pixel>::Subpixel>,
 {
-    filter(image, kernel, |channel, acc| *channel = S::clamp(acc))
+    filter(image, kernel, S::clamp)
+}
+#[cfg(feature = "rayon")]
+#[doc = generate_parallel_doc_comment!("filter_clamped")]
+pub fn filter_clamped_parallel<P, K, S>(
+    image: &Image<P>,
+    kernel: Kernel<K>,
+) -> Image<ChannelMap<P, S>>
+where
+    P: Sync,
+    P::Subpixel: Send + Sync,
+    <P as WithChannel<S>>::Pixel: Send + Sync,
+    S: Clamp<K> + Primitive + Send + Sync,
+    P: WithChannel<S>,
+    K: Num + Copy + Send + Sync + From<P::Subpixel>,
+{
+    filter_parallel(image, kernel, |x| S::clamp(x))
 }
 
 /// Returns horizontal correlations between an image and a 1d kernel.
@@ -441,6 +520,12 @@ where
 pub fn laplacian_filter(image: &GrayImage) -> Image<Luma<i16>> {
     filter_clamped(image, kernel::FOUR_LAPLACIAN_3X3)
 }
+#[must_use = "the function does not modify the original image"]
+#[cfg(feature = "rayon")]
+#[doc = generate_parallel_doc_comment!("laplacian_filter")]
+pub fn laplacian_filter_parallel(image: &GrayImage) -> Image<Luma<i16>> {
+    filter_clamped_parallel(image, kernel::FOUR_LAPLACIAN_3X3)
+}
 
 #[cfg(test)]
 mod tests {
@@ -578,7 +663,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 {
@@ -691,6 +776,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() {
@@ -706,7 +816,7 @@ mod tests {
 
         let k = &[1u8, 2u8];
         let kernel = Kernel::new(k, 2, 1);
-        let filtered = filter(&image, kernel, |c, a| *c = a);
+        let filtered = filter(&image, kernel, |x| x);
 
         let expected = gray_image!(
              9,  7;
@@ -721,7 +831,7 @@ mod tests {
 
         let k = &[2u8];
         let kernel = Kernel::new(k, 1, 1);
-        let filtered = filter(&image, kernel, |c, a| *c = a);
+        let filtered = filter(&image, kernel, |x| x);
 
         let expected = gray_image!();
         assert_pixels_eq!(filtered, expected);
@@ -740,8 +850,7 @@ mod tests {
             0.1, 0.2, 0.1
         ];
         let kernel = Kernel::new(k, 3, 3);
-        let filtered: Image<Luma<u8>> =
-            filter(&image, kernel, |c, a| *c = <u8 as Clamp<f32>>::clamp(a));
+        let filtered: Image<Luma<u8>> = filter(&image, kernel, <u8 as Clamp<f32>>::clamp);
 
         let expected = gray_image!(
             11,  7;
@@ -819,8 +928,8 @@ mod proptests {
             ker in arbitrary_image::<Luma<f32>>(1..20, 1..20),
         ) {
             let kernel = Kernel::new(&ker, ker.width(), ker.height());
-            let out: Image<Luma<f32>> = filter(&img, kernel, |dst, src| {
-                *dst = src;
+            let out: Image<Luma<f32>> = filter(&img, kernel, |x| {
+                x
             });
             assert_eq!(out.dimensions(), img.dimensions());
         }
@@ -919,6 +1028,22 @@ 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: Image<Luma<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},
@@ -12,6 +14,13 @@ pub fn sharpen3x3(image: &GrayImage) -> GrayImage {
     let identity_minus_laplacian = Kernel::new(&[0, -1, 0, -1, 5, -1, 0, -1, 0], 3, 3);
     filter_clamped(image, identity_minus_laplacian)
 }
+#[must_use = "the function does not modify the original image"]
+#[cfg(feature = "rayon")]
+#[doc = generate_parallel_doc_comment!("sharpen3x3")]
+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.
 ///

src/geometric_transformations.rs

@@ -432,7 +432,7 @@ pub fn warp_into<P>(
 /// # Examples
 /// Applying a wave pattern.
 /// ```
-/// use image::{ImageBuffer, Luma};
+/// use image::Luma;
 /// use imageproc::utils::gray_bench_image;
 /// use imageproc::geometric_transformations::*;
 ///