use image::{self, GenericImageView};
use rayon::prelude::*;
use std::f32::consts::*;
use std::*;
const TAU: f32 = PI * 2.0;
#[inline(always)]
fn clamp<T: PartialOrd>(f: T, lo: T, hi: T) -> T {
debug_assert!(lo < hi);
if f > hi {
hi
} else if f < lo {
lo
} else {
f
}
}
#[derive(Clone)]
pub struct Detection {
edges: Vec<Vec<Edge>>,
}
impl Detection {
pub fn width(&self) -> usize {
self.edges.len()
}
pub fn height(&self) -> usize {
self.edges[0].len()
}
pub fn interpolate(&self, x: f32, y: f32) -> Edge {
let ax = clamp(x.floor() as isize, 0, self.width() as isize - 1) as usize;
let ay = clamp(y.floor() as isize, 0, self.height() as isize - 1) as usize;
let bx = clamp(x.ceil() as isize, 0, self.width() as isize - 1) as usize;
let by = clamp(y.ceil() as isize, 0, self.height() as isize - 1) as usize;
let e1 = self.edges[ax][ay];
let e2 = self.edges[bx][ay];
let e3 = self.edges[ax][by];
let e4 = self.edges[bx][by];
let nx = (x.fract() + 1.0).fract();
let ny = (y.fract() + 1.0).fract();
let x1 = Edge {
magnitude: e1.magnitude * (1.0 - nx) + e2.magnitude * nx,
vec_x: e1.vec_x * (1.0 - nx) + e2.vec_x * nx,
vec_y: e1.vec_y * (1.0 - nx) + e2.vec_y * nx,
};
let x2 = Edge {
magnitude: e3.magnitude * (1.0 - nx) + e4.magnitude * nx,
vec_x: e3.vec_x * (1.0 - nx) + e4.vec_x * nx,
vec_y: e3.vec_y * (1.0 - nx) + e4.vec_y * nx,
};
Edge {
magnitude: x1.magnitude * (1.0 - ny) + x2.magnitude * ny,
vec_x: x1.vec_x * (1.0 - ny) + x2.vec_x * ny,
vec_y: x1.vec_y * (1.0 - ny) + x2.vec_y * ny,
}
}
pub fn as_image(&self) -> image::DynamicImage {
let img = image::RgbImage::from_fn(self.width() as u32, self.height() as u32, |x, y| {
let (h, s, v) = {
let edge = &self[(x as usize, y as usize)];
((edge.angle() + TAU) % TAU, 1.0, edge.magnitude())
};
let (r, g, b) = {
let c = v * s;
let x = c * (1.0 - ((h / FRAC_PI_3) % 2.0 - 1.0).abs());
let m = v - c;
let (r, g, b) = match h {
h if h < FRAC_PI_3 => (c, x, 0.0),
h if h < FRAC_PI_3 * 2.0 => (x, c, 0.0),
h if h < PI => (0.0, c, x),
h if h < PI + FRAC_PI_3 => (0.0, x, c),
h if h < PI + FRAC_PI_3 * 2.0 => (x, 0.0, c),
h if h < TAU => (c, 0.0, x),
_ => unreachable!(),
};
(r + m, g + m, b + m)
};
image::Rgb {
data: [
(r * 255.0).round() as u8,
(g * 255.0).round() as u8,
(b * 255.0).round() as u8,
],
}
});
image::DynamicImage::ImageRgb8(img)
}
}
impl ops::Index<usize> for Detection {
type Output = Edge;
fn index(&self, index: usize) -> &Self::Output {
let x = index % self.width();
let y = index / self.height();
&self.edges[x][y]
}
}
impl ops::Index<(usize, usize)> for Detection {
type Output = Edge;
fn index(&self, index: (usize, usize)) -> &Self::Output {
&self.edges[index.0][index.1]
}
}
#[derive(Copy, Clone, Debug)]
pub struct Edge {
vec_x: f32,
vec_y: f32,
magnitude: f32,
}
impl Edge {
fn new(vec_x: f32, vec_y: f32) -> Edge {
let vec_x = FRAC_1_SQRT_2 * clamp(vec_x, -1.0, 1.0);
let vec_y = FRAC_1_SQRT_2 * clamp(vec_y, -1.0, 1.0);
let magnitude = f32::hypot(vec_x, vec_y);
debug_assert!(0.0 <= magnitude && magnitude <= 1.0);
let frac_1_mag = if magnitude != 0.0 {
magnitude.recip()
} else {
1.0
};
Edge {
vec_x: vec_x * frac_1_mag,
vec_y: vec_y * frac_1_mag,
magnitude,
}
}
pub fn angle(&self) -> f32 {
f32::atan2(self.vec_y, self.vec_x)
}
pub fn dir(&self) -> (f32, f32) {
(self.vec_x * self.magnitude(), self.vec_y * self.magnitude())
}
pub fn dir_norm(&self) -> (f32, f32) {
(self.vec_x, self.vec_y)
}
pub fn magnitude(&self) -> f32 {
self.magnitude
}
}
pub fn canny<T: Into<image::GrayImage>>(
image: T,
sigma: f32,
strong_threshold: f32,
weak_threshold: f32,
) -> Detection {
let gs_image = image.into();
assert!(gs_image.width() > 0);
assert!(gs_image.height() > 0);
let edges = detect_edges(&gs_image, sigma);
let edges = minmax_suppression(&Detection { edges }, weak_threshold);
let edges = hysteresis(&edges, strong_threshold, weak_threshold);
Detection { edges }
}
fn filter_kernel(sigma: f32) -> (usize, Vec<(f32, f32)>) {
let size = (sigma * 10.0).round() as usize;
let mul_2_sigma_2 = 2.0 * sigma.powi(2);
let kernel = (0..size)
.flat_map(|y| {
(0..size).map(move |x| {
let (xf, yf) = (x as f32 - size as f32 / 2.0, y as f32 - size as f32 / 2.0);
let g = (-(xf.powi(2) + yf.powi(2)) / mul_2_sigma_2).exp() / mul_2_sigma_2;
(xf * g, yf * g)
})
})
.collect();
(size, kernel)
}
fn neighbour_pos_delta(theta: f32) -> (i32, i32) {
let neighbours = [
(1, 0),
(1, 1),
(0, 1),
(-1, 1),
(-1, 0),
(-1, -1),
(0, -1),
(1, -1),
];
let n = ((theta + TAU) % TAU) / TAU;
let i = (n * 8.0).round() as usize % 8;
neighbours[i]
}
fn detect_edges(image: &image::GrayImage, sigma: f32) -> Vec<Vec<Edge>> {
let (width, height) = (image.width() as i32, image.height() as i32);
let (ksize, g_kernel) = filter_kernel(sigma);
let ks = ksize as i32;
(0..width)
.into_par_iter()
.map(|g_ix| {
let ix = g_ix;
let kernel = &g_kernel;
(0..height)
.into_par_iter()
.map(move |iy| {
let mut sum_x = 0.0;
let mut sum_y = 0.0;
for kyi in 0..ks {
let ky = kyi - ks / 2;
for kxi in 0..ks {
let kx = kxi - ks / 2;
let k = unsafe {
let i = (kyi * ks + kxi) as usize;
debug_assert!(i < kernel.len());
kernel.get_unchecked(i)
};
let pix = unsafe {
let x = clamp(ix + kx, 0, width - 1);
let y = clamp(iy + ky, 0, height - 1);
f32::from(image.unsafe_get_pixel(x as u32, y as u32).data[0])
};
sum_x += pix * k.0;
sum_y += pix * k.1;
}
}
Edge::new(sum_x / 255.0, sum_y / 255.0)
})
.collect()
})
.collect()
}
fn minmax_suppression(edges: &Detection, weak_threshold: f32) -> Vec<Vec<Edge>> {
let (width, height) = (edges.edges.len(), edges.edges[0].len());
(0..width)
.into_par_iter()
.map(|x| {
(0..height)
.into_par_iter()
.map(|y| {
let edge = edges.edges[x][y];
if edge.magnitude < weak_threshold {
return Edge::new(0.0, 0.0);
}
let truncate = |f: f32| (f * 1e5).round() * 1e-6;
let mut select = 0;
let mut select_flip_bit = 1;
let directions = [1.0, -1.0];
let mut distances = [0i32; 2];
let mut seek_positions = [(x as f32, y as f32); 2];
let mut seek_magnitudes = [truncate(edge.magnitude); 2];
while (distances[0] - distances[1]).abs() <= 1 {
let seek_pos = &mut seek_positions[select];
let seek_magnitude = &mut seek_magnitudes[select];
let direction = directions[select];
seek_pos.0 += edge.dir_norm().0 * direction;
seek_pos.1 += edge.dir_norm().1 * direction;
let interpolated_magnitude =
truncate(edges.interpolate(seek_pos.0, seek_pos.1).magnitude());
let trunc_edge_magnitude = truncate(edge.magnitude);
let end =
interpolated_magnitude < trunc_edge_magnitude
|| *seek_magnitude > trunc_edge_magnitude && interpolated_magnitude < *seek_magnitude;
*seek_magnitude = interpolated_magnitude;
distances[select] += 1;
select ^= select_flip_bit;
if end {
if select_flip_bit == 0 {
break;
}
select_flip_bit = 0;
}
}
let is_apex =
distances[0] == distances[1]
|| (distances[0] - distances[1] == 1 && ((1.0 - edge.vec_x.abs()).abs() < 1e-5 || (1.0 - edge.vec_y.abs()).abs() < 1e-5));
if is_apex {
edge
} else {
Edge::new(0.0, 0.0)
}
})
.collect()
})
.collect()
}
fn hysteresis(
edges: &Vec<Vec<Edge>>,
strong_threshold: f32,
weak_threshold: f32,
) -> Vec<Vec<Edge>> {
assert!(0.0 < strong_threshold && strong_threshold < 1.0);
assert!(0.0 < weak_threshold && weak_threshold < 1.0);
assert!(weak_threshold < strong_threshold);
let (width, height) = (edges.len(), edges.first().unwrap().len());
let mut edges_out: Vec<Vec<Edge>> = vec![vec![Edge::new(0.0, 0.0); height]; width];
for x in 0..width {
for y in 0..height {
if edges[x][y].magnitude < strong_threshold
|| edges_out[x][y].magnitude >= strong_threshold
{
continue;
}
for side in &[0.0, PI] {
let mut current_pos = (x, y);
loop {
let edge = edges[current_pos.0][current_pos.1];
edges_out[current_pos.0][current_pos.1] = edge;
let (nb_pos, nb_magnitude) = [FRAC_PI_4, 0.0, -FRAC_PI_4]
.into_iter()
.map(|bearing| {
neighbour_pos_delta(edge.angle() + FRAC_PI_2 + side + bearing)
})
.filter_map(|(nb_dx, nb_dy)| {
let nb_x = current_pos.0 as i32 + nb_dx;
let nb_y = current_pos.1 as i32 + nb_dy;
if 0 <= nb_x && nb_x < width as i32 && 0 <= nb_y && nb_y < height as i32
{
let nb = (nb_x as usize, nb_y as usize);
Some((nb, edges[nb.0][nb.1].magnitude))
} else {
None
}
})
.fold(((0, 0), 0.0), |(max_pos, max_mag), (pos, mag)| {
if mag > max_mag {
(pos, mag)
} else {
(max_pos, max_mag)
}
});
if nb_magnitude < weak_threshold
|| edges_out[nb_pos.0][nb_pos.1].magnitude > weak_threshold
{
break;
}
current_pos = nb_pos;
}
}
}
}
edges_out
}
#[cfg(test)]
mod tests {
use super::*;
fn edges_to_image(edges: &Vec<Vec<Edge>>) -> image::GrayImage {
let (width, height) = (edges.len(), edges.first().unwrap().len());
let mut image =
image::GrayImage::from_pixel(width as u32, height as u32, image::Luma { data: [0] });
for x in 0..width {
for y in 0..height {
let edge = edges[x][y];
*image.get_pixel_mut(x as u32, y as u32) = image::Luma {
data: [(edge.magnitude * 255.0).round() as u8],
};
}
}
image
}
fn canny_output_stages<T: AsRef<str>>(
path: T,
sigma: f32,
strong_threshold: f32,
weak_threshold: f32,
) -> Detection {
let path = path.as_ref();
let image = image::open(path).unwrap();
let edges = detect_edges(&image.to_luma(), sigma);
let intermediage_d = Detection { edges };
intermediage_d
.as_image()
.save(format!("{}.0-vectors.png", path))
.unwrap();
edges_to_image(&intermediage_d.edges)
.save(format!("{}.1-edges.png", path))
.unwrap();
let edges = minmax_suppression(&intermediage_d, weak_threshold);
edges_to_image(&edges)
.save(format!("{}.2-minmax.png", path))
.unwrap();
let edges = hysteresis(&edges, strong_threshold, weak_threshold);
edges_to_image(&edges)
.save(format!("{}.3-hysteresis.png", path))
.unwrap();
let detection = Detection { edges };
detection
.as_image()
.save(format!("{}.4-result.png", path))
.unwrap();
detection
}
#[test]
fn neighbour_pos_delta_from_theta() {
let neighbours = [
(1, 0),
(1, 1),
(0, 1),
(-1, 1),
(-1, 0),
(-1, -1),
(0, -1),
(1, -1),
];
for nb in neighbours.iter() {
let d = neighbour_pos_delta(f32::atan2(nb.1 as f32, nb.0 as f32));
assert_eq!(*nb, d);
}
}
#[test]
fn edge_new() {
let e = Edge::new(1.0, 0.0);
assert!(1.0 - 1e-6 < e.vec_x && e.vec_x < 1.0 + 1e-6);
assert!(-1e-5 < e.vec_y && e.vec_y < 1e-6);
let e = Edge::new(1.0, 1.0);
assert!(FRAC_1_SQRT_2 - 1e-5 < e.vec_x && e.vec_x < FRAC_1_SQRT_2 + 1e-6);
assert!(FRAC_1_SQRT_2 - 1e-5 < e.vec_y && e.vec_y < FRAC_1_SQRT_2 + 1e-6);
assert!(1.0 - 1e-6 < e.magnitude && e.magnitude < 1.0 + 1e-6);
}
#[test]
fn detection_interpolate() {
let dummy = |magnitude| Edge {
magnitude,
vec_x: 0.0,
vec_y: 0.0,
};
let d = Detection {
edges: vec![vec![dummy(2.0), dummy(8.0)], vec![dummy(4.0), dummy(16.0)]],
};
assert!((d.interpolate(0.0, 0.0).magnitude() - 2.0).abs() <= 1e-6);
assert!((d.interpolate(1.0, 0.0).magnitude() - 4.0).abs() <= 1e-6);
assert!((d.interpolate(0.0, 1.0).magnitude() - 8.0).abs() <= 1e-6);
assert!((d.interpolate(1.0, 1.0).magnitude() - 16.0).abs() <= 1e-6);
assert!((d.interpolate(0.5, 0.0).magnitude() - 3.0).abs() <= 1e-6);
assert!((d.interpolate(0.0, 0.5).magnitude() - 5.0).abs() <= 1e-6);
assert!((d.interpolate(0.5, 1.0).magnitude() - 12.0).abs() <= 1e-6);
assert!((d.interpolate(1.0, 0.5).magnitude() - 10.0).abs() <= 1e-6);
}
#[test]
fn kernel_integral_in_bounds() {
for sigma_i in 1..200 {
let sigma = sigma_i as f32 / 10.0;
let (ksize, kernel) = filter_kernel(sigma);
assert!(ksize.pow(2) == kernel.len());
let mut sum_x = 0.0;
let mut sum_y = 0.0;
for (gx, gy) in kernel {
sum_x += gx;
sum_y += gy;
}
println!(
"sum = ({}, {}), sigma = {}, kernel_size = {}",
sum_x, sum_y, sigma, ksize
);
assert!(-0.0001 < sum_x && sum_x <= 0.0001);
assert!(-0.0001 < sum_y && sum_y <= 0.0001);
}
}
fn detect_vertical_line(detection: &Detection) -> usize {
let mut line_x = None;
for x in 0..detection.width() {
if detection.edges[x][detection.height() / 2].magnitude > 0.5 {
if line_x.is_some() {
panic!("the line is thicker than 1px");
}
line_x = Some(x)
}
}
let line_x = line_x.expect("no line detected");
let middle = detection.width() / 2;
assert!(middle - 1 <= line_x && line_x <= middle);
for y in 0..detection.height() {
let edge = detection.edges[line_x][y];
assert!(edge.magnitude > 0.0);
}
for x in 0..detection.width() {
if x == line_x {
continue;
}
for y in 0..detection.height() {
assert!(detection.edges[x][y].magnitude == 0.0);
}
}
line_x
}
#[test]
fn detect_vertical_line_simple() {
let d = canny_output_stages("testdata/line-simple.png", 1.2, 0.4, 0.05);
let x = detect_vertical_line(&d);
for y in 0..d.height() {
assert!(d.edges[x][y].angle().abs() < 1e-5);
}
}
#[test]
fn detect_vertical_line_fuzzy() {
let d = canny_output_stages("testdata/line-fuzzy.png", 2.0, 0.4, 0.05);
let x = detect_vertical_line(&d);
for y in 0..d.height() {
assert!(d.edges[x][y].angle().abs() < 0.01);
}
}
#[test]
fn detect_vertical_line_weakening() {
let d = canny_output_stages("testdata/line-weakening.png", 1.2, 0.7, 0.05);
detect_vertical_line(&d);
}
}
#[cfg(all(test, feature = "unstable"))]
mod benchmarks {
extern crate test;
use super::*;
static IMG_PATH: &str = "testdata/circle.png";
#[bench]
fn bench_filter_kernel_low_sigma(b: &mut test::Bencher) {
b.iter(|| filter_kernel(1.2));
}
#[bench]
fn bench_filter_kernel_high_sigma(b: &mut test::Bencher) {
b.iter(|| filter_kernel(5.0));
}
#[bench]
fn bench_detect_edges_low_sigma(b: &mut test::Bencher) {
let image = image::open(IMG_PATH).unwrap().to_luma();
b.iter(|| detect_edges(&image, 1.2));
}
#[bench]
fn bench_detect_edges_high_sigma(b: &mut test::Bencher) {
let image = image::open(IMG_PATH).unwrap().to_luma();
b.iter(|| detect_edges(&image, 5.0));
}
#[bench]
fn bench_minmax_suppression_low_sigma(b: &mut test::Bencher) {
let image = image::open(IMG_PATH).unwrap().to_luma();
let edges = Detection {
edges: detect_edges(&image, 1.2),
};
b.iter(|| minmax_suppression(&edges, 0.01));
}
#[bench]
fn bench_minmax_suppression_high_sigma(b: &mut test::Bencher) {
let image = image::open(IMG_PATH).unwrap().to_luma();
let edges = Detection {
edges: detect_edges(&image, 5.0),
};
b.iter(|| minmax_suppression(&edges, 0.01));
}
#[bench]
fn bench_hysteresis(b: &mut test::Bencher) {
let image = image::open(IMG_PATH).unwrap().to_luma();
let edges = Detection {
edges: detect_edges(&image, 1.2),
};
let edges = minmax_suppression(&edges, 0.1);
b.iter(|| hysteresis(&edges, 0.4, 0.1));
}
}