Puppet Warp: the mesh follows the object outline (#3105) (#3117)

The mesh was every square grid cell touching the layer's opaque pixels:
a staircase up to a cell outside the shape, and shapes closer than about
a cell shared cells, so pinning one dragged the other along.

build_mesh now traces the outline of the opaque region grown (or shrunk)
by Expansion from a signed distance field, simplifies it (Douglas-Peucker,
corners kept), splits its sides to the Density spacing, fills the inside
with a triangular lattice and triangulates the points (Bowyer-Watson
Delaunay, outline sides recovered by midpoint splits), keeping the
triangles inside the outline. Separate shapes get separate meshes, so an
unpinned one stays put, as in Photoshop. Vertices are numbered in reverse
Cuthill-McKee order to keep the solver's banded factor narrow, and the
density's spacing counts the grown region so a large Expansion on a small
layer no longer builds tens of thousands of vertices.

Scorecard TOOL-214-6 done; docs/gaps.md G7 updated.

Co-authored-by: Claude Opus 5.5 <[email protected]>
Co-authored-by: Brandon Thomas <[email protected]>

crates/algo/src/puppet.rs

@@ -1,10 +1,12 @@
 //! Edit › Puppet Warp: pins on a triangle mesh over the layer's opaque region, deformed
 //! as-rigidly-as-possible.
 //!
-//! * **Mesh.** The layer's coverage (alpha above a threshold, grown or shrunk by `expansion`) is
-//!   rasterized onto a coarse block grid; a square grid at the density's spacing keeps every
-//!   cell that touches covered blocks, split into two triangles (alternating diagonals). A
-//!   grid-clipped mesh is simple, robust and uniform, which is all ARAP needs.
+//! * **Mesh.** The layer's coverage (alpha above a threshold) becomes a signed distance field,
+//!   offset by `expansion` to grow or shrink the region. A square grid at the density's spacing,
+//!   split into triangles, is fitted to that outline (2D isosurface stuffing): vertices near the
+//!   outline snap onto it and the triangles it crosses are cut along it, so the mesh edge follows
+//!   the shape, and separate shapes get separate meshes. Inside, the mesh stays a uniform grid,
+//!   which is all ARAP needs.
 //! * **Deformation.** The per-triangle ARAP energy of Liu et al. 2008 (the triangle form of
 //!   Sorkine & Alexa 2007), `E = Σ_t Σ_(i,j)∈t w_ij |(v_i − v_j) − R_t (x_i − x_j)|²` with
 //!   cotangent weights, minimized by alternating a local step (best rotation per triangle from
@@ -18,6 +20,8 @@
 //! * **Rendering.** Triangles are rasterized with [`crate::warp::warp_triangles`] (tile-parallel,
 //!   premultiplied, bilinear/bicubic), back to front by the depth of their nearest pin.
 
+use std::collections::{HashMap, HashSet};
+
 use photocraft_geom::Rect;
 use photocraft_raster::Surface;
 use serde::{Deserialize, Serialize};
@@ -125,120 +129,619 @@ pub struct PuppetMesh {
 /// Coverage threshold: alpha above this counts as part of the puppet.
 const ALPHA_MIN: f32 = 0.05;
 
+/// The outline field's longest side, in blocks: how finely the mesh edge follows the outline.
+const FIELD_BLOCKS: f64 = 1024.0;
+
+/// The puppet's region as a signed distance field over a block grid: negative inside the
+/// opaque area grown (or shrunk) by the expansion, positive outside, in px.
+struct Field {
+    x0: f64,
+    y0: f64,
+    /// Block size (px).
+    g: f64,
+    bw: usize,
+    bh: usize,
+    /// Signed distance from each block centre to the opaque outline (px, negative inside).
+    sd: Vec<f32>,
+    expansion: f64,
+}
+
+impl Field {
+    fn new(src: &Surface, bounds: Rect, expansion: f64) -> Field {
+        let long = f64::from(bounds.width().max(bounds.height()));
+        let grow = expansion.max(0.0).ceil();
+        let g = ((long + 2.0 * grow + 8.0) / FIELD_BLOCKS).ceil().max(1.0) as i32;
+        let area = bounds.inflate(grow as i32 + 2 * g);
+        let bw = ((area.width() as i32 + g - 1) / g).max(1) as usize;
+        let bh = ((area.height() as i32 + g - 1) / g).max(1) as usize;
+        let mut cov = vec![false; bw * bh];
+        let fmt = src.format();
+        let n = fmt.channels();
+        let block = |bx: usize, by: usize| {
+            let (x, y) = (area.x0 + bx as i32 * g, area.y0 + by as i32 * g);
+            Rect::new(x, y, x + g, y + g)
+        };
+        if !fmt.alpha {
+            for (i, c) in cov.iter_mut().enumerate() {
+                *c = !block(i % bw, i / bw).intersect(&bounds).is_empty();
+            }
+        } else {
+            // One block row per task (rows are independent).
+            let scan = |by: usize, out: &mut [bool]| {
+                let r = Rect::new(area.x0, area.y0 + by as i32 * g, area.x0 + bw as i32 * g, area.y0 + (by as i32 + 1) * g).intersect(&bounds);
+                if r.is_empty() || !src.has_tiles_in(r) {
+                    return;
+                }
+                let row = src.read_region(r);
+                let rw = r.width() as usize;
+                for (i, p) in row.chunks_exact(n).enumerate() {
+                    if p.last().is_some_and(|&a| a > ALPHA_MIN) {
+                        let x = r.x0 + (i % rw) as i32;
+                        if let Some(o) = out.get_mut(((x - area.x0) / g) as usize) {
+                            *o = true;
+                        }
+                    }
+                }
+            };
+            #[cfg(not(target_arch = "wasm32"))]
+            {
+                use rayon::prelude::*;
+                cov.par_chunks_mut(bw).enumerate().for_each(|(by, out)| scan(by, out));
+            }
+            #[cfg(target_arch = "wasm32")]
+            cov.chunks_mut(bw).enumerate().for_each(|(by, out)| scan(by, out));
+        }
+        let outside = edt(&cov, true, bw, bh);
+        let inside = edt(&cov, false, bw, bh);
+        let gf = f64::from(g);
+        let sd = cov
+            .iter()
+            .zip(outside.iter().zip(&inside))
+            .map(|(&c, (&o, &i))| if c { (-(i.sqrt() - 0.5) * gf) as f32 } else { ((o.sqrt() - 0.5) * gf) as f32 })
+            .collect();
+        Field { x0: f64::from(area.x0), y0: f64::from(area.y0), g: gf, bw, bh, sd, expansion }
+    }
+
+    /// The field at a block centre (exact, no interpolation).
+    fn at(&self, bx: usize, by: usize) -> f64 {
+        self.sd.get(by * self.bw + bx).map_or(f64::MAX, |&d| f64::from(d) - self.expansion)
+    }
+
+    /// The field at `p` (bilinear between block centres; clamped to the grid, which ends
+    /// outside the region).
+    fn value(&self, p: [f64; 2]) -> f64 {
+        let fx = ((p[0] - self.x0) / self.g - 0.5).clamp(0.0, (self.bw - 1) as f64);
+        let fy = ((p[1] - self.y0) / self.g - 0.5).clamp(0.0, (self.bh - 1) as f64);
+        let (ix, iy) = (fx.floor() as usize, fy.floor() as usize);
+        let (ix1, iy1) = ((ix + 1).min(self.bw - 1), (iy + 1).min(self.bh - 1));
+        let (tx, ty) = (fx - ix as f64, fy - iy as f64);
+        let top = self.at(ix, iy) * (1.0 - tx) + self.at(ix1, iy) * tx;
+        let bottom = self.at(ix, iy1) * (1.0 - tx) + self.at(ix1, iy1) * tx;
+        top * (1.0 - ty) + bottom * ty
+    }
+
+    /// The outline (where the field crosses zero) as closed loops, traced by marching squares
+    /// over the block centres. The field's border blocks lie outside, so every loop closes.
+    fn outline(&self) -> Vec<Vec<[f64; 2]>> {
+        let (bw, bh) = (self.bw, self.bh);
+        if bw < 2 || bh < 2 {
+            return Vec::new();
+        }
+        let inside = |x: usize, y: usize| self.at(x, y) < 0.0;
+        // A crossing is named by its block edge: 2·block for the edge to the right neighbour,
+        // 2·block + 1 for the edge to the one below.
+        let id = |x: usize, y: usize, down: bool| 2 * (y * bw + x) + usize::from(down);
+        let point = |id: usize| {
+            let b = id / 2;
+            let (x, y) = (b % bw, b / bw);
+            let (x2, y2) = if id.is_multiple_of(2) { (x + 1, y) } else { (x, y + 1) };
+            let (fa, fb) = (self.at(x, y), self.at(x2, y2));
+            let t = if (fa - fb).abs() > 1e-12 { (fa / (fa - fb)).clamp(0.0, 1.0) } else { 0.5 };
+            let (dx, dy) = ((x2 - x) as f64, (y2 - y) as f64);
+            [self.x0 + (x as f64 + 0.5 + dx * t) * self.g, self.y0 + (y as f64 + 0.5 + dy * t) * self.g]
+        };
+        let mut links: HashMap<usize, [usize; 2]> = HashMap::new();
+        let mut link = |a: usize, b: usize| {
+            for (p, q) in [(a, b), (b, a)] {
+                let e = links.entry(p).or_insert([usize::MAX; 2]);
+                if e[0] == usize::MAX {
+                    e[0] = q;
+                } else {
+                    e[1] = q;
+                }
+            }
+        };
+        for y in 0..bh - 1 {
+            for x in 0..bw - 1 {
+                let c = [inside(x, y), inside(x + 1, y), inside(x + 1, y + 1), inside(x, y + 1)];
+                if c.iter().all(|&v| v == c[0]) {
+                    continue;
+                }
+                // Top, right, bottom, left.
+                let edges = [id(x, y, false), id(x + 1, y, true), id(x, y + 1, false), id(x, y, true)];
+                let cross = [c[0] != c[1], c[1] != c[2], c[3] != c[2], c[0] != c[3]];
+                let ks: Vec<usize> = (0..4).filter(|&k| cross[k]).collect();
+                if let [a, b] = ks[..] {
+                    link(edges[a], edges[b]);
+                } else if ks.len() == 4 {
+                    // A saddle: the centre decides which diagonal pair is joined.
+                    let centre = self.at(x, y) + self.at(x + 1, y) + self.at(x + 1, y + 1) + self.at(x, y + 1) < 0.0;
+                    if c[0] == centre {
+                        link(edges[0], edges[1]);
+                        link(edges[2], edges[3]);
+                    } else {
+                        link(edges[3], edges[0]);
+                        link(edges[1], edges[2]);
+                    }
+                }
+            }
+        }
+        let mut keys: Vec<usize> = links.keys().copied().collect();
+        keys.sort_unstable();
+        let mut seen: HashSet<usize> = HashSet::with_capacity(keys.len());
+        let mut loops = Vec::new();
+        for &start in &keys {
+            if seen.contains(&start) {
+                continue;
+            }
+            let mut lp = Vec::new();
+            let (mut prev, mut cur) = (usize::MAX, start);
+            loop {
+                seen.insert(cur);
+                lp.push(point(cur));
+                let Some(&[a, b]) = links.get(&cur) else { break };
+                let next = if a != prev { a } else { b };
+                if next == usize::MAX || seen.contains(&next) {
+                    break;
+                }
+                (prev, cur) = (cur, next);
+            }
+            if lp.len() >= 3 {
+                loops.push(lp);
+            }
+        }
+        loops
+    }
+}
+
+/// Squared Euclidean distance (in blocks) from every block to the nearest block whose `mask`
+/// equals `site` (Felzenszwalb & Huttenlocher's separable lower-envelope transform). About
+/// 1e20 when there is none.
+fn edt(mask: &[bool], site: bool, w: usize, h: usize) -> Vec<f64> {
+    const FAR: f64 = 1e20;
+    let mut d: Vec<f64> = mask.iter().map(|&m| if m == site { 0.0 } else { FAR }).collect();
+    let mut f = vec![0.0; w.max(h)];
+    let mut out = vec![0.0; w.max(h)];
+    let mut v = vec![0usize; w.max(h)];
+    let mut z = vec![0.0; w.max(h) + 1];
+    for x in 0..w {
+        for (y, fy) in f[..h].iter_mut().enumerate() {
+            *fy = d[y * w + x];
+        }
+        edt_1d(&f[..h], &mut out[..h], &mut v, &mut z);
+        for (y, o) in out[..h].iter().enumerate() {
+            d[y * w + x] = *o;
+        }
+    }
+    for row in d.chunks_mut(w) {
+        f[..w].copy_from_slice(row);
+        edt_1d(&f[..w], row, &mut v, &mut z);
+    }
+    d
+}
+
+/// One-dimensional squared distance transform of the sampled function `f`.
+fn edt_1d(f: &[f64], out: &mut [f64], v: &mut [usize], z: &mut [f64]) {
+    let n = f.len();
+    if n == 0 {
+        return;
+    }
+    let mut k = 0usize;
+    v[0] = 0;
+    z[0] = f64::NEG_INFINITY;
+    z[1] = f64::INFINITY;
+    let meet = |q: usize, p: usize| {
+        let (qf, pf) = (q as f64, p as f64);
+        ((f[q] + qf * qf) - (f[p] + pf * pf)) / (2.0 * (qf - pf))
+    };
+    for q in 1..n {
+        let mut s = meet(q, v[k]);
+        while k > 0 && s <= z[k] {
+            k -= 1;
+            s = meet(q, v[k]);
+        }
+        k += 1;
+        v[k] = q;
+        z[k] = s;
+        z[k + 1] = f64::INFINITY;
+    }
+    k = 0;
+    for (q, o) in out.iter_mut().enumerate() {
+        let qf = q as f64;
+        while z[k + 1] < qf {
+            k += 1;
+        }
+        let p = v[k];
+        *o = (qf - p as f64).powi(2) + f[p];
+    }
+}
+
+/// Distance from `p` to the segment `a`–`b`.
+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 };
+    dist2(p, [a[0] + dx * t, a[1] + dy * t]).sqrt()
+}
+
+/// Simplifies a closed loop (Douglas–Peucker, within `tol` px), which keeps its corners.
+fn simplify(lp: &[[f64; 2]], tol: f64) -> Vec<[f64; 2]> {
+    let n = lp.len();
+    if n < 4 {
+        return lp.to_vec();
+    }
+    let far = (1..n).max_by(|&a, &b| dist2(lp[0], lp[a]).total_cmp(&dist2(lp[0], lp[b]))).unwrap_or(n / 2);
+    let mut keep = vec![false; n];
+    keep[0] = true;
+    keep[far] = true;
+    // Spans (a, b) of the loop; b == n stands for point 0.
+    let mut stack = vec![(0, far), (far, n)];
+    while let Some((a, b)) = stack.pop() {
+        let (pa, pb) = (lp[a], lp[b % n]);
+        let mut best = (0.0, 0);
+        for (i, &p) in lp.iter().enumerate().take(b).skip(a + 1) {
+            let d = seg_dist(p, pa, pb);
+            if d > best.0 {
+                best = (d, i);
+            }
+        }
+        if best.0 > tol {
+            keep[best.1] = true;
+            stack.push((a, best.1));
+            stack.push((best.1, b));
+        }
+    }
+    let out: Vec<[f64; 2]> = (0..n).filter(|&i| keep[i]).map(|i| lp[i]).collect();
+    if out.len() >= 3 {
+        return out;
+    }
+    // A speck thinner than the tolerance: keep a few of its points.
+    (0..n).step_by((n / 4).max(1)).map(|i| lp[i]).collect()
+}
+
+/// Splits the sides of a closed loop longer than `step` into equal parts.
+fn densify(lp: &[[f64; 2]], step: f64) -> Vec<[f64; 2]> {
+    let mut out = Vec::with_capacity(lp.len());
+    for (k, &a) in lp.iter().enumerate() {
+        let b = lp[(k + 1) % lp.len()];
+        let parts = (dist2(a, b).sqrt() / step).ceil().clamp(1.0, 1e4) as usize;
+        for i in 0..parts {
+            out.push(lerp(a, b, i as f64 / parts as f64));
+        }
+    }
+    out
+}
+
+fn lerp(a: [f64; 2], b: [f64; 2], t: f64) -> [f64; 2] {
+    [a[0] + (b[0] - a[0]) * t, a[1] + (b[1] - a[1]) * t]
+}
+
+fn orient(a: [f64; 2], b: [f64; 2], c: [f64; 2]) -> f64 {
+    (b[0] - a[0]) * (c[1] - a[1]) - (b[1] - a[1]) * (c[0] - a[0])
+}
+
+/// Positive when `d` lies inside the circumcircle of the positively oriented `a`, `b`, `c`.
+fn in_circle(a: [f64; 2], b: [f64; 2], c: [f64; 2], d: [f64; 2]) -> f64 {
+    let (ax, ay) = (a[0] - d[0], a[1] - d[1]);
+    let (bx, by) = (b[0] - d[0], b[1] - d[1]);
+    let (cx, cy) = (c[0] - d[0], c[1] - d[1]);
+    (ax * ax + ay * ay) * (bx * cy - cx * by) - (bx * bx + by * by) * (ax * cy - cx * ay) + (cx * cx + cy * cy) * (ax * by - bx * ay)
+}
+
+const NONE: usize = usize::MAX;
+
+/// A triangle of the Delaunay triangulation: vertices (positively oriented) and the neighbour
+/// across the edge opposite each vertex.
+#[derive(Clone, Copy)]
+struct DTri {
+    v: [usize; 3],
+    n: [usize; 3],
+    alive: bool,
+}
+
+/// Incremental Delaunay triangulation (Bowyer–Watson), started from a triangle around the area.
+/// Vertices 0..3 are that triangle's.
+struct Delaunay {
+    pts: Vec<[f64; 2]>,
+    tris: Vec<DTri>,
+    stamp: Vec<u32>,
+    pass: u32,
+    last: usize,
+}
+
+impl Delaunay {
+    fn new(lo: [f64; 2], hi: [f64; 2]) -> Delaunay {
+        let (cx, cy) = (0.5 * (lo[0] + hi[0]), 0.5 * (lo[1] + hi[1]));
+        let d = (hi[0] - lo[0]).max(hi[1] - lo[1]).max(1.0) * 4.0;
+        let mut pts = vec![[cx - 2.0 * d, cy - d], [cx + 2.0 * d, cy - d], [cx, cy + 2.0 * d]];
+        if orient(pts[0], pts[1], pts[2]) < 0.0 {
+            pts.swap(1, 2);
+        }
+        Delaunay { pts, tris: vec![DTri { v: [0, 1, 2], n: [NONE; 3], alive: true }], stamp: vec![0], pass: 0, last: 0 }
+    }
+
+    /// The triangle holding `p`: a visibility walk from the last one made.
+    fn locate(&self, p: [f64; 2]) -> Option<usize> {
+        let mut t = if self.tris.get(self.last).is_some_and(|t| t.alive) { self.last } else { self.tris.iter().position(|t| t.alive)? };
+        for _ in 0..self.tris.len() + 3 {
+            let tri = self.tris.get(t)?;
+            let step = (0..3).find(|&k| orient(self.pts[tri.v[(k + 1) % 3]], self.pts[tri.v[(k + 2) % 3]], p) < 0.0);
+            match step {
+                None => return Some(t),
+                Some(k) if tri.n[k] != NONE => t = tri.n[k],
+                Some(_) => return None,
+            }
+        }
+        self.tris.iter().position(|tri| tri.alive && (0..3).all(|k| orient(self.pts[tri.v[(k + 1) % 3]], self.pts[tri.v[(k + 2) % 3]], p) >= 0.0))
+    }
+
+    /// Inserts `p`; returns its vertex index (an existing vertex's when `p` coincides with it).
+    fn insert(&mut self, p: [f64; 2]) -> Option<usize> {
+        let t0 = self.locate(p)?;
+        if let Some(&v) = self.tris[t0].v.iter().find(|&&v| dist2(self.pts[v], p) < 1e-12) {
+            return Some(v);
+        }
+        self.pass = self.pass.wrapping_add(1);
+        let pass = self.pass;
+        // The cavity: every triangle whose circumcircle holds `p`, grown from `t0`.
+        let mut bad = vec![t0];
+        self.stamp[t0] = pass;
+        let mut i = 0;
+        while let Some(&t) = bad.get(i) {
+            i += 1;
+            for k in 0..3 {
+                let nb = self.tris[t].n[k];
+                if nb == NONE || self.stamp[nb] == pass {
+                    continue;
+                }
+                let [a, b, c] = self.tris[nb].v.map(|v| self.pts[v]);
+                if in_circle(a, b, c, p) > 0.0 {
+                    self.stamp[nb] = pass;
+                    bad.push(nb);
+                }
+            }
+        }
+        let mut rim: Vec<(usize, usize, usize)> = Vec::new();
+        for &t in &bad {
+            let tri = self.tris[t];
+            for k in 0..3 {
+                let nb = tri.n[k];
+                if nb == NONE || self.stamp[nb] != pass {
+                    rim.push((tri.v[(k + 1) % 3], tri.v[(k + 2) % 3], nb));
+                }
+            }
+            self.tris[t].alive = false;
+        }
+        let pi = self.pts.len();
+        self.pts.push(p);
+        let base = self.tris.len();
+        for (i, &(a, b, o)) in rim.iter().enumerate() {
+            self.tris.push(DTri { v: [a, b, pi], n: [NONE, NONE, o], alive: true });
+            self.stamp.push(0);
+            if let Some(ot) = self.tris.get_mut(o) {
+                for k in 0..3 {
+                    let (x, y) = (ot.v[(k + 1) % 3], ot.v[(k + 2) % 3]);
+                    if (x == b && y == a) || (x == a && y == b) {
+                        ot.n[k] = base + i;
+                    }
+                }
+            }
+        }
+        for (i, &(a, b, _)) in rim.iter().enumerate() {
+            let across_b = rim.iter().position(|r| r.0 == b).map_or(NONE, |j| base + j);
+            let across_a = rim.iter().position(|r| r.1 == a).map_or(NONE, |j| base + j);
+            self.tris[base + i].n[0] = across_b;
+            self.tris[base + i].n[1] = across_a;
+        }
+        self.last = base;
+        Some(pi)
+    }
+
+    fn edges(&self) -> HashSet<(usize, usize)> {
+        let mut e = HashSet::new();
+        for t in self.tris.iter().filter(|t| t.alive) {
+            for k in 0..3 {
+                let (a, b) = (t.v[k], t.v[(k + 1) % 3]);
+                e.insert((a.min(b), a.max(b)));
+            }
+        }
+        e
+    }
+}
+
+/// Even-odd point-in-polygon test against closed outline segments, bucketed by row.
+struct Inside<'a> {
+    pts: &'a [[f64; 2]],
+    segs: &'a [(usize, usize)],
+    y0: f64,
+    row: f64,
+    rows: Vec<Vec<usize>>,
+}
+
+impl<'a> Inside<'a> {
+    fn new(pts: &'a [[f64; 2]], segs: &'a [(usize, usize)], row: f64) -> Self {
+        let y0 = segs.iter().flat_map(|&(a, b)| [pts[a][1], pts[b][1]]).fold(f64::MAX, f64::min);
+        let y1 = segs.iter().flat_map(|&(a, b)| [pts[a][1], pts[b][1]]).fold(f64::MIN, f64::max);
+        let n = if y1 >= y0 { ((y1 - y0) / row).floor().clamp(0.0, 1e5) as usize + 1 } else { 0 };
+        let mut rows = vec![Vec::new(); n];
+        for (i, &(a, b)) in segs.iter().enumerate() {
+            let (lo, hi) = (pts[a][1].min(pts[b][1]), pts[a][1].max(pts[b][1]));
+            let (r0, r1) = (((lo - y0) / row).floor() as usize, ((hi - y0) / row).floor() as usize);
+            for r in rows.iter_mut().take(r1 + 1).skip(r0) {
+                r.push(i);
+            }
+        }
+        Inside { pts, segs, y0, row, rows }
+    }
+
+    fn contains(&self, p: [f64; 2]) -> bool {
+        if p[1] < self.y0 {
+            return false;
+        }
+        let Some(row) = self.rows.get(((p[1] - self.y0) / self.row).floor() as usize) else { return false };
+        let mut odd = false;
+        for &i in row {
+            let (a, b) = (self.pts[self.segs[i].0], self.pts[self.segs[i].1]);
+            if (a[1] > p[1]) != (b[1] > p[1]) && p[0] < a[0] + (p[1] - a[1]) / (b[1] - a[1]) * (b[0] - a[0]) {
+                odd = !odd;
+            }
+        }
+        odd
+    }
+}
+
+/// Most outline points the mesh takes; a busier outline is simplified more coarsely.
+const MAX_OUTLINE_POINTS: usize = 20_000;
+
 /// Builds the mesh over the covered part of `src` inside `bounds`.
