diff --git a/nurbs/src/sampled_surface.rs b/nurbs/src/sampled_surface.rs index f3e02e4..c35da05 100644 --- a/nurbs/src/sampled_surface.rs +++ b/nurbs/src/sampled_surface.rs @@ -26,6 +26,13 @@ pub struct ProjectionScratch { } 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 { @@ -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, diff --git a/step/src/parse.rs b/step/src/parse.rs index 325bc4c..912ac29 100644 --- a/step/src/parse.rs +++ b/step/src/parse.rs @@ -119,7 +119,10 @@ impl<'a, T: Parse<'a>, const CAP: usize> Parse<'a> for ArrayVec { 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::::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(); diff --git a/triangulate/src/colored_mesh.rs b/triangulate/src/colored_mesh.rs index d1396b6..6e4a0f4 100644 --- a/triangulate/src/colored_mesh.rs +++ b/triangulate/src/colored_mesh.rs @@ -124,10 +124,12 @@ pub fn group_mesh_by_color(mesh: &Mesh) -> Result { bucket.indices.extend([a, b, c]); } - let mut submeshes = buckets.into_values().collect::>(); - // 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::>(); + // 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(), + }) } // --------------------------------------------------------------------------- @@ -135,7 +137,7 @@ pub fn group_mesh_by_color(mesh: &Mesh) -> Result { // --------------------------------------------------------------------------- /// 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::>(); + 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 { diff --git a/triangulate/src/surface.rs b/triangulate/src/surface.rs index 68e4fac..ac215d6 100644 --- a/triangulate/src/surface.rs +++ b/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 { + 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 { + 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 { let dropped = axis.iamax(); let scale = axis[dropped].abs(); diff --git a/triangulate/src/triangulate.rs b/triangulate/src/triangulate.rs index ea41292..fd504cf 100644 --- a/triangulate/src/triangulate.rs +++ b/triangulate/src/triangulate.rs @@ -62,7 +62,7 @@ fn collect_shape_instances<'a>( s: &'a StepFile, rep_instances: &HashMap, Vec>, shape_rep_relationship: &HashMap, Vec>>, -) -> HashMap, ShapeGeometry> { +) -> Vec<(RepresentationItem<'a>, ShapeGeometry)> { let mut todo: Vec<_> = rep_instances .iter() .flat_map(|(rep, mats)| mats.iter().copied().map(move |mat| (*rep, mat))) @@ -126,6 +126,15 @@ fn collect_shape_instances<'a>( }); } + // The maps above are visited in a random order. Mesh shapes, and their + // copies, in a fixed one so the output does not change from run to run. + let mut to_mesh: Vec<_> = to_mesh.into_iter().collect(); + to_mesh.sort_unstable_by_key(|(id, _)| id.0); + for (_, shape) in &mut to_mesh { + shape.instances.sort_unstable_by(|a, b| { + a.iter().map(|v| v.to_bits()).cmp(b.iter().map(|v| v.to_bits())) + }); + } to_mesh } @@ -839,13 +848,12 @@ fn shell( }; let v_start = mesh.verts.len(); let t_start = mesh.triangles.len(); - let mut scratch = Mesh::default(); - let mut edge_samples = HashMap::new(); - for face in faces { - scratch.verts.clear(); - scratch.triangles.clear(); - stats.num_faces += 1; - match advanced_face( + // Sample every edge once up front, so faces sharing an edge share its + // points exactly and the faces can then be meshed independently. + let edge_samples = sample_edges(s, faces, tolerance); + let mesh_face = |face: &Face| { + let mut scratch = Mesh::default(); + let result = advanced_face( s, *face, &mut scratch, @@ -853,8 +861,19 @@ fn shell( default_color, uncertainty, tolerance, - &mut edge_samples, - ) { + &edge_samples, + ); + (result, scratch) + }; + #[cfg(feature = "rayon")] + let meshed: Vec<_> = faces.par_iter().with_max_len(1).map(mesh_face).collect(); + #[cfg(not(feature = "rayon"))] + let meshed: Vec<_> = faces.iter().map(mesh_face).collect(); + // Results are merged in face order, so output does not depend on + // scheduling. + for (face, (result, mut scratch)) in faces.iter().zip(meshed) { + stats.num_faces += 1; + match result { Ok(()) => mesh.append(&mut scratch), Err(err) => { let surface_id = match &s[*face] { @@ -906,7 +925,7 @@ fn advanced_face( default_color: DVec3, uncertainty: f64, tolerance: f64, - edge_samples: &mut HashMap>, + edge_samples: &HashMap>, ) -> Result<(), Error> { // Closed shells may legally reference either ADVANCED_FACE or the more // general FACE_SURFACE; OCCT/KiCad translate both through FaceSurface. @@ -1605,7 +1624,7 @@ fn face_bound( b: FaceBound, edge_uses: &mut HashMap, tolerance: f64, - edge_samples: &mut HashMap>, + edge_samples: &HashMap>, ) -> Result<(Vec, usize), Error> { let (bound, orientation) = match &s[b] { Entity::FaceBound(b) => (b.bound, b.orientation), @@ -1636,8 +1655,46 @@ fn face_bound( } } +/// The 3D samples of every edge bounding `faces`, by edge curve id. A shell +/// has one STEP file and tolerance, so an edge's samples are the same for +/// every face that uses it; traversal and surface charts remain +/// face-specific. Edges that cannot be sampled are left out. +fn sample_edges(s: &StepFile, faces: &[Face], tolerance: f64) -> HashMap> { + let mut edges: Vec = faces + .iter() + .filter_map(|f| match &s[*f] { + Entity::AdvancedFace(face) => Some(&face.bounds[..]), + Entity::FaceSurface(face) => Some(&face.bounds[..]), + _ => None, + }) + .flatten() + .filter_map(|b| match &s[*b] { + Entity::FaceBound(b) => Some(b.bound), + Entity::FaceOuterBound(b) => Some(b.bound), + _ => None, + }) + .filter_map(|bound| match &s[bound] { + Entity::EdgeLoop(e) => Some(&e.edge_list[..]), + _ => None, + }) + .flatten() + .filter_map(|e| s.entity(*e).map(|edge| edge.edge_element.0)) + .collect(); + edges.sort_unstable(); + edges.dedup(); + let sample = |&id: &usize| { + edge_curve(s, EdgeCurve::new(id), true, tolerance) + .ok() + .map(|points| (id, points)) + }; + #[cfg(feature = "rayon")] + return edges.par_iter().filter_map(sample).collect(); + #[cfg(not(feature = "rayon"))] + return edges.iter().filter_map(sample).collect(); +} + fn edge_loop(s: &StepFile, edge_list: &[OrientedEdge], tolerance: f64, - edge_samples: &mut HashMap>, + edge_samples: &HashMap>, ) -> Result, Error> { let mut out = Vec::new(); for (i, e) in edge_list.iter().enumerate() { @@ -1649,12 +1706,15 @@ fn edge_loop(s: &StepFile, edge_list: &[OrientedEdge], tolerance: f64, let edge = s .entity(*e) .ok_or(Error::InvalidStepEntity("OrientedEdge"))?; - // A shell has one STEP file and tolerance. Cache only edge-owned 3D - // samples; traversal and surface charts remain face-specific. - let points = match edge_samples.entry(edge.edge_element.0) { - std::collections::hash_map::Entry::Occupied(entry) => entry.into_mut(), - std::collections::hash_map::Entry::Vacant(entry) => - entry.insert(edge_curve(s, edge.edge_element.cast(), true, tolerance)?), + // Edges that could not be sampled up front fail here, for each face + // that uses them. + let sampled; + let points = match edge_samples.get(&edge.edge_element.0) { + Some(points) => points, + None => { + sampled = edge_curve(s, edge.edge_element.cast(), true, tolerance)?; + &sampled + } }; if edge.orientation { out.extend_from_slice(points); @@ -2089,7 +2149,7 @@ mod tests { let step = StepFile::parse(&flat).unwrap(); let mut mesh = Mesh::default(); let color = DVec3::new(0.2, 0.4, 0.6); - advanced_face(&step, Id::new(72), &mut mesh, &HashMap::new(), color, 0., 0.01, &mut HashMap::new()).unwrap(); + advanced_face(&step, Id::new(72), &mut mesh, &HashMap::new(), color, 0., 0.01, &HashMap::new()).unwrap(); assert!(mesh.verts.iter().any(|v| v.pos == DVec3::new(1., 1., 0.))); assert!(mesh .verts @@ -2397,6 +2457,54 @@ mod tests { assert!(edges.is_empty()); } + #[test] + fn whole_sphere_bounded_by_a_pole_is_meshed() { + // A BGA solder ball: the whole sphere, bounded only by a vertex + // loop at the pole where it meets the package. + for (sense, outward) in [(".T.", 1.), (".F.", -1.)] { + let text = format!("ISO-10303-21;HEADER;ENDSEC;DATA; + #1=CARTESIAN_POINT('',(1.,2.,0.115)); + #2=DIRECTION('',(0.,0.,1.)); + #3=DIRECTION('',(1.,0.,0.)); + #4=AXIS2_PLACEMENT_3D('',#1,#2,#3); + #5=SPHERICAL_SURFACE('',#4,0.115); + #6=CARTESIAN_POINT('',(1.,2.,0.)); + #7=VERTEX_POINT('',#6); + #8=VERTEX_LOOP('',#7); + #9=FACE_BOUND('',#8,.T.); + #10=ADVANCED_FACE('',(#9),#5,{sense}); + ENDSEC;END-ISO-10303-21;"); + let flat = StepFile::strip_flatten(text.as_bytes()).unwrap(); + let step = StepFile::parse(&flat).unwrap(); + let tolerance = 0.01; + let mut mesh = Mesh::default(); + advanced_face(&step, Id::new(10), &mut mesh, &HashMap::new(), DVec3::zeros(), 0., tolerance, &HashMap::new()).unwrap(); + let center = DVec3::new(1., 2., 0.115); + for v in &mesh.verts { + assert!(((v.pos - center).norm() - 0.115).abs() < 1e-12); + assert!(((v.pos - center) / 0.115 * outward - v.norm).norm() < 1e-9); + } + // Closed: every edge is used once in each direction. + let mut edges = HashMap::<_, i32>::new(); + let mut area = 0.; + for t in &mesh.triangles { + let [a, b, c] = [t.verts.x, t.verts.y, t.verts.z].map(|i| mesh.verts[i as usize].pos); + let p = (a + b + c) / 3.; + assert!(0.115 - (p - center).norm() < tolerance); + let n = (b - a).cross(&(c - a)); + assert!(n.dot(&(p - center)) * outward > 0.); + area += n.norm() / 2.; + for (u, v) in [(t.verts.x, t.verts.y), (t.verts.y, t.verts.z), (t.verts.z, t.verts.x)] { + *edges.entry((u.min(v), u.max(v))).or_default() += if u < v { 1 } else { -1 }; + } + } + assert_eq!(edges.len() * 2, mesh.triangles.len() * 3); + assert!(edges.values().all(|&d| d == 0)); + let sphere = 4. * std::f64::consts::PI * 0.115f64.powi(2); + assert!((area - sphere).abs() < 0.05 * sphere, "{area} vs {sphere}"); + } + } + #[test] fn hemisphere_refines_interior_curvature_to_the_physical_budget() { let text = b"ISO-10303-21;HEADER;ENDSEC;DATA; @@ -2419,7 +2527,7 @@ mod tests { let mut previous = 0; for tolerance in [0.04,0.01] { let mut mesh = Mesh::default(); - advanced_face(&step, Id::new(13), &mut mesh, &HashMap::new(), DVec3::zeros(), 0., tolerance, &mut HashMap::new()).unwrap(); + advanced_face(&step, Id::new(13), &mut mesh, &HashMap::new(), DVec3::zeros(), 0., tolerance, &HashMap::new()).unwrap(); assert!(mesh.triangles.len() > previous); previous = mesh.triangles.len(); for t in &mesh.triangles { @@ -2488,12 +2596,14 @@ mod tests { let mut backward = edge_curve(&step, Id::new(id), false, 0.01).unwrap(); backward.reverse(); assert_eq!(forward, backward); - let mut samples = HashMap::new(); - let mut cached_backward = edge_loop(&step, &[Id::new(backward_id)], 0.01, &mut samples).unwrap(); - cached_backward.reverse(); - assert_eq!(forward, cached_backward); - assert_eq!(forward, edge_loop(&step, &[Id::new(forward_id)], 0.01, &mut samples).unwrap()); - assert_eq!(samples.len(), 1); + // Shared samples and sampling on demand agree. + let shared = HashMap::from([(id, forward.clone())]); + for samples in [&shared, &HashMap::new()] { + let mut cached_backward = edge_loop(&step, &[Id::new(backward_id)], 0.01, samples).unwrap(); + cached_backward.reverse(); + assert_eq!(forward, cached_backward); + assert_eq!(forward, edge_loop(&step, &[Id::new(forward_id)], 0.01, samples).unwrap()); + } } } } @@ -2528,8 +2638,44 @@ mod tests { .cloned() .collect(); let shapes = collect_shape_instances(&step, &roots, &HashMap::new()); - assert_eq!(shapes[&Id::new(10)].uncertainty, 2e-7); - assert!((shapes[&Id::new(11)].uncertainty - 2e-10).abs() < 1e-25); + // Shapes come in entity order. + assert_eq!(shapes.iter().map(|(id, _)| id.0).collect::>(), [10, 11]); + assert_eq!(shapes[0].1.uncertainty, 2e-7); + assert!((shapes[1].1.uncertainty - 2e-10).abs() < 1e-25); + } + + #[test] + fn shapes_and_their_copies_come_in_a_fixed_order() { + // Shape #10 is used by both representations, so it collects copies + // from each in whatever order the representation map is visited. + let text = b"ISO-10303-21;HEADER;ENDSEC;DATA; + #4=REPRESENTATION_CONTEXT('',''); + #7=SHAPE_REPRESENTATION('',(#10),#4); + #8=SHAPE_REPRESENTATION('',(#11,#10),#4); + #9=CLOSED_SHELL('',()); + #10=MANIFOLD_SOLID_BREP('',#9); + #11=MANIFOLD_SOLID_BREP('',#9); + ENDSEC;END-ISO-10303-21;"; + let flat = StepFile::strip_flatten(text).unwrap(); + let step = StepFile::parse(&flat).unwrap(); + let shift = |x: f64| glm::translation(&DVec3::new(x, 0., 0.)); + let order = || { + // Each new map hashes with new random keys. + let roots = HashMap::from([ + (Id::new(7), vec![shift(1.), shift(2.)]), + (Id::new(8), vec![shift(3.)]), + ]); + collect_shape_instances(&step, &roots, &HashMap::new()) + .iter() + .map(|(id, shape)| (id.0, shape.instances.iter().map(|m| m[(0, 3)]).collect::>())) + .collect::>() + }; + let first = order(); + assert_eq!(first.len(), 2); + assert_eq!(first.iter().map(|(_, copies)| copies.len()).sum::(), 4); + for _ in 0..16 { + assert_eq!(order(), first); + } } #[test]