initial implementation of fast blur

src/imageops/fast_blur.rs

@@ -0,0 +1,224 @@
+use std::{
+    cmp::{max, min},
+    ops::Deref,
+};
+
+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, Container>(image_buffer: &mut ImageBuffer<P, Container>, sigma: f32)
+where
+    Container: Deref<Target = [P::Subpixel]> + AsMut<[P::Subpixel]>,
+{
+    let (width, height) = image_buffer.dimensions();
+    let samples = image_buffer.as_flat_samples_mut().samples;
+
+    let num_passes = 3;
+
+    let boxes = boxes_for_gauss(sigma, num_passes);
+
+    for pass in 0..num_passes {
+        let radius = boxes[pass];
+        let mut horizontally_blurred = my_fast_horizontal_blur::<P::Subpixel>(
+            samples,
+            width as usize,
+            height as usize,
+            radius,
+            P::CHANNEL_COUNT as usize,
+        );
+        let totally_blurred = my_fast_vertical_blur::<P::Subpixel>(
+            &mut horizontally_blurred,
+            width as usize,
+            height as usize,
+            radius,
+            P::CHANNEL_COUNT as usize,
+        );
+        samples.copy_from_slice(&totally_blurred);
+    }
+}
+
+fn boxes_for_gauss(sigma: f32, n: usize) -> Vec<usize> {
+    let w_ideal = f32::sqrt((12.0 * sigma * sigma / (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 = (12.0 * sigma * sigma
+        - (n as f32) * (w_l) * (w_l)
+        - 4.0 * (n as f32) * (w_l)
+        - 3.0 * (n as f32))
+        / (-4.0 * (w_l) - 4.0);
+    let m = f32::round(m_ideal) as usize;
+
+    let mut box_sizes: Vec<usize> = vec![];
+    for i in 0..n {
+        box_sizes.push(if i < m { w_l as usize } else { w_u as usize });
+    }
+    return box_sizes;
+}
+
+fn channel_idx(channel: usize, idx: usize) -> usize {
+    3 * idx + channel
+}
+
+fn my_fast_horizontal_blur<P: Primitive>(
+    samples: &mut [P],
+    width: usize,
+    height: usize,
+    r: usize,
+    channel_num: usize,
+) -> Vec<P> {
+    let channel_size = width * height;
+
+    let mut out_samples: Vec<P> = vec![P::from(0).unwrap(); channel_size * 3];
+    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 i in 0..height {
+        for channel in 0..channel_num {
+            vals[channel] = {
+                let mut val = 0.0;
+                for x in (-(r as isize))..(r + 1) as isize {
+                    val += extended_f(samples, width, height, x, i as isize, channel)
+                        .to_f32()
+                        .unwrap_or(0.0);
+                }
+                val
+            };
+        }
+
+        for j in 0..width {
+            for channel in 0..channel_num {
+                let val = vals[channel] / (2.0 * r as f32 + 1.0);
+                let val = if val < min_value {
+                    min_value
+                } else if val > max_value {
+                    max_value
+                } else {
+                    val
+                };
+                let val = P::from(val).unwrap();
+
+                out_samples[channel_idx(channel, i * width + j)] = val;
+                vals[channel] = vals[channel]
+                    - extended_f(
+                        samples,
+                        width,
+                        height,
+                        j as isize - r as isize,
+                        i as isize,
+                        channel,
+                    )
+                    .to_f32()
+                    .unwrap_or(0.0)
+                    + extended_f(
+                        samples,
+                        width,
+                        height,
+                        { j + r + 1 } as isize,
+                        i as isize,
+                        channel,
+                    )
+                    .to_f32()
+                    .unwrap_or(0.0)
+            }
+        }
+    }
+
+    out_samples
+}
+
+fn my_fast_vertical_blur<P: Primitive>(
+    samples: &mut [P],
+    width: usize,
+    height: usize,
+    r: usize,
+    channel_num: usize,
+) -> Vec<P> {
+    let channel_size = width * height;
+
+    let mut out_samples: Vec<P> = 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 j in 0..width {
+        for channel in 0..channel_num {
+            vals[channel] = {
+                let mut val = 0.0;
+                for y in (-(r as isize))..(r + 1) as isize {
+                    val += extended_f(samples, width, height, j as isize, y, channel)
+                        .to_f32()
+                        .unwrap();
+                }
+                val
+            };
+        }
+
+        for i in 0..height {
+            for channel in 0..channel_num {
+                let val = vals[channel] / (2.0 * (r as f32) + 1.0);
+                let val = if val < min_value {
+                    min_value
+                } else if val > max_value {
+                    max_value
+                } else {
+                    val
+                };
+                let valc = P::from(val);
+                let val = match valc {
+                    Some(a) => a,
+                    None => {
+                        dbg!(val);
+                        panic!()
+                    }
+                };
+
+                out_samples[channel_idx(channel, i * width + j)] = val;
+                vals[channel] = vals[channel]
+                    - extended_f(
+                        samples,
+                        width,
+                        height,
+                        j as isize,
+                        i as isize - r as isize,
+                        channel,
+                    )
+                    .to_f32()
+                    .unwrap()
+                    + extended_f(
+                        samples,
+                        width,
+                        height,
+                        j as isize,
+                        i as isize + (r + 1) as isize,
+                        channel,
+                    )
+                    .to_f32()
+                    .unwrap();
+            }
+        }
+    }
+
+    out_samples
+}
+
+fn extended_f<P: Primitive>(
+    samples: &mut [P],
+    width: usize,
+    height: usize,
+    x: isize,
+    y: isize,
+    channel: usize,
+) -> P {
+    let x = min(width as isize - 1, max(0, x)) as usize;
+    let y = min(height as isize - 1, max(0, y)) as usize;
+    samples[channel_idx(channel, y * width + x)]
+}

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>(