Implement Box blur fast filter that could approximate gaussian filter

crates/kornia-imgproc/src/filter/kernels.rs

@@ -61,6 +61,38 @@ 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 +124,15 @@ 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]);
+    }
 }

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::{kernels, separable_filter, fast_horizontal_filter};
 
 /// Blur an image using a box blur filter
 ///
@@ -79,3 +79,107 @@ 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();
+    for (half_kernel_x_size, half_kernel_y_size) in half_kernel_x_sizes.iter().zip(half_kernel_y_sizes.iter()) {
+        let mut transposed = Image::<f32, C>::from_size_val(transposed_size, 0.0)?;
+
+        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))?;
+
+        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))?;
+
+        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(())
+    }
+}

crates/kornia-imgproc/src/filter/separable_filter.rs

@@ -88,6 +88,72 @@ pub fn separable_filter<const C: usize>(
     Ok(())
 }
 
+
+/// Apply a fast filter horizontally, take advantage of property where all
+/// weights are equal to filter an image at O(n) time
+/// 
+/// # 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];
+    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;
+            let mut kernel_pix_offset = pix_offset;
+            for _ in 0..half_kernel_x_size {
+                kernel_pix_offset += C;
+                row_acc[ch] += src_data[kernel_pix_offset];
+            }
+            leftmost_pixel[ch] = *source_pixel;
+            rightmost_pixel[ch] = src_data[pix_offset+((src.cols()-1)*C)];
+        } 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 +200,58 @@ mod tests {
 
         Ok(())
     }
+
+
+    #[test]
+    fn test_fast_horizontal_filter() -> Result<(), ImageError> {
+        let size = ImageSize {
+            width: 5,
+            height: 5,
+        };
+
+        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)?;
+
+        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)?;
+
+        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(())
+    }
 }