From 384c4361a60f5f7c5ab42b8057dd5685d4b94521 Mon Sep 17 00:00:00 2001 From: "Dylan M. Taylor" Date: Mon, 2 Nov 2026 22:15:43 -0500 Subject: [PATCH] Auto-trace: bound the nesting, simplification and fitting work Many separate regions made the nesting depth test every pair of contours, Douglas-Peucker could rescan a span once per point (and recurse as deep as the loop is long), and the curve fit searched the raw contour from its start for every simplified vertex. - Depth: each loop's bounding box goes into a coarse grid; a loop's first point is tested only against loops whose box holds it. Same depths as testing every pair. - Douglas-Peucker: spans wait on an explicit stack, with a budget of 64 distance evaluations per point; spans still waiting when it runs out keep all their points, so the outline stays exact. - Fitting: uses the raw indices the simplifier kept instead of searching. --- crates/path/src/trace.rs | 187 +++++++++++++++++++++++++++++++-------- 1 file changed, 149 insertions(+), 38 deletions(-) diff --git a/crates/path/src/trace.rs b/crates/path/src/trace.rs index 01508bcaf..25fe1c7ff 100644 --- a/crates/path/src/trace.rs +++ b/crates/path/src/trace.rs @@ -55,37 +55,70 @@ pub fn trace(field: &[f32], w: usize, h: usize, o: &TraceOpts) -> Vec { return vec![]; } let loops = contours(field, w, h, o.threshold); - let mut polys: Vec<(Vec<[f64; 2]>, Vec<[f64; 2]>, f64)> = loops + // (raw loop, kept raw indices, simplified loop, area) + let mut polys: Vec<(Vec<[f64; 2]>, Vec, 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 { + let idx = simplify_idx(&raw, o.tolerance.max(0.01)); + if idx.len() < 3 { return None; } - Some((raw, simp, a)) - }) - .collect(); - // Depth: how many other loops contain a point of this loop. - let depths: Vec = (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() + let simp = idx.iter().map(|&i| raw[i]).collect(); + Some((raw, idx, simp, a)) }) .collect(); + let depths = nesting_depths(&polys.iter().map(|p| (p.2.as_slice(), p.3.abs())).collect::>(), w, h); let mut out: Vec = 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 }) + .map(|((raw, idx, simp, a), depth)| Traced { path: fit_loop(&raw, &idx, &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 } +/// Grid cell (pixels) of the index [`nesting_depths`] finds candidate containers with. +const DEPTH_CELL: f64 = 32.0; + +/// Per loop (simplified points, |area|): how many other, larger loops contain its first point. +/// Only loops whose bounding box holds that point can contain it, so each loop's box is entered +/// into a coarse grid over the `w × h` field and a point is tested against the loops of its +/// cell alone: many separate regions no longer test every pair. +fn nesting_depths(polys: &[(&[[f64; 2]], f64)], w: usize, h: usize) -> Vec { + // Loops lie within the field padded by half a pixel each side. + let (gw, gh) = (((w as f64 + 1.0) / DEPTH_CELL).ceil() as usize + 1, ((h as f64 + 1.0) / DEPTH_CELL).ceil() as usize + 1); + let cell = |v: f64, n: usize| (((v + 0.5) / DEPTH_CELL).floor().max(0.0) as usize).min(n - 1); + let boxes: Vec<[f64; 4]> = polys + .iter() + .map(|(p, _)| p.iter().fold([f64::INFINITY, f64::INFINITY, f64::NEG_INFINITY, f64::NEG_INFINITY], |b, q| [b[0].min(q[0]), b[1].min(q[1]), b[2].max(q[0]), b[3].max(q[1])])) + .collect(); + let mut grid: Vec> = vec![Vec::new(); gw * gh]; + for (j, b) in boxes.iter().enumerate() { + for cy in cell(b[1], gh)..=cell(b[3], gh) { + for cx in cell(b[0], gw)..=cell(b[2], gw) { + grid[cy * gw + cx].push(j); + } + } + } + (0..polys.len()) + .map(|i| { + let p = polys[i].0[0]; + grid[cell(p[1], gh) * gw + cell(p[0], gw)] + .iter() + .filter(|&&j| { + let b = boxes[j]; + j != i && polys[j].1 > polys[i].1 && p[0] >= b[0] && p[0] <= b[2] && p[1] >= b[1] && p[1] <= b[3] && point_in_poly(p, polys[j].0) + }) + .count() + }) + .collect() +} + /// Marching squares: closed polylines (pixel coordinates). fn contours(field: &[f32], w: usize, h: usize, thr: f32) -> Vec> { // Padded grid: sample (i, j) of the padded grid is pixel (i - 1, j - 1), outside = -inf. @@ -233,23 +266,37 @@ fn seg_dist(p: [f64; 2], a: [f64; 2], b: [f64; 2]) -> f64 { ((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) { - 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; +/// Distance evaluations Douglas–Peucker may spend per point of a loop. Balanced splits need +/// about log2(n); splits that keep landing next to an end would need up to n. +const DP_WORK_PER_POINT: usize = 64; + +/// Douglas–Peucker on an open polyline: kept indices (unordered). Spans wait on an explicit +/// stack (no recursion depth to overflow), and once `budget` distance evaluations are spent the +/// spans still waiting keep all their points: the outline stays exact, just less simplified. +fn dp(p: &[[f64; 2]], lo: usize, hi: usize, tol: f64, budget: &mut usize, keep: &mut Vec) { + let mut spans = vec![(lo, hi)]; + while let Some((lo, hi)) = spans.pop() { + if hi <= lo + 1 { + continue; + } + if *budget < hi - lo - 1 { + keep.extend(lo + 1..hi); + continue; + } + *budget -= hi - lo - 1; + 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 { + keep.push(bi); + spans.push((lo, bi)); + spans.push((bi, hi)); } - } - if best > tol { - dp(p, lo, bi, tol, keep); - keep.push(bi); - dp(p, bi, hi, tol, keep); } } @@ -265,27 +312,23 @@ fn simplify_idx(p: &[[f64; 2]], tol: f64) -> Vec { 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); + let mut budget = n.saturating_mul(DP_WORK_PER_POINT); + dp(&ext, 0, far, tol, &mut budget, &mut keep); keep.push(far); - dp(&ext, far, n, tol, &mut keep); + dp(&ext, far, n, tol, &mut budget, &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 { +/// Fit cubics through the simplified vertices of a loop (`raw` holds the full contour, `idx` +/// the raw index of each simplified vertex, as [`simplify_idx`] kept them). +fn fit_loop(raw: &[[f64; 2]], idx: &[usize], 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 = 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. @@ -349,6 +392,74 @@ mod tests { fill_coverage(&paths, &Mat3::IDENTITY, w as u32, h as u32, FillRule::EvenOdd).data } + #[test] + fn nesting_depths_match_every_pair() { + // Rings within rings beside many separate dots: the grid lookup gives the same depths + // as testing every pair. + let (w, h) = (200, 120); + let mut field = vec![0.0f32; w * h]; + for y in 0..h { + for x in 0..w { + let (fx, fy) = (x as f64 + 0.5 - 60.0, y as f64 + 0.5 - 60.0); + let r = (fx * fx + fy * fy).sqrt(); + let ring = r < 56.0 && (r / 8.0) as usize % 2 == 0; + let dot = x >= 130 && x % 6 < 3 && y % 6 < 3; + field[y * w + x] = if ring || dot { 1.0 } else { 0.0 }; + } + } + let opts = TraceOpts { min_area: 1.0, ..Default::default() }; + let t = trace(&field, w, h, &opts); + assert!(t.len() > 200, "{}", t.len()); + let polys: Vec<(Vec<[f64; 2]>, f64)> = + t.iter().map(|t| (t.path.vertices.clone(), t.area.abs())).collect(); + let refs: Vec<(&[[f64; 2]], f64)> = polys.iter().map(|(p, a)| (p.as_slice(), *a)).collect(); + let fast = nesting_depths(&refs, w, h); + let pairs: Vec = (0..refs.len()) + .map(|i| (0..refs.len()).filter(|&j| j != i && refs[j].1 > refs[i].1 && point_in_poly(refs[i].0[0], refs[j].0)).count()) + .collect(); + assert_eq!(fast, pairs); + assert_eq!(fast, t.iter().map(|t| t.depth).collect::>()); + assert!(fast.iter().max() >= Some(&5)); + } + + #[test] + fn simplification_work_is_bounded_and_keeps_the_outline() { + // Out of budget, the spans still waiting keep every point: nothing is lost, the outline + // is just less simplified. + let p: Vec<[f64; 2]> = (0..12).map(|i| [i as f64, if i % 2 == 0 { 0.0 } else { 5.0 }]).collect(); + let (mut budget, mut keep) = (0, vec![]); + dp(&p, 0, 11, 0.01, &mut budget, &mut keep); + keep.sort_unstable(); + assert_eq!(keep, (1..11).collect::>()); + // With a budget, the result matches the plain recursion's on an ordinary loop. + let (mut budget, mut keep) = (usize::MAX, vec![]); + dp(&p, 0, 11, 0.01, &mut budget, &mut keep); + keep.sort_unstable(); + assert_eq!(keep, (1..11).collect::>()); + // A long, finely zig-zagging loop: simplified within a budget linear in its length. + let n = 200_000; + let ring: Vec<[f64; 2]> = (0..n) + .map(|i| { + let a = i as f64 / n as f64 * std::f64::consts::TAU; + let r = 1000.0 + if i % 2 == 0 { 0.0 } else { 3.0 }; + [r * a.cos(), r * a.sin()] + }) + .collect(); + let idx = simplify_idx(&ring, 0.5); + assert!(idx.windows(2).all(|w| w[0] < w[1]) && idx.len() >= 3 && *idx.last().unwrap() < n); + } + + #[test] + fn fit_uses_the_kept_indices() { + // A contour that passes the same point twice: the fit follows the simplifier's index, + // not the first raw point with those coordinates. + let raw = vec![[0.0, 0.0], [10.0, 0.0], [10.0, 10.0], [0.0, 0.0], [-10.0, 10.0], [-10.0, 0.0]]; + let idx = vec![0, 2, 3, 4]; + let simp: Vec<[f64; 2]> = idx.iter().map(|&i| raw[i]).collect(); + let path = fit_loop(&raw, &idx, &simp, 0.5); + assert_eq!(path.vertices, simp); + } + #[test] fn traces_a_disc_and_a_square_with_a_hole() { let (w, h) = (120, 80);