git.lucas.co / cce-designer
graphic design tool
git clone https://git.lucas.co/cce-designer.git

commitb59e6068431cee294069ef942803070a273b1125
parent6ea4f7482f
authorLucas Galante <lsgalante12@gmail.com>
date2026-09-29 21:24
perf: a remesh costs a third of what it did, and settles

The flip pass keeps a valence table where it built and sorted a list
eight times an edge; the closest-point search begins at a quarter of a
cell and works out the winner's normal alone. Both are held bit for bit
to the passes as first written.

The flip pass refuses a flip whose new edge the next split would cut,
which ends the three passes undoing one another on the same edges at
every step, and a remesh that finds nothing to do hands its input back
untouched, so what is kept by topology downstream is kept.

Co-Authored-By: Claude Fable 5.1 <noreply@anthropic.com>

 CLAUDE.md      |  43 +++++++++++++---
 src/main.rs    | 104 +++++++++++++++++++++++++++++++++++++-
 src/remesh.rs  | 156 +++++++++++++++++++++++++++++++++++++++++++++++++++++----
 src/spatial.rs |  67 +++++++++++++++++++++++--
 4 files changed, 349 insertions(+), 21 deletions(-)

diff --git a/CLAUDE.md b/CLAUDE.md
index c6d2d4f..9136ba0 100644
--- a/CLAUDE.md
+++ b/CLAUDE.md
@@ -1557,12 +1557,43 @@ a node costs is what the solve saves without it. The file is read and
 never written, and the disk cache is off for the run. On the project it
 was written for (2026-09-29; a pull, a Surface detangle with every row on
 and a remesh, 162 points growing to 525): 46 ms a frame as saved, 24
