Tessellation: point-in-trim tests skip runs of segments the ray canno… (#1309)

Thanks @Ian-cmd-ops.

vendor/README.md

@@ -142,6 +142,11 @@ Changes:
   `analyzers::collision_test_hooks`, behind a `test-hooks` feature that only
   `solidcraft-kernel`'s tests turn on, and `crates/kernel/tests/mesh_collision_equivalence.rs`
   compares the two.
+- `tessellation/triangulation.rs`: the inside test of a face's boundary also skips runs of 32
+  segments whose box, grown by the tolerance, the test's ray cannot meet. Same answers. The
+  old test is kept as `tessellation::include_test_hooks`, behind a `test-hooks` feature that
+  only `solidcraft-kernel`'s tests turn on, and `crates/kernel/tests/tess_include_equivalence.rs`
+  compares the two.
 - `lib.rs`: compiler warnings no longer fail the build. Example, test and bench targets are
   dropped (their sources are not vendored).
 

vendor/truck-meshalgo/src/tessellation/mod.rs

@@ -509,3 +509,6 @@ impl<C: PolylineableCurve, S: RobustMeshableSurface> RobustMeshableShape
 }
 
 mod triangulation;
+#[cfg(feature = "test-hooks")]
+#[doc(hidden)]
+pub use triangulation::include_test_hooks;

vendor/truck-meshalgo/src/tessellation/triangulation.rs

@@ -495,9 +495,49 @@ fn get_mindiff(u: f64, u0: f64, up: f64) -> f64 {
 }
 
 /// SolidCraft: the second field holds each loop's parameter box, so `include` skips loops that
-/// cannot hold the point (a face with hundreds of holes meshed in quadratic time).
+/// cannot hold the point (a face with hundreds of holes meshed in quadratic time). The third
+/// holds, for each loop, the box of each run of [`RUN`] segments in a row, so `include` tests
+/// only the runs its ray can meet (a long trim loop is otherwise crossed segment by segment
+/// for every grid point and every triangle).
 #[derive(Debug, Default, Clone)]
-struct PolyBoundary(Vec<Vec<Point2>>, Vec<BoundingBox<Point2>>);
+struct PolyBoundary(Vec<Vec<Point2>>, Vec<BoundingBox<Point2>>, Vec<Vec<BoundingBox<Point2>>>);
+
+/// Segments per run in [`PolyBoundary`]'s third field.
+const RUN: usize = 32;
+
+/// The box of each run of [`RUN`] segments of a closed loop (segment `i` joins point `i` to
+/// point `i + 1`, the last back to the first).
+fn run_boxes(curve: &[Point2]) -> Vec<BoundingBox<Point2>> {
+    let n = curve.len();
+    (0..n.div_ceil(RUN))
+        .map(|k| (k * RUN..((k + 1) * RUN).min(n)).flat_map(|i| [curve[i], curve[(i + 1) % n]]).collect())
+        .collect()
+}
+
+/// Whether the ray from `c` along `r`, started `margin` behind `c`, can meet `bdb` grown by
+/// `margin` on every side. A segment the ray crosses (or passes within the tolerance of `c`)
+/// has its crossing point in its box, so a run whose box the ray misses changes nothing.
+fn ray_meets_box(c: Point2, r: Vector2, bdb: &BoundingBox<Point2>, margin: f64) -> bool {
+    let (min, max) = (bdb.min(), bdb.max());
+    let (mut t0, mut t1) = (-margin, f64::INFINITY);
+    for i in 0..2 {
+        let (lo, hi) = (min[i] - margin, max[i] + margin);
+        if r[i].abs() < 1e-300 {
+            if c[i] < lo || c[i] > hi {
+                return false;
+            }
+            continue;
+        }
+        let (a, b) = ((lo - c[i]) / r[i], (hi - c[i]) / r[i]);
+        let (a, b) = if a <= b { (a, b) } else { (b, a) };
+        t0 = t0.max(a);
+        t1 = t1.min(b);
+        if t0 > t1 {
+            return false;
+        }
+    }
+    true
+}
 
 fn normalize_range(curve: &mut Vec<Point2>, compidx: usize, (u0, u1): (f64, f64)) {
     let p = curve[0];
@@ -649,7 +689,8 @@ impl PolyBoundary {
             }
         }
         let boxes = closed.iter().map(|curve| curve.iter().collect()).collect();
