From 7e62163c26834965da708691c0773725562679e9 Mon Sep 17 00:00:00 2001
From: Camden Dixie O'Brien <cdo@wip.sh>
Date: Tue, 27 Oct 2026 19:09:27 +0000
Subject: [PATCH] Implement automatic double threshold for Canny filter

---
 src/filter/canny/double_threshold.rs | 85 ++++++++++++++++++++++++++++
 1 file changed, 85 insertions(+)
 create mode 100644 src/filter/canny/double_threshold.rs

diff --git a/src/filter/canny/double_threshold.rs b/src/filter/canny/double_threshold.rs
new file mode 100644
index 0000000..c7b313e
--- /dev/null
+++ b/src/filter/canny/double_threshold.rs
@@ -0,0 +1,85 @@
+use crate::image::*;
+
+#[derive(Clone, Copy, PartialEq)]
+pub enum Candidate {
+    None,
+    Weak,
+    Strong,
+}
+
+impl Default for Candidate {
+    fn default() -> Self {
+        Self::None
+    }
+}
+
+pub fn double_threshold(img: Image<f32>) -> Image<Candidate> {
+    let t = otsu::<20>(&img.pixels);
+    let mut result = Image::<Candidate>::with_extent(img.extent);
+    result.for_each(|x, y, p| {
+        let value = img.at(x, y);
+        *p = if value < t / 2.0 {
+            Candidate::None
+        } else if value < t {
+            Candidate::Weak
+        } else {
+            Candidate::Strong
+        }
+    });
+    result
+}
+
+fn otsu<const N: usize>(values: &Vec<f32>) -> f32 {
+    let mut hist = [0.0; N];
+    for x in values {
+        let bin = std::cmp::min(N - 1, (*x * (N as f32)) as usize);
+        hist[bin] += 1.0;
+    }
+    for i in 0..N {
+        hist[i] /= values.len() as f32;
+    }
+
+    let mut max = 0.0;
+    let mut res = 0;
+    for t in 0..N {
+        let var = variance(t, hist);
+        if var > max {
+            res = t;
+            max = var;
+        }
+    }
+
+    res as f32 / N as f32
+}
+
+fn variance<const N: usize>(t: usize, hist: [f32; N]) -> f32 {
+    let omega_0 = omega(0, t, hist);
+    let omega_1 = omega(t, N, hist);
+
+    let mu_0 = mu(0, t, hist, omega_0);
+    let mu_1 = mu(t, N, hist, omega_1);
+
+    omega_0 * omega_1 * (mu_0 - mu_1).powf(2.0)
+}
+
+fn omega<const N: usize>(from: usize, to: usize, hist: [f32; N]) -> f32 {
+    let mut sum = 0.0;
+    for i in from..to {
+        sum += hist[i];
+    }
+    sum
+}
+
+fn mu<const N: usize>(
+    from: usize,
+    to: usize,
+    hist: [f32; N],
+    omega: f32,
+) -> f32 {
+    let mut sum = 0.0;
+    for i in from..to {
+        let bin_centre = (i as f32 + 0.5) / N as f32;
+        sum += bin_centre * hist[i];
+    }
+    sum / omega
+}