-without the detangle, 6 without the remesh. The remesh is most of it
-twice over: its own passes are some 20 ms at 525 points — the flip pass
-asks a point's valence eight times an edge, each a list built and sorted,
-and the projection's closest-point query is 6 µs — and it splits and
-collapses the same 93 edges at every step of a mesh that has stopped
-moving, so the topology the detangle keeps its lists by is new each step.
+without the detangle, 6 without the remesh. The remesh was most of it
+twice over, and three changes the same day brought the solve to 15 ms a
+frame, 11 once the mesh is at rest, where what is left is the detangle's
+own work at 524 points:
+
+- **The flip pass keeps a valence table** (`remesh::flip_pass`), counted
+  once and kept in step with the flips. It had asked `tris_of` for a
+  point's valence eight times an edge, each a list gathered, sorted and
+  counted: 9 ms of a remesh at 525 points.
+- **The closest-point search begins at a quarter of a cell**
+  (`TriGrid::closest`) and works a normal out for the winner alone. A
+  query on the surface, which is what a remesh's projection asks, is
+  answered from the cell it is in; it was 6 µs a query from the
+  twenty-seven cells about it, 12 ms a remesh. Where a search begins does
+  not change what it finds, since it ends only on a hit nearer than the
+  box searched is wide.
+- **A remesh settles.** The flip pass refuses a flip whose new edge the
+  next split would cut, as the collapse pass always refused its own.
+  Without the rule a long edge between two thin triangles was split and
+  its midpoint collapsed into a corner — the edge turned to its short
+  diagonal — and the flip pass, judging by valence alone, turned it back:
+  93 edges split, collapsed and flipped at every step of a mesh that had
+  stopped moving. And an iteration that finds nothing to split, collapse
+  or flip, with no relaxation asked, ends the remesh; on the first the
+  INPUT is handed back as it came, primitives and their order untouched,
+  so the topology the detangle keeps its lists by is the same from one
+  step to the next. The projection's grid is built when first wanted.
+
+The first two are how a remesh is run and not what it does: the passes as
+first written are kept under `cfg(test)` (`flip_pass_reference`,
+`TriGrid::closest_reference`, `remesh_reference`) and
+`the_remesh_matches_its_reference` holds the mesh to them bit for bit,
+step after step. The third changes what a remesh makes, where a flip
+would have made an overlong edge; `a_remesh_settles_and_then_leaves_the_mesh_alone`
+is its test, and fails without the rule ("still changing 162 edges after
+20 rounds"). `remesh::last_changes` is the count the test and the profile
+read.
 
 ### Simulation checkpoints
 
diff --git a/src/main.rs b/src/main.rs
index bed27a7..660b748 100644
--- a/src/main.rs
+++ b/src/main.rs
@@ -10573,6 +10573,105 @@ mod tests {
         assert!(coarsened.num_points() < dense.num_points());
     }
 
+    /// The same mesh in every particular: points, identities, primitives,
+    /// attributes.
+    fn same_mesh(a: &Detail, b: &Detail) -> bool {
+        a.positions() == b.positions()
+            && a.ids() == b.ids()
+            && a.num_prims() == b.num_prims()
+            && (0..a.num_prims()).all(|p| a.prim_points(p) == b.prim_points(p))
+            && a.points().names() == b.points().names()
+            && a.points().names().iter().all(|n| (0..a.num_points()).all(|p| a.points().value(n, p) == b.points().value(n, p)))
+    }
+
+    /// The faster flip pass and the faster closest-point search are HOW a
+    /// remesh is run and not what it does: against the passes as they were
+    /// first written, the same mesh bit for bit, step after step of a
+    /// surface being pulled about — and the same hit for a query on the
+    /// surface, off it and far from it.
+    #[test]
+    fn the_remesh_matches_its_reference() {
+        use crate::remesh::remesh_reference;
+        let fixtures = [
+            (sphere_detail(Vec3::ZERO, 1.0, 6, 8), 0.2),
+            (sphere_detail(Vec3::ZERO, 1.0, 24, 32), 0.5),
+            (sphere_detail(Vec3::new(0.3, -0.2, 1.0), 0.5, 10, 14), 0.1),
+        ];
+        for (start, target) in fixtures {
+            for relax in [0.0, 0.5] {
+                let settings = Settings { target, iterations: 3, relax, ..Default::default() };
+                let (mut a, mut b) = (start.clone(), start.clone());
+                for step in 0..6 {
+                    // A pull between remeshes, as a simulation makes one.
+                    for d in [&mut a, &mut b] {
+                        for p in 0..d.num_points() {
+                            let at = d.pos(p);
+                            let pull = Vec3::new(0.03, 0.0, 0.0) * (at.y * 3.0 + step as f32).sin();
+                            d.set_pos(p, at + pull);
+                        }
+                    }
+                    a = remesh(&a, settings);
+                    b = remesh_reference(&b, settings);
+                    assert!(same_mesh(&a, &b), "target {target}, relax {relax}, step {step}: {} points against {}", a.num_points(), b.num_points());
+                }
+                assert!(a.num_points() != start.num_points() || target == 0.1, "the fixture remeshes");
+            }
+        }
+
+        let sphere = sphere_detail(Vec3::ZERO, 1.0, 12, 16);
+        let grid = crate::spatial::TriGrid::build(&sphere);
+        let mut asked = 0;
+        for i in 0..4000 {
+            // On the surface, near it, inside it, far outside it.
+            let f = i as f32;
+            let dir = Vec3::new((f * 0.37).sin(), (f * 0.73).cos(), (f * 1.31).sin()).normalize_or_zero();
+            let p = dir * [1.0, 0.98, 1.05, 0.3, 0.0, 4.0, 40.0][i % 7];
+            let (new, old) = (grid.closest(p).unwrap(), grid.closest_reference(p).unwrap());
+            assert_eq!((new.point, new.distance, new.normal), (old.point, old.distance, old.normal), "query {p:?}");
+            asked += 1;
+        }
+        for p in 0..sphere.num_points() {
+            let (new, old) = (grid.closest(sphere.pos(p)).unwrap(), grid.closest_reference(sphere.pos(p)).unwrap());
+            assert_eq!((new.point, new.distance, new.normal), (old.point, old.distance, old.normal), "point {p}");
+        }
+        assert_eq!(asked, 4000);
+    }
+
+    /// A remesh SETTLES: run again on what it made, it comes to a mesh it
+    /// finds nothing to do to, and hands that back as it was given — the
+    /// primitives and their order untouched, which is what lets a step of
+    /// a simulation at rest keep what it knows of the topology. Until
+    /// 2026-09-29 the flip pass, judging by valence alone, turned back the
+    /// long edges the split and collapse had just turned, and the three
+    /// went round on the same edges for ever.
+    #[test]
+    fn a_remesh_settles_and_then_leaves_the_mesh_alone() {
+        for (start, target) in [
+            (sphere_detail(Vec3::ZERO, 1.0, 6, 8), 0.2),
+            (sphere_detail(Vec3::ZERO, 1.0, 24, 32), 0.5),
+            (sphere_detail(Vec3::ZERO, 0.5, 10, 14), 0.1),
+        ] {
+            let settings = Settings { target, iterations: 3, relax: 0.0, ..Default::default() };
+            let mut d = remesh(&start, settings);
+            assert!(crate::remesh::last_changes() > 0, "the fixture remeshes");
+            let mut rounds = 0;
+            loop {
+                let next = remesh(&d, settings);
+                if crate::remesh::last_changes() == 0 {
+                    assert!(same_mesh(&next, &d), "a remesh with nothing to do changed the mesh");
+                    break;
+                }
+                d = next;
+                rounds += 1;
+                assert!(rounds < 20, "target {target}: still changing {} edges after {rounds} rounds", crate::remesh::last_changes());
+            }
+            // What it settled on is still the mesh that was asked for.
+            let mean = mean_edge(&d);
+            assert!((mean - target).abs() < target * 0.5, "settled at {mean}, wanted about {target}");
+            assert!(d.is_closed(), "and still a closed surface");
+        }
+    }
+
     #[test]
     fn test_remesh_converges_rather_than_oscillating() {
         // The 4/3 and 4/5 thresholds exist so a split cannot produce edges the
@@ -11906,6 +12005,7 @@ mod tests {
             let mut cache = crate::geometry::SimCache::default();
             let mut times = Vec::new();
             let mut points = Vec::new();
+            let mut remeshed = Vec::new();
             for frame in 1..=frames {
                 let mut sim = crate::geometry::EvalSim::new(frame, 1, &mut cache);
                 let mut err = None;
@@ -11913,15 +12013,17 @@ mod tests {
                 let d = crate::geometry::generate_single_node_geometry_with_errors(root, node, &mut Vec::new(), &mut err, &mut sim);
                 times.push(t.elapsed().as_secs_f64() * 1000.0);
                 points.push(d.map_or(0, |d| d.num_points()));
+                remeshed.push(crate::remesh::last_changes());
             }
             let mean = times.iter().sum::<f64>() / times.len() as f64;
             let worst = times.iter().cloned().fold(0.0, f64::max);
             let at = |f: usize| times.get(f - 1).copied().unwrap_or(0.0);
             println!(
-                "{:>22}: {mean:7.2} ms a frame, worst {worst:7.2}; frame 2 {:.2}, 10 {:.2}, 30 {:.2}, last {:.2}; points {} -> {}",
+                "{:>22}: {mean:7.2} ms a frame, worst {worst:7.2}; frame 2 {:.2}, 10 {:.2}, 30 {:.2}, last {:.2}; points {} -> {}; edges remeshed at frame 2 {}, 30 {}, last {}",
                 way.map_or("as saved".to_string(), |n| format!("without {n}")),
                 at(2), at(10), at(30), at(frames as usize),
                 points.first().unwrap(), points.last().unwrap(),
+                remeshed.get(1).copied().unwrap_or(0), remeshed.get(29).copied().unwrap_or(0), remeshed.last().copied().unwrap_or(0),
             );
         }
     }
diff --git a/src/remesh.rs b/src/remesh.rs
index 946657d..88a67cc 100644
--- a/src/remesh.rs
+++ b/src/remesh.rs
@@ -11,7 +11,8 @@
 //!
 //! 1. **Split** every edge longer than 4/3 of the target length.
 //! 2. **Collapse** every edge shorter than 4/5 of it.
-//! 3. **Flip** edges that would bring their four points closer to valence 6.
+//! 3. **Flip** edges that would bring their four points closer to valence 6,
+//!    unless the new edge would be one the next split cuts.
 //! 4. **Relax** each point toward the centroid of its neighbours, with the
 //!    normal component removed so the pass smooths the triangulation without
 //!    moving the surface.
@@ -472,7 +473,75 @@ fn face_normal(m: &Mesh, tri: [u32; 3]) -> Vec3 {
 /// "Better" is total deviation from valence 6, which is the valence a regular
 /// triangulation of a plane has — the measure the paper uses, and the one that
 /// drives a mesh toward equilateral triangles.
-fn flip_pass(m: &mut Mesh) -> usize {
+fn flip_pass(m: &mut Mesh, target: f32) -> usize {
+    let long = target * 4.0 / 3.0;
+    // Each point's valence, counted once and kept in step with the flips.
+    // The first version asked `tris_of` for it — a list gathered, sorted
+    // and counted — eight times an edge, which at 525 points was nine
+    // milliseconds of a remesh that changed nothing.
+    let mut valence = vec![0i32; m.pos.len()];
+    for (t, tri) in m.tris.iter().enumerate() {
+        if m.dead_tri[t] {
+            continue;
+        }
+        for (i, &q) in tri.iter().enumerate() {
+            if !tri[..i].contains(&q) {
+                valence[q as usize] += 1;
+            }
+        }
+    }
+    let mut done = 0;
+    for (e, _) in m.edges() {
+        let tris = m.tris_on_edge(e[0], e[1]);
+        if tris.len() != 2 {
+            continue; // a boundary edge has nothing to flip into
+        }
+        let (t0, t1) = (tris[0], tris[1]);
+        let Some(&o0) = m.tris[t0].iter().find(|q| !e.contains(q)) else { continue };
+        let Some(&o1) = m.tris[t1].iter().find(|q| !e.contains(q)) else { continue };
+        if o0 == o1 {
+            continue;
+        }
+        let val = |p: u32| valence[p as usize];
+        let dev = |v: i32| (v - 6).abs();
+        let before = dev(val(e[0])) + dev(val(e[1])) + dev(val(o0)) + dev(val(o1));
+        let after = dev(val(e[0]) - 1) + dev(val(e[1]) - 1) + dev(val(o0) + 1) + dev(val(o1) + 1);
+        if after >= before {
+            continue;
+        }
+        // Refuse a flip whose new edge the next split would cut. A long edge
+        // between two thin triangles is split, and the midpoint collapsed
+        // into a corner — which is that edge turned to its short diagonal —
+        // and until 2026-09-29 this pass, judging by valence alone, turned
+        // it back: three passes undoing one another on the same edges at
+        // every step of a mesh that had stopped moving.
+        if (m.pos[o1 as usize] - m.pos[o0 as usize]).length() > long {
+            continue;
+        }
+        let n0 = face_normal(m, m.tris[t0]);
+        let (new0, new1) = ([o0, e[0], o1], [o1, e[1], o0]);
+        if face_normal(m, new0).dot(n0) <= 0.0 || face_normal(m, new1).dot(n0) <= 0.0 {
+            continue;
+        }
+        m.tris[t0] = new0;
+        m.tris[t1] = new1;
+        for &q in new0.iter().chain(new1.iter()) {
+            m.p2t[q as usize].push(t0);
+            m.p2t[q as usize].push(t1);
+        }
+        valence[e[0] as usize] -= 1;
+        valence[e[1] as usize] -= 1;
+        valence[o0 as usize] += 1;
+        valence[o1 as usize] += 1;
+        done += 1;
+    }
+    done
+}
+
+/// [`flip_pass`] as it was first written, kept to hold the faster one to.
+#[cfg(test)]
+fn flip_pass_reference(m: &mut Mesh, target: f32) -> usize {
+    let long = target * 4.0 / 3.0;
     let mut done = 0;
     for (e, _) in m.edges() {
         // Looked up now rather than taken from the snapshot, for the same
@@ -497,6 +566,15 @@ fn flip_pass(m: &mut Mesh) -> usize {
         if after >= before {
             continue;
         }
+        // Refuse a flip whose new edge the next split would cut. A long edge
+        // between two thin triangles is split, and the midpoint collapsed
+        // into a corner — which is that edge turned to its short diagonal —
+        // and until 2026-09-29 this pass, judging by valence alone, turned
+        // it back: three passes undoing one another on the same edges at
+        // every step of a mesh that had stopped moving.
+        if (m.pos[o1 as usize] - m.pos[o0 as usize]).length() > long {
+            continue;
+        }
         // Refuse a flip that would fold either new triangle against the
         // surface it came from.
         let n0 = face_normal(m, m.tris[t0]);
@@ -575,15 +653,20 @@ fn relax_pass(m: &mut Mesh, amount: f32) {
 /// first order: on anything curved the slide leaves the surface slightly, and
 /// the error compounds. Without this a sphere remeshed for fifty iterations is
 /// visibly smaller than the one it started as.
-fn project_pass(m: &mut Mesh, rest: &crate::spatial::TriGrid) {
+fn project_pass(m: &mut Mesh, rest: &crate::spatial::TriGrid, reference: bool) {
     if rest.is_empty() {
         return;
     }
+    let _ = reference;
     for p in 0..m.pos.len() {
         if m.dead_point[p] {
             continue;
         }
-        if let Some(hit) = rest.closest(m.pos[p]) {
+        #[cfg(test)]
+        let hit = if reference { rest.closest_reference(m.pos[p]) } else { rest.closest(m.pos[p]) };
+        #[cfg(not(test))]
+        let hit = rest.closest(m.pos[p]);
+        if let Some(hit) = hit {
             m.pos[p] = hit.point;
         }
     }
@@ -647,6 +730,28 @@ pub fn subdivide(input: &Detail, depth: usize) -> Detail {
 
 /// Remesh toward `settings.target` edge length.
 pub fn remesh(input: &Detail, settings: Settings) -> Detail {
+    remesh_by(input, settings, false)
+}
+
+/// [`remesh`] by the passes as they were first written: the same mesh, bit
+/// for bit, which is what makes the faster ones optimizations.
+#[cfg(test)]
+pub fn remesh_reference(input: &Detail, settings: Settings) -> Detail {
+    remesh_by(input, settings, true)
+}
+
+thread_local! {
+    /// How many edges the last remesh on this thread split, collapsed or
+    /// flipped: what a test and the profile read to see a mesh settle.
+    static LAST_CHANGES: std::cell::Cell<usize> = const { std::cell::Cell::new(0) };
+}
+
+/// The edges the last [`remesh`] on this thread split, collapsed or flipped.
+pub fn last_changes() -> usize {
+    LAST_CHANGES.with(|c| c.get())
+}
+
+fn remesh_by(input: &Detail, settings: Settings, reference: bool) -> Detail {
     if input.num_prims() == 0 || settings.target <= 0.0 {
         return input.clone();
     }
@@ -654,22 +759,51 @@ pub fn remesh(input: &Detail, settings: Settings) -> Detail {
     // the previous iteration's surface would chase the creep rather than
     // correct it, since each iteration's drift would become the next one's
     // idea of where the surface is.
-    let rest = settings.project.then(|| crate::spatial::TriGrid::build(input));
+    let mut rest: Option<crate::spatial::TriGrid> = None;
     let mut m = Mesh::from_detail(input);
-    for _ in 0..settings.iterations.min(20) {
+    let mut changed = 0;
+    LAST_CHANGES.with(|c| c.set(0));
+    for iteration in 0..settings.iterations.min(20) {
+        let mut done = 0;
         if settings.split {
-            split_pass(&mut m, settings.target);
+            done += split_pass(&mut m, settings.target);
         }
         if settings.collapse {
-            collapse_pass(&mut m, settings.target);
+            done += collapse_pass(&mut m, settings.target);
         }
         if settings.flip {
-            flip_pass(&mut m);
+            #[cfg(test)]
+            if reference {
+                done += flip_pass_reference(&mut m, settings.target);
+            } else {
+                done += flip_pass(&mut m, settings.target);
+            }
+            #[cfg(not(test))]
+            {
+                done += flip_pass(&mut m, settings.target);
+            }
+        }
+        // Nothing to split, collapse or flip, and nothing relaxed: the mesh
+        // is as this iteration found it, and every later one would find it
+        // the same. On the FIRST that is the input itself, handed back as it
+        // came — its primitives, their order and its polygons untouched, so
+        // what is kept by a mesh's topology downstream (the detangle's
+        // lists) is kept across a step that remeshed nothing.
+        if done == 0 && settings.relax <= 0.0 {
+            if iteration == 0 {
+                return input.clone();
+            }
+            break;
         }
+        changed += done;
         relax_pass(&mut m, settings.relax.clamp(0.0, 1.0));
-        if let Some(rest) = &rest {
-            project_pass(&mut m, rest);
+        if settings.project {
+            // Built when first wanted: a remesh that finds nothing to do
+            // never asks.
+            let rest = rest.get_or_insert_with(|| crate::spatial::TriGrid::build(input));
+            project_pass(&mut m, rest, reference);
         }
     }
+    LAST_CHANGES.with(|c| c.set(changed));
     m.into_detail()
 }
diff --git a/src/spatial.rs b/src/spatial.rs
index 2b913aa..02332a4 100644
--- a/src/spatial.rs
+++ b/src/spatial.rs
@@ -337,6 +337,70 @@ impl TriGrid {
                 })
                 .filter(|h| h.distance <= limit);
         }
+        // The nearest of the triangles named, and of those equally near the
+        // one with the lowest index — which is what the first version's
+        // `min_by` over a sorted list chose. Only the winner's normal is
+        // worked out: it was a cross product and a square root for every
+        // candidate, and a candidate is a triangle for every cell it is
+        // filed in.
+        let nearest = |among: &mut dyn Iterator<Item = u32>| -> Option<(f32, u32, Vec3)> {
+            let mut best: Option<(f32, u32, Vec3)> = None;
+            for i in among {
+                if best.is_some_and(|b| b.1 == i) {
+                    continue;
+                }
+                let t = self.tris[i as usize];
+                let q = closest_point_on_triangle(p, t[0], t[1], t[2]);
+                let d = (q - p).length();
+                let better = match best {
+                    None => true,
+                    Some((bd, bi, _)) => match d.partial_cmp(&bd) {
+                        Some(std::cmp::Ordering::Less) => true,
+                        Some(std::cmp::Ordering::Greater) => false,
+                        _ => i < bi,
+                    },
+                };
+                if better {
+                    best = Some((d, i, q));
+                }
+            }
+            best
+        };
+        let hit = |(distance, i, point): (f32, u32, Vec3)| {
+            let t = self.tris[i as usize];
+            Hit { point, distance, normal: (t[1] - t[0]).cross(t[2] - t[0]).normalize_or_zero() }
+        };
+
+        // From a quarter of a cell, where the first version began at a
+        // whole one: a query ON the surface — a remesh projecting its
+        // points, which is most of what is asked — is answered from the
+        // cell it is in and not from the twenty-seven about it. Where the
+        // search begins does not change what it finds: it ends only with a
+        // hit nearer than the box searched is wide, which nothing outside
+        // the box can beat or tie.
+        let mut reach = self.grid.cell * 0.25;
+        for _ in 0..14 {
+            let (a, b) = (self.grid.coord(p - Vec3::splat(reach)), self.grid.coord(p + Vec3::splat(reach)));
+            let mut among = (a[2]..=b[2]).flat_map(|z| (a[1]..=b[1]).map(move |y| (y, z))).flat_map(|(y, z)| {
+                (a[0]..=b[0]).flat_map(move |x| self.grid.buckets[self.grid.index([x, y, z])].iter().copied())
+            });
+            match nearest(&mut among) {
+                Some(h) if h.0 <= reach => return Some(hit(h)),
+                _ => reach *= 2.0,
+            }
+        }
+        // The grid is clamped to the geometry's bounds, so the doublings
+        // have gathered everything; whatever it found is the answer.
+        nearest(&mut (0..self.tris.len() as u32)).map(hit)
+    }
+
+    /// [`TriGrid::closest`] as it was first written, kept to hold the
+    /// faster one to: the same hit, bit for bit.
+    #[cfg(test)]
+    pub fn closest_reference(&self, p: Vec3) -> Option<Hit> {
+        if self.tris.is_empty() {
+            return None;
+        }
         let hit = |i: usize| {
             let t = self.tris[i];
             let q = closest_point_on_triangle(p, t[0], t[1], t[2]);
@@ -349,7 +413,6 @@ impl TriGrid {
         let nearer = |a: &Hit, b: &Hit| {
             a.distance.partial_cmp(&b.distance).unwrap_or(std::cmp::Ordering::Equal)
         };
-
         let mut reach = self.grid.cell;
         let mut scratch = Vec::new();
         for _ in 0..12 {
@@ -361,8 +424,6 @@ impl TriGrid {
                 _ => reach *= 2.0,
             }
         }
-        // The grid is clamped to the geometry's bounds, so twelve doublings
-        // have gathered everything; whatever it found is the answer.
         (0..self.tris.len()).map(hit).min_by(nearer)
     }
 }