diff --git a/crates/geom/src/lib.rs b/crates/geom/src/lib.rs index dc2cbbdd..77991ef6 100644 --- a/crates/geom/src/lib.rs +++ b/crates/geom/src/lib.rs @@ -11,7 +11,7 @@ mod plane; mod profile; mod vec; -pub use mesh::{Aabb3, Mesh, MeshMeasure}; +pub use mesh::{Aabb3, InsideIndex, Mesh, MeshMeasure}; pub use plane::Plane; pub use profile::{Loop2, MAX_BEZIER_DEGREE, Region2, Seg2, bezier}; pub use vec::{Vec2, Vec3}; diff --git a/crates/geom/src/mesh.rs b/crates/geom/src/mesh.rs index 8ad84338..32c3373a 100644 --- a/crates/geom/src/mesh.rs +++ b/crates/geom/src/mesh.rs @@ -227,3 +227,122 @@ impl Mesh { best } } + +/// Point-in-mesh queries in bulk: for each of `Mesh::contains`'s three ray directions, the +/// triangles binned on a grid in the plane square to it, so a ray tests only the triangles +/// whose shadow covers its cell. Answers exactly as `Mesh::contains` does. +pub struct InsideIndex<'a> { + mesh: &'a Mesh, + rays: Vec, +} + +struct RayGrid { + dir: Vec3, + u: Vec3, + v: Vec3, + min: (f64, f64), + cell: f64, + n: (usize, usize), + cells: Vec>, +} + +const RAY_DIRS: [Vec3; 3] = [Vec3::new(0.5773, 0.5774, 0.5776), Vec3::new(-0.6123, 0.3141, 0.7254), Vec3::new(0.2718, -0.8414, 0.4673)]; + +impl Mesh { + /// An index for many `contains` queries on this mesh. + pub fn inside_index(&self) -> InsideIndex<'_> { + let rays = RAY_DIRS + .iter() + .map(|d| { + let u = d.any_perp(); + let v = d.cross(u); + let tris: Vec<[(f64, f64); 3]> = self + .triangles + .iter() + .map(|t| match self.tri(t) { + Some([a, b, c]) => [(a.dot(u), a.dot(v)), (b.dot(u), b.dot(v)), (c.dot(u), c.dot(v))], + None => [(0.0, 0.0); 3], + }) + .collect(); + let (mut lo, mut hi) = ((f64::INFINITY, f64::INFINITY), (f64::NEG_INFINITY, f64::NEG_INFINITY)); + for t in &tris { + for p in t { + lo = (lo.0.min(p.0), lo.1.min(p.1)); + hi = (hi.0.max(p.0), hi.1.max(p.1)); + } + } + let span = (hi.0 - lo.0).max(hi.1 - lo.1).max(1e-9); + let side = ((tris.len() as f64).sqrt().ceil() as usize).clamp(1, 256); + let cell = span / side as f64 * (1.0 + 1e-9); + let n = (side + 1, side + 1); + let mut cells = vec![Vec::new(); n.0 * n.1]; + if lo.0.is_finite() { + for (i, t) in tris.iter().enumerate() { + let (x0, x1) = (t.iter().map(|p| p.0).fold(f64::INFINITY, f64::min), t.iter().map(|p| p.0).fold(f64::NEG_INFINITY, f64::max)); + let (y0, y1) = (t.iter().map(|p| p.1).fold(f64::INFINITY, f64::min), t.iter().map(|p| p.1).fold(f64::NEG_INFINITY, f64::max)); + let c = |x: f64, o: f64, m: usize| (((x - o) / cell).floor().max(0.0) as usize).min(m - 1); + for gx in c(x0, lo.0, n.0)..=c(x1, lo.0, n.0) { + for gy in c(y0, lo.1, n.1)..=c(y1, lo.1, n.1) { + if let Some(slot) = cells.get_mut(gx * n.1 + gy) { + slot.push(i as u32); + } + } + } + } + } + RayGrid { dir: *d, u, v, min: lo, cell, n, cells } + }) + .collect(); + InsideIndex { mesh: self, rays } + } +} + +impl InsideIndex<'_> { + fn crossings(&self, g: &RayGrid, o: Vec3) -> usize { + let (x, y) = (o.dot(g.u), o.dot(g.v)); + if !(g.min.0.is_finite()) { + return 0; + } + let (fx, fy) = ((x - g.min.0) / g.cell, (y - g.min.1) / g.cell); + if fx < 0.0 || fy < 0.0 { + return 0; + } + let (gx, gy) = (fx.floor() as usize, fy.floor() as usize); + if gx >= g.n.0 || gy >= g.n.1 { + return 0; + } + let Some(list) = g.cells.get(gx * g.n.1 + gy) else { return 0 }; + let dir = g.dir; + let mut n = 0; + for &i in list { + let Some(t) = self.mesh.triangles.get(i as usize) else { continue }; + let Some([a, b, c]) = self.mesh.tri(t) else { continue }; + let (e1, e2) = (b - a, c - a); + let p = dir.cross(e2); + let det = e1.dot(p); + if det.abs() < 1e-14 { + continue; + } + let inv = 1.0 / det; + let s = o - a; + let u = s.dot(p) * inv; + if !(0.0..1.0).contains(&u) { + continue; + } + let q = s.cross(e1); + let v = dir.dot(q) * inv; + if v < 0.0 || u + v >= 1.0 { + continue; + } + if e2.dot(q) * inv > 0.0 { + n += 1; + } + } + n + } + + /// Same answer as `Mesh::contains`. + pub fn contains(&self, p: Vec3) -> bool { + self.rays.iter().filter(|g| self.crossings(g, p) % 2 == 1).count() >= 2 + } +} diff --git a/crates/kernel/src/ops.rs b/crates/kernel/src/ops.rs index 3e34a144..2f97db93 100644 --- a/crates/kernel/src/ops.rs +++ b/crates/kernel/src/ops.rs @@ -49,19 +49,20 @@ fn mismatch(a: &solvecraft_geom::Mesh, b: &solvecraft_geom::Mesh, result: &solve if lo.x < hi.x && lo.y < hi.y && lo.z < hi.z { boxes.push((solvecraft_geom::Aabb3 { min: lo, max: hi }, 900)); } + let (ia_idx, ib_idx, ir_idx) = (a.inside_index(), b.inside_index(), result.inside_index()); let (mut n, mut bad) = (0usize, 0usize); for (bx, count) in boxes { let s = bx.size(); for i in 1..=count { let p = bx.min + Vec3::new(s.x * halton(i, 2), s.y * halton(i, 3), s.z * halton(i, 5)); - let (ia, ib) = (a.contains(p), b.contains(p)); + let (ia, ib) = (ia_idx.contains(p), ib_idx.contains(p)); let want = match op { BoolOp::Union => ia || ib, BoolOp::Cut => ia && !ib, BoolOp::Intersect => ia && ib, }; n += 1; - if result.contains(p) != want { + if ir_idx.contains(p) != want { bad += 1; } } @@ -75,23 +76,24 @@ fn overlap(a: &solvecraft_geom::Mesh, b: &solvecraft_geom::Mesh) -> Option<(f64, if a.triangles.len() + b.triangles.len() > 80_000 { return None; } - let frac = |inner: &solvecraft_geom::Mesh, outer: &solvecraft_geom::Mesh| -> Option { + let (ia, ib) = (a.inside_index(), b.inside_index()); + let frac = |inner: &solvecraft_geom::Mesh, inner_i: &solvecraft_geom::InsideIndex, outer_i: &solvecraft_geom::InsideIndex| -> Option { let bx = inner.bounds(); let s = bx.size(); let (mut n, mut hit) = (0usize, 0usize); for i in 1..=2000 { let p = bx.min + Vec3::new(s.x * halton(i, 2), s.y * halton(i, 3), s.z * halton(i, 5)); - if !inner.contains(p) { + if !inner_i.contains(p) { continue; } n += 1; - if outer.contains(p) { + if outer_i.contains(p) { hit += 1; } } (n >= 50).then(|| hit as f64 / n as f64) }; - Some((frac(b, a)?, frac(a, b)?)) + Some((frac(b, &ib, &ia)?, frac(a, &ia, &ib)?)) } /// Booleans that need no intersection: one body clear of the other, or inside it. These are @@ -163,22 +165,32 @@ pub fn boolean(a: &Body, b: &Body, op: BoolOp) -> Result> { a.require_brep("a boolean")?; b.require_brep("a boolean")?; let size = a.size().max(b.size()); - let (va, vb) = (volume(a), volume(b)); - let slack = 2e-3 * (va + vb) + 1e-9; let tol_m = size * 5e-4; let meshes = (a.tessellate(tol_m).ok(), b.tessellate(tol_m).ok()); - let plausible = |v: f64, r: &Body| { - let bounds = match op { - BoolOp::Union => v >= va.max(vb) - slack && v <= va + vb + slack, - BoolOp::Cut => v >= va - vb - slack && v <= va + slack, - BoolOp::Intersect => v <= va.min(vb) + slack, - }; - // Membership check against the operands (catches misclassified pieces). - let consistent = match (&meshes, r.tessellate(tol_m)) { - ((Some(ma), Some(mb)), Ok(mr)) => mismatch(ma, mb, &mr, op).is_none_or(|f| f < 0.004), + let (va, vb) = match &meshes { + (Some(ma), Some(mb)) => (ma.measure().volume, mb.measure().volume), + _ => (volume(a), volume(b)), + }; + let slack = 2e-3 * (va + vb) + 1e-9; + // Volume bounds, then membership against the operands (catches misclassified pieces), on + // one mesh of the result. Returns (plausible, volume). + let check = |r: &Body| -> (bool, f64) { + let Ok(mr) = r.tessellate(tol_m) else { return (false, f64::NAN) }; + let v = mr.measure().volume; + let bounds = v > 0.0 + && match op { + BoolOp::Union => v >= va.max(vb) - slack && v <= va + vb + slack, + BoolOp::Cut => v >= va - vb - slack && v <= va + slack, + BoolOp::Intersect => v <= va.min(vb) + slack, + }; + if !bounds { + return (false, v); + } + let consistent = match &meshes { + (Some(ma), Some(mb)) => mismatch(ma, mb, &mr, op).is_none_or(|f| f < 0.004), _ => true, }; - bounds && consistent + (consistent, v) }; // One body clear of or inside the other: no intersection needed. if let (Some(ma), Some(mb)) = &meshes @@ -212,8 +224,8 @@ pub fn boolean(a: &Body, b: &Body, op: BoolOp) -> Result> { && b.solid.boundaries().len() == 1 && let Some(Ok(body)) = crate::coplanar::glue(a, b) { - let v = volume(&body); - if v > 0.0 && plausible(v, &body) { + let (ok, _) = check(&body); + if ok { return Ok(Some(body)); } } @@ -222,8 +234,8 @@ pub fn boolean(a: &Body, b: &Body, op: BoolOp) -> Result> { let r = crate::coplanar::boolean_apart(a, b, op); APART.with(|c| c.set(false)); if let Some(Ok(Some(body))) = r { - let v = volume(&body); - if v > 0.0 && plausible(v, &body) { + let (ok, _) = check(&body); + if ok { return Ok(Some(body)); } } @@ -231,7 +243,7 @@ pub fn boolean(a: &Body, b: &Body, op: BoolOp) -> Result> { let mut empty_votes = 0; for (attempt, j) in JITTER.iter().enumerate() { let shift = mt::Vector3::new(j[0], j[1], j[2]) * (size * 0.01); - for k in [5e-4, 2e-3] { + for k in [2e-3, 5e-4] { let tol = (size * k).max(1e-5); let r = guard("boolean", || { let (mut sa, mut sb) = (a.deep_copy(), b.deep_copy()); @@ -257,8 +269,8 @@ pub fn boolean(a: &Body, b: &Body, op: BoolOp) -> Result> { } Ok(s) => match Body::new(crate::heal::heal(s, size)) { Ok(body) => { - let v = volume(&body); - if v > 0.0 && plausible(v, &body) { + let (ok, v) = check(&body); + if ok { return Ok(Some(body)); } last = format!("implausible result volume {v:.4}"); @@ -279,8 +291,8 @@ pub fn boolean(a: &Body, b: &Body, op: BoolOp) -> Result> { APART.with(|c| c.set(false)); match r { Some(Ok(Some(body))) => { - let v = volume(&body); - if v > 0.0 && plausible(v, &body) { + let (ok, v) = check(&body); + if ok { return Ok(Some(body)); } last = format!("{last}; faces pushed apart gave an implausible volume {v:.4}"); @@ -294,8 +306,8 @@ pub fn boolean(a: &Body, b: &Body, op: BoolOp) -> Result> { // polygon fallback. match crate::polybool::planar_boolean(a, b, op) { Ok(Some(body)) => { - let v = volume(&body); - if v > 0.0 && plausible(v, &body) { + let (ok, v) = check(&body); + if ok { return Ok(Some(body)); } last = format!("planar fallback gave an implausible volume {v:.4}"); @@ -313,8 +325,8 @@ pub fn boolean(a: &Body, b: &Body, op: BoolOp) -> Result> { APART.with(|c| c.set(false)); match r { Ok(Some(body)) => { - let v = volume(&body); - if v > 0.0 && plausible(v, &body) { + let (ok, v) = check(&body); + if ok { return Ok(Some(body)); } last = format!("{last}; the nudged tool gave an implausible volume {v:.4}"); diff --git a/crates/kernel/src/topo.rs b/crates/kernel/src/topo.rs index 960a5759..0d1d3ee3 100644 --- a/crates/kernel/src/topo.rs +++ b/crates/kernel/src/topo.rs @@ -411,7 +411,7 @@ pub fn merged_topology(b: &Body, mesh: &Mesh) -> Result { match (sa.search_nearest_parameter(q, None, 50), sb.search_nearest_parameter(q, None, 50)) { (Some((u, v)), Some((s, t))) => { let (x, y) = (sa.normal(u, v), sb.normal(s, t)); - x.x * y.x + x.y * y.y + x.z * y.z > 1.0 - 1e-4 + x.x * y.x + x.y * y.y + x.z * y.z > 1.0 - 1e-3 } _ => false, }