Create stats module

LICENSE.md

@@ -0,0 +1,108 @@
+# Komorebi License
+
+Version 2.0.0
+
+## Acceptance
+
+In order to get any license under these terms, you must agree
+to them as both strict obligations and conditions to all
+your licenses.
+
+## Copyright License
+
+The licensor grants you a copyright license for the software
+to do everything you might do with the software that would
+otherwise infringe the licensor's copyright in it for any
+permitted purpose. However, you may only distribute the source
+code of the software according to the [Distribution License](
+#distribution-license), you may only make changes according
+to the [Changes License](#changes-license), and you may not
+otherwise distribute the software or new works based on the
+software.
+
+## Distribution License
+
+The licensor grants you an additional copyright license to
+distribute copies of the source code of the software. Your
+license to distribute covers distributing the source code of
+the software with changes permitted by the [Changes License](
+#changes-license).
+
+## Changes License
+
+The licensor grants you an additional copyright license to
+make changes for any permitted purpose.
+
+## Patent License
+
+The licensor grants you a patent license for the software that
+covers patent claims the licensor can license, or becomes able
+to license, that you would infringe by using the software.
+
+## Personal Uses
+
+Personal use for research, experiment, and testing for
+the benefit of public knowledge, personal study, private
+entertainment, hobby projects, amateur pursuits, or religious
+observance, without any anticipated commercial application,
+is use for a permitted purpose.
+
+## Fair Use
+
+You may have "fair use" rights for the software under the
+law. These terms do not limit them.
+
+## No Other Rights
+
+These terms do not allow you to sublicense or transfer any of
+your licenses to anyone else, or prevent the licensor from
+granting licenses to anyone else.  These terms do not imply
+any other licenses.
+
+## Patent Defense
+
+If you make any written claim that the software infringes or
+contributes to infringement of any patent, your patent license
+for the software granted under these terms ends immediately. If
+your company makes such a claim, your patent license ends
+immediately for work on behalf of your company.
+
+## Violations
+
+The first time you are notified in writing that you have
+violated any of these terms, or done anything with the software
+not covered by your licenses, your licenses can nonetheless
+continue if you come into full compliance with these terms,
+and take practical steps to correct past violations, within
+32 days of receiving notice.  Otherwise, all your licenses
+end immediately.
+
+## No Liability
+
+***As far as the law allows, the software comes as is, without
+any warranty or condition, and the licensor will not be liable
+to you for any damages arising out of these terms or the use
+or nature of the software, under any kind of legal claim.***
+
+## Definitions
+
+The **licensor** is the individual or entity offering these
+terms, and the **software** is the software the licensor makes
+available under these terms.
+
+**You** refers to the individual or entity agreeing to these
+terms.
+
+**Your company** is any legal entity, sole proprietorship,
+or other kind of organization that you work for, plus all
+organizations that have control over, are under the control of,
+or are under common control with that organization.  **Control**
+means ownership of substantially all the assets of an entity,
+or the power to direct its management and policies by vote,
+contract, or otherwise.  Control can be direct or indirect.
+
+**Your licenses** are all the licenses granted to you for the
+software under these terms.
+
+**Use** means anything you do with the software requiring one
+of your licenses.

README.md

@@ -0,0 +1,31 @@
+# Cam's Computer Vision Library
+
+A basic library of (not remotely production-ready) computer vision
+functions.
+
+## Building
+
+CCVL is written in Rust and uses Cargo as its build system.  The Rust
+toolchain (including Cargo) is most easily installed with
+[Rustup](https://rustup.rs/).
+
+Given that you have a toolchain installed, the library can be built as
+usual for a Cargo project:
+
+```
+cargo build
+```
+
+You can also run one of the example binaries to see it in action:
+
+```
+cargo run --example canny foo.ff foo-edges.ff
+```
+
+Only the [farbfeld](https://tools.suckless.org/farbfeld/) image format
+is supported for input and output at the moment.
+
+## License
+
+Copyright © Camden Dixie O'Brien, available under the Komorebi 2.0.0
+license.  See [LICENSE.md](LICENSE.md) for details.

examples/contours.rs

@@ -0,0 +1,38 @@
+mod common;
+
+use ccvl::kernel::Kernel;
+use ccvl::*;
+
+fn draw(
+    img: &mut image::Image<image::SRGBA>,
+    path: &Vec<(usize, usize)>,
+    colour: image::SRGBA,
+) {
+    let n = path.len();
+    for i in 0..n {
+        let (x0, y0) = path[i];
+        let (x1, y1) = path[(i + 1) % n];
+        draw::line(x0, y0, x1, y1, img, colour);
+    }
+}
+
+fn approx(a: f32, b: f32, epsilon: f32) -> bool {
+    (a - b).abs() <= epsilon
+}
+
+fn main() -> std::process::ExitCode {
+    common::process_image(|mut img| {
+        let greyscale = img.map(|p| image::Luminance::from(p));
+        let blurred = kernel::gaussian::<5>(1.2).convolve_sat(&greyscale);
+        let thresholded = filter::threshold_auto(blurred);
+        let contours =
+            contour::find(&thresholded, image::Luminance { y: 0.0 });
+        for contour in contours {
+            draw(&mut img, &contour.vertices, match contour.kind {
+				contour::Kind::Outer => colours::GREEN,
+				contour::Kind::Inner => colours::BLUE,
+			});
+        }
+        img
+    })
+}

examples/hough.rs

@@ -0,0 +1,51 @@
+mod common;
+
+use ccvl::*;
+
+fn main() -> std::process::ExitCode {
+    common::process_image(|mut img| {
+        let greyscale = img.map(|p| image::Luminance::from(p));
+        let edges = filter::canny(greyscale);
+        for line in hough::lines(10.0, edges) {
+			let ((x0, y0), (x1, y1)) = endpoints(line, img.extent);
+			draw::line(x0, y0, x1, y1, &mut img, colours::RED);
+        }
+        img
+    })
+}
+
+fn endpoints(
+    l: hough::Line,
+    lim: image::Extent,
+) -> ((usize, usize), (usize, usize)) {
+    if l.theta.abs() < 1e-3 {
+        ((l.rho as usize, 0), (l.rho as usize, lim.y - 1))
+    } else if l.theta - std::f32::consts::PI / 2.0 < 1e-3 {
+        ((0, l.rho as usize), (lim.x - 1, l.rho as usize))
+    } else {
+        let w = lim.x as f32 - 1e-3;
+        let h = lim.y as f32 - 1e-3;
+
+        let ros = l.rho / l.theta.sin();
+        let roc = l.rho / l.theta.cos();
+        let t = l.theta.tan();
+
+        let candidates: [(f32, f32); 4] =
+            [(roc - h * t, h), (0.0, ros), (roc, 0.0), (w, ros - w / t)];
+        let res = candidates
+            .into_iter()
+            .filter(|p| in_bounds(p, lim))
+            .map(|(x, y)| (x as usize, y as usize))
+            .collect::<Vec<_>>();
+        if res.len() != 2 {
+            panic!("Unexpected number of endpoints for {:?}", l);
+        }
+        (res[0], res[1])
+    }
+}
+
+fn in_bounds(q: &(f32, f32), lim: image::Extent) -> bool {
+    let x_in = 0.0 <= q.0 && q.0 < lim.x as f32;
+    let y_in = 0.0 <= q.1 && q.1 < lim.y as f32;
+    x_in && y_in
+}

src/colours.rs

@@ -0,0 +1,21 @@
+use crate::*;
+
+macro_rules! hex {
+	( $value:literal ) => {
+		image::SRGBA {
+			r: (($value >> 16) & 0xff) as f32 / 255.0,
+			g: (($value >> 8) & 0xff) as f32 / 255.0,
+			b: ($value & 0xff) as f32 / 255.0,
+			a: 1.0,
+		}
+	}
+}
+
+pub const BLACK: image::SRGBA = hex!(0x000000);
+pub const RED: image::SRGBA = hex!(0xff0000);
+pub const GREEN: image::SRGBA = hex!(0x00ff00);
+pub const BLUE: image::SRGBA = hex!(0x0000ff);
+pub const CYAN: image::SRGBA = hex!(0x00ffff);
+pub const MAGENTA: image::SRGBA = hex!(0xff00ff);
+pub const YELLOW: image::SRGBA = hex!(0xffff00);
+pub const WHITE: image::SRGBA = hex!(0xffffff);

src/contour.rs

@@ -0,0 +1,437 @@
+use crate::*;
+
+#[derive(Clone, Copy, Debug, PartialEq)]
+pub enum Kind {
+    Outer,
+    Inner,
+}
+
+#[derive(Debug, PartialEq)]
+pub struct Contour {
+    pub vertices: Vec<(usize, usize)>,
+    pub children: Vec<usize>,
+    pub parent: Option<usize>,
+    pub kind: Kind,
+}
+
+#[derive(Clone, Copy, Debug, PartialEq)]
+enum Pixel {
+    Foreground,
+    Background,
+    TracedUp(usize),
+    TracedDown(usize),
+}
+
+#[derive(Clone, Copy, Debug, PartialEq)]
+enum Dir {
+    N,
+    NE,
+    E,
+    SE,
+    S,
+    SW,
+    W,
+    NW,
+}
+
+trait Coordinate {
+    fn prev_x(&self) -> Self;
+    fn next_x(&self) -> Self;
+    fn prev_y(&self) -> Self;
+    fn next_y(&self) -> Self;
+
+    fn in_bounds(&self, bounds: image::Extent) -> bool;
+    fn get<P: Copy>(&self, img: &image::Image<P>) -> Option<P>;
+    fn set<P>(&self, img: &mut image::Image<P>, value: P);
+}
+
+pub fn find<P>(img: &image::Image<P>, fg: P) -> Vec<Contour>
+where
+    P: PartialEq + Copy + Send + Sync,
+{
+    let mut map = img.map(|p| {
+        if p == fg {
+            Pixel::Foreground
+        } else {
+            Pixel::Background
+        }
+    });
+
+    let mut contours = Vec::new();
+    for y in 0..map.extent.y {
+        let mut parent = None;
+        for x in 0..map.extent.x {
+            let q = (x, y);
+            if foreground(&map, q) && background(&map, q.prev_x()) {
+                add(q, Kind::Outer, parent, &mut contours, &mut map);
+            } else if (foreground(&map, q) || traced_up(&map, q).is_some())
+                && background(&map, q.next_x())
+            {
+                add(q, Kind::Inner, parent, &mut contours, &mut map);
+            } else if let Some(label) = traced_up(&map, q) {
+                parent = Some(label);
+            } else if let Some(label) = traced_down(&map, q) {
+                parent = contours[label].parent;
+            }
+        }
+    }
+    contours
+}
+
+fn add(
+    start: (usize, usize),
+    kind: Kind,
+    parent: Option<usize>,
+    contours: &mut Vec<Contour>,
+    map: &mut image::Image<Pixel>,
+) {
+    let label = contours.len();
+    let dir = kind.find_start_dir(start, map);
+    let vertices = trace(map, start, dir, label);
+    contours.push(Contour::new(parent, kind, vertices));
+    if let Some(index) = parent {
+        contours[index].add_child(label);
+    }
+}
+
+fn trace(
+    img: &mut image::Image<Pixel>,
+    start_q: (usize, usize),
+    start_dir: Dir,
+    label: usize,
+) -> Vec<(usize, usize)> {
+    let mut q = start_q;
+    let mut dir = start_dir;
+
+    let mut vertices = Vec::new();
+    loop {
+        if foreground(img, q) {
+            q.set(img, Pixel::TracedUp(label));
+        }
+
+        let prev_dir = dir;
+        dir = dir.flip().next();
+        q = loop {
+            let qp = dir.from(q);
+            if !background(img, qp) {
+                if dir != prev_dir {
+                    vertices.push(q)
+                }
+                break qp;
+            }
+            if dir == Dir::E {
+                q.set(img, Pixel::TracedDown(label));
+            }
+            dir = dir.next();
+        };
+
+        if q == start_q && dir == start_dir || !q.in_bounds(img.extent) {
+            break;
+        }
+    }
+
+    // for y in 0..img.extent.y {
+    //     print!("  ");
+    //     for x in 0..img.extent.x {
+    //         if let Some(label) = traced_up(img, (x, y)) {
+    //             print!("{:1}", label);
+    //         } else if let Some(_) = traced_down(img, (x, y)) {
+    //             print!("-")
+    //         } else if foreground(img, (x, y)) {
+    //             print!(".");
+    //         } else if background(img, (x, y)) {
+    //             print!(" ");
+    //         }
+    //     }
+    //     println!();
+    // }
+
+    vertices
+}
+
+fn foreground(img: &image::Image<Pixel>, q: (usize, usize)) -> bool {
+    match q.get(img) {
+        Some(Pixel::Foreground) => true,
+        _ => false,
+    }
+}
+
+fn background(img: &image::Image<Pixel>, q: (usize, usize)) -> bool {
+    match q.get(img) {
+        Some(Pixel::Background) => true,
+        _ => false,
+    }
+}
+
+fn traced_up(img: &image::Image<Pixel>, q: (usize, usize)) -> Option<usize> {
+    match q.get(img) {
+        Some(Pixel::TracedUp(label)) => Some(label),
+        _ => None,
+    }
+}
+
+fn traced_down(
+    img: &image::Image<Pixel>,
+    q: (usize, usize),
+) -> Option<usize> {
+    match q.get(img) {
+        Some(Pixel::TracedDown(label)) => Some(label),
+        _ => None,
+    }
+}
+
+impl Kind {
+    fn find_start_dir(
+        self,
+        q: (usize, usize),
+        img: &image::Image<Pixel>,
+    ) -> Dir {
+        let mut dir = match self {
+            Kind::Outer => Dir::W,
+            Kind::Inner => Dir::E,
+        };
+        while background(img, dir.from(q)) {
+            dir = dir.prev()
+        }
+        dir.flip()
+    }
+}
+
+impl Dir {
+    fn next(self) -> Self {
+        match self {
+            Self::N => Self::NE,
+            Self::NE => Self::E,
+            Self::E => Self::SE,
+            Self::SE => Self::S,
+            Self::S => Self::SW,
+            Self::SW => Self::W,
+            Self::W => Self::NW,
+            Self::NW => Self::N,
+        }
+    }
+
+    fn prev(self) -> Self {
+        match self {
+            Self::N => Self::NW,
+            Self::NE => Self::N,
+            Self::E => Self::NE,
+            Self::SE => Self::E,
+            Self::S => Self::SE,
+            Self::SW => Self::S,
+            Self::W => Self::SW,
+            Self::NW => Self::W,
+        }
+    }
+
+    fn flip(self) -> Self {
+        match self {
+            Self::N => Self::S,
+            Self::NE => Self::SW,
+            Self::E => Self::W,
+            Self::SE => Self::NW,
+            Self::S => Self::N,
+            Self::SW => Self::NE,
+            Self::W => Self::E,
+            Self::NW => Self::SE,
+        }
+    }
+
+    fn from<C: Coordinate>(&self, q: C) -> C {
+        match self {
+            Self::N => q.prev_y(),
+            Self::NE => q.prev_y().next_x(),
+            Self::E => q.next_x(),
+            Self::SE => q.next_y().next_x(),
+            Self::S => q.next_y(),
+            Self::SW => q.next_y().prev_x(),
+            Self::W => q.prev_x(),
+            Self::NW => q.prev_y().prev_x(),
+        }
+    }
+}
+
+impl Coordinate for (usize, usize) {
+    fn prev_x(&self) -> Self {
+        (self.0.wrapping_sub(1), self.1)
+    }
+
+    fn next_x(&self) -> Self {
+        (self.0 + 1, self.1)
+    }
+
+    fn prev_y(&self) -> Self {
+        (self.0, self.1.wrapping_sub(1))
+    }
+
+    fn next_y(&self) -> Self {
+        (self.0, self.1 + 1)
+    }
+
+    fn in_bounds(&self, bounds: image::Extent) -> bool {
+        self.0 < bounds.x && self.1 < bounds.y
+    }
+
+    fn get<P: Copy>(&self, img: &image::Image<P>) -> Option<P> {
+        if self.in_bounds(img.extent) {
+            Some(img.at(self.0, self.1))
+        } else {
+            None
+        }
+    }
+
+    fn set<P>(&self, img: &mut image::Image<P>, value: P) {
+        if self.in_bounds(img.extent) {
+            *img.at_mut(self.0, self.1) = value;
+        }
+    }
+}
+
+impl Contour {
+    fn new(
+        parent: Option<usize>,
+        kind: Kind,
+        vertices: Vec<(usize, usize)>,
+    ) -> Self {
+        Self {
+            vertices,
+            children: Vec::new(),
+            parent,
+            kind,
+        }
+    }
+
+    fn add_child(&mut self, child: usize) {
+        self.children.push(child);
+    }
+}
+
+impl Default for Pixel {
+    fn default() -> Self {
+        Self::Background
+    }
+}
+
+#[cfg(test)]
+#[rustfmt::skip]
+mod tests {
+	use super::*;
+
+	#[test]
+	fn simple_block() {
+		let img = image::Image::<usize> {
+			extent: image::Extent { x: 4, y: 4 },
+			pixels: vec![
+				0, 0, 0, 0,
+				0, 1, 1, 0,
+				0, 1, 1, 0,
+				0, 0, 0, 0,
+			]
+		};
+
+		let contours = find(&img, 1);
+		assert_eq!(contours.len(), 1);
+
+		let contour = &contours[0];
+		assert_eq!(contour.kind, Kind::Outer);
+		assert_eq!(contour.vertices, vec![
+			(1, 1), (2, 1), (2, 2), (1, 2)
+		]);
+	}
+
+	#[test]
+	fn simple_block_with_hole() {
+		let img = image::Image::<usize> {
+			extent: image::Extent { x: 8, y: 8 },
+			pixels: vec![
+				0, 0, 0, 0, 0, 0, 0, 0,
+				0, 1, 1, 1, 1, 1, 1, 0,
+				0, 1, 1, 1, 1, 1, 1, 0,
+				0, 1, 1, 0, 0, 1, 1, 0,
+				0, 1, 1, 0, 0, 1, 1, 0,
+				0, 1, 1, 1, 1, 1, 1, 0,
+				0, 1, 1, 1, 1, 1, 1, 0,
+				0, 0, 0, 0, 0, 0, 0, 0,
+			]
+		};
+
+		let contours = find(&img, 1);
+		assert_eq!(contours.len(), 2);
+
+		let contour = &contours[0];
+		assert_eq!(contour.parent, None);
+		assert_eq!(contour.kind, Kind::Outer);
+		assert_eq!(contour.vertices, vec![
+			(1, 1), (6, 1), (6, 6), (1, 6)
+		]);
+		assert_eq!(contour.children, vec![ 1 ]);
+
+		let contour = &contours[1];
+		assert_eq!(contour.parent, Some(0));
+		assert_eq!(contour.kind, Kind::Inner);
+		assert_eq!(contour.vertices, vec![
+			(2, 3), (2, 4), (3, 5), (4, 5),
+			(5, 4), (5, 3), (4, 2), (3, 2),
+		]);
+	}
+
+	#[test]
+	fn diagonally_connected_block_with_hole() {
+		let img = image::Image::<usize> {
+			extent: image::Extent { x: 8, y: 8 },
+			pixels: vec![
+				0, 0, 0, 0, 0, 0, 0, 0,
+				0, 0, 0, 1, 1, 1, 1, 0,
+				0, 0, 0, 1, 1, 1, 1, 0,
+				0, 1, 1, 0, 0, 1, 1, 0,
+				0, 1, 1, 0, 0, 1, 1, 0,
+				0, 0, 0, 1, 1, 1, 1, 0,
+				0, 0, 0, 1, 1, 1, 1, 0,
+				0, 0, 0, 0, 0, 0, 0, 0,
+			]
+		};
+
+		let contours = find(&img, 1);
+		assert_eq!(contours.len(), 2);
+
+		let contour = &contours[0];
+		assert_eq!(contour.parent, None);
+		assert_eq!(contour.kind, Kind::Outer);
+		assert_eq!(contour.vertices, vec![
+			(3, 1), (6, 1), (6, 6), (3, 6), (3, 5),
+			(2, 4), (1, 4), (1, 3), (2, 3), (3, 2),
+		]);
+		assert_eq!(contour.children, vec![ 1 ]);
+
+		let contour = &contours[1];
+		assert_eq!(contour.parent, Some(0));
+		assert_eq!(contour.kind, Kind::Inner);
+		assert_eq!(contour.vertices, vec![
+			(2, 3), (2, 4), (3, 5), (4, 5),
+			(5, 4), (5, 3), (4, 2), (3, 2),
+		]);
+	}
+
+	#[test]
+	fn contour_passes_through_start() {
+		let img = image::Image::<usize> {
+			extent: image::Extent { x: 5, y: 5 },
+			pixels: vec![
+				0, 0, 0, 0, 0,
+				0, 0, 1, 1, 0,
+				0, 1, 0, 0, 0,
+				0, 1, 0, 0, 0,
+				0, 0, 0, 0, 0,
+			]
+		};
+
+		let contours = find(&img, 1);
+		assert_eq!(contours.len(), 1);
+
+		let contour = &contours[0];
+		assert_eq!(contour.kind, Kind::Outer);
+		assert_eq!(contour.vertices, vec![
+			(2, 1), (3, 1), (2, 1), (1, 2), (1, 3), (1, 2)
+		]);
+	}
+}

src/draw.rs

@@ -0,0 +1,43 @@
+use crate::*;
+
+pub fn line<P: Copy>(
+    x0: usize,
+    y0: usize,
+    x1: usize,
+    y1: usize,
+    img: &mut image::Image<P>,
+    colour: P,
+) {
+    let dx = x0.abs_diff(x1) as isize;
+    let dy = -(y0.abs_diff(y1) as isize);
+
+    let sx = if x0 < x1 { 1 } else { -1 };
+    let sy = if y0 < y1 { 1 } else { -1 };
+
+    let mut x = x0;
+    let mut y = y0;
+    let mut error = dx + dy;
+
+    loop {
+        if x >= img.extent.x || y >= img.extent.y {
+            return;
+        }
+
+        *img.at_mut(x, y) = colour;
+        let e2 = 2 * error;
+        if e2 >= dy {
+            if x == x1 {
+                break;
+            }
+            error += dy;
+            x = (x as isize + sx) as usize;
+        }
+        if e2 <= dx {
+            if y == y1 {
+                break;
+            }
+            error += dx;
+            y = (y as isize + sy) as usize;
+        }
+    }
+}

src/hough.rs

@@ -40,8 +40,8 @@ pub fn lines(
 ) -> 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 mu = stats::mean((&hough.pixels).into_iter().map(|p| *p));
+    let sigma = stats::stddev((&hough.pixels).into_iter().map(|p| *p), mu);
 
     let mut lines = Vec::<Line>::new();
     for y in 0..hough.extent.y {
@@ -59,29 +59,9 @@ pub fn lines(
             }
         }
     }
-	lines
+    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()
+    x as f32 * std::f32::consts::PI / THETA_N as f32
 }

src/lib.rs

@@ -1,4 +1,8 @@
+pub mod colours;
+pub mod contour;
+pub mod draw;
 pub mod filter;
 pub mod hough;
 pub mod image;
 pub mod kernel;
+pub mod stats;

src/stats.rs

@@ -0,0 +1,19 @@
+pub 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
+}
+
+pub 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()
+}