Implement Box blur fast filter that could approximate gaussian filter (#223)

light-le · edgarriba · web-flow · commit 0a357b5b6411 · 2025-02-06T10:29:11.000+01:00
* implement box_blur_fast_kernels_1d with test

* implement fast_horizontal_filter with test

* implement box_blur_fast ops with tests

* Apply docstring suggestions from code review

Co-authored-by: Edgar Riba &lt;edgar.riba@gmail.com&gt;

* add more tests for box_blur_fast_kernels_1d

* some fixes from suggestions

* optimize calculation within loop and fmt fix

* add bench for fast_box_blur

* fmt fixes

* simplify header doc

---------

Co-authored-by: Edgar Riba &lt;edgar.riba@gmail.com&gt;
diff --git a/crates/kornia-imgproc/benches/bench_filters.rs b/crates/kornia-imgproc/benches/bench_filters.rs
@@ -1,7 +1,7 @@
 use criterion::{black_box, criterion_group, criterion_main, BenchmarkId, Criterion};
 
 use kornia_image::Image;
-use kornia_imgproc::filter::gaussian_blur;
+use kornia_imgproc::filter::{box_blur_fast, gaussian_blur};
 
 use image::RgbImage;
 use imageproc::filter::gaussian_blur_f32;
@@ -51,6 +51,15 @@ fn bench_filters(c: &mut Criterion) {
                     b.iter(|| black_box(gaussian_blur_f32(&rgb_image, sigma)))
                 },
             );
+
+            group.bench_with_input(
+                BenchmarkId::new("box_blur_fast", &parameter_string),
+                &(&image, &output),
+                |b, i| {
+                    let (src, mut dst) = (i.0, i.1.clone());
+                    b.iter(|| black_box(box_blur_fast(src, &mut dst, (1.5, 1.5))))
+                },
+            );
         }
     }
 
diff --git a/crates/kornia-imgproc/src/filter/kernels.rs b/crates/kornia-imgproc/src/filter/kernels.rs
@@ -61,6 +61,37 @@ pub fn sobel_kernel_1d(kernel_size: usize) -> (Vec<f32>, Vec<f32>) {
     (kernel_x, kernel_y)
 }
 
+/// Create list of optimized box blur kernels based on gaussian sigma
+///
+/// https://www.peterkovesi.com/papers/FastGaussianSmoothing.pdf
+/// # Arguments
+///
+/// * `sigma` = The sigma of the gaussian kernel
+/// * `kernels` = The number of times the box blur kernels would be applied, ideally from 3-5
+///
+/// # Returns
+///
+/// A kernels-sized vector of the kernels.
+pub fn box_blur_fast_kernels_1d(sigma: f32, kernels: u8) -> Vec<usize> {
+    let n = kernels as f32;
+    let ideal_size = (12.0 * sigma * sigma / n + 1.0).sqrt();
+    let mut size_l = ideal_size.floor();
+    size_l -= if size_l % 2.0 == 0.0 { 1.0 } else { 0.0 };
+    let size_u = size_l + 2.0;
+
+    let ideal_m = (12.0 * sigma * sigma - n * size_l * size_l - 4.0 * n * size_l - 3.0 * n)
+        / (-4.0 * size_l - 4.0);
+    let mut boxes = Vec::new();
+    for i in 0..kernels {
+        if i < ideal_m.round() as u8 {
+            boxes.push(size_l as usize);
+        } else {
+            boxes.push(size_u as usize);
+        }
+    }
+    boxes
+}
+
 #[cfg(test)]
 mod tests {
     use super::*;
@@ -92,4 +123,14 @@ mod tests {
             assert_eq!(k, expected[i]);
         }
     }
+
+    #[test]
+    fn test_box_blur_fast_kernels_1d() {
+        assert_eq!(box_blur_fast_kernels_1d(0.5, 3), vec![1, 1, 1]);
+        assert_eq!(box_blur_fast_kernels_1d(0.5, 4), vec![1, 1, 1, 1]);
+        assert_eq!(box_blur_fast_kernels_1d(0.5, 5), vec![1, 1, 1, 1, 1]);
+
+        assert_eq!(box_blur_fast_kernels_1d(1.0, 3), vec![1, 1, 3]);
+        assert_eq!(box_blur_fast_kernels_1d(1.0, 5), vec![1, 1, 1, 1, 3]);
+    }
 }
diff --git a/crates/kornia-imgproc/src/filter/ops.rs b/crates/kornia-imgproc/src/filter/ops.rs
@@ -1,6 +1,6 @@
-use kornia_image::{Image, ImageError};
+use kornia_image::{Image, ImageError, ImageSize};
 
