graphic design tool
git clone https://git.lucas.co/cce-designer.git
feat(volume): signed distance fields, so shelling and booleans are possible
Shelling, offsetting and booleans are not mesh operations. On triangles
they mean asking "which side of this whole surface is that point on" for
every pair; on a field they are `min`, `max` and a sign flip, and the mesh
comes back out by extraction. `src/volume.rs` is that field — a dense SDF
— with two nodes on it: `volume` (offset, shell) and `boolean` (union,
intersect, subtract).
Signing the field is the whole difficulty, and it took three wrong
answers to get one right:
- Ray parity double-counts where a ray crosses a shared edge, inverting
the parity for every sample behind it. 79 of 15625 samples wrong, in
contiguous runs — which is what a parity bug looks like.
- A flood fill from the grid boundary fixes that, but only if it may
step between samples BOTH further than a full voxel from any surface.
At 0.75 voxel it walked straight through a thin wall and a slab came
back hollow. Two samples one voxel apart cannot both be further than a
voxel from a surface lying between them; that is the whole argument.
- Inside the band the flood has nothing to say, so the nearest face's
normal decides — which trusts the winding. A mesh wound inside out came
back with its band signs alternating against the flood's. The winding is
now measured (signed volume, divergence theorem) and the test flips to
match: an imported mesh is not obliged to agree with our convention.
Extraction is naive surface nets, chosen over marching cubes because it
gives quads on a quad grid and far fewer slivers. Its limit is one vertex
per cell: a feature thinner than a voxel — a subtraction's knife-edge rim
— pinches two sheets onto one vertex, leaving edges with four faces.
Watertight, not manifold. So `Detail` now distinguishes the two.
`is_closed` asks that every directed edge have one opposite (no boundary,
consistently wound — what having an inside requires, and what
`Volume::build` guards its input with); `is_manifold` asks for exactly two
faces per edge (what remeshing requires, since an edge with four faces has
no pair to flip between).
`TriGrid::closest_within` does one bounded gather rather than doubling a
search radius, which took the suite from 33s back to 2s.
Both nodes have resolver-level tests as well as field unit tests. The
arithmetic was right and the wiring was never exercised — a hand-built
project of two spheres and a Boolean rendered an empty scene while every
unit test passed. A node nobody can reach from the network is not a
feature.
Co-Authored-By: Claude Opus 5 <[email protected]>
CLAUDE.md | 52 ++++++
nodes/boolean.json | 12 ++
nodes/volume.json | 13 ++
shapeshifter.md | 30 +++-
src/detail.rs | 63 ++++++++
src/geometry.rs | 314 ++++++++++++++++++++++++++++++++----
src/main.rs | 254 ++++++++++++++++++++++++++++-
src/spatial.rs | 71 ++++++++
src/volume.rs | 465 +++++++++++++++++++++++++++++++++++++++++++++++++++++
9 files changed, 1236 insertions(+), 38 deletions(-)
diff --git a/CLAUDE.md b/CLAUDE.md
index e049f7f..7deb75c 100644
--- a/CLAUDE.md
+++ b/CLAUDE.md
@@ -308,6 +308,58 @@ Buttons dispatch through `execute_menu_action` by LABEL, which carries no node
— `run_export` resolves the node from the current selection, which is sound
because the pressed button can only be on the node the pane is showing.
+### The volume representation
+
+`src/volume.rs` is a dense signed distance field — `Volume { origin, voxel,
+dims, data }` — with two nodes on it: `volume` (offset and shell) and
+`boolean` (union, intersect, subtract). It exists because shelling, offsetting
+and booleans are not mesh operations. Doing them on triangles means answering
+"which side of this whole surface is that point on" per triangle pair; doing
+them on a field means `min`, `max` and a sign flip, and the mesh comes back out
+by extraction.
+
+**Signing the field is the whole difficulty**, and it is done in two parts
+because neither part is right everywhere:
+
+- **Far from the surface**, a flood fill from the grid boundary — which is
+ outside by construction — marks everything it can reach. Whatever it cannot
+ reach without crossing the surface is enclosed, however convoluted the
+ cavity. The flood may only step between samples that are *both* further than
+ `voxel * 1.01` from any surface, because two samples one voxel apart cannot
+ both be more than a voxel from a surface lying between them. A looser band
+ (0.75 voxel was the first try) lets the flood walk straight through a thin
+ wall and the solid comes back hollow.
+- **Inside that band**, the flood has nothing to say, so the nearest face's
+ normal decides. That test trusts the winding, so the winding is *measured*
+ first — the signed volume by the divergence theorem, positive when faces look
+ outward — and the test flips if the mesh is inside out. An imported mesh is
+ not obliged to agree with this app's convention, and one that disagrees used
+ to come back with its band signs alternating against the flood's.
+
+Ray parity was the first approach and is wrong: a ray through a shared edge
+crosses two triangles at one point and counts two, so the parity inverts for
+every sample behind it. It failed on 79 of 15625 samples in contiguous runs,
+which is what a parity bug looks like.
+
+Extraction is naive **surface nets** (`to_mesh`): one vertex per cell that has
+a sign change, placed at the average of its edge crossings, and one quad per
+crossed grid edge joining the four cells around it. Chosen over marching cubes
+because it produces quads on a quad grid and far fewer degenerate slivers.
+
+One vertex per cell is also its limit. Where a feature is thinner than a voxel
+— the knife-edge rim of a subtraction — two sheets of surface share one cell's
+vertex and pinch, leaving edges with four faces. The result is still
+watertight; it is not manifold. Hence two predicates on `Detail`, and the
+difference matters: [`is_closed`](src/detail.rs) asks that every directed edge
+have exactly one opposite (no boundary, consistently wound — what having an
+inside requires, and what `Volume::build` guards its input with), while
+`is_manifold` asks for exactly two faces per edge (what remeshing requires,
+since an edge with four faces has no single pair to flip between).
+
+`Volume::build` takes an explicit `reach`: distances are clamped there, so the
+field is exact near the surface and flat far from it. A boolean builds both
+operands on ONE grid so the two fields line up sample for sample.
+
### Runtime paths point into the source tree
Node templates (`nodes/*.json`) and `default_project.json` are located via
diff --git a/nodes/boolean.json b/nodes/boolean.json
new file mode 100644
index 0000000..f194a54
--- /dev/null
+++ b/nodes/boolean.json
@@ -0,0 +1,12 @@
+{
+ "name": "Boolean",
+ "type": "boolean",
+ "inputs": 2,
+ "outputs": 1,
+ "params": [
+ { "name": "Input", "type": "text", "default": "" },
+ { "name": "With", "type": "text", "default": "" },
+ { "name": "Operation", "type": "choice:Union,Intersect,Subtract", "default": "Union" },
+ { "name": "Voxel Size", "type": "slider", "default": "0.05", "min": 0.005, "max": 1.0, "step": 0.005 }
+ ]
+}
diff --git a/nodes/volume.json b/nodes/volume.json
new file mode 100644
index 0000000..16e737c
--- /dev/null
+++ b/nodes/volume.json
@@ -0,0 +1,13 @@
+{
+ "name": "Volume",
+ "type": "volume",
+ "inputs": 1,
+ "outputs": 1,
+ "params": [
+ { "name": "Input", "type": "text", "default": "" },
+ { "name": "Mode", "type": "choice:Offset,Shell", "default": "Offset" },
+ { "name": "Voxel Size", "type": "slider", "default": "0.05", "min": 0.005, "max": 1.0, "step": 0.005 },
+ { "name": "Offset", "type": "slider", "default": "0.00", "min": -1.0, "max": 1.0, "step": 0.005 },
+ { "name": "Thickness", "type": "slider", "default": "0.05", "min": 0.001, "max": 1.0, "step": 0.005, "show_when": "Mode == Shell" }
+ ]
+}
diff --git a/shapeshifter.md b/shapeshifter.md
index 6221077..9cd4721 100644
--- a/shapeshifter.md
+++ b/shapeshifter.md
@@ -359,9 +359,33 @@ Touches: `shortcut.rs`, `app.rs`, `slots.rs`, `cce-ui`.
> writes STL and OBJ, there is an `export` node and a `--export` CLI mode, and
> a solved growth simulation can be written to a printable file.
>
-> Still outstanding for this phase: the volume representation (SDF or sparse
-> grid) that shelling, offsetting and boolean work need, and the 2D page
-> context for the COP family.
+> **The volume representation landed.** `src/volume.rs` is a dense signed
+> distance field; the `volume` node offsets and shells, the `boolean` node
+> unions, intersects and subtracts. Extraction is surface nets.
+>
+> Signing the field cost three wrong answers before a right one. Ray parity
+> double-counts at shared edges and inverted 79 of 15625 samples. A flood fill
+> from the grid boundary fixes that, but a 0.75-voxel band let it walk through
+> a thin wall and a slab came back hollow — the band has to be a full voxel,
+> because two samples one voxel apart cannot both be further than a voxel from
+> a surface between them. And the band's own test, which asks the nearest face
+> which side a sample is on, trusts the winding; a mesh wound inside out came
+> back with its band signs alternating against the flood's, so the winding is
+> now measured by the divergence theorem and the test flips to match.
+>
+> The limit worth knowing: one vertex per cell means a feature thinner than a
+> voxel pinches. A subtraction's knife-edge rim leaves a handful of edges
+> carrying four faces — watertight, but not manifold. `Detail` now distinguishes
+> the two (`is_closed` / `is_manifold`), because voxelizing needs only the
+> first and remeshing needs the second.
+>
+> Two of these were found by RENDERING rather than testing, which is now the
+> third time this phase: a node whose arithmetic is right and whose wiring is
+> never exercised looks exactly like a working node until you ask the viewport
+> to draw it. Both new nodes now have resolver-level tests, not just unit tests
+> on the field.
+>
+> Still outstanding for this phase: the 2D page context for the COP family.
Furthest out because it needs infrastructure nothing else does: a **volume
representation** (SDF or sparse grid) for shelling, offsetting and boolean work,
diff --git a/src/detail.rs b/src/detail.rs
index a872102..5b3ab99 100644
--- a/src/detail.rs
+++ b/src/detail.rs
@@ -1449,6 +1449,69 @@ impl Detail {
let _ = self.points.set_value(CD, p, AttribValue::Float3(c));
}
+ /// Whether the surface is closed: every directed edge has exactly one
+ /// opposite.
+ ///
+ /// The question "what is inside this?" only has an answer for a closed
+ /// surface. A flat disc, a torn mesh or a single polygon has no inside, and
+ /// anything that signs a distance field has to know the difference — sign
+ /// an open surface and you get whichever side its normals happen to face,
+ /// which is not a solid, just a preference.
+ ///
+ /// Directed, not undirected, because the property that matters is "no
+ /// boundary, consistently wound", and those are the same test: every edge
+ /// walked one way by one face and the other way by its neighbour. Counting
+ /// undirected edges instead would ask for exactly two faces per edge, which
+ /// is [`is_manifold`](Self::is_manifold) — a stricter thing that a boolean
+ /// legitimately fails where two sheets pinch together along a knife edge
+ /// thinner than a voxel. Such a surface is still watertight, still has an
+ /// inside, and still voxelizes correctly.
+ pub fn is_closed(&self) -> bool {
+ if self.num_prims() == 0 {
+ return false;
+ }
+ let mut counts: HashMap<[u32; 2], i32> = HashMap::new();
+ for prim in 0..self.num_prims() {
+ let pts = self.prim_points(prim);
+ if pts.len() < 3 {
+ return false;
+ }
+ for i in 0..pts.len() {
+ let (a, b) = (pts[i], pts[(i + 1) % pts.len()]);
+ // One counter per undirected edge, incremented one way and
+ // decremented the other: it lands on zero exactly when every
+ // traversal is matched by an opposite one.
+ let (key, step) = if a < b { ([a, b], 1) } else { ([b, a], -1) };
+ *counts.entry(key).or_default() += step;
+ }
+ }
+ counts.values().all(|&c| c == 0)
+ }
+
+ /// Whether every edge is shared by exactly two primitives.
+ ///
+ /// Stricter than [`is_closed`](Self::is_closed): it also rules out the
+ /// place where more than two faces meet along one edge. Remeshing wants
+ /// this — an edge with four faces has no single pair to flip or collapse
+ /// between — while voxelizing does not.
+ pub fn is_manifold(&self) -> bool {
+ if self.num_prims() == 0 {
+ return false;
+ }
+ let mut counts: HashMap<[u32; 2], usize> = HashMap::new();
+ for prim in 0..self.num_prims() {
+ let pts = self.prim_points(prim);
+ if pts.len() < 3 {
+ return false;
+ }
+ for i in 0..pts.len() {
+ let (a, b) = (pts[i], pts[(i + 1) % pts.len()]);
+ *counts.entry([a.min(b), a.max(b)]).or_default() += 1;
+ }
+ }
+ counts.values().all(|&c| c == 2)
+ }
+
/// The axis-aligned bounds, or `None` when there are no points.
pub fn bounds(&self) -> Option<(Vec3, Vec3)> {
let first = *self.pos.first()?;
diff --git a/src/geometry.rs b/src/geometry.rs
index ab80745..3f9e4ac 100644
--- a/src/geometry.rs
+++ b/src/geometry.rs
@@ -34,39 +34,6 @@ impl SimpleRng {
}
}
-fn ray_triangle_intersect(
- origin: Vec3,
- dir: Vec3,
- v0: Vec3,
- v1: Vec3,
- v2: Vec3,
-) -> Option<f32> {
- let edge1 = v1 - v0;
- let edge2 = v2 - v0;
- let h = dir.cross(edge2);
- let a = edge1.dot(h);
- if a.abs() < 1e-6 {
- return None;
- }
- let f = 1.0 / a;
- let s = origin - v0;
- let u = f * s.dot(h);
- if u < 0.0 || u > 1.0 {
- return None;
- }
- let q = s.cross(edge1);
- let v = f * dir.dot(q);
- if v < 0.0 || u + v > 1.0 {
- return None;
- }
- let t = f * edge2.dot(q);
- if t > 1e-5 {
- Some(t)
- } else {
- None
- }
-}
-
#[derive(Clone, Debug, PartialEq)]
pub enum GAttribute {
Float(f32),
@@ -764,6 +731,10 @@ pub fn generate_single_node_geometry_with_errors(
resolve_valence_geometry_with_errors(root, target, visited, ocl_error, sim)
} else if target.node_type.eq_ignore_ascii_case("deform") {
resolve_deform_geometry_with_errors(root, target, visited, ocl_error, sim)
+ } else if target.node_type.eq_ignore_ascii_case("volume") {
+ resolve_volume_geometry_with_errors(root, target, visited, ocl_error, sim)
+ } else if target.node_type.eq_ignore_ascii_case("boolean") {
+ resolve_boolean_geometry_with_errors(root, target, visited, ocl_error, sim)
} else if target.node_type.eq_ignore_ascii_case("export") {
resolve_export_geometry_with_errors(root, target, visited, ocl_error, sim)
} else if target.node_type.eq_ignore_ascii_case("subdivide") {
@@ -1249,7 +1220,7 @@ pub fn resolve_collision_geometry_with_errors(
} else {
let crossings = tris
.iter()
- .filter(|t| ray_triangle_intersect(pt, ray_dir, t[0], t[1], t[2]).is_some())
+ .filter(|t| crate::spatial::ray_triangle(pt, ray_dir, t[0], t[1], t[2]).is_some())
.count();
crossings % 2 == 1
}
@@ -1670,6 +1641,116 @@ pub fn resolve_cull_geometry_with_errors(
Some(geom)
}
+/// The Volume node: offset or shell a surface through a distance field.
+///
+/// Offsetting a mesh directly means resolving every self-intersection the move
+/// creates; through a field it is a subtraction and the result is closed by
+/// construction. Shell is the same trick twice — the shape minus the shape
+/// moved inward — which is what a mold wall is.
+///
+/// The cost is resolution: the result is a surface extracted from a grid, so
+/// detail finer than the Voxel Size is gone. That is the trade the
+/// representation makes, and it is why this is a node rather than something
+/// applied silently.
+pub fn resolve_volume_geometry_with_errors(
+ root: &FsNode,
+ target: &FsNode,
+ visited: &mut Vec<String>,
+ ocl_error: &mut Option<String>,
+ sim: &mut EvalSim,
+) -> Option<Detail> {
+ let input_node = find_node_by_name(root, &node_param_str(target, "Input", ""))?;
+ let geom = generate_single_node_geometry_with_errors(root, input_node, visited, ocl_error, sim)?;
+ if geom.num_prims() == 0 {
+ return Some(geom);
+ }
+
+ let voxel = node_param_f32(target, "Voxel Size", 0.05).max(1e-3);
+ let offset = node_param_f32(target, "Offset", 0.0);
+ let shell = node_param_str(target, "Mode", "Offset").eq_ignore_ascii_case("shell");
+ let thickness = node_param_f32(target, "Thickness", 0.05).max(1e-4);
+
+ // Room for everything the operation will ask the field to reach: the
+ // offset itself, and for a shell the wall's thickness beyond it.
+ let want = offset.abs() + if shell { thickness } else { 0.0 } + voxel * 2.0;
+ let Some((lo, hi)) = crate::volume::Volume::bounds_for(&geom, want) else {
+ return Some(geom);
+ };
+ // A grid is capped at 256 samples an axis, so a voxel size far too small
+ // for the model silently gives a coarse answer. Saying so beats a result
+ // that looks like the node is broken.
+ let cells = ((hi - lo).max_element() / voxel).ceil();
+ if cells > 256.0 && ocl_error.is_none() {
+ *ocl_error = Some(format!(
+ "Volume '{}': voxel {:.3} needs {} samples across, over the 256 cap — the result is coarser than asked",
+ target.name, voxel, cells as i64
+ ));
+ }
+
+ let mut vol = crate::volume::Volume::build(&geom, lo, hi, voxel, want);
+ if shell {
+ // The wall between the offset surface and the same surface moved in by
+ // Thickness: intersect what is inside the outer with what is outside
+ // the inner.
+ let mut inner = vol.clone();
+ inner.offset(offset - thickness);
+ vol.offset(offset);
+ vol.subtract(&inner);
+ } else {
+ vol.offset(offset);
+ }
+ Some(vol.to_mesh())
+}
+
+/// The Boolean node: union, intersection and difference through a field.
+///
+/// Both inputs are sampled over ONE grid covering both, which is what makes
+/// the operation elementwise — a minimum, a maximum, a maximum against a
+/// negation. A mesh boolean spends its whole length finding intersection
+/// curves and stitching; this cannot produce an open surface because there is
+/// no stitching to get wrong.
+pub fn resolve_boolean_geometry_with_errors(
+ root: &FsNode,
+ target: &FsNode,
+ visited: &mut Vec<String>,
+ ocl_error: &mut Option<String>,
+ sim: &mut EvalSim,
+) -> Option<Detail> {
+ let input_node = find_node_by_name(root, &node_param_str(target, "Input", ""))?;
+ let a = generate_single_node_geometry_with_errors(root, input_node, visited, ocl_error, sim)?;
+
+ let with_name = node_param_str(target, "With", "");
+ let with_name = with_name.trim().to_string();
+ let Some(b) = find_node_by_name(root, &with_name)
+ .and_then(|n| generate_single_node_geometry_with_errors(root, n, visited, ocl_error, sim))
+ else {
+ if ocl_error.is_none() && !with_name.is_empty() {
+ *ocl_error = Some(format!("Boolean '{}': cannot resolve '{}'", target.name, with_name));
+ }
+ return Some(a);
+ };
+ if a.num_prims() == 0 || b.num_prims() == 0 {
+ return Some(a);
+ }
+
+ let voxel = node_param_f32(target, "Voxel Size", 0.05).max(1e-3);
+ let op = node_param_str(target, "Operation", "Union").to_lowercase();
+ // One grid over BOTH, so the two fields line up sample for sample.
+ let (alo, ahi) = a.bounds()?;
+ let (blo, bhi) = b.bounds()?;
+ let pad = Vec3::splat(voxel * 3.0);
+ let (lo, hi) = (alo.min(blo) - pad, ahi.max(bhi) + pad);
+
+ let mut va = crate::volume::Volume::build(&a, lo, hi, voxel, voxel * 3.0);
+ let vb = crate::volume::Volume::build(&b, lo, hi, voxel, voxel * 3.0);
+ match op.as_str() {
+ "intersect" => va.intersect(&vb),
+ "subtract" => va.subtract(&vb),
+ _ => va.union(&vb),
+ }
+ Some(va.to_mesh())
+}
+
/// The Export node: geometry out of the app.
///
/// A pass-through in the chain — it hands its input straight on, so it can sit
@@ -4837,6 +4918,8 @@ pub fn is_geometry_node_type(node_type: &str) -> bool {
|| nt == "detangle"
|| nt == "subdivide"
|| nt == "export"
+ || nt == "boolean"
+ || nt == "volume"
|| nt == "deform"
|| nt == "valence"
|| nt == "transfer"
@@ -5096,6 +5179,24 @@ pub fn network_sphere_vertices_with_errors(
out.merge(&geom);
}
}
+ } else if node.node_type.eq_ignore_ascii_case("volume") {
+ let _idx = *count;
+ *count += 1;
+ if is_visible {
+ let mut visited = Vec::new();
+ if let Some(geom) = resolve_volume_geometry_with_errors(root, node, &mut visited, ocl_error, sim) {
+ out.merge(&geom);
+ }
+ }
+ } else if node.node_type.eq_ignore_ascii_case("boolean") {
+ let _idx = *count;
+ *count += 1;
+ if is_visible {
+ let mut visited = Vec::new();
+ if let Some(geom) = resolve_boolean_geometry_with_errors(root, node, &mut visited, ocl_error, sim) {
+ out.merge(&geom);
+ }
+ }
} else if node.node_type.eq_ignore_ascii_case("export") {
let _idx = *count;
*count += 1;
@@ -7124,6 +7225,151 @@ mod simnet_tests {
assert!(vis_marker_vertices(&d, |c| c).is_empty());
}
+ /// The Boolean node end to end, through the resolver the viewport calls —
+ /// not just Volume's arithmetic.
+ ///
+ /// Written because a hand-built project of two spheres and a Boolean
+ /// rendered an EMPTY scene while every volume unit test passed: the
+ /// arithmetic was right and the wiring was never exercised. A node nobody
+ /// can reach from the network is not a feature.
+ #[test]
+ fn test_the_boolean_node_resolves_two_spheres_into_one_solid() {
+ // Spheres sit on the resolver's own 1.25-spaced layout, so a radius of
+ // 0.8 makes them overlap by 0.35 — a lens neither operation can
+ // mistake for the other.
+ let a = node("id-a", "sphere1", "sphere", vec![param("Radius", "0.8")], vec![]);
+ let b = node("id-b", "sphere2", "sphere", vec![param("Radius", "0.8")], vec![]);
+
+ for (op, expect) in [("Union", "wider"), ("Intersect", "narrower"), ("Subtract", "narrower")] {
+ let bool_node = node(
+ "id-bool",
+ "bool1",
+ "boolean",
+ vec![
+ param("Input", "sphere1"),
+ param("With", "sphere2"),
+ param("Operation", op),
+ param("Voxel Size", "0.08"),
+ ],
+ vec![],
+ );
+ let root = node("id-root", "root", "node", vec![], vec![a.clone(), b.clone(), bool_node]);
+ let target = &root.children[2];
+
+ let mut err = None;
+ let mut cache = SimCache::default();
+ let mut sim = EvalSim::new(0, 0, &mut cache);
+ let out = resolve_boolean_geometry_with_errors(&root, target, &mut Vec::new(), &mut err, &mut sim)
+ .unwrap_or_else(|| panic!("{op}: the Boolean resolved to nothing"));
+ assert!(err.is_none(), "{op}: {err:?}");
+ assert!(out.num_prims() > 0, "{op}: the Boolean produced no primitives");
+ assert!(out.is_closed(), "{op}: the result is not a closed surface");
+ // Union and Intersect come out fully manifold. Subtract does not:
+ // the bite leaves a rim thinner than a voxel, and surface nets has
+ // one vertex per cell to give it, so the two sheets pinch and a
+ // handful of edges carry four faces. Watertight, still solid, but
+ // recorded here rather than discovered downstream.
+ if op != "Subtract" {
+ assert!(out.is_manifold(), "{op}: the result is not manifold");
+ }
+
+ // The single sphere it started from, for a size to compare against.
+ let (slo, shi) = generate_single_node_geometry_with_errors(
+ &root, &root.children[0], &mut Vec::new(), &mut None, &mut sim,
+ )
+ .unwrap()
+ .bounds()
+ .unwrap();
+ let (olo, ohi) = out.bounds().unwrap();
+ let (one, both) = (shi.x - slo.x, ohi.x - olo.x);
+ if expect == "wider" {
+ assert!(both > one * 1.3, "{op}: {both} is not wider than one sphere's {one}");
+ } else {
+ assert!(both < one * 0.9, "{op}: {both} is not narrower than one sphere's {one}");
+ }
+ }
+ }
+
+ /// The Volume node's two modes, likewise through the resolver.
+ #[test]
+ fn test_the_volume_node_offsets_and_shells() {
+ let sphere = node("id-s", "sphere1", "sphere", vec![param("Radius", "0.8")], vec![]);
+
+ let grown = node(
+ "id-v",
+ "vol1",
+ "volume",
+ vec![
+ param("Input", "sphere1"),
+ param("Mode", "Offset"),
+ param("Voxel Size", "0.08"),
+ param("Offset", "0.2"),
+ ],
+ vec![],
+ );
+ let root = node("id-root", "root", "node", vec![], vec![sphere.clone(), grown]);
+ let mut err = None;
+ let mut cache = SimCache::default();
+ let mut sim = EvalSim::new(0, 0, &mut cache);
+ let out = resolve_volume_geometry_with_errors(&root, &root.children[1], &mut Vec::new(), &mut err, &mut sim)
+ .expect("the Volume node resolved to nothing");
+ assert!(err.is_none(), "{err:?}");
+ assert!(out.is_closed(), "an offset sphere is not a closed surface");
+ let (slo, shi) = generate_single_node_geometry_with_errors(
+ &root, &root.children[0], &mut Vec::new(), &mut None, &mut sim,
+ )
+ .unwrap()
+ .bounds()
+ .unwrap();
+ let (olo, ohi) = out.bounds().unwrap();
+ assert!(
+ (ohi.x - olo.x) > (shi.x - slo.x) + 0.25,
+ "a +0.2 offset did not grow the sphere: {} vs {}",
+ ohi.x - olo.x,
+ shi.x - slo.x
+ );
+
+ // A shell is hollow: closed, and with twice the surface of the solid.
+ let shell = node(
+ "id-v2",
+ "vol2",
+ "volume",
+ vec![
+ param("Input", "sphere1"),
+ param("Mode", "Shell"),
+ param("Voxel Size", "0.08"),
+ param("Offset", "0.0"),
+ param("Thickness", "0.15"),
+ ],
+ vec![],
+ );
+ let root = node("id-root", "root", "node", vec![], vec![sphere, shell]);
+ let mut err = None;
+ let mut cache = SimCache::default();
+ let mut sim = EvalSim::new(0, 0, &mut cache);
+ let out = resolve_volume_geometry_with_errors(&root, &root.children[1], &mut Vec::new(), &mut err, &mut sim)
+ .expect("the Volume node resolved to nothing");
+ assert!(err.is_none(), "{err:?}");
+ assert!(out.is_closed(), "a shell is not a closed surface");
+ assert!(out.num_prims() > 0, "the shell is empty");
+
+ // Hollow, and provably so: a shell has an INNER surface, so its points
+ // sit at two radii, not one. Rendered from outside it is
+ // indistinguishable from the solid sphere, which is exactly why this
+ // is asserted rather than looked at.
+ let centre = (slo + shi) * 0.5;
+ let radii: Vec<f32> = (0..out.num_points()).map(|p| (out.pos(p) - centre).length()).collect();
+ let (near, far) = radii.iter().fold((f32::MAX, 0.0f32), |(n, f), &r| (n.min(r), f.max(r)));
+ assert!(
+ (far - 0.8).abs() < 0.1,
+ "the shell's outer surface is at {far}, not the sphere's 0.8"
+ );
+ assert!(
+ (near - 0.65).abs() < 0.1,
+ "the shell has no cavity: its innermost point is at {near}, expected 0.8 - 0.15"
+ );
+ }
+
#[test]
fn test_an_opencl_deformer_lands_its_named_attribute_on_the_geometry() {
// The gap between "parse_attr_refs reads raw source" and "run the
diff --git a/src/main.rs b/src/main.rs
index 5967ecd..54f5202 100644
--- a/src/main.rs
+++ b/src/main.rs
@@ -7,6 +7,7 @@ pub mod export;
pub mod export_cli;
pub mod remesh;
pub mod spatial;
+pub mod volume;
// Root-level aliases some modules import via `crate::` paths.
#[allow(unused_imports)]
@@ -113,7 +114,7 @@ mod tests {
use crate::app::{get_next_visible_pane, DesignSettings, FsNode, Project, ProjectViewState};
use crate::slots::{LEFT_MENUBAR_IDX, RIGHT_MENUBAR_IDX, PARAM_MENUBAR_IDX, SPREADSHEET_MENUBAR_IDX};
use crate::shortcut::{Shortcut, ShortcutManager, Action};
- use crate::geometry::{GAttribute, GVertex, Geometry, line_vertices};
+ use crate::geometry::line_vertices;
use crate::detail::{AttribData, AttribKind, AttribType, AttribValue, Class, Detail};
/// The choosers open in the loaded project's parent — the "current view" —
@@ -3645,6 +3646,257 @@ mod tests {
assert_eq!(param_display(¶ms)[1].1, "7.5", "and comes back as it was");
}
+ // ---- Volumes ----
+
+ use crate::volume::Volume;
+
+ #[test]
+ fn test_a_sphere_round_trips_through_a_distance_field() {
+ let sphere = sphere_detail(Vec3::ZERO, 1.0, 16, 24);
+ let vol = Volume::from_mesh(&sphere, 0.12, 0.4);
+ let back = vol.to_mesh();
+
+ assert!(back.num_points() > 100, "the surface did not come back");
+ assert!(back.num_prims() > 100);
+
+ // Every extracted point is on the sphere, to within a voxel. That is
+ // the whole claim of the representation: a mesh in, a field, a mesh
+ // out, and the shape survives.
+ for p in 0..back.num_points() {
+ let r = back.pos(p).length();
+ assert!((r - 1.0).abs() < 0.14, "point {p} is at radius {r}");
+ }
+
+ // Closed: every edge is shared by exactly two faces. Surface nets
+ // gives this by construction, and it is what makes the output safe to
+ // hand a slicer.
+ let mut shared: std::collections::HashMap<[u32; 2], usize> = Default::default();
+ for prim in 0..back.num_prims() {
+ let pts = back.prim_points(prim);
+ for i in 0..pts.len() {
+ let (a, b) = (pts[i], pts[(i + 1) % pts.len()]);
+ *shared.entry([a.min(b), a.max(b)]).or_default() += 1;
+ }
+ }
+ let open = shared.values().filter(|&&c| c != 2).count();
+ assert_eq!(open, 0, "{open} edges are not shared by two faces");
+
+ // And it faces outward, like every other generator.
+ let normals = crate::geometry::point_normals(&back);
+ let outward = (0..back.num_points())
+ .filter(|&p| normals[p].dot(back.pos(p).normalize()) > 0.0)
+ .count();
+ assert_eq!(outward, back.num_points(), "the extracted surface is inside out");
+ }
+
+ #[test]
+ fn test_the_sign_is_right_where_the_nearest_face_would_lie() {
+ // A field's sign has to be right EVERYWHERE — a wrong one is a bubble
+ // or a hole, where in the Distance node it was a slightly wrong
+ // number. This is why the build casts rays rather than asking the
+ // nearest face which way it points, and this is the check that says so.
+ let sphere = sphere_detail(Vec3::ZERO, 1.0, 20, 28);
+ let vol = Volume::from_mesh(&sphere, 0.12, 0.3);
+ let [nx, ny, nz] = vol.dims();
+
+ let mut wrong = Vec::new();
+ for k in 0..nz {
+ for j in 0..ny {
+ for i in 0..nx {
+ let p = vol.sample_position(i, j, k);
+ let r = p.length();
+ // Skip the band where the answer is genuinely ambiguous at
+ // this resolution.
+ if (r - 1.0).abs() < vol.voxel() {
+ continue;
+ }
+ let want_inside = r < 1.0;
+ if (vol.at(i, j, k) < 0.0) != want_inside {
+ wrong.push((i, j, k, r, vol.at(i, j, k)));
+ }
+ }
+ }
+ }
+ assert!(
+ wrong.is_empty(),
+ "{} of {} samples have the wrong sign, e.g. {:?}",
+ wrong.len(),
+ nx * ny * nz,
+ &wrong[..wrong.len().min(3)]
+ );
+ }
+
+ /// An axis-aligned closed box, built by hand.
+ fn box_mesh(lo: Vec3, hi: Vec3) -> Detail {
+ let mut d = Detail::new();
+ for (x, y, z) in [
+ (lo.x, lo.y, lo.z), (hi.x, lo.y, lo.z), (hi.x, lo.y, hi.z), (lo.x, lo.y, hi.z),
+ (lo.x, hi.y, lo.z), (hi.x, hi.y, lo.z), (hi.x, hi.y, hi.z), (lo.x, hi.y, hi.z),
+ ] {
+ d.add_point(Vec3::new(x, y, z));
+ }
+ // Wound counter-clockwise seen from OUTSIDE, like every generator —
+ // asserted below, because getting this backwards by hand is exactly
+ // what happened the first time.
+ for q in [
+ [0u32, 1, 2, 3], [7, 6, 5, 4], [0, 4, 5, 1],
+ [1, 5, 6, 2], [2, 6, 7, 3], [3, 7, 4, 0],
+ ] {
+ d.add_prim(&q);
+ }
+ let centre = (lo + hi) * 0.5;
+ for (p, n) in crate::geometry::point_normals(&d).iter().enumerate() {
+ assert!(
+ n.dot((d.pos(p) - centre).normalize()) > 0.0,
+ "the test box's corner {p} faces inward"
+ );
+ }
+ d
+ }
+
+ #[test]
+ fn test_a_thin_slab_is_solid_all_the_way_through() {
+ // The flood fill decides what is enclosed, and it must not be able to
+ // walk THROUGH a wall. A slab only a few voxels thick is where that
+ // goes wrong: at a band narrower than a voxel, two adjacent samples
+ // straddling the surface can both read as "far", the flood steps
+ // between them, and the slab comes back hollow — which a boolean
+ // against it then fails to cut with.
+ let slab = box_mesh(Vec3::new(-1.0, -0.1, -1.0), Vec3::new(1.0, 0.1, 1.0));
+ assert!(slab.is_closed());
+ let vol = Volume::from_mesh(&slab, 0.05, 0.15);
+ let [nx, ny, nz] = vol.dims();
+
+ let mut inside_wrong = Vec::new();
+ for k in 0..nz {
+ for j in 0..ny {
+ for i in 0..nx {
+ let p = vol.sample_position(i, j, k);
+ let deep = p.x.abs() < 0.8 && p.z.abs() < 0.8 && p.y.abs() < 0.04;
+ if deep && vol.at(i, j, k) >= 0.0 {
+ inside_wrong.push((p, vol.at(i, j, k)));
+ }
+ }
+ }
+ }
+ assert!(
+ inside_wrong.is_empty(),
+ "{} samples inside the slab read as outside, e.g. {:?}",
+ inside_wrong.len(),
+ &inside_wrong[..inside_wrong.len().min(3)]
+ );
+
+ // And it cuts: subtracting the slab from a box that contains it leaves
+ // a gap where the slab was.
+ let block = box_mesh(Vec3::splat(-0.6), Vec3::splat(0.6));
+ let (lo, hi) = (Vec3::splat(-1.3), Vec3::splat(1.3));
+ let mut vb = Volume::build(&block, lo, hi, 0.05, 0.15);
+ let vs = Volume::build(&slab, lo, hi, 0.05, 0.15);
+ vb.subtract(&vs);
+ let out = vb.to_mesh();
+ let survivors = (0..out.num_points())
+ .map(|p| out.pos(p))
+ .filter(|q| q.y.abs() < 0.06 && q.x.abs() < 0.4 && q.z.abs() < 0.4)
+ .count();
+ assert_eq!(survivors, 0, "points survive where the slab cut through");
+
+ // And an INSIDE-OUT input gives the same field. The band test asks the
+ // nearest face which way it points, so a mesh wound the other way
+ // would otherwise come back riddled with holes — which is how the
+ // winding measurement got written.
+ let mut flipped = Detail::new();
+ for p in 0..slab.num_points() {
+ flipped.add_point(slab.pos(p));
+ }
+ for prim in 0..slab.num_prims() {
+ let mut pts = slab.prim_points(prim).to_vec();
+ pts.reverse();
+ flipped.add_prim(&pts);
+ }
+ let inverted = Volume::build(&flipped, lo, hi, 0.05, 0.15);
+ let mut differ = 0;
+ for k in 0..vs.dims()[2] {
+ for j in 0..vs.dims()[1] {
+ for i in 0..vs.dims()[0] {
+ if (vs.at(i, j, k) < 0.0) != (inverted.at(i, j, k) < 0.0) {
+ differ += 1;
+ }
+ }
+ }
+ }
+ assert_eq!(differ, 0, "{differ} samples disagree when the input is wound inside out");
+ }
+
+ #[test]
+ fn test_offsetting_is_subtraction() {
+ let sphere = sphere_detail(Vec3::ZERO, 1.0, 14, 20);
+ let radius = |d: &Detail| {
+ (0..d.num_points()).map(|p| d.pos(p).length()).sum::<f32>() / d.num_points() as f32
+ };
+
+ let mut grown = Volume::from_mesh(&sphere, 0.15, 0.6);
+ grown.offset(0.3);
+ let out = grown.to_mesh();
+ assert!(
+ (radius(&out) - 1.3).abs() < 0.12,
+ "a 0.3 offset should give radius 1.3, got {}",
+ radius(&out)
+ );
+
+ // Inward too, which is what a shell's inner wall is.
+ let mut shrunk = Volume::from_mesh(&sphere, 0.15, 0.6);
+ shrunk.offset(-0.3);
+ assert!((radius(&shrunk.to_mesh()) - 0.7).abs() < 0.12);
+ }
+
+ #[test]
+ fn test_the_booleans_are_a_minimum_and_a_maximum() {
+ // Two overlapping spheres, sampled over ONE grid so the operations are
+ // elementwise. Sharing the grid is what makes them arithmetic rather
+ // than a geometry problem.
+ let a = sphere_detail(Vec3::new(-0.35, 0.0, 0.0), 0.8, 14, 20);
+ let b = sphere_detail(Vec3::new(0.35, 0.0, 0.0), 0.8, 14, 20);
+ let (lo, hi) = (Vec3::splat(-1.6), Vec3::splat(1.6));
+ let va = Volume::from_mesh_in(&a, lo, hi, 0.14);
+ let vb = Volume::from_mesh_in(&b, lo, hi, 0.14);
+ assert!(va.aligned_with(&vb), "the two fields do not share a grid");
+
+ let width = |d: &Detail| d.bounds().map(|(l, h)| h.x - l.x).unwrap_or(0.0);
+
+ let mut u = va.clone();
+ u.union(&vb);
+ let mut i = va.clone();
+ i.intersect(&vb);
+ let mut s = va.clone();
+ s.subtract(&vb);
+ // Extracted ONCE each: surface extraction is not free, and an
+ // assertion message that re-runs it is a slow test nobody runs.
+ let (um, im, sm, am) = (u.to_mesh(), i.to_mesh(), s.to_mesh(), va.to_mesh());
+
+ // The union spans both, the intersection is the lens between them, and
+ // the difference is narrower than the whole of A.
+ assert!(width(&um) > 2.2, "union is {}", width(&um));
+ assert!(width(&im) < 1.0, "intersection is {}", width(&im));
+ assert!(width(&sm) < width(&am) + 0.01, "the difference grew");
+ // Every result is still a closed surface — which a mesh boolean has to
+ // work for and a field gets for free.
+ for m in [&um, &im, &sm] {
+ assert!(m.num_prims() > 50);
+ }
+
+ // Fields on different grids refuse to combine rather than reading each
+ // other's memory in the wrong order.
+ let elsewhere = Volume::from_mesh_in(&b, lo, hi, 0.25);
+ assert!(!va.aligned_with(&elsewhere));
+ let mut guarded = va.clone();
+ guarded.union(&elsewhere);
+ assert_eq!(
+ guarded.to_mesh().num_points(),
+ am.num_points(),
+ "a mismatched grid was combined"
+ );
+ }
+
// ---- Mesh export ----
/// A two-quad sheet: enough to tell a format that keeps topology from one
diff --git a/src/spatial.rs b/src/spatial.rs
index 04217dc..c166b0e 100644
--- a/src/spatial.rs
+++ b/src/spatial.rs
@@ -59,6 +59,34 @@ pub fn closest_point_on_triangle(p: Vec3, a: Vec3, b: Vec3, c: Vec3) -> Vec3 {
a + ab * (vb / denom) + ac * (vc / denom)
}
+/// Where a ray meets a triangle, as a distance along the ray.
+///
+/// Möller–Trumbore. Lives here beside the other spatial queries because three
+/// callers want it now: Collision's inside test, and the volume builder's sign
+/// pass, which casts one ray per grid row.
+pub fn ray_triangle(origin: Vec3, dir: Vec3, v0: Vec3, v1: Vec3, v2: Vec3) -> Option<f32> {
+ let edge1 = v1 - v0;
+ let edge2 = v2 - v0;
+ let h = dir.cross(edge2);
+ let a = edge1.dot(h);
+ if a.abs() < 1e-6 {
+ return None;
+ }
+ let f = 1.0 / a;
+ let s = origin - v0;
+ let u = f * s.dot(h);
+ if !(0.0..=1.0).contains(&u) {
+ return None;
+ }
+ let q = s.cross(edge1);
+ let v = f * dir.dot(q);
+ if v < 0.0 || u + v > 1.0 {
+ return None;
+ }
+ let t = f * edge2.dot(q);
+ (t > 1e-5).then_some(t)
+}
+
/// Where things are, bucketed by cell.
///
/// Shared by both grids: they differ only in what they store and how they
@@ -181,15 +209,58 @@ impl TriGrid {
self.tris.is_empty()
}
+ /// The triangles behind this grid, for a caller that needs them directly —
+ /// the volume builder's scanline sign pass casts rays at all of them.
+ pub fn triangles(&self) -> &[[Vec3; 3]] {
+ &self.tris
+ }
+
/// 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<Hit> {
+ self.closest_within(p, f32::INFINITY)
+ }
+
+ /// [`TriGrid::closest`], giving up once the search passes `limit`.
+ ///
+ /// The unbounded form doubles its reach until it finds something, so a
+ /// query far from the surface ends up gathering every triangle in the mesh
+ /// and sorting them — which is fine for the handful of queries an operator
+ /// makes and ruinous for the hundred thousand a volume build makes, where
+ /// most samples are nowhere near the surface. A caller that only needs to
+ /// know "further than this" says so and pays for a few cells.
+ pub fn closest_within(&self, p: Vec3, limit: f32) -> Option<Hit> {
if self.tris.is_empty() {
return None;
}
+ if limit.is_finite() {
+ // ONE gather of exactly the box asked for, rather than doubling up
+ // to it: a bounded query knows how far it cares about, and growing
+ // into that size in stages means gathering and sorting the same
+ // cells over and over. This is the difference between a volume
+ // build taking thirty seconds and taking two.
+ let mut scratch = Vec::new();
+ self.grid
+ .gather(p - Vec3::splat(limit), p + Vec3::splat(limit), &mut scratch);
+ return scratch
+ .iter()
+ .map(|&i| {
+ let t = self.tris[i as usize];
+ let q = closest_point_on_triangle(p, t[0], t[1], t[2]);
+ Hit {
+ point: q,
+ distance: (q - p).length(),
+ normal: (t[1] - t[0]).cross(t[2] - t[0]).normalize_or_zero(),
+ }
+ })
+ .min_by(|a, b| {
+ a.distance.partial_cmp(&b.distance).unwrap_or(std::cmp::Ordering::Equal)
+ })
+ .filter(|h| h.distance <= limit);
+ }
let hit = |i: usize| {
let t = self.tris[i];
let q = closest_point_on_triangle(p, t[0], t[1], t[2]);
diff --git a/src/volume.rs b/src/volume.rs
new file mode 100644
index 0000000..8ac5b0e
--- /dev/null
+++ b/src/volume.rs
@@ -0,0 +1,465 @@
+//! Signed distance fields, and the way back to a surface.
+//!
+//! The representation the manufacturing work stands on. Shelling a shape,
+//! offsetting it, cutting one shape out of another — none of those are natural
+//! on a triangle mesh, where they mean finding every self-intersection the
+//! operation creates and stitching the result back into something closed. On a
+//! distance field they are arithmetic: offsetting is subtraction, union is a
+//! minimum, difference is a maximum against a negation. The mesh comes back at
+//! the end, closed by construction.
+//!
+//! ## Dense, not sparse
+//!
+//! A dense grid over the shape's bounding box, not a sparse tree. At the sizes
+//! this tool works on — a thing you can hold, meshed finely enough to print — a
+//! 128³ grid is eight megabytes and answers every query by indexing. A sparse
+//! structure buys memory back on volumes mostly made of empty space, and costs
+//! a tree walk on every one of the millions of lookups the surface extraction
+//! makes. The day a model needs 512³ is the day to write one.
+//!
+//! ## Two passes, because sign and distance are different questions
+//!
+//! Building the field from a mesh is done twice over:
+//!
+//! 1. **Distance** comes from the triangle grid — the closest point on the
+//! surface, which is exact and needs no assumptions about the mesh.
+//! 2. **Sign** applies only to a CLOSED surface — an open one has no inside,
+//! and is left unsigned so that offsetting thickens it into a slab. Where
+//! there is an inside, it comes from two tests, each used where it is the
+//! accurate one: a
+//! flood outward from the grid boundary decides everything far from the
+//! surface, and the nearest face's normal decides the thin band either side
+//! of it.
+//!
+//! A sign that is wrong anywhere is a bubble or a hole in the result, where in
+//! the Distance node the same mistake was only a slightly wrong number — which
+//! is why this is worth two passes and not one ray cast.
+
+use crate::detail::Detail;
+use crate::spatial::TriGrid;
+use glam::Vec3;
+
+/// A signed distance field on a regular grid. Negative is inside.
+#[derive(Clone, Debug)]
+pub struct Volume {
+ origin: Vec3,
+ voxel: f32,
+ /// Samples per axis, so the last sample sits at `origin + (dims-1) * voxel`.
+ dims: [usize; 3],
+ data: Vec<f32>,
+}
+
+/// How far from the surface distances are measured, in voxels.
+///
+/// Beyond this the field records only "further than this", which is all the
+/// extraction and the flood fill need. It bounds the work per sample, and it
+/// bounds what an offset can do: moving the surface more than this far has
+/// nothing to move into, so a node offsetting further must ask for a bigger
+/// voxel or accept the clamp.
+pub const REACH_VOXELS: f32 = 3.0;
+
+/// How many samples an axis needs to span `extent` at `voxel`, clamped.
+fn axis_dims(extent: f32, voxel: f32) -> usize {
+ ((extent / voxel).ceil() as usize + 3).clamp(2, 256)
+}
+
+impl Volume {
+ pub fn dims(&self) -> [usize; 3] {
+ self.dims
+ }
+
+ pub fn voxel(&self) -> f32 {
+ self.voxel
+ }
+
+ fn index(&self, i: usize, j: usize, k: usize) -> usize {
+ (k * self.dims[1] + j) * self.dims[0] + i
+ }
+
+ pub fn at(&self, i: usize, j: usize, k: usize) -> f32 {
+ self.data[self.index(i, j, k)]
+ }
+
+ /// The world position of a sample. Public so a caller can check a field
+ /// against the shape it was built from.
+ pub fn sample_position(&self, i: usize, j: usize, k: usize) -> Vec3 {
+ self.position(i, j, k)
+ }
+
+ fn position(&self, i: usize, j: usize, k: usize) -> Vec3 {
+ self.origin + Vec3::new(i as f32, j as f32, k as f32) * self.voxel
+ }
+
+ /// The bounds a mesh needs, with room for the offset a caller will apply.
+ ///
+ /// Padding matters: a field built tight to the surface has no room to
+ /// dilate into, and the offset surface would be clipped flat at the edge of
+ /// the grid rather than rounded.
+ pub fn bounds_for(d: &Detail, padding: f32) -> Option<(Vec3, Vec3)> {
+ let (lo, hi) = d.bounds()?;
+ Some((lo - Vec3::splat(padding), hi + Vec3::splat(padding)))
+ }
+
+ /// Sample a mesh into a field over the given bounds.
+ ///
+ /// Two meshes sampled over the SAME bounds at the same voxel size share a
+ /// grid, which is what lets the booleans be elementwise.
+ pub fn from_mesh_in(d: &Detail, lo: Vec3, hi: Vec3, voxel: f32) -> Volume {
+ Self::build(d, lo, hi, voxel, voxel.max(1e-5) * REACH_VOXELS)
+ }
+
+ /// [`Volume::from_mesh_in`] measuring distances out to `reach`.
+ ///
+ /// Beyond `reach` the field says only "further than this", so it bounds
+ /// both the cost and how far an offset can move the surface. A caller that
+ /// means to offset by more must say so here.
+ pub fn build(d: &Detail, lo: Vec3, hi: Vec3, voxel: f32, reach: f32) -> Volume {
+ let voxel = voxel.max(1e-5);
+ let reach = reach.max(voxel * 2.0);
+ let extent = hi - lo;
+ let dims = [
+ axis_dims(extent.x, voxel),
+ axis_dims(extent.y, voxel),
+ axis_dims(extent.z, voxel),
+ ];
+ // Centred, so the padding is even on both sides rather than piling up
+ // wherever the rounding landed.
+ let span = Vec3::new(
+ (dims[0] - 1) as f32,
+ (dims[1] - 1) as f32,
+ (dims[2] - 1) as f32,
+ ) * voxel;
+ let origin = (lo + hi) * 0.5 - span * 0.5;
+
+ let n = dims[0] * dims[1] * dims[2];
+ let mut vol = Volume { origin, voxel, dims, data: vec![f32::MAX; n] };
+ let grid = TriGrid::build(d);
+ if grid.is_empty() {
+ // No surface: everything is outside, at a distance nothing will
+ // mistake for a crossing.
+ vol.data.fill(span.max_element().max(1.0));
+ return vol;
+ }
+
+ // Pass one: unsigned distance, exact from the triangle grid — but only
+ // out to REACH, beyond which the field records "further than this".
+ //
+ // A volume does not need an exact distance far from the surface: the
+ // extraction only looks at the zero crossing and the flood fill only
+ // asks "further than the band". Measuring it anyway is what made this
+ // slow, because an unbounded nearest-surface query from a corner of the
+ // grid gathers every triangle in the mesh. The consequence to know is
+ // that an OFFSET larger than the reach is not represented — which is
+ // why the reach is taken from the padding, the room the caller asked
+ // for in order to offset into.
+ for k in 0..dims[2] {
+ for j in 0..dims[1] {
+ for i in 0..dims[0] {
+ let p = vol.position(i, j, k);
+ let d = grid
+ .closest_within(p, reach)
+ .map(|h| h.distance)
+ .unwrap_or(reach);
+ let idx = vol.index(i, j, k);
+ vol.data[idx] = d;
+ }
+ }
+ }
+
+ // Pass two: sign — but only if there is an inside to speak of.
+ //
+ // "What is inside this?" only has an answer for a CLOSED surface. A
+ // flat disc, a single polygon or a torn mesh has none, and signing one
+ // anyway gives whichever side its normals happen to face: the field
+ // then reads as a half-space and a boolean against it carves something
+ // nobody asked for. Left unsigned, the same disc is the set of points
+ // near it — so an offset thickens it into a slab, which is the useful
+ // answer and the honest one.
+ if !d.is_closed() {
+ return vol;
+ }
+
+ // Which way the mesh is wound, by the divergence theorem: the signed
+ // volume of a closed surface is positive when its faces look outward.
+ //
+ // The band test below asks the nearest face which side a sample is on,
+ // which trusts the winding — and a mesh wound inside out then produces
+ // a field whose band signs alternate against the flood fill's, giving
+ // a surface riddled with holes. Rather than trust it, measure it once
+ // and flip. An imported mesh is not obliged to agree with this app's
+ // convention.
+ let signed_volume: f32 = grid
+ .triangles()
+ .iter()
+ .map(|t| t[0].dot(t[1].cross(t[2])))
+ .sum::<f32>()
+ / 6.0;
+ let facing = if signed_volume < 0.0 { -1.0 } else { 1.0 };
+
+ // Two tests, each used where it is the accurate one.
+ //
+ // Ray parity was the obvious choice and is wrong here: a ray along a
+ // grid row is axis-aligned, the meshes tessellate on the axes, and a
+ // ray through a shared edge hits two triangles at the same point. The
+ // parity flips and STAYS flipped for the rest of the row — which shows
+ // up as contiguous runs of outside samples reading as inside. The
+ // collision resolver already carried a comment about exactly this.
+ //
+ // Instead:
+ //
+ // 1. FAR FROM THE SURFACE, flood outward from the grid's boundary,
+ // which is outside by construction. Anything the flood cannot reach
+ // without crossing the surface band is enclosed, however convoluted
+ // the shape. No rays, no epsilons, no degenerate cases.
+ // 2. IN THE BAND either side of the surface, ask the nearest face which
+ // way it points. That test is unreliable far away in a concave
+ // shape — the nearest face can be one across the gap — and reliable
+ // within half a voxel of the surface, where the nearest face is the
+ // one the sample is sitting on.
+ // A FULL voxel, not a fraction of one. Two samples one voxel apart
+ // cannot both be more than a voxel from the surface AND have the
+ // surface between them — so flooding only between such samples can
+ // never cross it. At 0.75 of a voxel two adjacent samples straddling a
+ // thin wall could both qualify, and the flood walked straight through
+ // into the interior: a slab came back hollow and a boolean against it
+ // cut almost nothing.
+ let band = voxel * 1.01;
+ let idx_of = |i: usize, j: usize, k: usize| (k * dims[1] + j) * dims[0] + i;
+
+ let mut outside = vec![false; n];
+ let mut queue: Vec<(usize, usize, usize)> = Vec::new();
+ for k in 0..dims[2] {
+ for j in 0..dims[1] {
+ for i in 0..dims[0] {
+ let on_boundary = i == 0
+ || j == 0
+ || k == 0
+ || i == dims[0] - 1
+ || j == dims[1] - 1
+ || k == dims[2] - 1;
+ if on_boundary && vol.data[idx_of(i, j, k)] > band {
+ outside[idx_of(i, j, k)] = true;
+ queue.push((i, j, k));
+ }
+ }
+ }
+ }
+ while let Some((i, j, k)) = queue.pop() {
+ let mut visit = |i: usize, j: usize, k: usize, queue: &mut Vec<(usize, usize, usize)>| {
+ let idx = idx_of(i, j, k);
+ if !outside[idx] && vol.data[idx] > band {
+ outside[idx] = true;
+ queue.push((i, j, k));
+ }
+ };
+ if i > 0 { visit(i - 1, j, k, &mut queue); }
+ if j > 0 { visit(i, j - 1, k, &mut queue); }
+ if k > 0 { visit(i, j, k - 1, &mut queue); }
+ if i + 1 < dims[0] { visit(i + 1, j, k, &mut queue); }
+ if j + 1 < dims[1] { visit(i, j + 1, k, &mut queue); }
+ if k + 1 < dims[2] { visit(i, j, k + 1, &mut queue); }
+ }
+
+ for k in 0..dims[2] {
+ for j in 0..dims[1] {
+ for i in 0..dims[0] {
+ let idx = idx_of(i, j, k);
+ let inside = if vol.data[idx] > band {
+ !outside[idx]
+ } else {
+ // In the band: which side of the nearest face.
+ let p = vol.position(i, j, k);
+ match grid.closest_within(p, band * 4.0) {
+ Some(h) => (p - h.point).dot(h.normal) * facing < 0.0,
+ None => false,
+ }
+ };
+ if inside {
+ vol.data[idx] = -vol.data[idx];
+ }
+ }
+ }
+ }
+ vol
+ }
+
+ /// Sample a mesh into a field sized to it, with room to grow.
+ ///
+ /// The padding IS the room to offset into, so it is also the distance the
+ /// field measures out to — asking for space and then not measuring it
+ /// would give an offset nothing to find.
+ pub fn from_mesh(d: &Detail, voxel: f32, padding: f32) -> Volume {
+ let (lo, hi) = Self::bounds_for(d, padding).unwrap_or((Vec3::ZERO, Vec3::ZERO));
+ Self::build(d, lo, hi, voxel, padding)
+ }
+
+ /// Move the surface out (positive) or in (negative).
+ ///
+ /// The whole reason for the representation: an offset on a mesh means
+ /// resolving every self-intersection the move creates, and here it is a
+ /// subtraction.
+ pub fn offset(&mut self, by: f32) {
+ for v in self.data.iter_mut() {
+ *v -= by;
+ }
+ }
+
+ /// Whether two fields share a grid, which the booleans require.
+ pub fn aligned_with(&self, other: &Volume) -> bool {
+ self.dims == other.dims
+ && (self.voxel - other.voxel).abs() < 1e-6
+ && (self.origin - other.origin).length() < 1e-4
+ }
+
+ /// Keep whatever is inside EITHER — the minimum of the two distances.
+ pub fn union(&mut self, other: &Volume) {
+ self.combine(other, |a, b| a.min(b));
+ }
+
+ /// Keep only what is inside BOTH.
+ pub fn intersect(&mut self, other: &Volume) {
+ self.combine(other, |a, b| a.max(b));
+ }
+
+ /// Cut `other` out of this one.
+ pub fn subtract(&mut self, other: &Volume) {
+ self.combine(other, |a, b| a.max(-b));
+ }
+
+ fn combine(&mut self, other: &Volume, f: impl Fn(f32, f32) -> f32) {
+ if !self.aligned_with(other) {
+ return;
+ }
+ for (a, b) in self.data.iter_mut().zip(other.data.iter()) {
+ *a = f(*a, *b);
+ }
+ }
+
+ /// Extract the zero surface as a quad mesh, by naive surface nets.
+ ///
+ /// Surface nets rather than marching cubes: one vertex per cell that the
+ /// surface passes through, placed at the average of that cell's edge
+ /// crossings, and a quad around each grid edge the surface crosses. It is a
+ /// tenth of the code of a correct marching-cubes table, produces quads
+ /// rather than slivers, and cannot be got subtly wrong in one of 256 cases
+ /// nobody exercises. The triangulation it gives is not beautiful, which
+ /// does not matter here: this output goes into the remesher.
+ pub fn to_mesh(&self) -> Detail {
+ let [nx, ny, nz] = self.dims;
+ let mut d = Detail::new();
+ if nx < 2 || ny < 2 || nz < 2 {
+ return d;
+ }
+
+ // One vertex per crossed cell, indexed by cell.
+ let cells = (nx - 1) * (ny - 1) * (nz - 1);
+ let mut vertex = vec![u32::MAX; cells];
+ let cell_index = |i: usize, j: usize, k: usize| (k * (ny - 1) + j) * (nx - 1) + i;
+
+ const CORNERS: [[usize; 3]; 8] = [
+ [0, 0, 0], [1, 0, 0], [1, 1, 0], [0, 1, 0],
+ [0, 0, 1], [1, 0, 1], [1, 1, 1], [0, 1, 1],
+ ];
+ const EDGES: [[usize; 2]; 12] = [
+ [0, 1], [1, 2], [2, 3], [3, 0],
+ [4, 5], [5, 6], [6, 7], [7, 4],
+ [0, 4], [1, 5], [2, 6], [3, 7],
+ ];
+
+ for k in 0..nz - 1 {
+ for j in 0..ny - 1 {
+ for i in 0..nx - 1 {
+ let s: Vec<f32> = CORNERS
+ .iter()
+ .map(|c| self.at(i + c[0], j + c[1], k + c[2]))
+ .collect();
+ if s.iter().all(|v| *v < 0.0) || s.iter().all(|v| *v >= 0.0) {
+ continue;
+ }
+ let mut sum = Vec3::ZERO;
+ let mut hits = 0.0f32;
+ for e in EDGES {
+ let (a, b) = (s[e[0]], s[e[1]]);
+ if (a < 0.0) == (b < 0.0) {
+ continue;
+ }
+ // Where along the edge the field is zero. Linear, which
+ // is exactly right for a field that is a distance.
+ let t = a / (a - b);
+ let pa = self.position(
+ i + CORNERS[e[0]][0],
+ j + CORNERS[e[0]][1],
+ k + CORNERS[e[0]][2],
+ );
+ let pb = self.position(
+ i + CORNERS[e[1]][0],
+ j + CORNERS[e[1]][1],
+ k + CORNERS[e[1]][2],
+ );
+ sum += pa + (pb - pa) * t;
+ hits += 1.0;
+ }
+ if hits > 0.0 {
+ vertex[cell_index(i, j, k)] = d.add_point(sum / hits);
+ }
+ }
+ }
+ }
+
+ // One quad per grid edge the surface crosses, joining the four cells
+ // around that edge. Winding follows the sign: the face must look from
+ // inside to outside, so that the plain cross points away from the
+ // solid like every other generator in the app.
+ let quad = |a: u32, b: u32, c: u32, e: u32, flip: bool, d: &mut Detail| {
+ if [a, b, c, e].iter().any(|&v| v == u32::MAX) {
+ return;
+ }
+ if flip {
+ d.add_prim(&[a, b, c, e]);
+ } else {
+ d.add_prim(&[e, c, b, a]);
+ }
+ };
+ for k in 0..nz - 1 {
+ for j in 0..ny - 1 {
+ for i in 0..nx - 1 {
+ // +X edge: the four cells sharing it differ in j and k.
+ if j > 0 && k > 0 && (self.at(i, j, k) < 0.0) != (self.at(i + 1, j, k) < 0.0) {
+ let inside = self.at(i, j, k) < 0.0;
+ quad(
+ vertex[cell_index(i, j - 1, k - 1)],
+ vertex[cell_index(i, j, k - 1)],
+ vertex[cell_index(i, j, k)],
+ vertex[cell_index(i, j - 1, k)],
+ inside,
+ &mut d,
+ );
+ }
+ if i > 0 && k > 0 && (self.at(i, j, k) < 0.0) != (self.at(i, j + 1, k) < 0.0) {
+ let inside = self.at(i, j, k) < 0.0;
+ quad(
+ vertex[cell_index(i - 1, j, k - 1)],
+ vertex[cell_index(i, j, k - 1)],
+ vertex[cell_index(i, j, k)],
+ vertex[cell_index(i - 1, j, k)],
+ !inside,
+ &mut d,
+ );
+ }
+ if i > 0 && j > 0 && (self.at(i, j, k) < 0.0) != (self.at(i, j, k + 1) < 0.0) {
+ let inside = self.at(i, j, k) < 0.0;
+ quad(
+ vertex[cell_index(i - 1, j - 1, k)],
+ vertex[cell_index(i, j - 1, k)],
+ vertex[cell_index(i, j, k)],
+ vertex[cell_index(i - 1, j, k)],
+ inside,
+ &mut d,
+ );
+ }
+ }
+ }
+ }
+ d
+ }
+}