-        Self(closed, boxes)
+        let runs = closed.iter().map(|curve| run_boxes(curve)).collect();
+        Self(closed, boxes, runs)
     }
 
     /// whether `c` is included in the domain with boundary = `self`.
@@ -663,11 +704,22 @@ impl PolyBoundary {
             let (min, max) = (bdb.min(), bdb.max());
             !(c.x < min.x - margin || c.x > max.x + margin || c.y < min.y - margin || c.y > max.y + margin)
         };
+        // Wide enough for every crossing the test below can find: within the tolerance of `c`
+        // (the boundary case) and anywhere along the ray, with room for rounding.
+        let scale = c.x.abs().max(c.y.abs()).max(1.0);
+        let ray_margin = 4.0 * TOLERANCE + 1e-9 * scale;
         self.0
             .iter()
             .zip(&self.1)
-            .filter(|(_, bdb)| near(bdb))
-            .flat_map(|(vec, _)| vec.iter().circular_tuple_windows())
+            .zip(&self.2)
+            .filter(|((_, bdb), _)| near(bdb))
+            .flat_map(|((vec, _), runs)| {
+                let n = vec.len();
+                runs.iter()
+                    .enumerate()
+                    .filter(move |(_, rb)| ray_meets_box(c, r, rb, ray_margin))
+                    .flat_map(move |(k, _)| (k * RUN..((k + 1) * RUN).min(n)).map(move |i| (&vec[i], &vec[(i + 1) % n])))
+            })
             .try_fold(0_i32, move |counter, (p0, p1)| {
                 let a = p0 - c;
                 let b = p1 - c;
@@ -1192,3 +1244,125 @@ fn par_bench() {
     });
     println!("{}ms", instant.elapsed().as_millis());
 }
