diff --git a/examples/hough.rs b/examples/hough.rs
new file mode 100644
index 0000000..7be8ad0
--- /dev/null
+++ b/examples/hough.rs
@@ -0,0 +1,114 @@
+mod common;
+
+use ccvl::filter::*;
+use ccvl::hough;
+use ccvl::image::*;
+
+const LINE: SRGBA = SRGBA {
+    r: 1.0,
+    g: 0.0,
+    b: 0.0,
+    a: 1.0,
+};
+
+fn main() -> std::process::ExitCode {
+    common::process_image(|mut img| {
+        let greyscale = img.map(|p| Luminance::from(p));
+        let edges = canny(greyscale);
+        for line in hough::lines(10.0, edges) {
+            draw(line, &mut img);
+        }
+        img
+    })
+}
+
+fn in_bounds(q: &(f32, f32), lim: 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
+}
+
+fn endpoints(
+    l: hough::Line,
+    lim: Extent,
+) -> ((usize, usize), (usize, usize)) {
+    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 draw(l: hough::Line, img: &mut Image<SRGBA>) {
+    if l.theta.abs() < 1e-3 {
+        // vertical line
+        let x = l.rho as usize;
+        for y in 0..img.extent.y {
+            if x >= img.extent.x || y >= img.extent.y {
+                return;
+            }
+            *img.at_mut(x, y) = LINE;
+        }
+    } else if (l.theta - std::f32::consts::PI / 2.0) < 1e-3 {
+        // horizontal line
+        let y = l.rho as usize;
+        for x in 0..img.extent.x {
+            if x >= img.extent.x || y >= img.extent.y {
+                return;
+            }
+            *img.at_mut(x, y) = LINE;
+        }
+    } else {
+        bresenham(endpoints(l, img.extent), img);
+    }
+}
+
+fn bresenham(qs: ((usize, usize), (usize, usize)), img: &mut Image<SRGBA>) {
+    let ((x0, y0), (x1, y1)) = qs;
+
+    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) = LINE;
+        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;
+        }
+    }
+}