+///
+/// The outline of the opaque region grown by `expansion` is traced, simplified (corners kept)
+/// and split into sides no longer than the density's spacing; a triangular lattice at that
+/// spacing fills the inside, at least half a step from the outline. The mesh is the
+/// Delaunay triangulation of those points with the outline's sides recovered (split until
+/// each is an edge), minus the triangles outside the outline.
 pub fn build_mesh(src: &Surface, bounds: Rect, density: PuppetDensity, expansion: f64) -> PuppetMesh {
     if bounds.is_empty() {
         return PuppetMesh::default();
     }
-    let long = f64::from(bounds.width().max(bounds.height()));
+    let e = if expansion.is_finite() { expansion.clamp(-200.0, 200.0) } else { 0.0 };
+    // The density counts steps across the region the mesh covers (grown by the expansion).
+    let long = f64::from(bounds.width().max(bounds.height())) + 2.0 * e.max(0.0);
     let spacing = (long / density.cells()).max(2.0);
-    let e = expansion.clamp(-200.0, 200.0);
-    // Coverage on a block grid (block = a quarter cell or less).
-    let g = ((spacing / 4.0).floor() as i32).max(1);
-    let pad = e.max(0.0).ceil() as i32 + g;
-    let area = bounds.inflate(pad);
-    let bw = (area.width() as i32 + g - 1) / g;
-    let bh = (area.height() as i32 + g - 1) / g;
-    let (bw, bh) = (bw as usize, bh as usize);
-    let mut cov = vec![false; bw * bh];
-    let fmt = src.format();
-    let n = fmt.channels();
-    if !fmt.alpha {
-        for by in 0..bh {
-            for bx in 0..bw {
-                let r = Rect::new(area.x0 + bx as i32 * g, area.y0 + by as i32 * g, area.x0 + (bx as i32 + 1) * g, area.y0 + (by as i32 + 1) * g);
-                cov[by * bw + bx] = !r.intersect(&bounds).is_empty();
-            }
-        }
-    } else {
-        // One block row per task (rows are independent).
-        let scan = |by: usize, out: &mut [bool]| {
-            let r = Rect::new(area.x0, area.y0 + by as i32 * g, area.x0 + bw as i32 * g, area.y0 + (by as i32 + 1) * g).intersect(&bounds);
-            if r.is_empty() || !src.has_tiles_in(r) {
-                return;
-            }
-            let row = src.read_region(r);
-            let rw = r.width() as usize;
-            for (i, p) in row.chunks_exact(n).enumerate() {
-                if p[n - 1] > ALPHA_MIN {
-                    let x = r.x0 + (i % rw) as i32;
-                    out[((x - area.x0) / g) as usize] = true;
-                }
-            }
-        };
-        #[cfg(not(target_arch = "wasm32"))]
-        {
-            use rayon::prelude::*;
-            cov.par_chunks_mut(bw).enumerate().for_each(|(by, out)| scan(by, out));
-        }
-        #[cfg(target_arch = "wasm32")]
-        cov.chunks_mut(bw).enumerate().for_each(|(by, out)| scan(by, out));
+    let field = Field::new(src, bounds, e);
+    let loops = field.outline();
+    let empty = PuppetMesh { spacing, ..Default::default() };
+    if loops.is_empty() {
+        return empty;
     }
-    // Grow / shrink by the expansion (round stamp in block units).
-    let eb = (e.abs() / f64::from(g)).round() as i32;
-    if eb > 0 {
-        let src_cov = cov.clone();
-        let grow = e > 0.0;
-        for by in 0..bh as i32 {
-            for bx in 0..bw as i32 {
-                let mut hit = !grow;
-                'stamp: for dy in -eb..=eb {
-                    for dx in -eb..=eb {
-                        if dx * dx + dy * dy > eb * eb {
-                            continue;
-                        }
-                        let (x, y) = (bx + dx, by + dy);
-                        let v = x >= 0 && y >= 0 && x < bw as i32 && y < bh as i32 && src_cov[y as usize * bw + x as usize];
-                        if grow && v {
-                            hit = true;
-                            break 'stamp;
-                        }
-                        if !grow && !v {
-                            hit = false;
-                            break 'stamp;
-                        }
-                    }
-                }
-                cov[by as usize * bw + bx as usize] = hit;
+    let mut tol = (0.5 * field.g).max(0.5);
+    let rings = loop {
+        let rings: Vec<Vec<[f64; 2]>> = loops.iter().map(|l| densify(&simplify(l, tol), spacing)).collect();
+        if rings.iter().map(Vec::len).sum::<usize>() <= MAX_OUTLINE_POINTS || tol > spacing {
+            break rings;
+        }
+        tol *= 2.0;
+    };
+    let (lo, hi) = ([field.x0, field.y0], [field.x0 + field.bw as f64 * field.g, field.y0 + field.bh as f64 * field.g]);
+    let mut dl = Delaunay::new(lo, hi);
+    let mut segs: Vec<(usize, usize)> = Vec::new();
+    for ring in &rings {
+        let idx: Vec<usize> = ring.iter().filter_map(|&p| dl.insert(p)).collect();
+        for (k, &a) in idx.iter().enumerate() {
+            let b = idx[(k + 1) % idx.len()];
+            if a != b {
+                segs.push((a, b));
             }
         }
     }
-    // Grid cells touching covered blocks.
-    let ox = f64::from(area.x0);
-    let oy = f64::from(area.y0);
-    let cw = (f64::from(bw as i32 * g) / spacing).ceil() as usize;
-    let ch = (f64::from(bh as i32 * g) / spacing).ceil() as usize;
-    let mut index = vec![usize::MAX; (cw + 1) * (ch + 1)];
-    let mut mesh = PuppetMesh { spacing, ..Default::default() };
-    for cj in 0..ch {
-        for ci in 0..cw {
-            let bx0 = ((ci as f64 * spacing) / f64::from(g)).floor() as usize;
-            let by0 = ((cj as f64 * spacing) / f64::from(g)).floor() as usize;
-            let bx1 = (((ci + 1) as f64 * spacing) / f64::from(g)).ceil().min(bw as f64) as usize;
-            let by1 = (((cj + 1) as f64 * spacing) / f64::from(g)).ceil().min(bh as f64) as usize;
-            let any = (by0..by1).any(|y| (bx0..bx1).any(|x| cov[y * bw + x]));
-            if !any {
+    // The inside: a triangular lattice, kept half a step clear of the outline.
+    let rise = spacing * 3f64.sqrt() / 2.0;
+    let rows = ((hi[1] - lo[1]) / rise).ceil() as usize;
+    let cols = ((hi[0] - lo[0]) / spacing).ceil() as usize;
+    for j in 0..=rows {
+        for i in 0..=cols {
+            let p = [lo[0] + (i as f64 + 0.5 * (j % 2) as f64) * spacing, lo[1] + j as f64 * rise];
+            if field.value(p) < -0.5 * spacing {
+                dl.insert(p);
+            }
+        }
+    }
+    // Recover the outline's sides: split a side the triangulation lacks at its midpoint.
+    let mut budget = segs.len() * 2;
+    for _ in 0..8 {
+        let edges = dl.edges();
+        let mut next = Vec::with_capacity(segs.len());
+        let mut split = false;
+        for &(a, b) in &segs {
+            if budget == 0 || edges.contains(&(a.min(b), a.max(b))) {
+                next.push((a, b));
                 continue;
             }
-            let mut v = |i: usize, j: usize| {
-                let k = j * (cw + 1) + i;
-                if index[k] == usize::MAX {
-                    index[k] = mesh.verts.len();
-                    mesh.verts.push([ox + i as f64 * spacing, oy + j as f64 * spacing]);
+            match dl.insert(lerp(dl.pts[a], dl.pts[b], 0.5)) {
+                Some(c) if c != a && c != b => {
+                    next.push((a, c));
+                    next.push((c, b));
+                    budget -= 1;
+                    split = true;
                 }
-                index[k]
-            };
-            let (a, b, c, d) = (v(ci, cj), v(ci + 1, cj), v(ci + 1, cj + 1), v(ci, cj + 1));
-            if (ci + cj) % 2 == 0 {
-                mesh.tris.push([a, b, c]);
-                mesh.tris.push([a, c, d]);
-            } else {
-                mesh.tris.push([a, b, d]);
-                mesh.tris.push([b, c, d]);
+                _ => next.push((a, b)),
+            }
+        }
+        segs = next;
+        if !split {
+            break;
+        }
+    }
+    let inside = Inside::new(&dl.pts, &segs, spacing);
+    let mut tris: Vec<[usize; 3]> = Vec::new();
+    for t in dl.tris.iter().filter(|t| t.alive && t.v.iter().all(|&v| v >= 3)) {
+        let [a, b, c] = t.v.map(|v| dl.pts[v]);
+        if orient(a, b, c) > 1e-9 && inside.contains([(a[0] + b[0] + c[0]) / 3.0, (a[1] + b[1] + c[1]) / 3.0]) {
+            tris.push(t.v);
+        }
+    }
+    if tris.is_empty() {
+        return empty;
+    }
+    // Keep the vertices in use, numbered for a narrow band in the solver's matrix.
+    let pos = &dl.pts;
+    let order = band_order(pos, &tris);
+    let mut remap = vec![usize::MAX; pos.len()];
+    for (new, &old) in order.iter().enumerate() {
+        remap[old] = new;
+    }
+    PuppetMesh { verts: order.iter().map(|&i| pos[i]).collect(), tris: tris.iter().map(|t| t.map(|i| remap[i])).collect(), spacing }
+}
+
+/// The vertices `tris` use, in reverse Cuthill–McKee order (breadth first from the topmost
+/// vertex of each piece, lower degrees first), which keeps the solver's banded factor small.
+fn band_order(pos: &[[f64; 2]], tris: &[[usize; 3]]) -> Vec<usize> {
+    let mut adj: Vec<Vec<usize>> = vec![Vec::new(); pos.len()];
+    for t in tris {
+        for k in 0..3 {
+            let (a, b) = (t[k], t[(k + 1) % 3]);
+            adj[a].push(b);
+            adj[b].push(a);
+        }
+    }
+    for a in &mut adj {
+        a.sort_unstable();
+        a.dedup();
+    }
+    let mut starts: Vec<usize> = (0..pos.len()).filter(|&i| !adj[i].is_empty()).collect();
+    starts.sort_by(|&a, &b| pos[a][1].total_cmp(&pos[b][1]).then(pos[a][0].total_cmp(&pos[b][0])));
+    let mut seen = vec![false; pos.len()];
+    let mut order = Vec::with_capacity(starts.len());
+    for s in starts {
+        if seen[s] {
+            continue;
+        }
+        seen[s] = true;
+        let first = order.len();
+        order.push(s);
+        let mut i = first;
+        while let Some(&v) = order.get(i) {
+            i += 1;
+            let mut next: Vec<usize> = adj[v].iter().copied().filter(|&w| !seen[w]).collect();
+            next.sort_by_key(|&w| adj[w].len());
+            for w in next {
+                seen[w] = true;
+                order.push(w);
             }
         }
     }
-    mesh
+    order.reverse();
+    order
 }
 
 /// A compressed sparse row matrix (symmetric, used by CG).
@@ -330,9 +833,9 @@ impl Csr {
     }
 }
 
-/// Banded Cholesky factor `A = L Lᵀ` (grid meshes number their vertices row by row, so the
-/// system's bandwidth is about two mesh rows): factored once per pin set, then each global step
-/// is two triangular solves.
+/// Banded Cholesky factor `A = L Lᵀ` (the mesh numbers its vertices in reverse Cuthill–McKee
+/// order, so the system's bandwidth is about one breadth-first front): factored once per pin
+/// set, then each global step is two triangular solves.
 struct BandChol {
     n: usize,
     p: usize,
@@ -833,6 +1336,43 @@ mod tests {
         assert!(more.tris.len() > m.tris.len() && m.tris.len() > fewer.tris.len());
     }
 
+    fn in_mesh(m: &PuppetMesh, p: [f64; 2]) -> bool {
+        m.tris.iter().any(|t| {
+            let [a, b, c] = t.map(|i| m.verts[i]);
+            let (d0, d1, d2) = (orient(a, b, p), orient(b, c, p), orient(c, a, p));
+            (d0 >= 0.0 && d1 >= 0.0 && d2 >= 0.0) || (d0 <= 0.0 && d1 <= 0.0 && d2 <= 0.0)
+        })
+    }
+
+    #[test]
+    fn mesh_follows_the_outline() {
+        // #3105: the mesh was every grid cell touching the shape, a staircase up to a cell
+        // outside it. It now follows the outline grown by the expansion.
+        let mut s = Surface::new(PixelFormat::new(ColorMode::Rgb, SampleType::U8, true));
+        for y in 40..160 {
+            for x in 40..160 {
+                if (f64::from(x) + 0.5 - 100.0).hypot(f64::from(y) + 0.5 - 100.0) <= 60.0 {
+                    s.write_pixel(x, y, &[0.3, 0.6, 0.9, 1.0]);
+                }
+            }
+        }
+        let m = build_mesh(&s, s.content_bounds(), PuppetDensity::Normal, 2.0);
+        let r = |v: &[f64; 2]| (v[0] - 100.0).hypot(v[1] - 100.0);
+        assert!(m.verts.iter().all(|v| r(v) <= 63.0), "a vertex outside the grown disc");
+        assert!(m.verts.iter().filter(|v| r(v) > 61.0).count() >= 40, "the outline carries vertices");
+        let area: f64 = m.tris.iter().map(|t| orient(m.verts[t[0]], m.verts[t[1]], m.verts[t[2]]).abs() / 2.0).sum();
+        let disc = std::f64::consts::PI * 62.0 * 62.0;
+        assert!((area / disc - 1.0).abs() < 0.03, "{area} vs {disc}");
+        // Sharp corners and a concave corner: every opaque pixel is still inside the mesh, even
+        // without an expansion, and the empty corner of the L is not.
+        let l = blob(SampleType::U8);
+        let m = build_mesh(&l, l.content_bounds(), PuppetDensity::Fewer, 0.0);
+        for (x, y) in (20..36).flat_map(|x| (10..70).map(move |y| (x, y))).chain((36..70).flat_map(|x| (56..70).map(move |y| (x, y)))) {
+            assert!(in_mesh(&m, [f64::from(x) + 0.5, f64::from(y) + 0.5]), "({x}, {y}) is outside the mesh");
+        }
+        assert!(!in_mesh(&m, [40.0, 50.0]) && !in_mesh(&m, [60.0, 20.0]));
+    }
+
     #[test]
     fn no_pins_or_unmoved_pins_are_identity() {
         for st in [SampleType::U8, SampleType::U16, SampleType::F32] {

crates/engine/src/distort_cmds/tests.rs

@@ -180,6 +180,42 @@ fn puppet_two_pins_bend_the_block() {
     assert!(after.content_bounds().y0 < 6, "{:?}", after.content_bounds());
 }
 
+#[test]
+fn puppet_mesh_follows_the_outline_so_separate_shapes_move_separately() {
+    // #3105: the mesh was every grid cell touching opaque pixels, so two bars 24 px apart (more
+    // than a mesh cell at Fewer Points) shared cells and a pin on one dragged the other along.
+    // The mesh now follows the outline grown by the expansion: two meshes, and the unpinned one
+    // stays put.
+    for depth in [8u64, 16, 32] {
+        let mut s = Session::new();
+        s.execute("file.new", json!({"width": 260, "height": 140, "depth": depth})).unwrap();
+        s.execute("layer.new.layer", json!({})).unwrap();
+        s.edit("paint", |doc, active| {
+            let surf = doc.layer_mut(active.unwrap()).unwrap().surface_mut().unwrap();
+            for (y0, y1) in [(20, 60), (84, 120)] {
+                for y in y0..y1 {
+                    for x in 10..250 {
+                        let v = ((x / 4 + y / 4) % 2) as f32;
+                        surf.write_pixel(x, y, &[v, 1.0 - v, (y as f32) / 140.0, 1.0]);
+                    }
+                }
+            }
+            Ok(())
+        })
+        .unwrap();
+        let before = active_surface(&s);
+        s.execute(PUPPET, json!({"pins": [{"src": [130, 40], "dst": [130, 30]}], "density": "fewer", "expansion": 2, "interpolation": "bilinear"})).unwrap();
+        let after = active_surface(&s).convert(before.format());
+        let lower = Rect::new(0, 70, 260, 140);
+        let still = worst(&after, &before, lower);
+        assert!(still <= 1.0 / 255.0, "{depth}: the unpinned bar stays ({still})");
+        assert!(
+            worst(&after, &photocraft_algo::resample::translate_surface(&before, 0, -10), Rect::new(0, 0, 260, 60)) <= 1.0 / 255.0,
+            "{depth}: the pinned bar moves"
+        );
+    }
+}
+
 #[test]
 fn perspective_identity_known_homography_and_mask() {
     for depth in [8u64, 16, 32] {

docs/gaps.md

@@ -68,8 +68,7 @@ what an agent can't finish alone.
 - **Missing:** Select and Mask workspace layout
   (TOOL-214-1, Refine Hair #2671; its view modes preview live since #1754, Refine Edge Brush #2972); Content-Aware Fill workspace with a paintable sampling area
   (TOOL-214-15); Vanishing Point plane editor (TOOL-214-14); Content-Aware Scale through the
-  transform box (TOOL-214-7); Puppet Warp mesh following the outline (TOOL-214-6); Content-Aware
-  Move's Transform On Drop (TOOL-213-4).
+  transform box (TOOL-214-7); Content-Aware Move's Transform On Drop (TOOL-213-4).
 - **Evidence:** scorecard rows above, all partial or missing.
 - **Impact:** these are the tools retouchers and compositors use daily; a slider dialog is not a
   substitute.

docs/scorecard.md

@@ -8,7 +8,7 @@ Generated by `cargo xtask scorecard` from `perf/budgets.toml`, `perf/baseline.js
 |---|---:|---:|---:|---|
 | [Performance](#performance) | 3 | 22 | 22 | 56 scenarios: 3 within budget, 22 over budget or failing, 9 tracked only, 17 not run, 5 not measurable yet (baseline 2026-10-05) |
 | [File compatibility](#file-compatibility) | 6 | 13 | 8 | io corpus: oracle ≥ 146/170, round trip ≥ 169/170; psd-tools corpus: oracle ≥ 236/309, round trip ≥ 308/309; photoshop oracles: oracle ≥ 133/256, round trip ≥ 258/256 |
-| [Tools](#tools) | 27 | 14 | 2 |  |
+| [Tools](#tools) | 28 | 14 | 1 |  |
 | [Workspace and UI](#workspace-and-ui) | 12 | 6 | 5 | settings that do nothing: 30 of 155 (target 0) |
 | [Type](#type) | 4 | 2 | 3 |  |
 | [Automation](#automation) | 3 | 2 | 6 | 631 of 631 menu items live (100.0%) |
@@ -150,7 +150,7 @@ Corpus floors (`crates/io/tests/corpus.rs`): the enforced minimums; the corpus t
 | TOOL-214-3 | Color Range test: Sampled Colors at Fuzziness 0 selects the exact match only | partial | [#214](https://github.com/storytold/photocraft/issues/214) | Fuzziness-0 tests cover midtones/shadows only (selection_cmds.rs:444-447) |
 | TOOL-214-4 | Move tool: box-select layers, Auto-Select hover outline, double-click fit, drag to a tab | partial | [#214](https://github.com/storytold/photocraft/issues/214) | Auto-Select hover outline (Preferences › Tools › Show Highlight on Rollover, #2834) and box-select layers from an empty spot (move_box.rs, layer.selectInRect, #2844) done; double-click fit and drag to a tab open |
 | TOOL-214-5 | Right-click menu during Free Transform (scale, rotate, skew, warp, flip, rotate 90/180) | done | [#214](https://github.com/storytold/photocraft/issues/214) | canvas_tool_menu.rs TRANSFORM_MENU: modes, Content-Aware Scale / Puppet Warp greyed, Rotate 180 / 90 CW / 90 CCW and Flip turn the open box (transform_tool::intercept, #2832) |
-| TOOL-214-6 | Puppet Warp mesh follows the object outline (constrained triangulation) | missing | [#214](https://github.com/storytold/photocraft/issues/214) | algo/puppet.rs:136 builds a block-grid mesh |
+| TOOL-214-6 | Puppet Warp mesh follows the object outline (constrained triangulation) | done | [#214](https://github.com/storytold/photocraft/issues/214) | algo/puppet.rs build_mesh: traced outline of the opaque region grown by Expansion, triangulated with a Density-spaced lattice (conforming Delaunay); separate shapes get separate meshes (#3105) |
 | TOOL-214-7 | Content-Aware Scale through the Free Transform box with a live proxy preview | done | [#214](https://github.com/storytold/photocraft/issues/214) | Edit > Content-Aware Scale opens a scale-only Free Transform box (cas_transform.rs) with Amount, Protect and Protect Skin Tones in the options bar; edit.contentAwareScale carves a 384 px copy in the background for the live preview, and Enter runs it with the box's size and x/y as one step (#3096) |
 | TOOL-214-8 | A measured ruler angle prefills Image Rotation > Arbitrary (0.01 deg) | done | [#214](https://github.com/storytold/photocraft/issues/214) | analysis_ui.rs: image.rotation.arbitrary opens seeded with the ruler angle and the direction that straightens it; arbitrary_rotation_starts_at_the_ruler_angle checks 5.7106 deg within 1e-3 and cw/ccw for three measured lines |
 | TOOL-214-9 | Crop: W x H x Resolution with units, saved crop presets, Straighten in the options bar | done | [#214](https://github.com/storytold/photocraft/issues/214) | Ratio dropdown and fields, Straighten (TOOL-214-21); W x H x Resolution with units, built-in size presets and Front Image (crop_size.rs, image.crop targetWidth/targetHeight/resolution, #2443); New / Delete Crop Preset save and delete crop tool presets (crop_presets.rs, tool.presets.* with tool crop, persisted with the preferences) |

scorecard/tools.toml

@@ -148,9 +148,9 @@ note = "canvas_tool_menu.rs TRANSFORM_MENU: modes, Content-Aware Scale / Puppet
 [[item]]
 id = "TOOL-214-6"
 target = "Puppet Warp mesh follows the object outline (constrained triangulation)"
-status = "missing"
+status = "done"
 issue = 214
-note = "algo/puppet.rs:136 builds a block-grid mesh"
+note = "algo/puppet.rs build_mesh: traced outline of the opaque region grown by Expansion, triangulated with a Density-spaced lattice (conforming Delaunay); separate shapes get separate meshes (#3105)"
 
 [[item]]
 id = "TOOL-214-7"