Refactor kernel module

src/kernel.rs

@@ -1,29 +1,20 @@
 use crate::*;
 
+pub mod gaussian;
+pub use gaussian::*;
+
 pub trait Kernel {
     fn convolve<P>(&self, img: &image::Image<P>) -> image::Image<P>
     where
         P: image::Pixel;
 }
 
-fn sum<const N: usize>(kernel: &[[f32; N]; N]) -> f32 {
-    let mut x = 0.0;
-    for i in 0..N {
-        for j in 0..N {
-            x += kernel[j][i];
-        }
-    }
-    x
-}
-
 impl<const N: usize> Kernel for [[f32; N]; N] {
     fn convolve<P>(&self, img: &image::Image<P>) -> image::Image<P>
     where
         P: image::Pixel,
     {
         let l = N / 2;
-        let mag = sum(self);
-
         let extent = image::Extent {
             x: img.extent.x - 2 * l,
             y: img.extent.y - 2 * l,
@@ -32,10 +23,9 @@ impl<const N: usize> Kernel for [[f32; N]; N] {
         result.for_each(|x, y, p| {
             for i in 0..N {
                 for j in 0..N {
-                    *p += img.at(x + i, y + j) * self[j][i];
+                    *p += img.at(x + j, y + i) * self[i][j];
                 }
             }
-            *p /= mag;
         });
         result
     }
@@ -83,24 +73,22 @@ mod tests {
             let input = image::Image {
                 extent: image::Extent { x: 5, y: 5 },
                 pixels: vec![
-					10.0, 10.0, 10.0, 10.0, 10.0,
+                    10.0, 10.0, 10.0, 10.0, 10.0,
-					10.0, p1,   p2,   p3,   10.0,
+                    10.0, p1,   p2,   p3,   10.0,
-					10.0, p4,   p5,   p6,   10.0,
+                    10.0, p4,   p5,   p6,   10.0,
-					10.0, p7,   p8,   p9,   10.0,
+                    10.0, p7,   p8,   p9,   10.0,
-					10.0, 10.0, 10.0, 10.0, 10.0,
+                    10.0, 10.0, 10.0, 10.0, 10.0,
-				],
+                ],
             };
             let kernel = [[k1, k2, k3], [k4, k5, k6], [k7, k8, k9]];
             let output = kernel.convolve(&input);
             prop_assert_eq!(output.extent.x, 3);
             prop_assert_eq!(output.extent.y, 3);
 
-			let mag = k1 + k2 + k3 + k4 + k5 + k6 + k7 + k8 + k9;
+            let exp = p1 * k1 + p2 * k2 + p3 * k3
-			let sum = p1 * k1 + p2 * k2 + p3 * k3
+                + p4 * k4 + p5 * k5 + p6 * k6
-				+ p4 * k4 + p5 * k5 + p6 * k6
+                + p7 * k7 + p8 * k8 + p9 * k9;
-				+ p7 * k7 + p8 * k8 + p9 * k9;
+            prop_assert!((output.at(1, 1) - exp).abs() < 1e-6);
-			let exp = sum / mag;
-			prop_assert!((output.at(1, 1) - exp).abs() < 1e-6);
         }
     }
 }

src/kernel/gaussian.rs

@@ -0,0 +1,25 @@
+pub fn gaussian<const N: usize>(sigma: f32) -> [[f32; N]; N] {
+	let mut kernel = [[0.0; N]; N];
+
+    let midpoint = (N / 2 + 1) as f32;
+	let v = sigma * sigma;
+	let k = 1.0 / (2.0 * std::f32::consts::PI * v).sqrt();
+	let mut sum = 0.0;
+	for y in 0..N {
+		for x in 0..N {
+			let dx = x as f32 - midpoint;
+			let dy = y as f32 - midpoint;
+			let s = (dx * dx + dy * dy).sqrt();
+			kernel[y][x] = k * (- s / (2.0 * v)).exp();
+			sum += kernel[y][x];
+		}
+	}
+
+	for y in 0..N {
+		for x in 0..N {
+			kernel[y][x] /= sum;
+		}
+	}
+
+	kernel
+}