418 lines
17 KiB
Rust
418 lines
17 KiB
Rust
//! Auto-trace: outlines of a scalar field (a layer's alpha, a colour channel or luminance) as
|
||
//! Bezier paths (Layer ▸ Auto-trace).
|
||
//!
|
||
//! 1. **Marching squares** over the samples at pixel centres, with linear interpolation of the
|
||
//! threshold crossing on each cell edge (saddles resolved by the cell's mean), on a field padded
|
||
//! with "outside" so every contour closes.
|
||
//! 2. **Douglas–Peucker** simplification of each closed loop to the Tolerance (pixels).
|
||
//! 3. **Bezier fitting**: each simplified span is fitted with one cubic (least squares with fixed
|
||
//! end tangents, Newton reparameterisation; [`crate::fit`]). The tangent at each vertex is the
|
||
//! chord direction at Corner Roundness 0 (straight segments, sharp corners) and the smooth
|
||
//! (Catmull–Rom) direction at 100 %, blended between.
|
||
//!
|
||
//! Loops smaller than Minimum Area are dropped. Each loop gets its nesting depth (0 = outermost),
|
||
//! so callers can add outer outlines and subtract holes.
|
||
|
||
use std::collections::HashMap;
|
||
|
||
use effectcraft_keyframe::ShapePath;
|
||
use kurbo::{Point, Vec2};
|
||
|
||
use crate::fit::{fit_single, unit};
|
||
|
||
/// Auto-trace options.
|
||
#[derive(Clone, Copy, Debug, PartialEq)]
|
||
pub struct TraceOpts {
|
||
/// Inside where the field is at least this (0…1).
|
||
pub threshold: f32,
|
||
/// Douglas–Peucker tolerance in pixels.
|
||
pub tolerance: f64,
|
||
/// Loops with a smaller area (square pixels) are dropped.
|
||
pub min_area: f64,
|
||
/// 0 = polygons, 1 = fully smooth.
|
||
pub corner_roundness: f64,
|
||
}
|
||
|
||
impl Default for TraceOpts {
|
||
fn default() -> Self {
|
||
TraceOpts { threshold: 0.5, tolerance: 1.0, min_area: 10.0, corner_roundness: 0.5 }
|
||
}
|
||
}
|
||
|
||
/// A traced outline.
|
||
#[derive(Clone, Debug, PartialEq)]
|
||
pub struct Traced {
|
||
pub path: ShapePath,
|
||
/// Number of other outlines that contain this one (even = filled, odd = hole).
|
||
pub depth: usize,
|
||
/// Signed area of the simplified polygon (pixels²).
|
||
pub area: f64,
|
||
}
|
||
|
||
/// Trace `field` (`w × h`, row-major, sample (x, y) at pixel centre (x + 0.5, y + 0.5)).
|
||
pub fn trace(field: &[f32], w: usize, h: usize, o: &TraceOpts) -> Vec<Traced> {
|
||
if w == 0 || h == 0 || field.len() < w * h {
|
||
return vec![];
|
||
}
|
||
let loops = contours(field, w, h, o.threshold);
|
||
let mut polys: Vec<(Vec<[f64; 2]>, Vec<[f64; 2]>, f64)> = loops
|
||
.into_iter()
|
||
.filter_map(|raw| {
|
||
let a = signed_area(&raw);
|
||
if a.abs() < o.min_area.max(0.0) || raw.len() < 3 {
|
||
return None;
|
||
}
|
||
let simp = simplify_closed(&raw, o.tolerance.max(0.01));
|
||
if simp.len() < 3 {
|
||
return None;
|
||
}
|
||
Some((raw, simp, a))
|
||
})
|
||
.collect();
|
||
// Depth: how many other loops contain a point of this loop.
|
||
let depths: Vec<usize> = (0..polys.len())
|
||
.map(|i| {
|
||
let p = polys[i].1[0];
|
||
(0..polys.len()).filter(|&j| j != i && polys[j].2.abs() > polys[i].2.abs() && point_in_poly(p, &polys[j].1)).count()
|
||
})
|
||
.collect();
|
||
let mut out: Vec<Traced> = polys
|
||
.drain(..)
|
||
.zip(depths)
|
||
.map(|((raw, simp, a), depth)| Traced { path: fit_loop(&raw, &simp, o.corner_roundness.clamp(0.0, 1.0)), depth, area: a })
|
||
.collect();
|
||
// Outer outlines first, then holes, islands…; larger first within a depth.
|
||
out.sort_by(|a, b| a.depth.cmp(&b.depth).then(b.area.abs().total_cmp(&a.area.abs())));
|
||
out
|
||
}
|
||
|
||
/// Marching squares: closed polylines (pixel coordinates).
|
||
fn contours(field: &[f32], w: usize, h: usize, thr: f32) -> Vec<Vec<[f64; 2]>> {
|
||
// Padded grid: sample (i, j) of the padded grid is pixel (i - 1, j - 1), outside = -inf.
|
||
let (gw, gh) = (w + 2, h + 2);
|
||
let val = |i: usize, j: usize| -> f32 { if i == 0 || j == 0 || i > w || j > h { f32::NEG_INFINITY } else { field[(j - 1) * w + (i - 1)] } };
|
||
let inside = |i: usize, j: usize| val(i, j) >= thr;
|
||
// Edge ids: horizontal edge (i, j)-(i+1, j) = 2 * (j * gw + i); vertical (i, j)-(i, j+1) = +1.
|
||
let pos = |i: usize, j: usize| [i as f64 - 0.5, j as f64 - 0.5];
|
||
let cross = |a: (usize, usize), b: (usize, usize)| -> [f64; 2] {
|
||
let (va, vb) = (val(a.0, a.1), val(b.0, b.1));
|
||
let t = if va.is_finite() && vb.is_finite() && (vb - va).abs() > 1e-12 {
|
||
((thr - va) / (vb - va)).clamp(0.0, 1.0) as f64
|
||
} else if va.is_finite() {
|
||
// b is padding: the crossing sits half a pixel out.
|
||
if va >= thr { 0.5 } else { 0.0 }
|
||
} else if vb >= thr {
|
||
0.5
|
||
} else {
|
||
1.0
|
||
};
|
||
let (pa, pb) = (pos(a.0, a.1), pos(b.0, b.1));
|
||
[pa[0] + (pb[0] - pa[0]) * t, pa[1] + (pb[1] - pa[1]) * t]
|
||
};
|
||
let hid = |i: usize, j: usize| 2 * (j * gw + i);
|
||
let vid = |i: usize, j: usize| 2 * (j * gw + i) + 1;
|
||
// Directed segments (from edge, to edge) keeping the inside on the left.
|
||
let mut next: HashMap<usize, usize> = HashMap::new();
|
||
let mut point: HashMap<usize, [f64; 2]> = HashMap::new();
|
||
for j in 0..gh - 1 {
|
||
for i in 0..gw - 1 {
|
||
// Corners: tl (i, j), tr (i+1, j), br (i+1, j+1), bl (i, j+1).
|
||
let (tl, tr, br, bl) = (inside(i, j), inside(i + 1, j), inside(i + 1, j + 1), inside(i, j + 1));
|
||
let code = (tl as u8) << 3 | (tr as u8) << 2 | (br as u8) << 1 | bl as u8;
|
||
if code == 0 || code == 15 {
|
||
continue;
|
||
}
|
||
let top = hid(i, j);
|
||
let bottom = hid(i, j + 1);
|
||
let left = vid(i, j);
|
||
let right = vid(i + 1, j);
|
||
let mut seg = |a: usize, b: usize| {
|
||
next.insert(a, b);
|
||
};
|
||
// Walk with the inside on the right in image coordinates (y down), which is
|
||
// counter-clockwise on screen for outer loops… orientation is fixed afterwards.
|
||
match code {
|
||
1 => seg(left, bottom),
|
||
2 => seg(bottom, right),
|
||
3 => seg(left, right),
|
||
4 => seg(right, top),
|
||
5 => {
|
||
// tr + bl inside (saddle).
|
||
let centre = (val(i, j) + val(i + 1, j) + val(i + 1, j + 1) + val(i, j + 1)) / 4.0;
|
||
if centre >= thr {
|
||
seg(left, top);
|
||
seg(right, bottom);
|
||
} else {
|
||
seg(left, bottom);
|
||
seg(right, top);
|
||
}
|
||
}
|
||
6 => seg(bottom, top),
|
||
7 => seg(left, top),
|
||
8 => seg(top, left),
|
||
9 => seg(top, bottom),
|
||
10 => {
|
||
// tl + br inside (saddle).
|
||
let centre = (val(i, j) + val(i + 1, j) + val(i + 1, j + 1) + val(i, j + 1)) / 4.0;
|
||
if centre >= thr {
|
||
seg(top, right);
|
||
seg(bottom, left);
|
||
} else {
|
||
seg(top, left);
|
||
seg(bottom, right);
|
||
}
|
||
}
|
||
11 => seg(top, right),
|
||
12 => seg(right, left),
|
||
13 => seg(right, bottom),
|
||
14 => seg(bottom, left),
|
||
_ => {}
|
||
}
|
||
for (e, a, b) in [(top, (i, j), (i + 1, j)), (bottom, (i, j + 1), (i + 1, j + 1)), (left, (i, j), (i, j + 1)), (right, (i + 1, j), (i + 1, j + 1))]
|
||
{
|
||
if inside(a.0, a.1) != inside(b.0, b.1) {
|
||
point.entry(e).or_insert_with(|| cross(a, b));
|
||
}
|
||
}
|
||
}
|
||
}
|
||
let mut out = vec![];
|
||
let mut keys: Vec<usize> = next.keys().copied().collect();
|
||
keys.sort_unstable();
|
||
for start in keys {
|
||
if !next.contains_key(&start) {
|
||
continue;
|
||
}
|
||
let mut poly = vec![];
|
||
let mut e = start;
|
||
while let Some(n) = next.remove(&e) {
|
||
if let Some(p) = point.get(&e) {
|
||
poly.push(*p);
|
||
}
|
||
e = n;
|
||
if e == start {
|
||
break;
|
||
}
|
||
}
|
||
if poly.len() >= 3 {
|
||
out.push(poly);
|
||
}
|
||
}
|
||
out
|
||
}
|
||
|
||
/// Shoelace area (positive = clockwise on screen with y down).
|
||
pub fn signed_area(p: &[[f64; 2]]) -> f64 {
|
||
let n = p.len();
|
||
let mut a = 0.0;
|
||
for i in 0..n {
|
||
let (q, r) = (p[i], p[(i + 1) % n]);
|
||
a += q[0] * r[1] - r[0] * q[1];
|
||
}
|
||
a * 0.5
|
||
}
|
||
|
||
fn point_in_poly(p: [f64; 2], poly: &[[f64; 2]]) -> bool {
|
||
let mut inside = false;
|
||
let n = poly.len();
|
||
let mut j = n - 1;
|
||
for i in 0..n {
|
||
let (a, b) = (poly[i], poly[j]);
|
||
if (a[1] > p[1]) != (b[1] > p[1]) && p[0] < (b[0] - a[0]) * (p[1] - a[1]) / (b[1] - a[1]) + a[0] {
|
||
inside = !inside;
|
||
}
|
||
j = i;
|
||
}
|
||
inside
|
||
}
|
||
|
||
fn seg_dist(p: [f64; 2], a: [f64; 2], b: [f64; 2]) -> f64 {
|
||
let (dx, dy) = (b[0] - a[0], b[1] - a[1]);
|
||
let l2 = dx * dx + dy * dy;
|
||
let t = if l2 > 0.0 { (((p[0] - a[0]) * dx + (p[1] - a[1]) * dy) / l2).clamp(0.0, 1.0) } else { 0.0 };
|
||
((p[0] - a[0] - t * dx).powi(2) + (p[1] - a[1] - t * dy).powi(2)).sqrt()
|
||
}
|
||
|
||
/// Douglas–Peucker on an open polyline: kept indices.
|
||
fn dp(p: &[[f64; 2]], lo: usize, hi: usize, tol: f64, keep: &mut Vec<usize>) {
|
||
if hi <= lo + 1 {
|
||
return;
|
||
}
|
||
let (mut best, mut bi) = (0.0, lo);
|
||
for (k, q) in p.iter().enumerate().take(hi).skip(lo + 1) {
|
||
let d = seg_dist(*q, p[lo], p[hi]);
|
||
if d > best {
|
||
best = d;
|
||
bi = k;
|
||
}
|
||
}
|
||
if best > tol {
|
||
dp(p, lo, bi, tol, keep);
|
||
keep.push(bi);
|
||
dp(p, bi, hi, tol, keep);
|
||
}
|
||
}
|
||
|
||
/// Douglas–Peucker on a closed loop: indices into `p` of the kept vertices, in order.
|
||
fn simplify_idx(p: &[[f64; 2]], tol: f64) -> Vec<usize> {
|
||
let n = p.len();
|
||
// Split at the two mutually farthest-ish points (first point and the point farthest from it).
|
||
let far = (0..n).max_by(|&a, &b| {
|
||
let d = |i: usize| (p[i][0] - p[0][0]).powi(2) + (p[i][1] - p[0][1]).powi(2);
|
||
d(a).total_cmp(&d(b))
|
||
});
|
||
let far = far.unwrap_or(0);
|
||
let mut ext: Vec<[f64; 2]> = p.to_vec();
|
||
ext.push(p[0]);
|
||
let mut keep = vec![0];
|
||
dp(&ext, 0, far, tol, &mut keep);
|
||
keep.push(far);
|
||
dp(&ext, far, n, tol, &mut keep);
|
||
keep.sort_unstable();
|
||
keep.dedup();
|
||
keep.retain(|&i| i < n);
|
||
keep
|
||
}
|
||
|
||
fn simplify_closed(p: &[[f64; 2]], tol: f64) -> Vec<[f64; 2]> {
|
||
simplify_idx(p, tol).into_iter().map(|i| p[i]).collect()
|
||
}
|
||
|
||
/// Fit cubics through the simplified vertices of a loop (`raw` holds the full contour).
|
||
fn fit_loop(raw: &[[f64; 2]], simp: &[[f64; 2]], round: f64) -> ShapePath {
|
||
let n = simp.len();
|
||
if round <= 1e-6 {
|
||
return ShapePath::polygon(simp, true);
|
||
}
|
||
// Raw index of each simplified vertex (they are raw points).
|
||
let idx: Vec<usize> = simp.iter().map(|s| raw.iter().position(|r| r == s).unwrap_or(0)).collect();
|
||
let pt = |a: [f64; 2]| Point::new(a[0], a[1]);
|
||
let v = |a: [f64; 2], b: [f64; 2]| Vec2::new(b[0] - a[0], b[1] - a[1]);
|
||
// Per vertex: (direction leaving, direction arriving), blended between chord and smooth.
|
||
let dirs: Vec<(Vec2, Vec2)> = (0..n)
|
||
.map(|i| {
|
||
let (prev, cur, nxt) = (simp[(i + n - 1) % n], simp[i], simp[(i + 1) % n]);
|
||
let smooth = unit(v(prev, nxt)).unwrap_or(Vec2::new(1.0, 0.0));
|
||
let out_chord = unit(v(cur, nxt)).unwrap_or(smooth);
|
||
let in_chord = unit(v(prev, cur)).unwrap_or(smooth);
|
||
let leave = unit(out_chord * (1.0 - round) + smooth * round).unwrap_or(out_chord);
|
||
let arrive = unit(in_chord * (1.0 - round) + smooth * round).unwrap_or(in_chord);
|
||
(leave, arrive)
|
||
})
|
||
.collect();
|
||
let mut path = ShapePath::polygon(simp, true);
|
||
for i in 0..n {
|
||
let j = (i + 1) % n;
|
||
let (a, b) = (idx[i], idx[j]);
|
||
let mut pts: Vec<Point> =
|
||
if b > a { raw[a..=b].iter().map(|q| pt(*q)).collect() } else { raw[a..].iter().chain(raw[..=b].iter()).map(|q| pt(*q)).collect() };
|
||
if pts.len() < 2 {
|
||
pts = vec![pt(simp[i]), pt(simp[j])];
|
||
}
|
||
let chord = pts[0].distance(pts[pts.len() - 1]);
|
||
let (t0, t1) = (dirs[i].0, dirs[j].1);
|
||
let (c, _, _) = if pts.len() >= 3 { fit_single(&pts, t0, t1) } else { (kurbo::CubicBez::new(pts[0], pts[0], pts[1], pts[1]), 0.0, 0) };
|
||
// Keep handles pointing along the tangents and within the chord.
|
||
let max = chord * 0.5 * round.max(0.34);
|
||
let l0 = (c.p1 - c.p0).dot(t0).clamp(0.0, max);
|
||
let l1 = (c.p3 - c.p2).dot(t1).clamp(0.0, max);
|
||
let (l0, l1) = if pts.len() < 3 { (chord / 3.0 * round, chord / 3.0 * round) } else { (l0, l1) };
|
||
path.out_tangents[i] = [t0.x * l0, t0.y * l0];
|
||
path.in_tangents[j] = [-t1.x * l1, -t1.y * l1];
|
||
}
|
||
path
|
||
}
|
||
|
||
#[cfg(test)]
|
||
mod tests {
|
||
use super::*;
|
||
use crate::{FillRule, fill_coverage, to_kurbo};
|
||
use effectcraft_geom::Mat3;
|
||
|
||
fn iou(a: &[bool], b: &[f32]) -> f64 {
|
||
let (mut i, mut u) = (0.0, 0.0);
|
||
for (x, y) in a.iter().zip(b) {
|
||
let y = *y >= 0.5;
|
||
if *x && y {
|
||
i += 1.0;
|
||
}
|
||
if *x || y {
|
||
u += 1.0;
|
||
}
|
||
}
|
||
i / u
|
||
}
|
||
|
||
/// Rasterise traced outlines with even-odd fill (as the Add / Subtract mask stack does).
|
||
fn raster(t: &[Traced], w: usize, h: usize) -> Vec<f32> {
|
||
let paths: Vec<_> = t.iter().map(|t| to_kurbo(&t.path)).collect();
|
||
fill_coverage(&paths, &Mat3::IDENTITY, w as u32, h as u32, FillRule::EvenOdd).data
|
||
}
|
||
|
||
#[test]
|
||
fn traces_a_disc_and_a_square_with_a_hole() {
|
||
let (w, h) = (120, 80);
|
||
let mut field = vec![0.0f32; w * h];
|
||
let mut truth = vec![false; w * h];
|
||
for y in 0..h {
|
||
for x in 0..w {
|
||
let (fx, fy) = (x as f64 + 0.5, y as f64 + 0.5);
|
||
let disc = (fx - 30.0).powi(2) + (fy - 40.0).powi(2) <= 20.0f64.powi(2);
|
||
let square = (60.0..110.0).contains(&fx) && (15.0..65.0).contains(&fy);
|
||
let hole = (75.0..95.0).contains(&fx) && (30.0..50.0).contains(&fy);
|
||
let on = disc || (square && !hole);
|
||
truth[y * w + x] = on;
|
||
field[y * w + x] = if on { 1.0 } else { 0.0 };
|
||
}
|
||
}
|
||
for round in [0.0, 0.5, 1.0] {
|
||
let t = trace(&field, w, h, &TraceOpts { corner_roundness: round, ..Default::default() });
|
||
assert_eq!(t.len(), 3, "disc, square, hole");
|
||
assert_eq!(t.iter().filter(|t| t.depth == 1).count(), 1, "one hole");
|
||
let r = raster(&t, w, h);
|
||
let score = iou(&truth, &r);
|
||
// Full roundness rounds the square's corners off as well.
|
||
let min = if round >= 1.0 { 0.88 } else { 0.96 };
|
||
assert!(score > min, "IoU {score} at roundness {round}");
|
||
}
|
||
// A disc stays a disc when smooth.
|
||
let disc: Vec<f32> =
|
||
(0..w * h).map(|i| if ((i % w) as f64 + 0.5 - 30.0).powi(2) + ((i / w) as f64 + 0.5 - 40.0).powi(2) <= 400.0 { 1.0 } else { 0.0 }).collect();
|
||
let t = trace(&disc, w, h, &TraceOpts { corner_roundness: 1.0, tolerance: 1.5, ..Default::default() });
|
||
let truth_disc: Vec<bool> = disc.iter().map(|v| *v > 0.5).collect();
|
||
let score = iou(&truth_disc, &raster(&t, w, h));
|
||
assert!(score > 0.97, "disc IoU {score}");
|
||
assert!(t[0].path.vertices.len() < 40, "simplified to {}", t[0].path.vertices.len());
|
||
// The square traces to 4 corners at zero roundness.
|
||
let t = trace(&field, w, h, &TraceOpts { corner_roundness: 0.0, tolerance: 1.0, ..Default::default() });
|
||
let sq = t.iter().find(|t| t.depth == 0 && t.path.vertices.iter().all(|v| v[0] > 55.0)).expect("square");
|
||
assert!(sq.path.vertices.len() <= 8, "{:?}", sq.path.vertices);
|
||
}
|
||
|
||
#[test]
|
||
fn min_area_and_threshold() {
|
||
let (w, h) = (40, 40);
|
||
let mut field = vec![0.0f32; w * h];
|
||
// A 2×2 speck and a 20×20 block at 0.6.
|
||
for y in 2..4 {
|
||
for x in 2..4 {
|
||
field[y * w + x] = 1.0;
|
||
}
|
||
}
|
||
for y in 10..30 {
|
||
for x in 10..30 {
|
||
field[y * w + x] = 0.6;
|
||
}
|
||
}
|
||
assert_eq!(trace(&field, w, h, &TraceOpts { min_area: 10.0, ..Default::default() }).len(), 1);
|
||
assert_eq!(trace(&field, w, h, &TraceOpts { min_area: 1.0, ..Default::default() }).len(), 2);
|
||
// Above 0.6 only the speck remains (and it is too small).
|
||
assert!(trace(&field, w, h, &TraceOpts { threshold: 0.7, ..Default::default() }).is_empty());
|
||
// Edge-touching shapes close along the border.
|
||
let full = vec![1.0f32; 100];
|
||
let t = trace(&full, 10, 10, &TraceOpts::default());
|
||
assert_eq!(t.len(), 1);
|
||
assert!((t[0].area.abs() - 100.0).abs() < 2.0, "{}", t[0].area);
|
||
}
|
||
}
|