-use super::{kernels, separable_filter};
+use super::{fast_horizontal_filter, kernels, separable_filter};
 
 /// Blur an image using a box blur filter
 ///
@@ -79,3 +79,108 @@ pub fn sobel<const C: usize>(
 
     Ok(())
 }
+
+/// Blur an image using a box blur filter multiple times to achieve a near gaussian blur
+///
+/// # Arguments
+///
+/// * `src` - The source image with shape (H, W, C).
+/// * `dst` - The destination image with shape (H, W, C).
+/// * `kernel_size` - The size of the kernel (kernel_x, kernel_y).
+/// * `sigma` - The sigma of the gaussian kernel, xy-ordered.
+///
+/// PRECONDITION: `src` and `dst` must have the same shape.
+pub fn box_blur_fast<const C: usize>(
+    src: &Image<f32, C>,
+    dst: &mut Image<f32, C>,
+    sigma: (f32, f32),
+) -> Result<(), ImageError> {
+    let half_kernel_x_sizes = kernels::box_blur_fast_kernels_1d(sigma.0, 3);
+    let half_kernel_y_sizes = kernels::box_blur_fast_kernels_1d(sigma.1, 3);
+
+    let transposed_size = ImageSize {
+        width: src.size().height,
+        height: src.size().width,
+    };
+
+    let mut input_img = src.clone();
+    let mut transposed = Image::<f32, C>::from_size_val(transposed_size, 0.0)?;
+
+    for (half_kernel_x_size, half_kernel_y_size) in
+        half_kernel_x_sizes.iter().zip(half_kernel_y_sizes.iter())
+    {
+        fast_horizontal_filter(&input_img, &mut transposed, *half_kernel_x_size)?;
+        fast_horizontal_filter(&transposed, dst, *half_kernel_y_size)?;
+
+        input_img = dst.clone();
+    }
+
+    Ok(())
+}
+
+#[cfg(test)]
+mod tests {
+    use super::*;
+
+    #[test]
+    fn test_box_blur_fast() -> Result<(), ImageError> {
+        let size = ImageSize {
+            width: 5,
+            height: 5,
+        };
+
+        #[rustfmt::skip]
+        let img = Image::new(
+            size,
+            (0..25).map(|x| x as f32).collect(),
+        )?;
+        let mut dst = Image::<_, 1>::from_size_val(size, 0.0)?;
+
+        box_blur_fast(&img, &mut dst, (0.5, 0.5))?;
+
+        #[rustfmt::skip]
+        assert_eq!(
+            dst.as_slice(),
+            &[
+                4.444444, 4.9259257, 5.7037034, 6.4814816, 6.962963,
+                6.851851, 7.3333335, 8.111111, 8.888889, 9.370372,
+                10.740741, 11.222222, 12.0, 12.777779, 13.259262,
+                14.629628, 15.111112, 15.888888, 16.666666, 17.14815,
+                17.037035, 17.518518, 18.296295, 19.074074, 19.555555,
+            ],
+        );
+
+        Ok(())
+    }
+
+    #[test]
+    fn test_gaussian_blur() -> Result<(), ImageError> {
+        let size = ImageSize {
+            width: 5,
+            height: 5,
+        };
+
+        #[rustfmt::skip]
+        let img = Image::new(
+            size,
+            (0..25).map(|x| x as f32).collect(),
+        )?;
+
+        let mut dst = Image::<_, 1>::from_size_val(size, 0.0)?;
+
+        gaussian_blur(&img, &mut dst, (3, 3), (0.5, 0.5))?;
+
+        #[rustfmt::skip]
+        assert_eq!(
+            dst.as_slice(),
+            &[
+                0.57097936, 1.4260278, 2.3195207, 3.213014, 3.5739717,
+                4.5739717, 5.999999, 7.0, 7.999999, 7.9349294,
+                9.041435, 10.999999, 12.0, 12.999998, 12.402394,
+                13.5089, 15.999998, 17.0, 17.999996, 16.86986,
+                15.58594, 18.230816, 19.124311, 20.017801, 18.588936,
+            ]
+        );
+        Ok(())
+    }
+}
diff --git a/crates/kornia-imgproc/src/filter/separable_filter.rs b/crates/kornia-imgproc/src/filter/separable_filter.rs
@@ -88,6 +88,70 @@ pub fn separable_filter<const C: usize>(
     Ok(())
 }
 
