diff --git a/triangulate/src/surface.rs b/triangulate/src/surface.rs index b75b5b9..8cdc2f6 100644 --- a/triangulate/src/surface.rs +++ b/triangulate/src/surface.rs @@ -194,6 +194,20 @@ pub struct PreparedSurface<'a> { uncertainty: f64, } +/// A prepared surface that owns its geometry, so a face can pause meshing +/// between stages without re-preparing its chart. +pub struct OwnedPreparedSurface { + surface: Surface, + chart: FaceChart, + uncertainty: f64, +} + +impl OwnedPreparedSurface { + pub fn prepared(&self) -> PreparedSurface<'_> { + PreparedSurface { surface: &self.surface, chart: self.chart.clone(), uncertainty: self.uncertainty } + } +} + impl Surface { pub fn new_nurbs(surf: SampledSurface<4>) -> Self { if let Some(normal) = surf.surf.bilinear_plane_normal() { @@ -617,6 +631,13 @@ impl Surface { uncertainty, }) } + + pub fn prepare_owned(self, verts: &[Vertex], boundary_edges: &[(usize, usize)], same_sense: bool, + uncertainty: f64, has_seam: bool, + ) -> Result { + let PreparedSurface { chart, .. } = self.prepare(verts, boundary_edges, same_sense, uncertainty, has_seam)?; + Ok(OwnedPreparedSurface { surface: self, chart, uncertainty }) + } } impl PreparedSurface<'_> { @@ -980,7 +1001,7 @@ impl PreparedSurface<'_> { /// Interior refinement cannot fix a constraint drawn through the wrong /// surface region (for example, a polar diameter instead of a rim arc). pub fn refine_boundary(&self, pts: &mut Vec<(f64, f64)>, edges: &mut Vec<(usize, usize)>, - verts: &mut Vec, tolerance: f64, + verts: &mut Vec, tolerance: f64, splits: &mut Vec<[DVec3; 2]>, ) -> Result<(), Error> { let mut scratch = ProjectionScratch::default(); let mut i = 0; @@ -990,7 +1011,6 @@ impl PreparedSurface<'_> { let pb = DVec2::new(pts[b].0, pts[b].1); let va = verts[a]; let vb = verts[b]; - let edge = vb.pos - va.pos; let Some(start) = self.raise(pa) else { i += 1; continue; }; let Some(end) = self.raise(pb) else { i += 1; continue; }; // Measure chart distortion separately from source curve/surface @@ -1005,7 +1025,9 @@ impl PreparedSurface<'_> { }) }); if !needs_split { i += 1; continue; } - let pos = va.pos + edge*0.5; + // Symmetric in the endpoints, so that faces sharing the segment + // can split it at exactly the same point. + let pos = (va.pos + vb.pos)*0.5; if pos == va.pos || pos == vb.pos || pts.len() >= 1_000_000 { return Err(Error::InvalidGeometry("boundary chart approximation did not converge")); } @@ -1018,6 +1040,7 @@ impl PreparedSurface<'_> { verts.push(Vertex { pos, norm: DVec3::zeros(), color: va.color }); edges[i] = (a, mid); edges.push((mid, b)); + splits.push([va.pos, vb.pos]); } Ok(()) } @@ -2376,7 +2399,7 @@ mod tests { let mut edges = vec![(0,1)]; let prepared = surface.prepare(&verts,&edges,true,0.,false).unwrap(); let mut pts = prepared.lower_verts(&verts).unwrap(); - prepared.refine_boundary(&mut pts,&mut edges,&mut verts,0.01).unwrap(); + prepared.refine_boundary(&mut pts,&mut edges,&mut verts,0.01,&mut Vec::new()).unwrap(); assert_eq!(verts.len(),2); assert_eq!([verts[0].pos,verts[1].pos],[a,b]); } diff --git a/triangulate/src/triangulate.rs b/triangulate/src/triangulate.rs index fcaadb2..bfb800c 100644 --- a/triangulate/src/triangulate.rs +++ b/triangulate/src/triangulate.rs @@ -13,7 +13,7 @@ use crate::{ mesh, mesh::{Mesh, Triangle}, stats::{FailureKind, Stats, TessellationFailure}, - surface::Surface, + surface::{OwnedPreparedSurface, Surface}, Error, }; use nurbs::{KnotVector, NURBSSurface, SampledCurve, SampledSurface}; @@ -866,18 +866,15 @@ fn shell( // 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( + meshed_face( s, std::slice::from_ref(face), - &mut scratch, styled_item_colors, default_color, uncertainty, tolerance, &edge_samples, - ); - (result, scratch) + ) }; #[cfg(feature = "rayon")] let mut meshed: Vec<_> = faces.par_iter().with_max_len(1).map(mesh_face).collect(); @@ -886,7 +883,7 @@ fn shell( // Exporters can split one face's loops between faces on copies of the // same cylinder, so each part alone extends to infinity. Mesh them as one. let split_key = |i: usize| { - if !matches!(meshed[i].0, Err(Error::InvalidGeometry(UNBOUNDED_FACE))) { + if !matches!(meshed[i], Err(Error::InvalidGeometry(UNBOUNDED_FACE))) { return None; } let Entity::AdvancedFace(face) = &s[faces[i]] else { return None; }; @@ -897,25 +894,32 @@ fn shell( let keys: Vec<_> = (0..faces.len()).map(split_key).collect(); for i in 0..faces.len() { let Some(key) = keys[i] else { continue; }; - let group: Vec<_> = (i..faces.len()).filter(|&j| keys[j] == Some(key) && meshed[j].0.is_err()).collect(); + let group: Vec<_> = (i..faces.len()).filter(|&j| keys[j] == Some(key) && meshed[j].is_err()).collect(); if group.len() < 2 { continue; } let parts: Vec<_> = group.iter().map(|&j| faces[j]).collect(); - let mut scratch = Mesh::default(); - if advanced_face(s, &parts, &mut scratch, styled_item_colors, default_color, uncertainty, tolerance, &edge_samples).is_ok() { - meshed[i] = (Ok(()), scratch); + if let Ok(face) = meshed_face(s, &parts, styled_item_colors, default_color, uncertainty, tolerance, &edge_samples) { + meshed[i] = Ok(face); for &j in &group[1..] { - meshed[j] = (Ok(()), Mesh::default()); + meshed[j] = Ok(MeshedFace::default()); } } } + // Each face splits the boundary segments its own chart needs split. Split + // them in every face that shares them too, so faces meet without + // T-junctions. + let splits: FastSet = meshed.iter().flatten().flat_map(|face| &face.splits).copied().collect(); + #[cfg(feature = "rayon")] + meshed.par_iter_mut().flatten().for_each(|face| face.share_splits(&splits)); + #[cfg(not(feature = "rayon"))] + meshed.iter_mut().flatten().for_each(|face| face.share_splits(&splits)); // Results are merged in face order, so output does not depend on // scheduling. - for (face, (result, mut scratch)) in faces.iter().zip(meshed) { + for (face, result) in faces.iter().zip(meshed) { stats.num_faces += 1; match result { - Ok(()) => mesh.append(&mut scratch), + Ok(mut meshed) => mesh.append(&mut meshed.mesh), Err(err) => { let surface_id = match &s[*face] { Entity::AdvancedFace(f) => Some(f.face_geometry.0), @@ -958,6 +962,34 @@ fn shell( stats.num_shells += 1; } +/// Hashes keys of indices and exact float bits, which need no protection +/// against collisions, with one multiply-rotate per word: much cheaper than +/// SipHash in the per-edge and per-triangle maps below. +#[derive(Default)] +struct FastHasher(u64); + +impl std::hash::Hasher for FastHasher { + fn write(&mut self, bytes: &[u8]) { + for word in bytes.chunks(8) { + let mut padded = [0; 8]; + padded[..word.len()].copy_from_slice(word); + self.write_u64(u64::from_le_bytes(padded)); + } + } + fn write_u64(&mut self, word: u64) { + self.0 = (self.0.rotate_left(5) ^ word).wrapping_mul(0x517c_c1b7_2722_0a95); + } + fn write_usize(&mut self, word: usize) { + self.write_u64(word as u64); + } + fn finish(&self) -> u64 { + self.0 + } +} + +type FastMap = HashMap>; +type FastSet = HashSet>; + const UNBOUNDED_FACE: &str = "face boundary does not enclose a finite region"; fn advanced_face( @@ -970,6 +1002,99 @@ fn advanced_face( tolerance: f64, edge_samples: &HashMap>, ) -> Result<(), Error> { + debug_assert!(mesh.verts.is_empty() && mesh.triangles.is_empty()); + *mesh = meshed_face(s, faces, styled_item_colors, default_color, uncertainty, tolerance, edge_samples)?.mesh; + Ok(()) +} + +/// A boundary segment by its endpoints' bits, in either direction. +type Segment = [[u64; 3]; 2]; + +fn segment(a: DVec3, b: DVec3) -> Segment { + let [a, b] = [a, b].map(|p| [p.x.to_bits(), p.y.to_bits(), p.z.to_bits()]); + if a <= b { [a, b] } else { [b, a] } +} + +/// A meshed face, ready to take in the boundary splits of the faces it +/// shares edges with. +#[derive(Default)] +struct MeshedFace { + mesh: Mesh, + /// The boundary segments this face's own chart needed split. + splits: Vec, + /// Faces meshed in a chart keep it to place and orient new vertices. + boundary: Option, +} + +struct ChartBoundary { + surface: OwnedPreparedSurface, + same_sense: bool, + /// Chart coordinates of the mesh vertices. + uv: Vec<(f64, f64)>, + /// Mesh edges that border one triangle. + edges: Vec<(usize, usize)>, +} + +impl MeshedFace { + /// Split each boundary edge that is in `splits` at its midpoint, the + /// point the face that needed the split used, and fan its triangle in + /// two. The midpoint is on the edge, so the surface does not change. + /// This releases the chart, here in parallel rather than in the merge. + fn share_splits(&mut self, splits: &FastSet) { + let Some(mut boundary) = self.boundary.take() else { return; }; + let boundary = &mut boundary; + let Mesh { verts, triangles } = &mut self.mesh; + let split = |verts: &[mesh::Vertex], a: usize, b: usize| splits.contains(&segment(verts[a].pos, verts[b].pos)); + let mut pending: Vec<_> = boundary.edges.iter().copied().filter(|&(a, b)| split(verts, a, b)).collect(); + if pending.is_empty() { + return; + } + let edge_key = |a: usize, b: usize| (a.min(b), a.max(b)); + let mut owner: FastMap<_, _> = pending.iter().map(|&(a, b)| (edge_key(a, b), 0)).collect(); + for (t, triangle) in triangles.iter().enumerate() { + let v = [triangle.verts.x, triangle.verts.y, triangle.verts.z].map(|i| i as usize); + for k in 0..3 { + if let Some(owner) = owner.get_mut(&edge_key(v[k], v[(k + 1) % 3])) { + *owner = t; + } + } + } + let prepared = boundary.surface.prepared(); + let sign = if boundary.same_sense { 1. } else { -1. }; + while let Some((a, b)) = pending.pop() { + let t = owner[&edge_key(a, b)]; + let v = [triangles[t].verts.x, triangles[t].verts.y, triangles[t].verts.z].map(|i| i as usize); + let k = (0..3).find(|&k| edge_key(v[k], v[(k + 1) % 3]) == edge_key(a, b)).expect("owner has the edge"); + let (x, y, z) = (v[k], v[(k + 1) % 3], v[(k + 2) % 3]); + let m = verts.len(); + let pos = (verts[a].pos + verts[b].pos) * 0.5; + let (ua, ub) = (boundary.uv[a], boundary.uv[b]); + let uv = ((ua.0 + ub.0) * 0.5, (ua.1 + ub.1) * 0.5); + let norm = prepared.normal(pos, DVec2::new(uv.0, uv.1)) * sign; + verts.push(mesh::Vertex { pos, norm, color: verts[a].color }); + boundary.uv.push(uv); + triangles[t].verts = U32Vec3::new(x as u32, m as u32, z as u32); + owner.insert(edge_key(x, m), t); + // The new triangle takes over the edge y z. + for edge in [edge_key(m, y), edge_key(y, z)] { + owner.insert(edge, triangles.len()); + } + triangles.push(Triangle { verts: U32Vec3::new(m as u32, y as u32, z as u32) }); + pending.extend([(x, m), (m, y)].iter().copied().filter(|&(a, b)| split(verts, a, b))); + } + } +} + +fn meshed_face( + s: &StepFile, + faces: &[Face], + styled_item_colors: &HashMap, + default_color: DVec3, + uncertainty: f64, + tolerance: f64, + edge_samples: &HashMap>, +) -> Result { + let mut mesh = Mesh::default(); // Closed shells may legally reference either ADVANCED_FACE or the more // general FACE_SURFACE; OCCT/KiCad translate both through FaceSurface. // Faces meshed together use the bounds of all and the first one's surface. @@ -990,11 +1115,10 @@ fn advanced_face( .unwrap_or(default_color); info!("triangulating face {} (geometry {})", f.0, face_geometry.0); - debug_assert!(mesh.verts.is_empty() && mesh.triangles.is_empty()); - if let Some(&b) = bounds.first() { if let Entity::PolyLoop(_) = &s[bound_loop(s, b)?.0] { - return poly_loop_face(s, face_geometry, same_sense, &bounds, face_color, mesh); + poly_loop_face(s, face_geometry, same_sense, &bounds, face_color, &mut mesh)?; + return Ok(MeshedFace { mesh, ..Default::default() }); } } @@ -1067,7 +1191,7 @@ fn advanced_face( } if !bounds.is_empty() && num_pts == 0 { - return Ok(()); + return Ok(MeshedFace { mesh, ..Default::default() }); } // Swept surfaces use the actual trims to choose a finite NURBS domain. @@ -1085,14 +1209,14 @@ fn advanced_face( .all(|&(forward, reverse)| forward == reverse) { if let Some(full) = surf.untrimmed_mesh(face_color, same_sense, tolerance) { - *mesh = full; - return Ok(()); + return Ok(MeshedFace { mesh: full, ..Default::default() }); } } // Add curvature samples before constraint insertion. The CDT subdivides // constraints at existing vertices, including samples exactly on an edge. - let prepared = surf.prepare(&mesh.verts, &edges, same_sense, uncertainty, has_seam)?; + let surface = surf.prepare_owned(&mesh.verts, &edges, same_sense, uncertainty, has_seam)?; + let prepared = surface.prepared(); let mut pts = crate::timing::time("face:lower_verts", || prepared.lower_verts(&mesh.verts))?; prepared.continue_trims(&mut pts,&mesh.verts,&edges,tolerance); if !prepared.cut_periodic(&mut pts, &mut edges, &mut mesh.verts, tolerance)? { @@ -1105,7 +1229,8 @@ fn advanced_face( if had_boundary && edges.is_empty() { return Err(Error::InvalidGeometry("face boundary cancels completely")); } - prepared.refine_boundary(&mut pts, &mut edges, &mut mesh.verts, tolerance)?; + let mut splits = Vec::new(); + prepared.refine_boundary(&mut pts, &mut edges, &mut mesh.verts, tolerance, &mut splits)?; let mut constraints: Vec<_> = edges.iter().map(|&(a, b)| (a, b, true)).collect(); let bonus_points = pts.len(); crate::timing::time("face:add_steiner_points", || { @@ -1146,14 +1271,14 @@ fn advanced_face( })) }); let t = result.map_err(|_| Error::TriangulationPanic)??; - crate::timing::time("face:refine_and_normals", || { + let boundary = crate::timing::time("face:refine_and_normals", || { let mut triangles: Vec<_> = t.triangles().map(|(a,b,c)| [a,b,c]).collect(); let edge_key = |a: usize,b: usize| (a.min(b),a.max(b)); - let mut uses = HashMap::new(); + let mut uses = FastMap::default(); for &[a,b,c] in &triangles { for (a,b) in [(a,b),(b,c),(c,a)] { *uses.entry(edge_key(a,b)).or_insert(0) += 1; } } - let boundary: HashSet<_> = uses.into_iter().filter_map(|(edge,count)| (count == 1).then_some(edge)).collect(); + let boundary: FastSet<_> = uses.into_iter().filter_map(|(edge,count)| (count == 1).then_some(edge)).collect(); // Refine in spatial geometry, without rerunning a chart-metric CDT that // repeatedly reconnects long spatial edges on highly stretched charts. // Shared edge midpoints and a subdivision table keep the mesh conforming. @@ -1165,7 +1290,7 @@ fn advanced_face( ]; let mut surface_positions: Vec<_> = pts.iter().zip(&mesh.verts) .map(|(&(u,v), vertex)| prepared.raise(DVec2::new(u,v)).unwrap_or(vertex.pos)).collect(); - let children = |[a,b,c]: [usize;3], splits: &HashMap<(usize,usize),usize>| { + let children = |[a,b,c]: [usize;3], splits: &FastMap<(usize,usize),usize>| { let mut nodes = [a,b,c,0,0,0]; let mut mask = 0; for (i,(a,b)) in [(a,b),(b,c),(c,a)].iter().copied().enumerate() { @@ -1176,7 +1301,7 @@ fn advanced_face( // Keep red leaves, not their temporary green completion. Refining a // failing green child promotes its owner; repeated green-only splitting // can otherwise preserve diameter while its altitude tends to zero. - let mut splits = HashMap::new(); + let mut splits = FastMap::default(); loop { let mut marked = Vec::with_capacity(triangles.len()); for &[a,b,c] in &triangles { @@ -1264,7 +1389,15 @@ fn advanced_face( v.norm = prepared.normal(v.pos, DVec2::new(u, w)) * if same_sense { 1. } else { -1. }; } info!("face {} done", face_id); - Ok(()) + Ok(boundary) + })?; + // Sorted, so that splits apply in the same order on every run. + let mut edges: Vec<_> = boundary.into_iter().collect(); + edges.sort_unstable(); + Ok(MeshedFace { + splits: splits.into_iter().map(|[a, b]| segment(a, b)).collect(), + boundary: Some(ChartBoundary { surface, same_sense, uv: pts, edges }), + mesh, }) } @@ -2098,7 +2231,7 @@ fn resolve_crossing_edges( let bits = |v: f64| if v == 0.0 { 0 } else { v.to_bits() }; (bits(x), bits(y)) }; - let mut vertices = HashMap::new(); + let mut vertices = FastMap::default(); for (i, &p) in pts.iter().enumerate() { vertices.entry(key(p)).or_insert(i); } diff --git a/triangulate/tests/examples.rs b/triangulate/tests/examples.rs index 6affdd2..f3f5870 100644 --- a/triangulate/tests/examples.rs +++ b/triangulate/tests/examples.rs @@ -16,5 +16,26 @@ fn checked_in_models_tessellate_without_errors() { assert!(stats.num_faces > 0, "{} contained no faces", name); assert_eq!(stats.num_errors(), 0, "{} had tessellation errors", name); assert_eq!(stats.num_panics(), 0, "{} had tessellation panics", name); + + // The models are closed solids. Faces that share an edge share its + // points, so every edge borders exactly two triangles. + let mut ids = std::collections::HashMap::new(); + let id: Vec = mesh + .verts + .iter() + .map(|v| { + let next = ids.len(); + *ids.entry([v.pos.x, v.pos.y, v.pos.z].map(f64::to_bits)).or_insert(next) + }) + .collect(); + let mut uses = std::collections::HashMap::new(); + for t in &mesh.triangles { + let [a, b, c] = [t.verts.x, t.verts.y, t.verts.z].map(|i| id[i as usize]); + for (p, q) in [(a, b), (b, c), (c, a)] { + *uses.entry((p.min(q), p.max(q))).or_insert(0) += 1; + } + } + let open = uses.values().filter(|&&n| n != 2).count(); + assert_eq!(open, 0, "{} has {} edges without exactly two triangles", name, open); } }