graphic design tool
git clone https://git.lucas.co/cce-designer.git
feat(phase3): collision, and the spatial index three operators wanted
Detangle pushes a surface off itself: a point repulsion over Iterations
passes, Thickness in edge lengths, points within Rings of each other excluded
— the algorithm otls-audit.md records for developer_surface_detangle.
The ring exclusion is the whole trick. Every point is within a thickness of
its own neighbours by construction, because that is what an edge is, so a
naive repulsion blows the mesh apart from the inside. Excluding the
topological neighbourhood leaves exactly the pairs that are near in SPACE and
far across the SURFACE, which is what a self-intersection is. Thickness is in
edge lengths so the setting survives a remesh; an absolute distance would stop
separating the moment the mesh got finer.
Suture resolves against a second input and fuses what keeps touching. A point
within Distance Threshold is pushed back out and its contact counter goes up;
one that is clear has its counter reset, because contact has to be SUSTAINED
to count — that is the difference between two surfaces brushing past each
other and two surfaces growing into each other. Past Fusion Threshold, points
within Distance Threshold of each other become one.
Both run on src/spatial.rs, built once because three operators arrived wanting
it at the same moment: a uniform grid over triangles for "what surface point
is nearest", one over points for "what is within this radius". The third
consumer is the remesher's projection pass, which was the piece the last
commit left out — and it now holds a relaxed surface on the mesh it was given.
Three things the tests caught:
Suture could never fuse. Resolving a contact pushed the point to exactly the
threshold, and the next pass compared strictly, so every resolved contact read
as released and no counter could reach the fusion threshold. Contact is
unsustainable by construction unless the comparison is inclusive.
fuse_points left folded primitives behind. The point compaction only knows
about points that went away; a triangle whose corners merged onto each other
still references three live points and has to be dropped where the merge
happens.
And the projection test was measuring the wrong thing. Projection holds points
on the INPUT surface — a faceted sphere, whose edge midpoints are legitimately
inside the ideal one — so radius measures the discretization, not the drift.
It now measures distance from the input surface, where the difference is a
factor of forty.
NOT ported, and not guessed at: developer_surface_adapt (eight read parameters,
no description anywhere) and developer_surface_open (no parameters, no
description). Unlike Charge, whose parameter names carried a reading, these
carry nothing.
Co-Authored-By: Claude Opus 5 <[email protected]>
nodes/detangle.json | 13 +++
nodes/remesh.json | 56 +++++++++--
nodes/suture.json | 13 +++
shapeshifter.md | 18 +++-
src/detail.rs | 63 ++++++++++++
src/geometry.rs | 273 ++++++++++++++++++++++++++++++++++++++++++++++++++++
src/main.rs | 262 +++++++++++++++++++++++++++++++++++++++++++++++++
src/remesh.rs | 50 ++++++++--
src/spatial.rs | 246 ++++++++++++++++++++++++++++++++++++++++++++++
9 files changed, 978 insertions(+), 16 deletions(-)
diff --git a/nodes/detangle.json b/nodes/detangle.json
new file mode 100644
index 0000000..8562b59
--- /dev/null
+++ b/nodes/detangle.json
@@ -0,0 +1,13 @@
+{
+ "name": "Detangle",
+ "type": "detangle",
+ "inputs": 1,
+ "outputs": 1,
+ "params": [
+ { "name": "Input", "type": "text", "default": "" },
+ { "name": "Thickness", "type": "slider", "default": "1.00", "min": 0.0, "max": 8.0, "step": 0.05 },
+ { "name": "Rings", "type": "spinbox", "default": "2", "min": 0.0, "max": 6.0, "step": 1.0 },
+ { "name": "Iterations", "type": "spinbox", "default": "4", "min": 1.0, "max": 32.0, "step": 1.0 },
+ { "name": "Group", "type": "text", "default": "" }
+ ]
+}
diff --git a/nodes/remesh.json b/nodes/remesh.json
index 9b62608..d454ccd 100644
--- a/nodes/remesh.json
+++ b/nodes/remesh.json
@@ -4,12 +4,54 @@
"inputs": 1,
"outputs": 1,
"params": [
- { "name": "Input", "type": "text", "default": "" },
- { "name": "Target Length", "type": "slider", "default": "0.10", "min": 0.001, "max": 1.0, "step": 0.005 },
- { "name": "Iterations", "type": "spinbox", "default": "3", "min": 1.0, "max": 20.0, "step": 1.0 },
- { "name": "Relax", "type": "slider", "default": "0.50", "min": 0.0, "max": 1.0, "step": 0.01 },
- { "name": "Split", "type": "choice:true,false", "default": "true" },
- { "name": "Collapse", "type": "choice:true,false", "default": "true" },
- { "name": "Flip", "type": "choice:true,false", "default": "true" }
+ {
+ "name": "Input",
+ "type": "text",
+ "default": ""
+ },
+ {
+ "name": "Target Length",
+ "type": "slider",
+ "default": "0.10",
+ "min": 0.001,
+ "max": 1.0,
+ "step": 0.005
+ },
+ {
+ "name": "Iterations",
+ "type": "spinbox",
+ "default": "3",
+ "min": 1.0,
+ "max": 20.0,
+ "step": 1.0
+ },
+ {
+ "name": "Relax",
+ "type": "slider",
+ "default": "0.50",
+ "min": 0.0,
+ "max": 1.0,
+ "step": 0.01
+ },
+ {
+ "name": "Split",
+ "type": "choice:true,false",
+ "default": "true"
+ },
+ {
+ "name": "Collapse",
+ "type": "choice:true,false",
+ "default": "true"
+ },
+ {
+ "name": "Flip",
+ "type": "choice:true,false",
+ "default": "true"
+ },
+ {
+ "name": "Project",
+ "type": "choice:true,false",
+ "default": "true"
+ }
]
}
diff --git a/nodes/suture.json b/nodes/suture.json
new file mode 100644
index 0000000..394ae60
--- /dev/null
+++ b/nodes/suture.json
@@ -0,0 +1,13 @@
+{
+ "name": "Suture",
+ "type": "suture",
+ "inputs": 2,
+ "outputs": 1,
+ "params": [
+ { "name": "Input", "type": "text", "default": "" },
+ { "name": "Against", "type": "text", "default": "" },
+ { "name": "Distance Threshold", "type": "slider", "default": "0.05", "min": 0.0, "max": 1.0, "step": 0.005 },
+ { "name": "Fusion Threshold", "type": "spinbox", "default": "3", "min": 1.0, "max": 100.0, "step": 1.0 },
+ { "name": "Counter", "type": "text", "default": "contact" }
+ ]
+}
diff --git a/shapeshifter.md b/shapeshifter.md
index 0f7aa68..ee37425 100644
--- a/shapeshifter.md
+++ b/shapeshifter.md
@@ -231,9 +231,21 @@ Touches: `geometry.rs`, `app.rs`, `render.rs`, `playbar.rs`.
> sides, so nothing looked wrong. Develop is the first operator whose answer
> depends on it, and it grew the surface inward.
>
-> Outstanding: `adapt`, `open`, and collision (`detangle`, `suture`) — which
-> want the spatial index the remesher's missing surface-projection pass also
-> needs.
+> Collision is in. `detangle` is a point repulsion over Iterations passes with
+> Thickness in edge lengths and a Rings exclusion — the audit's own
+> description of the Shapeshifter algorithm. `suture` resolves against a second
+> input, counts SUSTAINED contact, and fuses points past Fusion Threshold
+> within Distance Threshold. Both, and the remesher's projection pass, run on
+> `src/spatial.rs` — a uniform grid built once for all three.
+>
+> **`adapt` and `open` are not ported and should not be guessed at.**
+> `developer_surface_adapt` has eight read parameters and no description
+> anywhere; `developer_surface_open` has none at all. Unlike Charge, whose
+> parameter names carried a reading, these carry nothing. Say what they do and
+> they are a short job each.
+>
+> Remaining: `subdivide`, which is the one documented Surface operator still
+> missing and is the easy one.
`Develop` is easy — displace along the normal by a development attribute.
**Remesh is the hard one**, and it is load-bearing: without topology that keeps
diff --git a/src/detail.rs b/src/detail.rs
index 3eefe78..a872102 100644
--- a/src/detail.rs
+++ b/src/detail.rs
@@ -1334,6 +1334,69 @@ impl Detail {
self.invalidate();
}
+ /// Merge points onto representatives: point `p` becomes `rep[p]`.
+ ///
+ /// A representative keeps its identity and its values — the same choice
+ /// the remesher's collapse makes, and for the same reason: one of the two
+ /// is a point the solver has been writing to, and the merge should cost
+ /// the simulation as little memory as it can. Primitives are rewired, and
+ /// one left with a repeated corner is dropped, because a triangle with two
+ /// corners in the same place is not a triangle.
+ ///
+ /// Chains are followed, so `rep` need not already be flat: a fuse that
+ /// pointed a at b and b at c leaves everything at c.
+ pub fn fuse_points(&mut self, rep: &[u32]) {
+ let n = self.num_points();
+ let root = |mut p: u32| {
+ // Bounded rather than trusting the map to be acyclic: a cycle in a
+ // caller's representative map would otherwise hang the app.
+ for _ in 0..n {
+ let next = rep.get(p as usize).copied().unwrap_or(p);
+ if next == p {
+ break;
+ }
+ p = next;
+ }
+ p
+ };
+ for v in self.vert_point.iter_mut() {
+ *v = root(*v);
+ }
+
+ // A primitive whose corners collapsed onto each other is not a
+ // primitive any more. Dropped here rather than left for the point
+ // compaction, which only knows about points that went away — these
+ // ones all still exist, they have just stopped being distinct.
+ let mut vert_point = Vec::with_capacity(self.vert_point.len());
+ let mut prim_start = vec![0u32];
+ let mut kept_prims: Vec<u32> = Vec::new();
+ let mut kept_verts: Vec<u32> = Vec::new();
+ for prim in 0..self.num_prims() {
+ let range = self.prim_verts(prim);
+ let pts = &self.vert_point[range.clone()];
+ let mut uniq = pts.to_vec();
+ uniq.sort_unstable();
+ uniq.dedup();
+ if uniq.len() < pts.len() || uniq.len() < 3 {
+ continue;
+ }
+ for v in range {
+ kept_verts.push(v as u32);
+ vert_point.push(self.vert_point[v]);
+ }
+ prim_start.push(vert_point.len() as u32);
+ kept_prims.push(prim as u32);
+ }
+ self.verts = self.verts.gather(&kept_verts);
+ self.prims = self.prims.gather(&kept_prims);
+ self.vert_point = vert_point;
+ self.prim_start = prim_start;
+
+ let keep: Vec<bool> = (0..n).map(|p| root(p as u32) as usize == p).collect();
+ self.invalidate();
+ self.keep_points(&keep);
+ }
+
/// Keep the points `keep` marks true, dropping the rest.
pub fn keep_points(&mut self, keep: &[bool]) {
let idx: Vec<u32> = (0..self.num_points() as u32)
diff --git a/src/geometry.rs b/src/geometry.rs
index e382ff6..ffe2195 100644
--- a/src/geometry.rs
+++ b/src/geometry.rs
@@ -740,6 +740,10 @@ pub fn generate_single_node_geometry_with_errors(
resolve_neighbour_geometry_with_errors(root, target, visited, ocl_error, sim)
} else if target.node_type.eq_ignore_ascii_case("time") {
resolve_time_geometry_with_errors(root, target, visited, ocl_error, sim)
+ } else if target.node_type.eq_ignore_ascii_case("detangle") {
+ resolve_detangle_geometry_with_errors(root, target, visited, ocl_error, sim)
+ } else if target.node_type.eq_ignore_ascii_case("suture") {
+ resolve_suture_geometry_with_errors(root, target, visited, ocl_error, sim)
} else if target.node_type.eq_ignore_ascii_case("remesh") {
resolve_remesh_geometry_with_errors(root, target, visited, ocl_error, sim)
} else if target.node_type.eq_ignore_ascii_case("develop") {
@@ -1292,6 +1296,254 @@ pub fn resolve_relax_geometry_with_errors(
Some(geom)
}
+/// The Detangle node: push a surface off itself.
+///
+/// Growth folds a surface into its own neighbourhood long before it looks
+/// wrong from outside, and a diffusion across a self-intersecting mesh reads
+/// neighbours that are topologically far away. This resolves it the way the
+/// plugin's Detangle does: a point repulsion over Iterations passes, with
+/// Thickness measured in edge lengths and points within Rings of each other
+/// excluded.
+///
+/// The ring exclusion is the whole trick. Every point is within a thickness of
+/// its own neighbours by construction — that is what an edge is — so a naive
+/// repulsion would blow the mesh apart from the inside. Excluding the
+/// topological neighbourhood leaves exactly the pairs that are near in SPACE
+/// but far across the SURFACE, which is what a self-intersection is.
+pub fn resolve_detangle_geometry_with_errors(
+ root: &FsNode,
+ target: &FsNode,
+ visited: &mut Vec<String>,
+ ocl_error: &mut Option<String>,
+ sim: &mut EvalSim,
+) -> Option<Detail> {
+ let input_name = node_param_str(target, "Input", "");
+ if input_name.is_empty() {
+ return None;
+ }
+ let input_node = find_node_by_name(root, &input_name)?;
+ let mut geom = generate_single_node_geometry_with_errors(root, input_node, visited, ocl_error, sim)?;
+ apply_detangle(&mut geom, target);
+ Some(geom)
+}
+
+pub(crate) fn apply_detangle(geom: &mut Detail, target: &FsNode) {
+ let n = geom.num_points();
+ if n == 0 || geom.num_prims() == 0 {
+ return;
+ }
+ // Thickness in EDGE LENGTHS, so the setting means the same thing before
+ // and after a remesh — an absolute distance would stop separating the
+ // moment the mesh got finer.
+ let edges = geom.edges().to_vec();
+ if edges.is_empty() {
+ return;
+ }
+ let mean_edge = edges
+ .iter()
+ .map(|e| (geom.pos(e[1] as usize) - geom.pos(e[0] as usize)).length())
+ .sum::<f32>()
+ / edges.len() as f32;
+ let thickness = node_param_f32(target, "Thickness", 1.0).max(0.0) * mean_edge;
+ if thickness <= 0.0 {
+ return;
+ }
+ let rings = node_param_f32(target, "Rings", 2.0).clamp(0.0, 6.0) as usize;
+ let iterations = node_param_f32(target, "Iterations", 4.0).clamp(1.0, 32.0) as usize;
+ let group = node_param_str(target, "Group", "");
+ let group = group.trim().to_string();
+ let movable: Vec<bool> = (0..n)
+ .map(|p| group.is_empty() || geom.points().in_group(&group, p))
+ .collect();
+
+ // The excluded neighbourhood, once: it is topology, and the pass does not
+ // change topology.
+ let excluded: Vec<Vec<u32>> = (0..n)
+ .map(|p| {
+ let mut seen = vec![p as u32];
+ let mut frontier = vec![p as u32];
+ for _ in 0..rings {
+ let mut next = Vec::new();
+ for &q in &frontier {
+ for &r in geom.point_neighbours(q as usize) {
+ if !seen.contains(&r) {
+ seen.push(r);
+ next.push(r);
+ }
+ }
+ }
+ if next.is_empty() {
+ break;
+ }
+ frontier = next;
+ }
+ seen.sort_unstable();
+ seen
+ })
+ .collect();
+
+ let mut pos: Vec<Vec3> = (0..n).map(|p| geom.pos(p)).collect();
+ let mut near = Vec::new();
+ for _ in 0..iterations {
+ let grid = crate::spatial::PointGrid::build(&pos, thickness);
+ // Gathered against the positions at the START of the pass and applied
+ // at the end, so the result does not depend on the order points are
+ // visited in — a Gauss-Seidel sweep here would make the same mesh
+ // untangle differently depending on how its points were numbered.
+ let mut push = vec![Vec3::ZERO; n];
+ for p in 0..n {
+ grid.within(pos[p], thickness, &mut near);
+ for &q in &near {
+ let q = q as usize;
+ if q == p || excluded[p].binary_search(&(q as u32)).is_ok() {
+ continue;
+ }
+ let d = pos[p] - pos[q];
+ let len = d.length();
+ if len >= thickness {
+ continue;
+ }
+ // Two coincident points have no direction to separate along;
+ // nudging along an arbitrary axis at least breaks the tie.
+ let dir = if len < 1e-9 {
+ Vec3::new((p % 7) as f32 - 3.0, (p % 5) as f32 - 2.0, 1.0).normalize_or_zero()
+ } else {
+ d / len
+ };
+ push[p] += dir * ((thickness - len) * 0.5);
+ }
+ }
+ for p in 0..n {
+ if movable[p] {
+ pos[p] += push[p];
+ }
+ }
+ }
+ for (p, v) in pos.iter().enumerate() {
+ geom.set_pos(p, *v);
+ }
+}
+
+/// The Suture node: resolve a surface against another, and fuse what keeps
+/// touching.
+///
+/// Two thresholds, as the plugin has them. A point closer to the `Against`
+/// geometry than Distance Threshold is in contact: it is pushed back out to
+/// that distance, and its contact counter goes up. A point that is NOT in
+/// contact has its counter reset to zero — contact has to be sustained to
+/// count, which is the difference between two surfaces brushing past each
+/// other and two surfaces growing into each other.
+///
+/// Once a point's counter passes Fusion Threshold, it fuses with any other
+/// such point within Distance Threshold: they become one point, keeping the
+/// lower index's identity and values. That is the operation the name is
+/// about — a surface that has been pressed against itself for long enough
+/// stops being two surfaces.
+pub fn resolve_suture_geometry_with_errors(
+ root: &FsNode,
+ target: &FsNode,
+ visited: &mut Vec<String>,
+ ocl_error: &mut Option<String>,
+ sim: &mut EvalSim,
+) -> Option<Detail> {
+ let input_name = node_param_str(target, "Input", "");
+ if input_name.is_empty() {
+ return None;
+ }
+ let input_node = find_node_by_name(root, &input_name)?;
+ let mut geom = generate_single_node_geometry_with_errors(root, input_node, visited, ocl_error, sim)?;
+
+ let against_name = node_param_str(target, "Against", "");
+ let against_name = against_name.trim().to_string();
+ let against = if against_name.is_empty() {
+ None
+ } else {
+ find_node_by_name(root, &against_name)
+ .and_then(|n| generate_single_node_geometry_with_errors(root, n, visited, ocl_error, sim))
+ };
+ apply_suture(&mut geom, against.as_ref(), target);
+ Some(geom)
+}
+
+pub(crate) fn apply_suture(geom: &mut Detail, against: Option<&Detail>, target: &FsNode) {
+ let n = geom.num_points();
+ if n == 0 {
+ return;
+ }
+ let distance = node_param_f32(target, "Distance Threshold", 0.05).max(0.0);
+ let fusion = node_param_f32(target, "Fusion Threshold", 3.0).max(1.0) as i32;
+ let counter = node_param_str(target, "Counter", "contact");
+ let counter = counter.trim().to_string();
+ if counter.is_empty() || distance <= 0.0 {
+ return;
+ }
+
+ // The counter is LIVE: it is the memory that makes sustained contact
+ // different from a brush, and a derivative one would reset every step and
+ // never reach the threshold.
+ geom.points_mut()
+ .get_or_create(&counter, AttribValue::Int(0));
+
+ let mut counts: Vec<i32> = (0..n)
+ .map(|p| geom.points().value(&counter, p).map(|v| v.as_f32() as i32).unwrap_or(0))
+ .collect();
+
+ if let Some(other) = against.filter(|o| o.num_prims() > 0) {
+ let grid = crate::spatial::TriGrid::build(other);
+ for p in 0..n {
+ let here = geom.pos(p);
+ let Some((closest, dist)) = grid.closest(here) else { continue };
+ // Inclusive, with room for the float error: a point resolved to
+ // exactly the threshold last step is STILL in contact this step.
+ // Comparing strictly would have every resolved contact read as
+ // released on the next pass, and no counter could ever reach the
+ // fusion threshold — contact would be unsustainable by
+ // construction.
+ if dist > distance * (1.0 + 1e-3) {
+ counts[p] = 0;
+ continue;
+ }
+ counts[p] += 1;
+ // Pushed back out along the line to the surface. A point sitting
+ // exactly on it has no such line, and is left where it is rather
+ // than shoved in an invented direction.
+ let away = here - closest;
+ if away.length_squared() > 1e-12 {
+ geom.set_pos(p, closest + away.normalize() * distance);
+ }
+ }
+ }
+
+ for p in 0..n {
+ let _ = geom
+ .points_mut()
+ .set_value(&counter, p, AttribValue::Int(counts[p]));
+ }
+
+ // Fuse the sustained contacts that are near each other. Lowest index wins,
+ // so the result does not depend on visit order.
+ let welded: Vec<usize> = (0..n).filter(|&p| counts[p] >= fusion).collect();
+ if welded.len() < 2 {
+ return;
+ }
+ let pos: Vec<Vec3> = welded.iter().map(|&p| geom.pos(p)).collect();
+ let grid = crate::spatial::PointGrid::build(&pos, distance);
+ let mut rep: Vec<u32> = (0..n as u32).collect();
+ let mut near = Vec::new();
+ for (i, &p) in welded.iter().enumerate() {
+ grid.within(pos[i], distance, &mut near);
+ for &j in &near {
+ let q = welded[j as usize];
+ if q > p {
+ // Only ever point a higher index at a lower one, so the map
+ // cannot contain a cycle.
+ rep[q] = rep[q].min(p as u32);
+ }
+ }
+ }
+ geom.fuse_points(&rep);
+}
+
/// The Remesh node: keep the triangulation proportional to the surface.
///
/// The counterweight to Develop. Growth pushes points apart and the triangles
@@ -1327,6 +1579,7 @@ pub(crate) fn remesh_settings(target: &FsNode) -> crate::remesh::Settings {
split: node_param_str(target, "Split", "true") == "true",
collapse: node_param_str(target, "Collapse", "true") == "true",
flip: node_param_str(target, "Flip", "true") == "true",
+ project: node_param_str(target, "Project", "true") == "true",
}
}
@@ -3684,6 +3937,8 @@ pub fn is_geometry_node_type(node_type: &str) -> bool {
|| nt == "visualize"
|| nt == "develop"
|| nt == "remesh"
+ || nt == "suture"
+ || nt == "detangle"
|| nt == "attribute"
|| nt == "simnet"
}
@@ -3831,6 +4086,24 @@ pub fn network_sphere_vertices_with_errors(
out.merge(&geom);
}
}
+ } else if node.node_type.eq_ignore_ascii_case("detangle") {
+ let _idx = *count;
+ *count += 1;
+ if is_visible {
+ let mut visited = Vec::new();
+ if let Some(geom) = resolve_detangle_geometry_with_errors(root, node, &mut visited, ocl_error, sim) {
+ out.merge(&geom);
+ }
+ }
+ } else if node.node_type.eq_ignore_ascii_case("suture") {
+ let _idx = *count;
+ *count += 1;
+ if is_visible {
+ let mut visited = Vec::new();
+ if let Some(geom) = resolve_suture_geometry_with_errors(root, node, &mut visited, ocl_error, sim) {
+ out.merge(&geom);
+ }
+ }
} else if node.node_type.eq_ignore_ascii_case("remesh") {
let _idx = *count;
*count += 1;
diff --git a/src/main.rs b/src/main.rs
index 6484064..bd943fb 100644
--- a/src/main.rs
+++ b/src/main.rs
@@ -4,6 +4,7 @@ pub mod application;
pub mod curve_tool;
pub mod detail;
pub mod remesh;
+pub mod spatial;
// Root-level aliases some modules import via `crate::` paths.
#[allow(unused_imports)]
@@ -3596,6 +3597,267 @@ mod tests {
assert_eq!(out.ids(), sphere.ids());
}
+ fn phase3_node(ty: &str, params: &[(&str, &str)]) -> FsNode {
+ FsNode {
+ id: format!("id-{ty}"),
+ name: format!("{ty} 1"),
+ node_type: ty.into(),
+ children: vec![],
+ params: params
+ .iter()
+ .map(|(name, default)| crate::app::ParamDef {
+ name: (*name).into(),
+ label: String::new(),
+ param_type: "text".into(),
+ default: (*default).into(),
+ options: vec![],
+ min: None,
+ max: None,
+ step: None,
+ })
+ .collect(),
+ geometry_visible: true,
+ position: (0.0, 0.0),
+ inputs: 1,
+ outputs: 1,
+ }
+ }
+
+ #[test]
+ fn test_the_projection_pass_stops_a_remeshed_surface_creeping() {
+ use crate::spatial::TriGrid;
+ // Tangential relaxation slides points within the surface, but "within"
+ // is only true to first order: on anything curved the slide leaves the
+ // surface a little, and the error compounds.
+ //
+ // Measured as distance from the INPUT SURFACE, not as radius. The
+ // projection holds points on the mesh it was given — which is a
+ // faceted sphere, whose edge midpoints are legitimately inside the
+ // ideal one. Radius would be measuring the discretization, not the
+ // drift.
+ let sphere = sphere_detail(Vec3::ZERO, 1.0, 10, 14);
+ let rest = TriGrid::build(&sphere);
+ let drift = |d: &Detail| {
+ (0..d.num_points())
+ .filter_map(|p| rest.closest(d.pos(p)).map(|(_, dist)| dist))
+ .sum::<f32>()
+ / d.num_points() as f32
+ };
+
+ let cfg = |project| Settings {
+ target: 0.25,
+ iterations: 16,
+ relax: 1.0,
+ project,
+ ..Default::default()
+ };
+ let drifted = remesh(&sphere, cfg(false));
+ let held = remesh(&sphere, cfg(true));
+
+ assert!(drift(&drifted) > 1e-3, "relaxation did not drift at all: {}", drift(&drifted));
+ assert!(
+ drift(&held) < drift(&drifted) / 4.0,
+ "projection barely helped: {} vs {}",
+ drift(&held),
+ drift(&drifted)
+ );
+ assert!(drift(&held) < 1e-4, "projected points are off the surface: {}", drift(&held));
+ }
+
+ #[test]
+ fn test_the_tri_grid_finds_the_nearest_surface_point() {
+ use crate::spatial::TriGrid;
+ let sphere = sphere_detail(Vec3::ZERO, 1.0, 12, 16);
+ let grid = TriGrid::build(&sphere);
+
+ // A point well outside: the closest surface point is along the ray to
+ // the centre, at about the radius.
+ let (q, d) = grid.closest(Vec3::new(3.0, 0.0, 0.0)).unwrap();
+ assert!((d - 2.0).abs() < 0.05, "distance {d}");
+ assert!(q.x > 0.9 && q.y.abs() < 0.2 && q.z.abs() < 0.2, "{q:?}");
+
+ // A point ON the surface finds itself.
+ let on = sphere.pos(20);
+ let (_, d) = grid.closest(on).unwrap();
+ assert!(d < 1e-3, "a point on the surface is {d} from it");
+
+ // The centre is a radius from everywhere, and the search still
+ // terminates — the expanding box has to grow several times to find
+ // anything at all.
+ let (_, d) = grid.closest(Vec3::ZERO).unwrap();
+ assert!((d - 1.0).abs() < 0.05, "from the centre: {d}");
+
+ assert!(TriGrid::build(&Detail::new()).closest(Vec3::ZERO).is_none());
+ }
+
+ #[test]
+ fn test_detangle_separates_what_is_near_in_space_but_far_across_the_surface() {
+ // Two sheets pressed closer together than the thickness. They share no
+ // topology, so every pair between them is a self-intersection in
+ // waiting; within each sheet, neighbours are closer than the thickness
+ // by construction and must NOT be pushed apart.
+ let mut d = Detail::new();
+ let step = 0.25;
+ let gap = 0.05;
+ let mut rows = Vec::new();
+ for sheet in 0..2 {
+ let mut row = Vec::new();
+ for i in 0..4 {
+ for j in 0..4 {
+ row.push(d.add_point(Vec3::new(
+ i as f32 * step,
+ sheet as f32 * gap,
+ j as f32 * step,
+ )));
+ }
+ }
+ rows.push(row);
+ }
+ for sheet in 0..2 {
+ for i in 0..3 {
+ for j in 0..3 {
+ let at = |a: usize, b: usize| rows[sheet][a * 4 + b];
+ d.add_prim(&[at(i, j), at(i + 1, j), at(i + 1, j + 1)]);
+ d.add_prim(&[at(i, j), at(i + 1, j + 1), at(i, j + 1)]);
+ }
+ }
+ }
+ let before_gap = d.pos(16).y - d.pos(0).y;
+ assert!((before_gap - gap).abs() < 1e-6);
+ let before_edge = (d.pos(1) - d.pos(0)).length();
+
+ let node = phase3_node(
+ "detangle",
+ &[("Thickness", "1.00"), ("Rings", "2"), ("Iterations", "6")],
+ );
+ crate::geometry::apply_detangle(&mut d, &node);
+
+ // The sheets moved apart.
+ let after_gap = d.pos(16).y - d.pos(0).y;
+ assert!(after_gap > before_gap * 2.0, "sheets did not separate: {after_gap}");
+ // But the mesh did not explode: points that are neighbours ACROSS the
+ // surface are within a thickness of each other by construction, and a
+ // repulsion that did not exclude them would blow every sheet apart
+ // from the inside.
+ let after_edge = (d.pos(1) - d.pos(0)).length();
+ assert!(
+ (after_edge / before_edge - 1.0).abs() < 0.35,
+ "the sheet stretched from {before_edge} to {after_edge}"
+ );
+ }
+
+ #[test]
+ fn test_detangle_is_independent_of_point_order() {
+ // Gathered against the start-of-pass positions and applied at the end,
+ // so the same tangle untangles the same way however its points are
+ // numbered. A Gauss-Seidel sweep would not.
+ let mut a = sphere_detail(Vec3::ZERO, 0.5, 6, 8);
+ let node = phase3_node("detangle", &[("Thickness", "2.00"), ("Rings", "1"), ("Iterations", "3")]);
+ let mut b = a.clone();
+ crate::geometry::apply_detangle(&mut a, &node);
+ crate::geometry::apply_detangle(&mut b, &node);
+ assert_eq!(a.positions(), b.positions());
+ }
+
+ #[test]
+ fn test_suture_counts_sustained_contact_before_it_fuses() {
+ // A grid sitting just above a collider it is in contact with.
+ let collider = {
+ let mut c = Detail::new();
+ let pts: Vec<u32> = [
+ Vec3::new(-1.0, 0.0, -1.0),
+ Vec3::new(1.0, 0.0, -1.0),
+ Vec3::new(1.0, 0.0, 1.0),
+ Vec3::new(-1.0, 0.0, 1.0),
+ ]
+ .iter()
+ .map(|&p| c.add_point(p))
+ .collect();
+ c.add_prim(&[pts[0], pts[1], pts[2]]);
+ c.add_prim(&[pts[0], pts[2], pts[3]]);
+ c
+ };
+ let mut sheet = Detail::new();
+ for i in 0..3 {
+ sheet.add_point(Vec3::new(i as f32 * 0.02, 0.01, 0.0));
+ }
+ sheet.add_point(Vec3::new(0.0, 5.0, 0.0)); // far away, never in contact
+
+ let node = phase3_node(
+ "suture",
+ &[("Distance Threshold", "0.10"), ("Fusion Threshold", "3"), ("Counter", "contact")],
+ );
+ let count = |d: &Detail, p: usize| d.points().value("contact", p).unwrap().as_f32() as i32;
+
+ // Contact accrues, and the contacting points are pushed out to the
+ // threshold.
+ crate::geometry::apply_suture(&mut sheet, Some(&collider), &node);
+ assert_eq!(count(&sheet, 0), 1);
+ assert!((sheet.pos(0).y - 0.10).abs() < 1e-3, "not pushed out: {}", sheet.pos(0).y);
+ // A point out of contact stays at zero — contact has to be SUSTAINED
+ // to count, which is the difference between brushing past and growing
+ // together.
+ assert_eq!(count(&sheet, 3), 0);
+
+ crate::geometry::apply_suture(&mut sheet, Some(&collider), &node);
+ assert_eq!(count(&sheet, 0), 2);
+ assert_eq!(sheet.num_points(), 4, "nothing fuses below the threshold");
+
+ // The third crossing takes them past Fusion Threshold, and the three
+ // contacting points — all within Distance Threshold of each other —
+ // become one. The far point is untouched.
+ crate::geometry::apply_suture(&mut sheet, Some(&collider), &node);
+ assert_eq!(sheet.num_points(), 2, "sustained contact did not fuse");
+ }
+
+ #[test]
+ fn test_suture_without_a_collider_changes_nothing_but_the_counter() {
+ let mut sheet = sphere_detail(Vec3::ZERO, 0.5, 4, 6);
+ let before = sheet.positions().to_vec();
+ let node = phase3_node("suture", &[("Distance Threshold", "0.10"), ("Counter", "contact")]);
+ crate::geometry::apply_suture(&mut sheet, None, &node);
+ assert_eq!(sheet.positions(), &before[..]);
+ // The counter exists so the chain downstream can read it either way.
+ assert!(sheet.points().has("contact"));
+ assert_eq!(sheet.points().kind("contact"), AttribKind::Live);
+ }
+
+ #[test]
+ fn test_fusing_points_keeps_the_representative_and_drops_folded_faces() {
+ let mut d = quad_grid();
+ let keep_id = d.id(0).unwrap();
+ d.points_mut().create("mass", AttribValue::Float(0.0));
+ d.points_mut().set_value("mass", 0, AttribValue::Float(7.0)).unwrap();
+ d.points_mut().set_value("mass", 1, AttribValue::Float(9.0)).unwrap();
+
+ // Point 1 merges onto point 0.
+ let mut rep: Vec<u32> = (0..9).collect();
+ rep[1] = 0;
+ d.fuse_points(&rep);
+
+ assert_eq!(d.num_points(), 8);
+ // The representative keeps its identity AND its values — the same
+ // choice the remesher's collapse makes, and for the same reason.
+ assert_eq!(d.id(0), Some(keep_id));
+ assert_eq!(d.points().value("mass", 0), Some(AttribValue::Float(7.0)));
+ // Every surviving primitive is still a real triangle or quad; the ones
+ // that ended up with a repeated corner are gone.
+ for prim in 0..d.num_prims() {
+ let pts = d.prim_points(prim);
+ let mut uniq = pts.to_vec();
+ uniq.sort_unstable();
+ uniq.dedup();
+ assert_eq!(uniq.len(), pts.len(), "prim {prim} folded onto itself");
+ }
+ // A chain resolves: a→b→c leaves everything at c.
+ let mut e = quad_grid();
+ let mut chain: Vec<u32> = (0..9).collect();
+ chain[2] = 1;
+ chain[1] = 0;
+ e.fuse_points(&chain);
+ assert_eq!(e.num_points(), 7);
+ }
+
// ---- Phase 2: the solver contract ----
#[test]
diff --git a/src/remesh.rs b/src/remesh.rs
index ae339a6..a622ecc 100644
--- a/src/remesh.rs
+++ b/src/remesh.rs
@@ -34,11 +34,11 @@
//! memory as it can.
//! - A **flip** and a **relax** change no attributes at all.
//!
-//! What it does NOT do is project back onto the input surface, which the
-//! paper's fifth pass does. Tangential relaxation alone lets a surface creep
-//! slightly over many iterations; the projection needs a spatial index over the
-//! original triangles, which is the same index Detangle and Suture will want,
-//! and is better built once for all three.
+//! The paper's fifth pass is here too: after relaxing, every point is pulled
+//! back onto the surface the remesh started from. Tangential relaxation alone
+//! lets a surface creep — each point slides a little, and over a few dozen
+//! iterations a sphere quietly shrinks — so the projection is what makes it
+//! safe to run a remesh every frame of a solve, which is the whole point.
use crate::detail::{AttribData, AttribValue, Detail, PointId};
use glam::Vec3;
@@ -56,11 +56,21 @@ pub struct Settings {
pub split: bool,
pub collapse: bool,
pub flip: bool,
+ /// Pull relaxed points back onto the input surface.
+ pub project: bool,
}
impl Default for Settings {
fn default() -> Self {
- Self { target: 0.1, iterations: 3, relax: 0.5, split: true, collapse: true, flip: true }
+ Self {
+ target: 0.1,
+ iterations: 3,
+ relax: 0.5,
+ split: true,
+ collapse: true,
+ flip: true,
+ project: true,
+ }
}
}
@@ -559,11 +569,36 @@ fn relax_pass(m: &mut Mesh, amount: f32) {
}
}
+/// Pull every point back onto the surface the remesh started from.
+///
+/// Relaxation slides points within the surface, but "within" is only true to
+/// 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) {
+ if rest.is_empty() {
+ return;
+ }
+ for p in 0..m.pos.len() {
+ if m.dead_point[p] {
+ continue;
+ }
+ if let Some((q, _)) = rest.closest(m.pos[p]) {
+ m.pos[p] = q;
+ }
+ }
+}
+
/// Remesh toward `settings.target` edge length.
pub fn remesh(input: &Detail, settings: Settings) -> Detail {
if input.num_prims() == 0 || settings.target <= 0.0 {
return input.clone();
}
+ // Built once from the INPUT and reused by every iteration: projecting onto
+ // 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 m = Mesh::from_detail(input);
for _ in 0..settings.iterations.min(20) {
if settings.split {
@@ -576,6 +611,9 @@ pub fn remesh(input: &Detail, settings: Settings) -> Detail {
flip_pass(&mut m);
}
relax_pass(&mut m, settings.relax.clamp(0.0, 1.0));
+ if let Some(rest) = &rest {
+ project_pass(&mut m, rest);
+ }
}
m.into_detail()
}
diff --git a/src/spatial.rs b/src/spatial.rs
new file mode 100644
index 0000000..31161ae
--- /dev/null
+++ b/src/spatial.rs
@@ -0,0 +1,246 @@
+//! Uniform grids for the "what is near this?" questions.
+//!
+//! Three operators want one of these and all three arrived at once, which is
+//! why it is built here rather than inside any of them: the remesher's
+//! projection pass asks for the closest point on a surface, Detangle asks
+//! which points are within a thickness, and Suture asks both.
+//!
+//! A uniform grid rather than a BVH because the geometry these run on is
+//! already near-uniform — that is what remeshing is for — and a grid sized to
+//! the mesh's own scale has no degenerate case on it. A surface with wildly
+//! varying triangle sizes would want a tree, and that is the day to write one.
+
+use crate::detail::Detail;
+use glam::Vec3;
+
+/// The closest point to `p` on triangle `(a, b, c)`.
+///
+/// The Voronoi-region walk from Ericson's *Real-Time Collision Detection*
+/// §5.1.5: test the three vertex regions, then the three edge regions, and
+/// what is left is the face interior.
+pub fn closest_point_on_triangle(p: Vec3, a: Vec3, b: Vec3, c: Vec3) -> Vec3 {
+ let (ab, ac, ap) = (b - a, c - a, p - a);
+ let (d1, d2) = (ab.dot(ap), ac.dot(ap));
+ if d1 <= 0.0 && d2 <= 0.0 {
+ return a;
+ }
+ let bp = p - b;
+ let (d3, d4) = (ab.dot(bp), ac.dot(bp));
+ if d3 >= 0.0 && d4 <= d3 {
+ return b;
+ }
+ let vc = d1 * d4 - d3 * d2;
+ if vc <= 0.0 && d1 >= 0.0 && d3 <= 0.0 {
+ let denom = d1 - d3;
+ let v = if denom.abs() < 1e-20 { 0.0 } else { d1 / denom };
+ return a + ab * v;
+ }
+ let cp = p - c;
+ let (d5, d6) = (ab.dot(cp), ac.dot(cp));
+ if d6 >= 0.0 && d5 <= d6 {
+ return c;
+ }
+ let vb = d5 * d2 - d1 * d6;
+ if vb <= 0.0 && d2 >= 0.0 && d6 <= 0.0 {
+ let denom = d2 - d6;
+ let w = if denom.abs() < 1e-20 { 0.0 } else { d2 / denom };
+ return a + ac * w;
+ }
+ let va = d3 * d6 - d5 * d4;
+ if va <= 0.0 && (d4 - d3) >= 0.0 && (d5 - d6) >= 0.0 {
+ let denom = (d4 - d3) + (d5 - d6);
+ let w = if denom.abs() < 1e-20 { 0.0 } else { (d4 - d3) / denom };
+ return b + (c - b) * w;
+ }
+ let denom = va + vb + vc;
+ if denom.abs() < 1e-20 {
+ return a;
+ }
+ a + ab * (vb / denom) + ac * (vc / denom)
+}
+
+/// Where things are, bucketed by cell.
+///
+/// Shared by both grids: they differ only in what they store and how they
+/// answer, not in how they divide space.
+struct Grid {
+ min: Vec3,
+ cell: f32,
+ dims: [i32; 3],
+ buckets: Vec<Vec<u32>>,
+}
+
+impl Grid {
+ /// A grid over `bounds` with cells about `cell` across, capped so that a
+ /// pathological request cannot ask for a billion buckets.
+ fn new(min: Vec3, max: Vec3, cell: f32) -> Grid {
+ let span = (max - min).max(Vec3::splat(1e-6));
+ let cell = cell.max(span.max_element() / 128.0).max(1e-6);
+ let dims = [
+ ((span.x / cell).ceil() as i32 + 1).clamp(1, 256),
+ ((span.y / cell).ceil() as i32 + 1).clamp(1, 256),
+ ((span.z / cell).ceil() as i32 + 1).clamp(1, 256),
+ ];
+ let n = (dims[0] * dims[1] * dims[2]) as usize;
+ Grid { min, cell, dims, buckets: vec![Vec::new(); n] }
+ }
+
+ fn coord(&self, p: Vec3) -> [i32; 3] {
+ let rel = (p - self.min) / self.cell;
+ [
+ (rel.x.floor() as i32).clamp(0, self.dims[0] - 1),
+ (rel.y.floor() as i32).clamp(0, self.dims[1] - 1),
+ (rel.z.floor() as i32).clamp(0, self.dims[2] - 1),
+ ]
+ }
+
+ fn index(&self, c: [i32; 3]) -> usize {
+ ((c[2] * self.dims[1] + c[1]) * self.dims[0] + c[0]) as usize
+ }
+
+ fn insert(&mut self, p: Vec3, id: u32) {
+ let i = self.index(self.coord(p));
+ self.buckets[i].push(id);
+ }
+
+ /// Everything in the cells overlapping the box, deduplicated.
+ fn gather(&self, lo: Vec3, hi: Vec3, out: &mut Vec<u32>) {
+ out.clear();
+ let (a, b) = (self.coord(lo), self.coord(hi));
+ for z in a[2]..=b[2] {
+ for y in a[1]..=b[1] {
+ for x in a[0]..=b[0] {
+ out.extend_from_slice(&self.buckets[self.index([x, y, z])]);
+ }
+ }
+ }
+ out.sort_unstable();
+ out.dedup();
+ }
+}
+
+/// A grid over a surface's triangles, for asking what the nearest surface
+/// point is.
+pub struct TriGrid {
+ tris: Vec<[Vec3; 3]>,
+ grid: Grid,
+}
+
+impl TriGrid {
+ pub fn build(d: &Detail) -> TriGrid {
+ let tris: Vec<[Vec3; 3]> = d
+ .triangulate(|pos, _| Vec3::from(pos))
+ .chunks_exact(3)
+ .map(|t| [t[0], t[1], t[2]])
+ .collect();
+ let (min, max) = d.bounds().unwrap_or((Vec3::ZERO, Vec3::ZERO));
+ // Cells about the size of a triangle: small enough that a cell holds
+ // few, large enough that one triangle spans few.
+ let mean = if tris.is_empty() {
+ 1.0
+ } else {
+ tris.iter()
+ .map(|t| (t[1] - t[0]).length().max((t[2] - t[0]).length()))
+ .sum::<f32>()
+ / tris.len() as f32
+ };
+ let mut grid = Grid::new(min, max, mean.max(1e-5));
+ // A triangle goes in every cell its bounding box touches, so a lookup
+ // that finds a cell finds every triangle passing through it.
+ for (i, t) in tris.iter().enumerate() {
+ let lo = t[0].min(t[1]).min(t[2]);
+ let hi = t[0].max(t[1]).max(t[2]);
+ let (a, b) = (grid.coord(lo), grid.coord(hi));
+ for z in a[2]..=b[2] {
+ for y in a[1]..=b[1] {
+ for x in a[0]..=b[0] {
+ let idx = grid.index([x, y, z]);
+ grid.buckets[idx].push(i as u32);
+ }
+ }
+ }
+ }
+ TriGrid { tris, grid }
+ }
+
+ pub fn is_empty(&self) -> bool {
+ self.tris.is_empty()
+ }
+
+ /// The closest point on the surface, and its distance.
+ ///
+ /// Searches an expanding box until the best hit is closer than the box is
+ /// wide — at which point nothing outside can beat it, because anything out
+ /// there is at least that far away.
+ pub fn closest(&self, p: Vec3) -> Option<(Vec3, f32)> {
+ if self.tris.is_empty() {
+ return None;
+ }
+ let mut reach = self.grid.cell;
+ let mut scratch = Vec::new();
+ for _ in 0..12 {
+ self.grid
+ .gather(p - Vec3::splat(reach), p + Vec3::splat(reach), &mut scratch);
+ let best = scratch
+ .iter()
+ .map(|&i| {
+ let t = self.tris[i as usize];
+ let q = closest_point_on_triangle(p, t[0], t[1], t[2]);
+ (q, (q - p).length())
+ })
+ .min_by(|a, b| a.1.partial_cmp(&b.1).unwrap_or(std::cmp::Ordering::Equal));
+ match best {
+ Some(hit) if hit.1 <= reach => return Some(hit),
+ _ => reach *= 2.0,
+ }
+ }
+ // The grid is clamped to the geometry's bounds, so twelve doublings
+ // have gathered everything; whatever it found is the answer.
+ let best = self
+ .tris
+ .iter()
+ .map(|t| {
+ let q = closest_point_on_triangle(p, t[0], t[1], t[2]);
+ (q, (q - p).length())
+ })
+ .min_by(|a, b| a.1.partial_cmp(&b.1).unwrap_or(std::cmp::Ordering::Equal));
+ best
+ }
+}
+
+/// A grid over a point set, for asking which points are near a place.
+pub struct PointGrid {
+ points: Vec<Vec3>,
+ grid: Grid,
+}
+
+impl PointGrid {
+ pub fn build(points: &[Vec3], cell: f32) -> PointGrid {
+ let (min, max) = points.iter().fold(
+ (Vec3::splat(f32::MAX), Vec3::splat(f32::MIN)),
+ |(lo, hi), &p| (lo.min(p), hi.max(p)),
+ );
+ let (min, max) = if points.is_empty() { (Vec3::ZERO, Vec3::ZERO) } else { (min, max) };
+ let mut grid = Grid::new(min, max, cell);
+ for (i, &p) in points.iter().enumerate() {
+ grid.insert(p, i as u32);
+ }
+ PointGrid { points: points.to_vec(), grid }
+ }
+
+ /// Indices of points within `radius` of `p`, excluding nothing — the
+ /// caller decides what does not count as a neighbour, because "not
+ /// itself" and "not topologically adjacent" are different questions.
+ pub fn within(&self, p: Vec3, radius: f32, out: &mut Vec<u32>) {
+ let mut scratch = Vec::new();
+ self.grid
+ .gather(p - Vec3::splat(radius), p + Vec3::splat(radius), &mut scratch);
+ let r2 = radius * radius;
+ out.clear();
+ out.extend(
+ scratch
+ .into_iter()
+ .filter(|&i| (self.points[i as usize] - p).length_squared() <= r2),
+ );
+ }
+}