Implement basic Hough transform

src/hough.rs

@@ -0,0 +1,87 @@
+use crate::*;
+
+const THETA_N: usize = 200;
+
+#[derive(Clone, Copy, Debug)]
+pub struct Line {
+    pub theta: f32,
+    pub rho: f32,
+    pub significance: f32,
+}
+
+pub fn transform(img: image::Image<image::Luminance>) -> image::Image<f32> {
+    let diag = ((img.extent.x.pow(2) + img.extent.y.pow(2)) as f32).sqrt();
+    let mut res = image::Image::<f32>::with_extent(image::Extent {
+        x: THETA_N,
+        y: diag as usize,
+    });
+    for y in 0..img.extent.y {
+        for x in 0..img.extent.x {
+            if img.at(x, y).y != 1.0 {
+                continue;
+            }
+            for a in 0..THETA_N {
+                let theta = theta_from(a);
+                let rho = x as f32 * theta.cos() + y as f32 * theta.sin();
+                if rho < 0.0 {
+                    continue;
+                }
+
+                *res.at_mut(a, rho as usize) += 1.0;
+            }
+        }
+    }
+    res
+}
+
+pub fn lines(
+    min_significance: f32,
+    img: image::Image<image::Luminance>,
+) -> Vec<Line> {
+    let hough = transform(img);
+
+    let mu = mean((&hough.pixels).into_iter().map(|p| *p));
+    let sigma = stddev((&hough.pixels).into_iter().map(|p| *p), mu);
+
+    let mut lines = Vec::<Line>::new();
+    for y in 0..hough.extent.y {
+        for x in 0..hough.extent.x {
+            let value = hough.at(x, y);
+            let significance = (value - mu) / sigma;
+            if significance >= min_significance {
+                let theta = theta_from(x);
+                let rho = y as f32;
+                lines.push(Line {
+                    theta,
+                    rho,
+                    significance,
+                });
+            }
+        }
+    }
+	lines
+}
+
+fn theta_from(x: usize) -> f32 {
+	x as f32 * std::f32::consts::PI / THETA_N as f32
+}
+
+fn mean(xs: impl Iterator<Item = f32>) -> f32 {
+    let mut n = 0;
+    let mut sum = 0.0;
+    for x in xs {
+        sum += x;
+        n += 1;
+    }
+    sum / n as f32
+}
+
+fn stddev(xs: impl Iterator<Item = f32>, mean: f32) -> f32 {
+    let mut n = 0;
+    let mut sum = 0.0;
+    for x in xs {
+        sum += (x - mean).powf(2.0);
+        n += 1;
+    }
+    (sum / (n - 1) as f32).sqrt()
+}

src/lib.rs

@@ -1,3 +1,4 @@
 pub mod filter;
+pub mod hough;
 pub mod image;
 pub mod kernel;