+/// Apply a fast filter horizontally using cumulative kernel
+///
+/// # Arguments
+///
+/// * `src` - The source image with shape (H, W, C).
+/// * `dst` - The destination image with transposed shape (W, H, C).
+/// * `half_kernel_x_size` - Half of the kernel at weight 1. The total size would be 2*this+1
+pub(crate) fn fast_horizontal_filter<const C: usize>(
+    src: &Image<f32, C>,
+    dst: &mut Image<f32, C>,
+    half_kernel_x_size: usize,
+) -> Result<(), ImageError> {
+    let src_data = src.as_slice();
+    let dst_data = dst.as_slice_mut();
+    let mut row_acc = [0.0; C];
+
+    let mut leftmost_pixel = [0.0; C];
+    let mut rightmost_pixel = [0.0; C];
+
+    let pixels_between_first_last_cols = (src.cols() - 1) * C;
+    let kernel_pix_offset_diffs: Vec<usize> =
+        (0..half_kernel_x_size).map(|p| (p + 1) * C).collect();
+    for (pix_offset, source_pixel) in src_data.iter().enumerate() {
+        let ch = pix_offset % C;
+        let rc = pix_offset / C;
+        let c = rc % src.cols();
+        let r = rc / src.cols();
+
+        let transposed_r = c;
+        let transposed_c = r;
+        let transposed_pix_offset = transposed_r * src.rows() * C + transposed_c * C + ch;
+
+        if c == 0 {
+            row_acc[ch] = *source_pixel * (half_kernel_x_size + 1) as f32;
+            for pix_diff in &kernel_pix_offset_diffs {
+                row_acc[ch] += src_data[pix_offset + pix_diff]
+            }
+            leftmost_pixel[ch] = *source_pixel;
+            rightmost_pixel[ch] = src_data[pix_offset + pixels_between_first_last_cols];
+        } else {
+            row_acc[ch] -= match c.checked_sub(half_kernel_x_size + 1) {
+                Some(_) => {
+                    let prv_leftmost_pix_offset = pix_offset - C * (half_kernel_x_size + 1);
+                    src_data[prv_leftmost_pix_offset]
+                }
+                None => leftmost_pixel[ch],
+            };
+
+            let rightmost_x = c + half_kernel_x_size;
+
+            row_acc[ch] += match rightmost_x {
+                x if x < src.cols() => {
+                    let rightmost_pix_offset = pix_offset + C * half_kernel_x_size;
+                    src_data[rightmost_pix_offset]
+                }
+                _ => rightmost_pixel[ch],
+            };
+        }
+        dst_data[transposed_pix_offset] = row_acc[ch] / (half_kernel_x_size * 2 + 1) as f32;
+    }
+
+    Ok(())
+}
+
 #[cfg(test)]
 mod tests {
     use super::*;
@@ -134,4 +198,60 @@ mod tests {
 
         Ok(())
     }
+
+    #[test]
+    fn test_fast_horizontal_filter() -> Result<(), ImageError> {
+        let size = ImageSize {
+            width: 5,
+            height: 5,
+        };
+
+        #[rustfmt::skip]
+        let img = Image::new(
+            size,
+            vec![
+                0.0, 0.0, 0.0, 0.0, 0.0,
+                0.0, 0.0, 0.0, 0.0, 0.0,
+                0.0, 0.0, 9.0, 0.0, 0.0,
+                0.0, 0.0, 0.0, 0.0, 0.0,
+                0.0, 0.0, 0.0, 0.0, 0.0,
+            ],
+        )?;
+
+        let mut transposed = Image::<_, 1>::from_size_val(size, 0.0)?;
+
+        fast_horizontal_filter(&img, &mut transposed, 1)?;
+
+        #[rustfmt::skip]
+        assert_eq!(
+            transposed.as_slice(),
+            &[
+                0.0, 0.0, 0.0, 0.0, 0.0,
+                0.0, 0.0, 3.0, 0.0, 0.0,
+                0.0, 0.0, 3.0, 0.0, 0.0,
+                0.0, 0.0, 3.0, 0.0, 0.0,
+                0.0, 0.0, 0.0, 0.0, 0.0,
+            ]
+        );
+
+        let mut dst = Image::<_, 1>::from_size_val(size, 0.0)?;
+
+        fast_horizontal_filter(&transposed, &mut dst, 1)?;
+
+        #[rustfmt::skip]
+        assert_eq!(
+            dst.as_slice(),
+            &[
+                0.0, 0.0, 0.0, 0.0, 0.0,
+                0.0, 1.0, 1.0, 1.0, 0.0,
+                0.0, 1.0, 1.0, 1.0, 0.0,
+                0.0, 1.0, 1.0, 1.0, 0.0,
+                0.0, 0.0, 0.0, 0.0, 0.0,
+            ]
+        );
+        let xsum = dst.as_slice().iter().sum::<f32>();
+        assert_eq!(xsum, 9.0);
+
+        Ok(())
+    }
 }
