Display tessellation plans (3/4): display meshes from plans, GPU evaluation in the viewport (#805)

Display meshes are built from tessellation plans, with GPU evaluation in the viewport.

crates/geom/src/mesh.rs

@@ -396,29 +396,59 @@ impl Mesh {
     }
 
     /// Closest intersection of a ray with the mesh: (distance along `dir`, triangle index).
+    ///
+    /// Watertight (Woop, Benthin & Wald, "Watertight Ray/Triangle Intersection", JCGT 2(1),
+    /// 2013): the vertices are moved into a frame where the ray runs along +z through the
+    /// origin, and each edge's side test is computed from those per-vertex values alone, so the
+    /// two triangles that share an edge get exactly opposite values there. A ray through a
+    /// shared edge or vertex hits at least one of them; it never slips between them.
     pub fn raycast(&self, o: Vec3, dir: Vec3) -> Option<(f64, usize)> {
+        let comp = |v: Vec3, k: usize| match k {
+            0 => v.x,
+            1 => v.y,
+            _ => v.z,
+        };
+        // The ray's dominant axis becomes z; x and y follow it cyclically, swapped to keep
+        // the frame right-handed when the ray runs towards −z.
+        let a = [dir.x.abs(), dir.y.abs(), dir.z.abs()];
+        let kz = if a[0] >= a[1] && a[0] >= a[2] {
+            0
+        } else if a[1] >= a[2] {
+            1
+        } else {
+            2
+        };
+        let dz = comp(dir, kz);
+        if !(dz.is_finite() && dz != 0.0) {
+            return None;
+        }
+        let (mut kx, mut ky) = ((kz + 1) % 3, (kz + 2) % 3);
+        if dz < 0.0 {
+            std::mem::swap(&mut kx, &mut ky);
+        }
+        let (sx, sy, sz) = (comp(dir, kx) / dz, comp(dir, ky) / dz, 1.0 / dz);
+        // A vertex in the ray's frame: (x, y) across the ray, z scaled along it.
+        let frame = |p: Vec3| {
+            let q = p - o;
+            let z = comp(q, kz);
+            (comp(q, kx) - sx * z, comp(q, ky) - sy * z, sz * z)
+        };
         let mut best: Option<(f64, usize)> = None;
         for (i, t) in self.triangles.iter().enumerate() {
-            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 {
+            let Some([pa, pb, pc]) = self.tri(t) else { continue };
+            let ((ax, ay, az), (bx, by, bz), (cx, cy, cz)) = (frame(pa), frame(pb), frame(pc));
+            let u = cx * by - cy * bx;
+            let v = ax * cy - ay * cx;
+            let w = bx * ay - by * ax;
+            if (u < 0.0 || v < 0.0 || w < 0.0) && (u > 0.0 || v > 0.0 || w > 0.0) {
                 continue;
             }
-            let inv = 1.0 / det;
-            let s = o - a;
-            let u = s.dot(p) * inv;
-            if !(0.0..=1.0).contains(&u) {
+            let det = u + v + w;
+            if det == 0.0 {
                 continue;
             }
-            let q = s.cross(e1);
-            let v = dir.dot(q) * inv;
-            if v < 0.0 || u + v > 1.0 {
-                continue;
-            }
-            let t = e2.dot(q) * inv;
-            if t > 1e-9 && best.is_none_or(|(bt, _)| t < bt) {
+            let t = (u * az + v * bz + w * cz) / det;
+            if t.is_finite() && t > 1e-9 && best.is_none_or(|(bt, _)| t < bt) {
                 best = Some((t, i));
             }
         }