Fix projection hang, parser panic, nondeterminism and sphere faces (ENG-3712)

nurbs/src/sampled_surface.rs

@@ -26,6 +26,13 @@ pub struct ProjectionScratch<const N: usize> {
 }
 
 const PROJECTION_TOL: f64 = 64. * f64::EPSILON;
+/// Patches one projection may examine before it settles for its incumbent.
+/// Ordinary queries examine at most a few hundred, and tens of thousands on
+/// the worst surfaces seen. When the incumbent is the pole of a collapsed
+/// edge, every patch along that edge reaches it: each subdivision halves
+/// their gap to the incumbent but doubles their number, so certifying it at
+/// f64 resolution would examine millions of patches for a single point.
+const MAX_PROJECTION_PATCHES: usize = 1 << 16;
 
 #[derive(Debug, Clone)]
 struct SurfaceCell<const N: usize> {
@@ -544,7 +551,14 @@ where
                     cell.spans));
             }
         }
+        let mut examined = 0;
         while let Some((patch, lo, hi, spans)) = queue.pop_front() {
+            // The incumbent is a feasible point no farther than any patch
+            // evaluated so far; keep it rather than refine without end.
+            examined += 1;
+            if examined > MAX_PROJECTION_PATCHES {
+                break;
+            }
             let mid = (lo + hi) * 0.5;
             let d = self.surf.derivs_in_span::<2>(mid, spans, p);
             let normal = d[1][0].cross(&d[0][1]);
@@ -682,6 +696,38 @@ mod tests {
         close((sampled.surf.point(uv)-target).norm_squared(), 0.02, 1e-12);
     }
 
+    #[test]
+    fn projection_near_a_collapsed_edge_terminates() {
+        // Face #62825 of a Samtec FMC connector (ASP-134486-01): the u = max
+        // row of control points collapses to one pole. Every patch along that
+        // edge reaches the pole, so certifying this projection at f64
+        // resolution examined millions of patches.
+        let rows = [
+            [(-5.589503327713219782, 3.795967870219560059, 24.88456541861770077, 1.),
+             (-5.632056464503211402, 3.715554243323587080, 24.90841966077580594, 0.8294535963452530636),
+             (-5.582626279250799683, 3.666787751982639953, 24.97187345677104986, 1.)],
+            [(-5.569374495857790208, 3.760876034061380047, 24.84803886198849909, 1.),
+             (-5.599499045440100353, 3.721231129362655210, 24.88421881863200369, 0.9379770243229920146),
+             (-5.582626266339880239, 3.666787755634880241, 24.90864600195275003, 1.)],
+            [(-5.570182848018389699, 3.711200168230984531, 24.80680506253315087, 1.),
+             (-5.574716636290549587, 3.689228228810422117, 24.82235407503383229, 1.000927775394292008),
+             (-5.582626295171389863, 3.666787741024605030, 24.83698483065439788, 1.)],
+            [(-5.582626325353579944, 3.666787764895184765, 24.76551111168974728, 1.); 3],
+        ];
+        let surf = NDBSplineSurface::new(true, true,
+            KnotVector::from_multiplicities(3, &[0.001418343725273965962, 0.1860978202366415091], &[4, 4]),
+            KnotVector::from_multiplicities(2, &[0., 0.09081581946402909888], &[3, 3]),
+            rows.iter().map(|row| row.iter().map(|&(x, y, z, w)|
+                nalgebra_glm::DVec4::new(x * w, y * w, z * w, w)).collect()).collect());
+        let sampled = SampledSurface::new(surf);
+        let target = DVec3::new(-5.58262627960619, 3.666787751487475, 24.765511111689747);
+        // The pole is the nearest point: the target lies just beyond it.
+        let uv = sampled.uv_from_point(target).unwrap();
+        let pole = rows[3][0];
+        let pole = DVec3::new(pole.0, pole.1, pole.2);
+        close((sampled.surf.point(uv) - target).norm(), (pole - target).norm(), 1e-12);
+    }
+
     #[test]
     fn single_cell_projection_resolves_close_branches_in_parameter_space() {
         let surf = NDBSplineSurface::new(true, true,

step/src/parse.rs

@@ -119,7 +119,10 @@ impl<'a, T: Parse<'a>, const CAP: usize> Parse<'a> for ArrayVec<T, CAP> {
                 e => e?,
             };
             s = s_;
-            out.push(o);
+            // Some exporters write more items than the bound allows, such
+            // as four-component directions. Keep the first ones, as OCCT
+            // does, rather than reject the file.
+            let _ = out.try_push(o);
         }
         let (s, _) = char(')')(s)?;
         Ok((s, out))
@@ -367,6 +370,13 @@ mod tests {
         assert!(<&str>::parse("'unterminated").is_err());
     }
 
+    #[test]
+    fn bounded_lists_keep_their_first_items() {
+        // Some exporters write four-component directions.
+        let (rest, ratios) = ArrayVec::<f64, 3>::parse("(0.,0.,1.,0.),next").unwrap();
+        assert_eq!((rest, ratios.as_slice()), (",next", &[0., 0., 1.][..]));
+    }
+
     #[test]
     fn test_parse_entity_decl() {
         parse_entity_decl(b"#3=SHAPE_DEFINITION_REPRESENTATION(#4,#10);").unwrap();

triangulate/src/colored_mesh.rs

@@ -124,18 +124,20 @@ pub fn group_mesh_by_color(mesh: &Mesh) -> Result<TessellatedMesh, String> {
         bucket.indices.extend([a, b, c]);
     }
 
-    let mut submeshes = buckets.into_values().collect::<Vec<_>>();
-    // Largest groups first for deterministic ordering.
-    submeshes.sort_by(|a, b| a.indices.len().cmp(&b.indices.len()).reverse());
-    Ok(TessellatedMesh { submeshes })
+    let mut submeshes = buckets.into_iter().collect::<Vec<_>>();
+    // Largest groups first. Map order is random, so colour breaks ties.
+    submeshes.sort_by_key(|(key, submesh)| (std::cmp::Reverse(submesh.indices.len()), *key));
+    Ok(TessellatedMesh {
+        submeshes: submeshes.into_iter().map(|(_, submesh)| submesh).collect(),
+    })
 }
 
 // ---------------------------------------------------------------------------
 // Internals
 // ---------------------------------------------------------------------------
 
 /// Quantised RGBA colour key for bucketing.
-#[derive(Clone, Copy, Debug, PartialEq, Eq, Hash)]
+#[derive(Clone, Copy, Debug, PartialEq, Eq, PartialOrd, Ord, Hash)]
 struct ColorKey(u8, u8, u8, u8);
 
 impl ColorKey {
@@ -261,6 +263,29 @@ mod tests {
         assert_eq!(green.indices, [0, 1, 2]);
     }
 
+    #[test]
+    fn orders_equal_sized_color_groups_the_same_every_time() {
+        let mut mesh = Mesh {
+            verts: (0..6)
+                .map(|i| vertex(DVec3::new(i as f64, (i % 3) as f64, 0.0), DVec3::new(0.0, 0.0, 1.0)))
+                .collect(),
+            triangles: vec![
+                Triangle { verts: U32Vec3::new(0, 1, 2) },
+                Triangle { verts: U32Vec3::new(3, 4, 5) },
+            ],
+        };
+        for (i, v) in mesh.verts.iter_mut().enumerate() {
+            v.color = if i < 3 { DVec3::x() } else { DVec3::y() };
+        }
+        // Each call buckets colours in a new map, with new random keys.
+        let colors = || group_mesh_by_color(&mesh).unwrap().submeshes.iter().map(|s| s.color).collect::<Vec<_>>();
+        let first = colors();
+        assert_eq!(first.len(), 2);
+        for _ in 0..16 {
+            assert_eq!(colors(), first);
+        }
+    }
+
     #[test]
     fn keeps_coincident_vertices_separate_across_hard_edges() {
         let mesh = Mesh {

triangulate/src/surface.rs

@@ -328,6 +328,9 @@ impl Surface {
     /// Tessellate a compact surface with no physical trim. Identify seam
     /// vertices by wrapped grid indices instead of cutting a planar polygon.
     pub fn untrimmed_mesh(&self, color: DVec3, same_sense: bool, tolerance: f64) -> Option<Mesh> {
+        if let Self::Sphere { location, radius } = self {
+            return Self::sphere_mesh(*location, *radius, color, same_sense, tolerance);
+        }
         let Self::Torus {
             mat,
             major_radius,
@@ -376,6 +379,57 @@ impl Surface {
         Some(mesh)
     }
 
+    /// A whole sphere, as a face bounded only by a vertex loop at a pole is:
+    /// rings of latitude between two pole vertices.
+    fn sphere_mesh(
+        center: DVec3,
+        radius: f64,
+        color: DVec3,
+        same_sense: bool,
+        tolerance: f64,
+    ) -> Option<Mesh> {
+        if !(radius > 0.) {
+            return None;
+        }
+        // As for the torus: a chord of angle a sags r*a²/8 off the sphere;
+        // allocate half the error budget to each angular direction.
+        let nu = (2. * PI * (radius / (2. * tolerance)).sqrt()).ceil().max(6.) as u32;
+        let nv = (nu / 2).max(3);
+        nu.checked_mul(nv)?;
+        let sense = if same_sense { 1. } else { -1. };
+        // The south pole, nv - 1 rings of nu vertices each, the north pole.
+        let mut normals = vec![DVec3::new(0., 0., -1.)];
+        for v in 1..nv {
+            let (sv, cv) = (PI * v as f64 / nv as f64).sin_cos();
+            normals.extend((0..nu).map(|u| {
+                let (su, cu) = (2. * PI * u as f64 / nu as f64).sin_cos();
+                DVec3::new(sv * cu, sv * su, -cv)
+            }));
+        }
+        normals.push(DVec3::new(0., 0., 1.));
+        let mut mesh = Mesh::default();
+        mesh.verts = normals
+            .iter()
+            .map(|&n| Vertex { pos: center + n * radius, norm: n * sense, color })
+            .collect();
+        let (south, north) = (0, mesh.verts.len() as u32 - 1);
+        let ring = |v: u32, u: u32| 1 + (v - 1) * nu + u % nu;
+        // Counter-clockwise seen from outside.
+        let mut push = |a: u32, b: u32, c: u32| {
+            let verts = if same_sense { [a, b, c] } else { [a, c, b] };
+            mesh.triangles.push(Triangle { verts: verts.into() });
+        };
+        for u in 0..nu {
+            push(south, ring(1, u + 1), ring(1, u));
+            for v in 1..nv - 1 {
+                push(ring(v, u), ring(v, u + 1), ring(v + 1, u + 1));
+                push(ring(v, u), ring(v + 1, u + 1), ring(v + 1, u));
+            }
+            push(north, ring(nv - 1, u), ring(nv - 1, u + 1));
+        }
+        Some(mesh)
+    }
+
     pub fn new_plane(axis: DVec3) -> Result<Self, Error> {
         let dropped = axis.iamax();
         let scale = axis[dropped].abs();