Mesh: a ray through an edge two triangles share crosses one of them

Mesh::ray_crossings and InsideIndex counted a triangle crossed by a
Möller-Trumbore test with half-open barycentric bounds (u in [0, 1), u + v < 1).
Those bounds depend on each triangle's own vertex order, so on an edge two
triangles share the ray was in both or in neither: the bracket example's hole
centre lies on a diagonal of the base's top triangulation and was counted 0 or 2
times instead of once (bracket_example_cuts_two_through_holes_and_reopens).

Both now use ray_crosses, a watertight test: each edge's side of the ray is the
signed volume dir · ((p - o) × (q - o)), which is exactly antisymmetric in
floating point, so the two triangles of an edge see it with opposite signs, and a
ray exactly on an edge's line is put on one side by the edge's endpoints in
lexicographic order. New geom test ray_crossings_count_a_shared_edge_once crosses
a grid at every inner edge, straight down and obliquely.

Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>

crates/geom/src/mesh.rs

@@ -600,33 +600,10 @@ impl Mesh {
         self.edges.extend(o.edges.iter().cloned());
     }
 
-    /// Number of triangles a ray crosses (for inside/outside tests on closed meshes).
+    /// Number of triangles a ray crosses (for inside/outside tests on closed meshes), each
+    /// crossing counted once ([`ray_crosses`]).
     pub fn ray_crossings(&self, o: Vec3, dir: Vec3) -> usize {
-        let mut n = 0;
-        for t in &self.triangles {
-            let Some([a, b, c]) = self.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
+        self.triangles.iter().filter_map(|t| self.tri(t)).filter(|&[a, b, c]| ray_crosses(a, b, c, o, dir)).count()
     }
 
     /// Point inside a closed mesh (majority vote of three skewed rays).
@@ -854,33 +831,11 @@ impl InsideIndex<'_> {
             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
+        list.iter()
+            .filter_map(|&i| self.mesh.triangles.get(i as usize))
+            .filter_map(|t| self.mesh.tri(t))
+            .filter(|&[a, b, c]| ray_crosses(a, b, c, o, g.dir))
+            .count()
     }
 
     /// Same answer as `Mesh::contains`.
@@ -888,3 +843,23 @@ impl InsideIndex<'_> {
         self.rays.iter().filter(|g| self.crossings(g, p) % 2 == 1).count() >= 2
     }
 }
+
+/// Whether the ray from `o` along `dir` crosses the triangle `a b c` in front of `o`. Watertight:
+/// a ray through an edge two triangles share crosses exactly one of them. Each edge's side of the
+/// ray is the signed volume `dir · ((p - o) × (q - o))`, exactly antisymmetric in floating point
+/// (swapping `p` and `q` negates every rounded term), so both triangles see one value with
+/// opposite signs; a ray exactly on the edge's line is put on one side by the edge's endpoints
+/// in lexicographic order, the same for both. A triangle the ray runs along is not crossed.
+fn ray_crosses(a: Vec3, b: Vec3, c: Vec3, o: Vec3, dir: Vec3) -> bool {
+    let side = |p: Vec3, q: Vec3| {
+        let v = dir.dot((p - o).cross(q - o));
+        if v != 0.0 { v > 0.0 } else { (p.x, p.y, p.z) < (q.x, q.y, q.z) }
+    };
+    let normal = (b - a).cross(c - a);
+    let along = dir.dot(normal);
+    if along == 0.0 || !along.is_finite() {
+        return false;
+    }
+    let s = side(a, b);
+    side(b, c) == s && side(c, a) == s && (a - o).dot(normal) / along > 0.0
+}

crates/geom/src/tests.rs

@@ -97,6 +97,44 @@ fn raycast_never_slips_between_triangles() {
     }
 }
 
+/// A ray through an edge that two triangles share crosses exactly one of them: the same grid,
+/// crossed at every inner edge's midpoint (and the middle of each cell, on its diagonal), straight
+/// down and obliquely. The bracket example's hole centre fell on such an edge of the base's top
+/// and was counted once for two faces.
+#[test]
+fn ray_crossings_count_a_shared_edge_once() {
+    let n = 6;
+    let positions: Vec<Vec3> =
+        (0..=n).flat_map(|j| (0..=n).map(move |i| Vec3::new(f64::from(i) * 40.0 / 6.0 + 100.0, f64::from(j) * 5.0, 20.0))).collect();
+    let id = |i: u32, j: u32| j * (n + 1) + i;
+    let mut triangles = Vec::new();
+    for j in 0..n {
+        for i in 0..n {
+            triangles.push([id(i, j), id(i + 1, j), id(i + 1, j + 1)]);
+            triangles.push([id(i, j), id(i + 1, j + 1), id(i, j + 1)]);
+        }
+    }
+    let tri_face = vec![0; triangles.len()];
+    let m = Mesh { positions, triangles, tri_face, ..Default::default() };
+    let dirs = [Vec3::new(-0.57735, 0.57735, -0.57735), Vec3::new(0.3, -0.7, -0.648), Vec3::new(-0.123456789, -0.2, -1.0), Vec3::new(0.0, 0.0, -1.0)];
+    for d in dirs {
+        for j in 1..2 * n {
+            for i in 1..2 * n {
+                if i % 2 == 0 && j % 2 == 0 {
+                    continue;
+                }
+                let p = Vec3::new(f64::from(i) * 20.0 / 6.0 + 100.0, f64::from(j) * 2.5, 20.0);
+                assert_eq!(m.ray_crossings(p - d * 50.0, d), 1, "{p:?} along {d:?}");
+                assert_eq!(m.ray_crossings(p + d * 50.0, d), 0, "{p:?} along {d:?}, from beyond");
+            }
+        }
+    }
+    // A closed cube: a ray through its top's diagonal crosses the top and the bottom.
+    let cube = cube_mesh(2.0);
+    assert_eq!(cube.ray_crossings(Vec3::new(0.5, 0.5, 5.0), Vec3::new(0.0, 0.0, -1.0)), 2);
+    assert_eq!(cube.ray_crossings(Vec3::new(0.5, 0.5, 1.0), Vec3::new(0.0, 0.0, -1.0)), 1);
+}
+
 #[test]
 fn planes() {
     assert_eq!(Plane::XY.normal(), Vec3::Z);