+
+/// `include` before the runs: kept to check the new one answers the same everywhere.
+///
+/// Not part of the API. Built only for tests, or with the `test-hooks` feature, so the
+/// SolidCraft kernel tests can run this check in CI (this crate's own tests do not run there).
+#[cfg(any(test, feature = "test-hooks"))]
+#[doc(hidden)]
+pub mod include_test_hooks {
+    use super::*;
+
+    fn old_include(b: &PolyBoundary, c: Point2) -> bool {
+        let t = 2.0 * std::f64::consts::PI * HashGen::hash1(c);
+        let r = Vector2::new(f64::cos(t), f64::sin(t));
+        let margin = 2.0 * TOLERANCE;
+        let near = |bdb: &BoundingBox<Point2>| {
+            let (min, max) = (bdb.min(), bdb.max());
+            !(c.x < min.x - margin || c.x > max.x + margin || c.y < min.y - margin || c.y > max.y + margin)
+        };
+        b.0.iter()
+            .zip(&b.1)
+            .filter(|(_, bdb)| near(bdb))
+            .flat_map(|(vec, _)| vec.iter().circular_tuple_windows())
+            .try_fold(0_i32, move |counter, (p0, p1): (&Point2, &Point2)| {
+                let a = p0 - c;
+                let b = p1 - c;
+                let s0 = r.x * a.y - r.y * a.x;
+                let s1 = r.x * b.y - r.y * b.x;
+                let s2 = a.x * b.y - a.y * b.x;
+                let x = s2 / (s1 - s0);
+                if x.so_small() && s0 * s1 < 0.0 {
+                    None
+                } else if x > 0.0 && s0 <= 0.0 && s1 > 0.0 {
+                    Some(counter + 1)
+                } else if x > 0.0 && s0 >= 0.0 && s1 < 0.0 {
+                    Some(counter - 1)
+                } else {
+                    Some(counter)
+                }
+            })
+            .map(|counter| counter > 0)
+            .unwrap_or(false)
+    }
+
+    struct Lcg(u64);
+    impl Lcg {
+        fn unit(&mut self) -> f64 {
+            self.0 = self.0.wrapping_mul(6364136223846793005).wrapping_add(1442695040888963407);
+            (self.0 >> 11) as f64 / (1u64 << 53) as f64
+        }
+    }
+
+    fn boundary(loops: Vec<Vec<Point2>>) -> PolyBoundary {
+        let boxes = loops.iter().map(|l| l.iter().collect()).collect();
+        let runs = loops.iter().map(|l| run_boxes(l)).collect();
+        PolyBoundary(loops, boxes, runs)
+    }
+
+    /// The point-in-trim test with run boxes answers as the old one did, on random loops
+    /// (some on a grid) probed at random, on vertices, on edges and just either side of them.
+    /// Returns how many points were checked.
+    pub fn include_answers_as_before(seeds: std::ops::Range<u64>) -> std::result::Result<usize, String> {
+        let mut checked = 0usize;
+        for seed in seeds {
+            let mut g = Lcg(seed * 104729 + 3);
+            // A wobbly outer loop and a few holes, up to 300 segments each, some on a grid.
+            let mut loops = Vec::new();
+            for l in 0..1 + (g.unit() * 4.0) as usize {
+                let n = 3 + (g.unit() * 300.0) as usize;
+                let (cx, cy, rad) = if l == 0 { (0.0, 0.0, 10.0) } else { (g.unit() * 8.0 - 4.0, g.unit() * 8.0 - 4.0, 0.5 + g.unit() * 2.0) };
+                let grid = g.unit() < 0.3;
+                let mut lp: Vec<Point2> = (0..n)
+                    .map(|i| {
+                        let a = 2.0 * std::f64::consts::PI * i as f64 / n as f64;
+                        let rr = rad * (0.7 + 0.3 * g.unit());
+                        let (x, y) = (cx + rr * a.cos(), cy + rr * a.sin());
+                        if grid { Point2::new((x * 4.0).round() / 4.0, (y * 4.0).round() / 4.0) } else { Point2::new(x, y) }
+                    })
+                    .collect();
+                if l > 0 {
+                    lp.reverse();
+                }
+                loops.push(lp);
+            }
+            let b = boundary(loops.clone());
+            let mut wrong = None;
+            let mut probe = |c: Point2| {
+                if wrong.is_none() && b.include(c) != old_include(&b, c) {
+                    wrong = Some(c);
+                }
+                checked += 1;
+            };
+            for _ in 0..400 {
+                probe(Point2::new(g.unit() * 24.0 - 12.0, g.unit() * 24.0 - 12.0));
+            }
+            for lp in &loops {
+                for i in 0..lp.len() {
+                    let (p, q) = (lp[i], lp[(i + 1) % lp.len()]);
+                    let t = g.unit();
+                    let on = p + (q - p) * t;
+                    let d = q - p;
+                    let nrm = Vector2::new(-d.y, d.x) * (1.0 / d.magnitude().max(1e-300));
+                    probe(p);
+                    probe(on);
+                    probe(on + nrm * (TOLERANCE * 0.5));
+                    probe(on + nrm * (TOLERANCE * 3.0));
+                    probe(on - nrm * (TOLERANCE * 3.0));
+                }
+            }
+            if let Some(c) = wrong {
+                return Err(format!("seed {seed}: point {c:?} answered differently"));
+            }
+        }
+        Ok(checked)
+    }
+
+    #[cfg(test)]
+    #[test]
+    fn hook_passes() {
+        let n = include_answers_as_before(0..120);
+        assert!(matches!(n, Ok(n) if n > 100_000), "{n:?}");
+    }
+}