diff --git a/examples/adjacency_demo.rs b/examples/adjacency_demo.rs index ab5f82c..3ef3587 100644 --- a/examples/adjacency_demo.rs +++ b/examples/adjacency_demo.rs @@ -145,10 +145,24 @@ fn main() { // Demonstrate adaptive refinement println!("\n7. Adaptive mesh refinement:"); let refined = tessellated.adaptive_refine(0.5, 2.0, 15.0); - let (refined_vertex_map, refined_adjacency_map) = refined.build_connectivity(); println!(" Original triangles: {}", tessellated.polygons.len()); println!(" After refinement: {}", refined.polygons.len()); + // Re-analyze connectivity to show changes + let (refined_vertex_map, refined_adjacency_map) = refined.build_connectivity(); + let refined_avg_valence = if !refined_adjacency_map.is_empty() { + refined_adjacency_map.values().map(|v| v.len()).sum::() as Real + / refined_adjacency_map.len() as Real + } else { + 0.0 + }; + + println!( + " Refined unique vertices: {}", + refined_vertex_map.vertex_count() + ); + println!(" Refined average valence: {:.2}", refined_avg_valence); + if refined.polygons.len() > tessellated.polygons.len() { println!(" ✓ Mesh was refined based on quality criteria"); } else { diff --git a/readme.md b/readme.md index 29dc742..c8126b6 100644 --- a/readme.md +++ b/readme.md @@ -501,6 +501,10 @@ and memory usage: - allocations should be kept to a minimum. Memory should be read-only if possible, clone if necessary, and offer the choice of transmut in place or create new copy via appropriate functions +- For performance-critical operations like BSP and SDF, we use a modular backend +system with dependency inversion. This allows switching between serial and parallel +implementations at compile time via the "parallel" feature flag, ensuring both +flexibility and high performance. ## Todo - when triangulating, detect T junctions with other polygons with shared edges, diff --git a/src/io/stl.rs b/src/io/stl.rs index 9459fb7..807ba5d 100644 --- a/src/io/stl.rs +++ b/src/io/stl.rs @@ -406,7 +406,6 @@ impl Sketch { } } - // // (C) Encode into a binary STL buffer // let mut cursor = Cursor::new(Vec::new()); diff --git a/src/lib.rs b/src/lib.rs index e528624..e0d9c83 100644 --- a/src/lib.rs +++ b/src/lib.rs @@ -3,7 +3,6 @@ //! //! ![Example CSG output][Example CSG output] #![cfg_attr(doc, doc = doc_image_embed::embed_image!("Example CSG output", "docs/csg.png"))] -//! //! # Features //! #### Default //! - **f64**: use f64 as Real diff --git a/src/main.rs b/src/main.rs index 12e8d4b..fceb367 100644 --- a/src/main.rs +++ b/src/main.rs @@ -252,7 +252,7 @@ fn main() { let ray_origin = Point3::new(0.0, 0.0, -5.0); let ray_dir = Vector3::new(0.0, 0.0, 1.0); // pointing along +Z let hits = cube.ray_intersections(&ray_origin, &ray_dir); - println!("Ray hits on the cube: {:?}", hits); + println!("Ray hits on the cube: {hits:?}"); } // 12) Polyhedron example (simple tetrahedron): @@ -297,9 +297,9 @@ fn main() { // 14) Mass properties (just printing them) let (mass, com, principal_frame) = cube.mass_properties(1.0); - println!("Cube mass = {}", mass); - println!("Cube center of mass = {:?}", com); - println!("Cube principal inertia local frame = {:?}", principal_frame); + println!("Cube mass = {mass}"); + println!("Cube center of mass = {com:?}"); + println!("Cube principal inertia local frame = {principal_frame:?}"); // 1) Create a cube from (-1,-1,-1) to (+1,+1,+1) // (By default, CSG::cube(None) is from -1..+1 if the "radius" is [1,1,1].) @@ -349,11 +349,11 @@ fn main() { ); } - //let poor_geometry_shape = moved_cube.difference(&sphere); + // let poor_geometry_shape = moved_cube.difference(&sphere); //#[cfg(feature = "earclip-io")] - //let retriangulated_shape = poor_geometry_shape.triangulate_earclip(); + // let retriangulated_shape = poor_geometry_shape.triangulate_earclip(); //#[cfg(all(feature = "earclip-io", feature = "stl-io"))] - //let _ = fs::write("stl/retriangulated.stl", retriangulated_shape.to_stl_binary("retriangulated").unwrap()); + // let _ = fs::write("stl/retriangulated.stl", retriangulated_shape.to_stl_binary("retriangulated").unwrap()); let sphere_test = Mesh::sphere(1.0, 16, 8, None); let cube_test = Mesh::cube(1.0, None); @@ -724,8 +724,8 @@ fn main() { let _ = fs::write("stl/octahedron.stl", oct.to_stl_ascii("octahedron")); } - //let dodec = CSG::dodecahedron(15.0, None); - //let _ = fs::write("stl/dodecahedron.stl", dodec.to_stl_ascii("")); + // let dodec = CSG::dodecahedron(15.0, None); + // let _ = fs::write("stl/dodecahedron.stl", dodec.to_stl_ascii("")); #[cfg(feature = "stl-io")] { @@ -1041,19 +1041,17 @@ fn main() { ); } - /* - let helical = CSG::helical_involute_gear( - 2.0, // module - 20, // z - 20.0, // pressure angle - 0.05, 0.02, 14, - 25.0, // face-width - 15.0, // helix angle β [deg] - 40, // axial slices (resolution of the twist) - None, - ); - let _ = fs::write("stl/helical.stl", helical.to_stl_ascii("helical")); - */ + // let helical = CSG::helical_involute_gear( + // 2.0, // module + // 20, // z + // 20.0, // pressure angle + // 0.05, 0.02, 14, + // 25.0, // face-width + // 15.0, // helix angle β [deg] + // 40, // axial slices (resolution of the twist) + // None, + // ); + // let _ = fs::write("stl/helical.stl", helical.to_stl_ascii("helical")); // Bézier curve demo #[cfg(feature = "stl-io")] @@ -1081,8 +1079,10 @@ fn main() { let bspline_ctrl = &[[0.0, 0.0], [1.0, 2.5], [3.0, 3.0], [5.0, 0.0], [6.0, -1.5]]; let bspline_2d = Sketch::bspline( bspline_ctrl, - /* degree p = */ 3, - /* seg/span */ 32, + // degree p = + 3, + // seg/span + 32, None, ); let _ = fs::write("stl/bspline_2d.stl", bspline_2d.to_stl_ascii("bspline_2d")); @@ -1092,8 +1092,8 @@ fn main() { println!("{:#?}", bezier_3d.to_bevy_mesh()); // a quick thickening just like the Bézier - //let bspline_3d = bspline_2d.extrude(0.25); - //let _ = fs::write( + // let bspline_3d = bspline_2d.extrude(0.25); + // let _ = fs::write( // "stl/bspline_extruded.stl", // bspline_3d.to_stl_ascii("bspline_extruded"), //); diff --git a/src/mesh/algorithm/convex_hull/mod.rs b/src/mesh/algorithm/convex_hull/mod.rs new file mode 100644 index 0000000..9ab8954 --- /dev/null +++ b/src/mesh/algorithm/convex_hull/mod.rs @@ -0,0 +1,16 @@ +//! Mesh convex hull algorithms. + +pub mod serial; +pub mod traits; + +#[cfg(feature = "parallel")] +pub mod parallel; + +// Re-export core types +pub use traits::ConvexHullOps; + +#[cfg(not(feature = "parallel"))] +pub use serial::SerialConvexHullOps; + +#[cfg(feature = "parallel")] +pub use parallel::ParallelConvexHullOps; diff --git a/src/mesh/algorithm/convex_hull/parallel.rs b/src/mesh/algorithm/convex_hull/parallel.rs new file mode 100644 index 0000000..521632e --- /dev/null +++ b/src/mesh/algorithm/convex_hull/parallel.rs @@ -0,0 +1,129 @@ +//! Parallel implementations of convex hull operations. + +use super::traits::ConvexHullOps; +use crate::float_types::Real; +use crate::mesh::{Mesh, polygon::Polygon, vertex::Vertex}; +use crate::traits::CSG; +use chull::ConvexHullWrapper; +use nalgebra::{Point3, Vector3}; +use rayon::prelude::*; +use std::fmt::Debug; + +/// Parallel implementation of `ConvexHullOps`. +pub struct ParallelConvexHullOps; + +impl ParallelConvexHullOps { + pub fn new() -> Self { + Self + } +} + +impl ConvexHullOps for ParallelConvexHullOps { + fn convex_hull(&self, mesh: &Mesh) -> Mesh { + // Gather all (x, y, z) coordinates from the polygons in parallel + let points: Vec> = mesh + .polygons + .par_iter() + .flat_map(|poly| poly.vertices.par_iter().map(|v| v.pos)) + .collect(); + + let points_for_hull: Vec> = + points.par_iter().map(|p| vec![p.x, p.y, p.z]).collect(); + + // Attempt to compute the convex hull using the robust wrapper + let hull = match ConvexHullWrapper::try_new(&points_for_hull, None) { + Ok(h) => h, + Err(_) => { + // Fallback to an empty CSG if hull generation fails + return Mesh::new(); + }, + }; + + let (verts, indices) = hull.vertices_indices(); + + // Reconstruct polygons as triangles + let mut polygons = Vec::new(); + for tri in indices.chunks(3) { + let v0 = &verts[tri[0]]; + let v1 = &verts[tri[1]]; + let v2 = &verts[tri[2]]; + let vv0 = Vertex::new(Point3::new(v0[0], v0[1], v0[2]), Vector3::zeros()); + let vv1 = Vertex::new(Point3::new(v1[0], v1[1], v1[2]), Vector3::zeros()); + let vv2 = Vertex::new(Point3::new(v2[0], v2[1], v2[2]), Vector3::zeros()); + polygons.push(Polygon::new(vec![vv0, vv1, vv2], None)); + } + + Mesh::from_polygons(&polygons, mesh.metadata.clone()) + } + + fn minkowski_sum(&self, mesh: &Mesh, other: &Mesh) -> Mesh { + // Collect all vertices (x, y, z) from self + let verts_a: Vec> = mesh + .polygons + .par_iter() + .flat_map(|poly| poly.vertices.par_iter().map(|v| v.pos)) + .collect(); + + // Collect all vertices from other + let verts_b: Vec> = other + .polygons + .par_iter() + .flat_map(|poly| poly.vertices.par_iter().map(|v| v.pos)) + .collect(); + + if verts_a.is_empty() || verts_b.is_empty() { + return Mesh::new(); + } + + // For Minkowski, add every point in A to every point in B + let sum_points: Vec<_> = verts_a + .par_iter() + .flat_map(|a| verts_b.par_iter().map(move |b| a + b.coords)) + .map(|v| vec![v.x, v.y, v.z]) + .collect(); + + // Early return if no points generated + if sum_points.is_empty() { + return Mesh::new(); + } + + // Compute convex hull with proper error handling + let hull = match ConvexHullWrapper::try_new(&sum_points, None) { + Ok(h) => h, + Err(_) => return Mesh::new(), // Robust fallback for degenerate cases + }; + let (verts, indices) = hull.vertices_indices(); + + // Reconstruct polygons with proper normal vector calculation + let polygons: Vec> = indices + .par_chunks_exact(3) + .filter_map(|tri| { + let v0 = &verts[tri[0]]; + let v1 = &verts[tri[1]]; + let v2 = &verts[tri[2]]; + + let p0 = Point3::new(v0[0], v0[1], v0[2]); + let p1 = Point3::new(v1[0], v1[1], v1[2]); + let p2 = Point3::new(v2[0], v2[1], v2[2]); + + // Calculate proper normal vector using cross product + let edge1 = p1 - p0; + let edge2 = p2 - p0; + let normal = edge1.cross(&edge2); + + // Filter out degenerate triangles + if normal.norm_squared() > Real::EPSILON { + let normalized_normal = normal.normalize(); + let vv0 = Vertex::new(p0, normalized_normal); + let vv1 = Vertex::new(p1, normalized_normal); + let vv2 = Vertex::new(p2, normalized_normal); + Some(Polygon::new(vec![vv0, vv1, vv2], None)) + } else { + None + } + }) + .collect(); + + Mesh::from_polygons(&polygons, mesh.metadata.clone()) + } +} diff --git a/src/mesh/convex_hull.rs b/src/mesh/algorithm/convex_hull/serial.rs similarity index 70% rename from src/mesh/convex_hull.rs rename to src/mesh/algorithm/convex_hull/serial.rs index 55046fc..ea9d4f4 100644 --- a/src/mesh/convex_hull.rs +++ b/src/mesh/algorithm/convex_hull/serial.rs @@ -1,28 +1,32 @@ -//! The [convex hull](https://en.wikipedia.org/wiki/Convex_hull) of a shape is the smallest convex set that contains it. -//! It may be visualized as the shape enclosed by a rubber band stretched around the subset. -//! -//! This is the set:\ -//! ![Pre-ConvexHull demo image][Pre-ConvexHull demo image] -//! -//! And this is the convex hull of that set:\ -//! ![ConvexHull demo image][ConvexHull demo image] -#![cfg_attr(doc, doc = doc_image_embed::embed_image!("Pre-ConvexHull demo image", "docs/convex_hull_before_nobackground.png"))] -#![cfg_attr(doc, doc = doc_image_embed::embed_image!("ConvexHull demo image", "docs/convex_hull_nobackground.png"))] +//! Serial implementations of convex hull operations. +use super::traits::ConvexHullOps; use crate::float_types::Real; -use crate::mesh::Mesh; -use crate::mesh::polygon::Polygon; -use crate::mesh::vertex::Vertex; +use crate::mesh::{Mesh, polygon::Polygon, vertex::Vertex}; use crate::traits::CSG; use chull::ConvexHullWrapper; use nalgebra::{Point3, Vector3}; use std::fmt::Debug; -impl Mesh { - /// Compute the convex hull of all vertices in this Mesh. - pub fn convex_hull(&self) -> Mesh { +/// Serial implementation of `ConvexHullOps`. +pub struct SerialConvexHullOps; + +impl Default for SerialConvexHullOps { + fn default() -> Self { + Self::new() + } +} + +impl SerialConvexHullOps { + pub const fn new() -> Self { + Self + } +} + +impl ConvexHullOps for SerialConvexHullOps { + fn convex_hull(&self, mesh: &Mesh) -> Mesh { // Gather all (x, y, z) coordinates from the polygons - let points: Vec> = self + let points: Vec> = mesh .polygons .iter() .flat_map(|poly| poly.vertices.iter().map(|v| v.pos)) @@ -54,19 +58,12 @@ impl Mesh { polygons.push(Polygon::new(vec![vv0, vv1, vv2], None)); } - Mesh::from_polygons(&polygons, self.metadata.clone()) + Mesh::from_polygons(&polygons, mesh.metadata.clone()) } - /// Compute the Minkowski sum: self ⊕ other - /// - /// **Mathematical Foundation**: For convex sets A and B, A ⊕ B = {a + b | a ∈ A, b ∈ B}. - /// By the Minkowski sum theorem, the convex hull of all pairwise vertex sums equals - /// the Minkowski sum of the convex hulls of A and B. - /// - /// **Algorithm**: O(|A| × |B|) vertex combinations followed by O(n log n) convex hull computation. - pub fn minkowski_sum(&self, other: &Mesh) -> Mesh { + fn minkowski_sum(&self, mesh: &Mesh, other: &Mesh) -> Mesh { // Collect all vertices (x, y, z) from self - let verts_a: Vec> = self + let verts_a: Vec> = mesh .polygons .iter() .flat_map(|poly| poly.vertices.iter().map(|v| v.pos)) @@ -132,6 +129,6 @@ impl Mesh { }) .collect(); - Mesh::from_polygons(&polygons, self.metadata.clone()) + Mesh::from_polygons(&polygons, mesh.metadata.clone()) } } diff --git a/src/mesh/algorithm/convex_hull/traits.rs b/src/mesh/algorithm/convex_hull/traits.rs new file mode 100644 index 0000000..9aa8f87 --- /dev/null +++ b/src/mesh/algorithm/convex_hull/traits.rs @@ -0,0 +1,13 @@ +//! Traits for convex hull operations. + +use crate::mesh::Mesh; +use std::fmt::Debug; + +/// Trait for convex hull and related operations. +pub trait ConvexHullOps { + /// Computes the convex hull of the mesh. + fn convex_hull(&self, mesh: &Mesh) -> Mesh; + + /// Computes the Minkowski sum of two meshes. + fn minkowski_sum(&self, mesh: &Mesh, other: &Mesh) -> Mesh; +} diff --git a/src/mesh/algorithm/mod.rs b/src/mesh/algorithm/mod.rs new file mode 100644 index 0000000..60c2f97 --- /dev/null +++ b/src/mesh/algorithm/mod.rs @@ -0,0 +1,4 @@ +//! Mesh algorithms with dependency inversion. + +pub mod convex_hull; +pub mod smoothing; diff --git a/src/mesh/algorithm/smoothing/mod.rs b/src/mesh/algorithm/smoothing/mod.rs new file mode 100644 index 0000000..1d1f9cb --- /dev/null +++ b/src/mesh/algorithm/smoothing/mod.rs @@ -0,0 +1,16 @@ +//! Mesh smoothing algorithms. + +pub mod serial; +pub mod traits; + +#[cfg(feature = "parallel")] +pub mod parallel; + +// Re-export core types +pub use traits::SmoothingOps; + +#[cfg(not(feature = "parallel"))] +pub use serial::SerialSmoothingOps; + +#[cfg(feature = "parallel")] +pub use parallel::ParallelSmoothingOps; diff --git a/src/mesh/algorithm/smoothing/parallel.rs b/src/mesh/algorithm/smoothing/parallel.rs new file mode 100644 index 0000000..c69c594 --- /dev/null +++ b/src/mesh/algorithm/smoothing/parallel.rs @@ -0,0 +1,406 @@ +//! Parallel implementations of mesh smoothing operations. + +use super::traits::SmoothingOps; +use crate::float_types::Real; +use crate::mesh::{Mesh, polygon::Polygon, vertex::Vertex}; +use nalgebra::{Point3, Vector3}; +use rayon::prelude::*; +use std::collections::HashMap; +use std::fmt::Debug; + +/// Parallel implementation of `SmoothingOps`. +pub struct ParallelSmoothingOps; + +impl ParallelSmoothingOps { + pub fn new() -> Self { + Self + } +} + +impl SmoothingOps for ParallelSmoothingOps { + fn laplacian_smooth( + &self, + mesh: &Mesh, + lambda: Real, + iterations: usize, + preserve_boundaries: bool, + ) -> Mesh { + let (vertex_map, adjacency) = mesh.build_connectivity(); + let mut smoothed_polygons = mesh.polygons.clone(); + + for iteration in 0..iterations { + // Build current vertex position mapping in parallel + let current_positions: HashMap> = smoothed_polygons + .par_iter() + .flat_map(|polygon| { + polygon.vertices.par_iter().filter_map(|vertex| { + for (pos, idx) in &vertex_map.position_to_index { + if (vertex.pos - pos).norm() < vertex_map.epsilon { + return Some((*idx, vertex.pos)); + } + } + None + }) + }) + .collect(); + + // Compute Laplacian for each vertex in parallel + let laplacian_updates: HashMap> = adjacency + .par_iter() + .map(|(&vertex_idx, neighbors)| { + if let Some(¤t_pos) = current_positions.get(&vertex_idx) { + // Check if this is a boundary vertex + if preserve_boundaries && neighbors.len() < 4 { + // Boundary vertex - skip smoothing + return (vertex_idx, current_pos); + } + + // Compute neighbor average + let (neighbor_sum, valid_neighbors) = neighbors + .par_iter() + .filter_map(|&neighbor_idx| { + current_positions.get(&neighbor_idx).map(|&p| (p.coords, 1)) + }) + .reduce( + || (Vector3::zeros(), 0), + |(sum_a, count_a), (coords_b, count_b)| { + (sum_a + coords_b, count_a + count_b) + }, + ); + + if valid_neighbors > 0 { + let neighbor_avg = + Point3::from(neighbor_sum / valid_neighbors as Real); + let laplacian = neighbor_avg - current_pos; + let new_pos = current_pos + laplacian * lambda; + (vertex_idx, new_pos) + } else { + (vertex_idx, current_pos) + } + } else { + (vertex_idx, Point3::origin()) // Should not happen if connectivity is correct + } + }) + .collect(); + + // Apply updates to mesh vertices in parallel + smoothed_polygons.par_iter_mut().for_each(|polygon| { + for vertex in &mut polygon.vertices { + // Find the global index for this vertex + for (pos, idx) in &vertex_map.position_to_index { + if (vertex.pos - pos).norm() < vertex_map.epsilon { + if let Some(&new_pos) = laplacian_updates.get(idx) { + vertex.pos = new_pos; + } + break; + } + } + } + // Recompute polygon plane and normals after smoothing + polygon.set_new_normal(); + }); + + // Progress feedback for long smoothing operations + if iterations > 10 && iteration % (iterations / 10) == 0 { + eprintln!( + "Smoothing progress: {}/{} iterations", + iteration + 1, + iterations + ); + } + } + + Mesh::from_polygons(&smoothed_polygons, mesh.metadata.clone()) + } + + fn taubin_smooth( + &self, + mesh: &Mesh, + lambda: Real, + mu: Real, + iterations: usize, + preserve_boundaries: bool, + ) -> Mesh { + let (vertex_map, adjacency) = mesh.build_connectivity(); + let mut smoothed_polygons = mesh.polygons.clone(); + + for _ in 0..iterations { + // --- Lambda (shrinking) pass --- + let current_positions: HashMap> = smoothed_polygons + .par_iter() + .flat_map(|polygon| { + polygon.vertices.par_iter().filter_map(|vertex| { + for (pos, idx) in &vertex_map.position_to_index { + if (vertex.pos - pos).norm() < vertex_map.epsilon { + return Some((*idx, vertex.pos)); + } + } + None + }) + }) + .collect(); + + let updates: HashMap> = adjacency + .par_iter() + .map(|(&vertex_idx, neighbors)| { + if let Some(¤t_pos) = current_positions.get(&vertex_idx) { + if preserve_boundaries && neighbors.len() < 4 { + return (vertex_idx, current_pos); + } + + let (neighbor_sum, valid_neighbors) = neighbors + .par_iter() + .filter_map(|&neighbor_idx| { + current_positions.get(&neighbor_idx).map(|&p| (p.coords, 1)) + }) + .reduce( + || (Vector3::zeros(), 0), + |(sum_a, count_a), (coords_b, count_b)| { + (sum_a + coords_b, count_a + count_b) + }, + ); + + if valid_neighbors > 0 { + let neighbor_avg = + Point3::from(neighbor_sum / valid_neighbors as Real); + let laplacian = neighbor_avg - current_pos; + (vertex_idx, current_pos + laplacian * lambda) + } else { + (vertex_idx, current_pos) + } + } else { + (vertex_idx, Point3::origin()) + } + }) + .collect(); + + smoothed_polygons.par_iter_mut().for_each(|polygon| { + for vertex in &mut polygon.vertices { + for (pos, idx) in &vertex_map.position_to_index { + if (vertex.pos - pos).norm() < vertex_map.epsilon { + if let Some(&new_pos) = updates.get(idx) { + vertex.pos = new_pos; + } + break; + } + } + } + }); + + // --- Mu (inflating) pass --- + let current_positions: HashMap> = smoothed_polygons + .par_iter() + .flat_map(|polygon| { + polygon.vertices.par_iter().filter_map(|vertex| { + for (pos, idx) in &vertex_map.position_to_index { + if (vertex.pos - pos).norm() < vertex_map.epsilon { + return Some((*idx, vertex.pos)); + } + } + None + }) + }) + .collect(); + + let updates: HashMap> = adjacency + .par_iter() + .map(|(&vertex_idx, neighbors)| { + if let Some(¤t_pos) = current_positions.get(&vertex_idx) { + if preserve_boundaries && neighbors.len() < 4 { + return (vertex_idx, current_pos); + } + + let (neighbor_sum, valid_neighbors) = neighbors + .par_iter() + .filter_map(|&neighbor_idx| { + current_positions.get(&neighbor_idx).map(|&p| (p.coords, 1)) + }) + .reduce( + || (Vector3::zeros(), 0), + |(sum_a, count_a), (coords_b, count_b)| { + (sum_a + coords_b, count_a + count_b) + }, + ); + + if valid_neighbors > 0 { + let neighbor_avg = + Point3::from(neighbor_sum / valid_neighbors as Real); + let laplacian = neighbor_avg - current_pos; + (vertex_idx, current_pos + laplacian * mu) + } else { + (vertex_idx, current_pos) + } + } else { + (vertex_idx, Point3::origin()) + } + }) + .collect(); + + smoothed_polygons.par_iter_mut().for_each(|polygon| { + for vertex in &mut polygon.vertices { + for (pos, idx) in &vertex_map.position_to_index { + if (vertex.pos - pos).norm() < vertex_map.epsilon { + if let Some(&new_pos) = updates.get(idx) { + vertex.pos = new_pos; + } + break; + } + } + } + }); + } + + // Final pass to recompute normals + smoothed_polygons.par_iter_mut().for_each(|polygon| { + polygon.set_new_normal(); + }); + + Mesh::from_polygons(&smoothed_polygons, mesh.metadata.clone()) + } + + fn adaptive_refine( + &self, + mesh: &Mesh, + quality_threshold: Real, + max_edge_length: Real, + curvature_threshold_deg: Real, + ) -> Mesh { + let qualities = mesh.analyze_triangle_quality(); + let (mut vertex_map, _adjacency) = mesh.build_connectivity(); + let mut polygon_map: HashMap> = HashMap::new(); + + for (poly_idx, poly) in mesh.polygons.iter().enumerate() { + for vertex in &poly.vertices { + let v_idx = vertex_map.get_or_create_index(vertex.pos); + polygon_map.entry(v_idx).or_default().push(poly_idx); + } + } + + // Pre-populate the vertex map with all edge vertices + let edge_vertices: Vec> = mesh + .polygons + .par_iter() + .flat_map(|polygon| { + polygon + .edges() + .collect::>() + .into_par_iter() + .flat_map(|edge| vec![edge.0.pos, edge.1.pos]) + .collect::>() + }) + .collect(); + + for pos in edge_vertices { + vertex_map.get_or_create_index(pos); + } + + let refined_polygons: Vec> = mesh + .polygons + .par_iter() + .enumerate() + .flat_map(|(i, polygon)| { + let mut should_refine = false; + + // Quality and edge length check + if i < qualities.len() { + let quality = &qualities[i]; + if quality.quality_score < quality_threshold + || Self::max_edge_length(&polygon.vertices) > max_edge_length + { + should_refine = true; + } + } + + // Curvature check + if !should_refine { + 'edge_loop: for edge in polygon.edges() { + let v1_idx = vertex_map.get_index(&edge.0.pos).unwrap(); + let v2_idx = vertex_map.get_index(&edge.1.pos).unwrap(); + + if let (Some(p1_indices), Some(p2_indices)) = + (polygon_map.get(&v1_idx), polygon_map.get(&v2_idx)) + { + for &p1_idx in p1_indices { + if p1_idx == i { + continue; + } + for &p2_idx in p2_indices { + if p1_idx == p2_idx { + let other_poly = &mesh.polygons[p1_idx]; + let angle = Mesh::dihedral_angle(polygon, other_poly); + if angle > curvature_threshold_deg.to_radians() { + should_refine = true; + break 'edge_loop; + } + } + } + } + } + } + } + + if should_refine { + let subdivided = polygon.subdivide_triangles(1.try_into().unwrap()); + subdivided + .into_par_iter() + .map(move |triangle| { + let vertices = triangle.to_vec(); + Polygon::new(vertices, polygon.metadata.clone()) + }) + .collect::>() + } else { + vec![polygon.clone()] + } + }) + .collect(); + + Mesh::from_polygons(&refined_polygons, mesh.metadata.clone()) + } + + fn remove_poor_triangles(&self, mesh: &Mesh, min_quality: Real) -> Mesh { + let qualities = mesh.analyze_triangle_quality(); + let filtered_polygons: Vec> = mesh + .polygons + .par_iter() + .enumerate() + .filter_map(|(i, polygon)| { + let keep_triangle = if i < qualities.len() { + let quality = &qualities[i]; + quality.quality_score >= min_quality + && quality.area > Real::EPSILON + && quality.min_angle > (5.0 as Real).to_radians() + && quality.aspect_ratio < 20.0 + } else { + true // Keep if we can't assess quality + }; + + if keep_triangle { + Some(polygon.clone()) + } else { + None + } + }) + .collect(); + + Mesh::from_polygons(&filtered_polygons, mesh.metadata.clone()) + } +} + +impl ParallelSmoothingOps { + /// Calculate maximum edge length in a polygon + pub fn max_edge_length(vertices: &[Vertex]) -> Real { + if vertices.len() < 2 { + return 0.0; + } + + vertices + .par_iter() + .enumerate() + .map(|(i, vertex)| { + let next_vertex = &vertices[(i + 1) % vertices.len()]; + (next_vertex.pos - vertex.pos).norm() + }) + .max_by(|a, b| a.partial_cmp(b).unwrap()) + .unwrap_or(0.0) + } +} diff --git a/src/mesh/smoothing.rs b/src/mesh/algorithm/smoothing/serial.rs similarity index 71% rename from src/mesh/smoothing.rs rename to src/mesh/algorithm/smoothing/serial.rs index d0f1ad2..3af5088 100644 --- a/src/mesh/smoothing.rs +++ b/src/mesh/algorithm/smoothing/serial.rs @@ -1,40 +1,37 @@ +//! Serial implementations of mesh smoothing operations. + +use super::traits::SmoothingOps; use crate::float_types::Real; -use crate::mesh::Mesh; -use crate::mesh::polygon::Polygon; -use crate::mesh::vertex::Vertex; +use crate::mesh::{Mesh, polygon::Polygon, vertex::Vertex}; use nalgebra::Point3; use std::collections::HashMap; use std::fmt::Debug; -impl Mesh { - /// **Mathematical Foundation: True Laplacian Mesh Smoothing with Global Connectivity** - /// - /// Implements proper discrete Laplacian smoothing using global mesh connectivity: - /// - /// ## **Discrete Laplacian Operator** - /// For each vertex v with neighbors N(v): - /// ```text - /// L(v) = (1/|N(v)|) · Σ(n∈N(v)) (n - v) - /// ``` - /// - /// ## **Global Connectivity Benefits** - /// - **Proper Neighborhoods**: Uses actual mesh connectivity, not just polygon edges - /// - **Uniform Weighting**: Each neighbor contributes equally to smoothing - /// - **Boundary Detection**: Automatically detects and preserves mesh boundaries - /// - **Volume Preservation**: Better volume preservation than local smoothing - /// - /// ## **Algorithm Improvements** - /// - **Epsilon-based Vertex Matching**: Robust floating-point coordinate handling - /// - **Manifold Preservation**: Ensures mesh topology is maintained - /// - **Feature Detection**: Can preserve sharp features based on neighbor count - pub fn laplacian_smooth( +/// Serial implementation of `SmoothingOps`. +pub struct SerialSmoothingOps; + +impl Default for SerialSmoothingOps { + fn default() -> Self { + Self::new() + } +} + +impl SerialSmoothingOps { + pub const fn new() -> Self { + Self + } +} + +impl SmoothingOps for SerialSmoothingOps { + fn laplacian_smooth( &self, + mesh: &Mesh, lambda: Real, iterations: usize, preserve_boundaries: bool, ) -> Mesh { - let (vertex_map, adjacency) = self.build_connectivity(); - let mut smoothed_polygons = self.polygons.clone(); + let (vertex_map, adjacency) = mesh.build_connectivity(); + let mut smoothed_polygons = mesh.polygons.clone(); for iteration in 0..iterations { // Build current vertex position mapping @@ -111,34 +108,19 @@ impl Mesh { } } - Mesh::from_polygons(&smoothed_polygons, self.metadata.clone()) + Mesh::from_polygons(&smoothed_polygons, mesh.metadata.clone()) } - /// **Mathematical Foundation: Taubin Mesh Smoothing** - /// - /// Implements Taubin's feature-preserving mesh smoothing algorithm, which reduces - /// shrinkage compared to standard Laplacian smoothing. - /// - /// ## **Taubin's Algorithm** - /// This method involves two steps per iteration: - /// 1. **Shrinking Step**: Apply standard Laplacian smoothing with a positive factor `lambda`. - /// `v' = v + λ * L(v)` - /// 2. **Inflating Step**: Apply a second Laplacian step with a negative factor `mu`. - /// `v'' = v' + μ * L(v')` - /// - /// Typically, `0 < λ < -μ`. A common choice is `mu = -λ / (1 - λ)`. - /// This combination effectively smooths the mesh while minimizing volume loss. - /// - /// Returns a new, smoothed CSG object. - pub fn taubin_smooth( + fn taubin_smooth( &self, + mesh: &Mesh, lambda: Real, mu: Real, iterations: usize, preserve_boundaries: bool, ) -> Mesh { - let (vertex_map, adjacency) = self.build_connectivity(); - let mut smoothed_polygons = self.polygons.clone(); + let (vertex_map, adjacency) = mesh.build_connectivity(); + let mut smoothed_polygons = mesh.polygons.clone(); for _ in 0..iterations { // --- Lambda (shrinking) pass --- @@ -249,44 +231,29 @@ impl Mesh { polygon.set_new_normal(); } - Mesh::from_polygons(&smoothed_polygons, self.metadata.clone()) + Mesh::from_polygons(&smoothed_polygons, mesh.metadata.clone()) } - /// **Mathematical Foundation: Adaptive Mesh Refinement** - /// - /// Intelligently refine mesh based on geometric criteria: - /// - /// ## **Refinement Criteria** - /// - **Quality threshold**: Refine triangles with quality score < threshold - /// - **Size variation**: Refine where edge lengths vary significantly - /// - **Curvature**: Refine high-curvature regions (based on normal variation) - /// - **Feature detection**: Preserve sharp edges and corners - /// - /// ## **Refinement Strategy** - /// 1. **Quality-based**: Subdivide poor-quality triangles - /// 2. **Size-based**: Subdivide triangles larger than target size - /// 3. **Curvature-based**: Subdivide where surface curves significantly - /// - /// This provides better mesh quality compared to uniform subdivision. - pub fn adaptive_refine( + fn adaptive_refine( &self, + mesh: &Mesh, quality_threshold: Real, max_edge_length: Real, curvature_threshold_deg: Real, ) -> Mesh { - let qualities = self.analyze_triangle_quality(); - let (mut vertex_map, _adjacency) = self.build_connectivity(); + let qualities = mesh.analyze_triangle_quality(); + let (mut vertex_map, _adjacency) = mesh.build_connectivity(); let mut refined_polygons = Vec::new(); let mut polygon_map: HashMap> = HashMap::new(); - for (poly_idx, poly) in self.polygons.iter().enumerate() { + for (poly_idx, poly) in mesh.polygons.iter().enumerate() { for vertex in &poly.vertices { let v_idx = vertex_map.get_or_create_index(vertex.pos); polygon_map.entry(v_idx).or_default().push(poly_idx); } } - for (i, polygon) in self.polygons.iter().enumerate() { + for (i, polygon) in mesh.polygons.iter().enumerate() { let mut should_refine = false; // Quality and edge length check @@ -314,8 +281,8 @@ impl Mesh { } for &p2_idx in p2_indices { if p1_idx == p2_idx { - let other_poly = &self.polygons[p1_idx]; - let angle = Self::dihedral_angle(polygon, other_poly); + let other_poly = &mesh.polygons[p1_idx]; + let angle = Mesh::dihedral_angle(polygon, other_poly); if angle > curvature_threshold_deg.to_radians() { should_refine = true; break 'edge_loop; @@ -338,45 +305,14 @@ impl Mesh { } } - Mesh::from_polygons(&refined_polygons, self.metadata.clone()) + Mesh::from_polygons(&refined_polygons, mesh.metadata.clone()) } - /// Calculate maximum edge length in a polygon - fn max_edge_length(vertices: &[Vertex]) -> Real { - if vertices.len() < 2 { - return 0.0; - } - - let mut max_length: Real = 0.0; - for i in 0..vertices.len() { - let j = (i + 1) % vertices.len(); - let edge_length = (vertices[j].pos - vertices[i].pos).norm(); - max_length = max_length.max(edge_length); - } - max_length - } - - /// **Mathematical Foundation: Feature-Preserving Mesh Optimization** - /// - /// Remove poor-quality triangles while preserving important geometric features: - /// - /// ## **Quality-Based Filtering** - /// Remove triangles that meet criteria: - /// - **Sliver triangles**: min_angle < threshold (typically 5°) - /// - **Needle triangles**: aspect_ratio > threshold (typically 20) - /// - **Small triangles**: area < threshold - /// - /// ## **Feature Preservation** - /// - **Sharp edges**: Preserve edges with large dihedral angles - /// - **Boundaries**: Maintain mesh boundaries - /// - **Topology**: Ensure mesh remains manifold - /// - /// Returns cleaned mesh with improved triangle quality. - pub fn remove_poor_triangles(&self, min_quality: Real) -> Mesh { - let qualities = self.analyze_triangle_quality(); + fn remove_poor_triangles(&self, mesh: &Mesh, min_quality: Real) -> Mesh { + let qualities = mesh.analyze_triangle_quality(); let mut filtered_polygons = Vec::new(); - for (i, polygon) in self.polygons.iter().enumerate() { + for (i, polygon) in mesh.polygons.iter().enumerate() { let keep_triangle = if i < qualities.len() { let quality = &qualities[i]; quality.quality_score >= min_quality @@ -392,6 +328,23 @@ impl Mesh { } } - Mesh::from_polygons(&filtered_polygons, self.metadata.clone()) + Mesh::from_polygons(&filtered_polygons, mesh.metadata.clone()) + } +} + +impl SerialSmoothingOps { + /// Calculate maximum edge length in a polygon + pub fn max_edge_length(vertices: &[Vertex]) -> Real { + if vertices.len() < 2 { + return 0.0; + } + + let mut max_length: Real = 0.0; + for i in 0..vertices.len() { + let j = (i + 1) % vertices.len(); + let edge_length = (vertices[j].pos - vertices[i].pos).norm(); + max_length = max_length.max(edge_length); + } + max_length } } diff --git a/src/mesh/algorithm/smoothing/traits.rs b/src/mesh/algorithm/smoothing/traits.rs new file mode 100644 index 0000000..1eec353 --- /dev/null +++ b/src/mesh/algorithm/smoothing/traits.rs @@ -0,0 +1,39 @@ +//! Traits for mesh smoothing operations. + +use crate::float_types::Real; +use crate::mesh::Mesh; +use std::fmt::Debug; + +/// Trait for mesh smoothing and refinement operations. +pub trait SmoothingOps { + /// Applies Laplacian smoothing to the mesh. + fn laplacian_smooth( + &self, + mesh: &Mesh, + lambda: Real, + iterations: usize, + preserve_boundaries: bool, + ) -> Mesh; + + /// Applies Taubin smoothing to the mesh. + fn taubin_smooth( + &self, + mesh: &Mesh, + lambda: Real, + mu: Real, + iterations: usize, + preserve_boundaries: bool, + ) -> Mesh; + + /// Adaptively refines the mesh based on geometric criteria. + fn adaptive_refine( + &self, + mesh: &Mesh, + quality_threshold: Real, + max_edge_length: Real, + curvature_threshold_deg: Real, + ) -> Mesh; + + /// Removes poor-quality triangles from the mesh. + fn remove_poor_triangles(&self, mesh: &Mesh, min_quality: Real) -> Mesh; +} diff --git a/src/mesh/bsp.rs b/src/mesh/bsp.rs deleted file mode 100644 index 0bfa4f5..0000000 --- a/src/mesh/bsp.rs +++ /dev/null @@ -1,338 +0,0 @@ -//! [BSP](https://en.wikipedia.org/wiki/Binary_space_partitioning) tree node structure and operations - -#[cfg(not(feature = "parallel"))] -use crate::float_types::EPSILON; - -#[cfg(not(feature = "parallel"))] -use crate::mesh::vertex::Vertex; - -use crate::float_types::Real; -use crate::mesh::plane::{BACK, COPLANAR, FRONT, Plane, SPANNING}; -use crate::mesh::polygon::Polygon; -use std::fmt::Debug; - -/// A [BSP](https://en.wikipedia.org/wiki/Binary_space_partitioning) tree node, containing polygons plus optional front/back subtrees -#[derive(Debug, Clone)] -pub struct Node { - /// Splitting plane for this node *or* **None** for a leaf that - /// only stores polygons. - pub plane: Option, - - /// Polygons in *front* half‑spaces. - pub front: Option>>, - - /// Polygons in *back* half‑spaces. - pub back: Option>>, - - /// Polygons that lie *exactly* on `plane` - /// (after the node has been built). - pub polygons: Vec>, -} - -impl Default for Node { - fn default() -> Self { - Self::new() - } -} - -impl Node { - /// Create a new empty BSP node - pub const fn new() -> Self { - Self { - plane: None, - front: None, - back: None, - polygons: Vec::new(), - } - } - - /// Creates a new BSP node from polygons - pub fn from_polygons(polygons: &[Polygon]) -> Self { - let mut node = Self::new(); - if !polygons.is_empty() { - node.build(polygons); - } - node - } - - /// Invert all polygons in the BSP tree - #[cfg(not(feature = "parallel"))] - pub fn invert(&mut self) { - // Flip all polygons and plane in this node - self.polygons.iter_mut().for_each(|p| p.flip()); - if let Some(ref mut plane) = self.plane { - plane.flip(); - } - - if let Some(ref mut front) = self.front { - front.invert(); - } - if let Some(ref mut back) = self.back { - back.invert(); - } - - std::mem::swap(&mut self.front, &mut self.back); - } - - pub fn pick_best_splitting_plane(&self, polygons: &[Polygon]) -> Plane { - const K_SPANS: Real = 8.0; // Weight for spanning polygons - const K_BALANCE: Real = 1.0; // Weight for front/back balance - - let mut best_plane = polygons[0].plane.clone(); - let mut best_score = Real::MAX; - - // Take a sample of polygons as candidate planes - let sample_size = polygons.len().min(20); - for p in polygons.iter().take(sample_size) { - let plane = &p.plane; - let mut num_front = 0; - let mut num_back = 0; - let mut num_spanning = 0; - - for poly in polygons { - match plane.classify_polygon(poly) { - COPLANAR => {}, // Not counted for balance - FRONT => num_front += 1, - BACK => num_back += 1, - SPANNING => num_spanning += 1, - _ => num_spanning += 1, // Treat any other combination as spanning - } - } - - let score = K_SPANS * num_spanning as Real - + K_BALANCE * ((num_front - num_back) as Real).abs(); - - if score < best_score { - best_score = score; - best_plane = plane.clone(); - } - } - best_plane - } - - /// Recursively remove all polygons in `polygons` that are inside this BSP tree - /// **Mathematical Foundation**: Uses plane classification to determine polygon visibility. - /// Polygons entirely in BACK half-space are clipped (removed). - /// **Algorithm**: O(n log d) where n is polygon count, d is tree depth. - #[cfg(not(feature = "parallel"))] - pub fn clip_polygons(&self, polygons: &[Polygon]) -> Vec> { - // If this node has no plane (i.e. it’s empty), just return - if self.plane.is_none() { - return polygons.to_vec(); - } - - let plane = self.plane.as_ref().unwrap(); - - // Pre-allocate for better performance - let mut front_polys = Vec::with_capacity(polygons.len()); - let mut back_polys = Vec::with_capacity(polygons.len()); - - // Optimized polygon splitting with iterator patterns - for polygon in polygons { - let (coplanar_front, coplanar_back, mut front_parts, mut back_parts) = - plane.split_polygon(polygon); - - // Efficient coplanar polygon classification using iterator chain - for coplanar_poly in coplanar_front.into_iter().chain(coplanar_back.into_iter()) { - if plane.orient_plane(&coplanar_poly.plane) == FRONT { - front_parts.push(coplanar_poly); - } else { - back_parts.push(coplanar_poly); - } - } - - front_polys.append(&mut front_parts); - back_polys.append(&mut back_parts); - } - - // Recursively clip with optimized pattern - let mut result = if let Some(front_node) = &self.front { - front_node.clip_polygons(&front_polys) - } else { - front_polys - }; - - if let Some(back_node) = &self.back { - result.extend(back_node.clip_polygons(&back_polys)); - } - - result - } - - /// Remove all polygons in this BSP tree that are inside the other BSP tree - #[cfg(not(feature = "parallel"))] - pub fn clip_to(&mut self, bsp: &Node) { - self.polygons = bsp.clip_polygons(&self.polygons); - if let Some(ref mut front) = self.front { - front.clip_to(bsp); - } - if let Some(ref mut back) = self.back { - back.clip_to(bsp); - } - } - - /// Return all polygons in this BSP tree using an iterative approach, - /// avoiding potential stack overflow of recursive approach - pub fn all_polygons(&self) -> Vec> { - let mut result = Vec::new(); - let mut stack = vec![self]; - - while let Some(node) = stack.pop() { - result.extend_from_slice(&node.polygons); - - // Use iterator to add child nodes more efficiently - stack.extend( - [&node.front, &node.back] - .iter() - .filter_map(|child| child.as_ref().map(|boxed| boxed.as_ref())), - ); - } - result - } - - /// Build a BSP tree from the given polygons - #[cfg(not(feature = "parallel"))] - pub fn build(&mut self, polygons: &[Polygon]) { - if polygons.is_empty() { - return; - } - - // Choose the best splitting plane using a heuristic if not already set. - if self.plane.is_none() { - self.plane = Some(self.pick_best_splitting_plane(polygons)); - } - let plane = self.plane.as_ref().unwrap(); - - // Pre-allocate with estimated capacity for better performance - let mut front = Vec::with_capacity(polygons.len() / 2); - let mut back = Vec::with_capacity(polygons.len() / 2); - - // Optimized polygon classification using iterator pattern - // **Mathematical Theorem**: Each polygon is classified relative to the splitting plane - for polygon in polygons { - let (coplanar_front, coplanar_back, mut front_parts, mut back_parts) = - plane.split_polygon(polygon); - - // Extend collections efficiently with iterator chains - self.polygons.extend(coplanar_front); - self.polygons.extend(coplanar_back); - front.append(&mut front_parts); - back.append(&mut back_parts); - } - - // Build child nodes using lazy initialization pattern for memory efficiency - if !front.is_empty() { - self.front - .get_or_insert_with(|| Box::new(Node::new())) - .build(&front); - } - - if !back.is_empty() { - self.back - .get_or_insert_with(|| Box::new(Node::new())) - .build(&back); - } - } - - /// Slices this BSP node with `slicing_plane`, returning: - /// - All polygons that are coplanar with the plane (within EPSILON), - /// - A list of line‐segment intersections (each a [Vertex; 2]) from polygons that span the plane. - #[cfg(not(feature = "parallel"))] - pub fn slice(&self, slicing_plane: &Plane) -> (Vec>, Vec<[Vertex; 2]>) { - let all_polys = self.all_polygons(); - - let mut coplanar_polygons = Vec::new(); - let mut intersection_edges = Vec::new(); - - for poly in &all_polys { - let vcount = poly.vertices.len(); - if vcount < 2 { - continue; // degenerate polygon => skip - } - - // Use iterator chain to compute vertex types more efficiently - let types: Vec<_> = poly - .vertices - .iter() - .map(|vertex| slicing_plane.orient_point(&vertex.pos)) - .collect(); - - let polygon_type = types.iter().fold(0, |acc, &vertex_type| acc | vertex_type); - - // Based on the combined classification of its vertices: - match polygon_type { - COPLANAR => { - // The entire polygon is in the plane, so push it to the coplanar list. - coplanar_polygons.push(poly.clone()); - }, - - FRONT | BACK => { - // Entirely on one side => no intersection. We skip it. - }, - - SPANNING => { - // The polygon crosses the plane. We'll gather the intersection points - // (the new vertices introduced on edges that cross the plane). - let crossing_points: Vec<_> = (0..vcount) - .filter_map(|i| { - let j = (i + 1) % vcount; - let ti = types[i]; - let tj = types[j]; - let vi = &poly.vertices[i]; - let vj = &poly.vertices[j]; - - if (ti | tj) == SPANNING { - let denom = slicing_plane.normal().dot(&(vj.pos - vi.pos)); - if denom.abs() > EPSILON { - let intersection = (slicing_plane.offset() - - slicing_plane.normal().dot(&vi.pos.coords)) - / denom; - Some(vi.interpolate(vj, intersection)) - } else { - None - } - } else { - None - } - }) - .collect(); - - // Convert crossing points to intersection edges - intersection_edges.extend( - crossing_points - .chunks_exact(2) - .map(|chunk| [chunk[0].clone(), chunk[1].clone()]), - ); - }, - - _ => { - // Shouldn't happen in a typical classification, but we can ignore - }, - } - } - - (coplanar_polygons, intersection_edges) - } -} - -#[cfg(test)] -mod tests { - use crate::mesh::bsp::Node; - use crate::mesh::polygon::Polygon; - use crate::mesh::vertex::Vertex; - use nalgebra::{Point3, Vector3}; - - #[test] - fn test_bsp_basic_functionality() { - let vertices = vec![ - Vertex::new(Point3::new(0.0, 0.0, 0.0), Vector3::new(0.0, 0.0, 1.0)), - Vertex::new(Point3::new(1.0, 0.0, 0.0), Vector3::new(0.0, 0.0, 1.0)), - Vertex::new(Point3::new(0.5, 1.0, 0.0), Vector3::new(0.0, 0.0, 1.0)), - ]; - let polygon: Polygon = Polygon::new(vertices, None); - let polygons = vec![polygon]; - - let node = Node::from_polygons(&polygons); - assert!(!node.all_polygons().is_empty()); - } -} diff --git a/src/mesh/bsp/mod.rs b/src/mesh/bsp/mod.rs new file mode 100644 index 0000000..cff3bb6 --- /dev/null +++ b/src/mesh/bsp/mod.rs @@ -0,0 +1,121 @@ +//! Binary Space Partitioning (BSP) tree implementation +//! +//! This module provides BSP tree operations with dependency inversion, +//! allowing for different algorithm implementations (serial/parallel). + +pub mod node; +pub mod traits; + +#[cfg(not(feature = "parallel"))] +pub mod serial; + +#[cfg(feature = "parallel")] +pub mod parallel; + +// Re-export core types for backward compatibility +pub use node::Node; +pub use traits::{BalancedSplittingStrategy, BspOps, SplittingPlaneStrategy}; + +#[cfg(not(feature = "parallel"))] +pub use serial::SerialBspOps; + +#[cfg(feature = "parallel")] +pub use parallel::ParallelBspOps; + +// Backward compatibility implementations on Node +use crate::mesh::plane::Plane; +use crate::mesh::polygon::Polygon; +use crate::mesh::vertex::Vertex; +use std::fmt::Debug; + +impl Node { + /// Creates a new BSP node from polygons + pub fn from_polygons(polygons: &[Polygon]) -> Self { + let mut node = Self::new(); + if !polygons.is_empty() { + node.build(polygons); + } + node + } + + /// Invert all polygons in the BSP tree + #[cfg(not(feature = "parallel"))] + pub fn invert(&mut self) { + let ops = SerialBspOps::new(); + ops.invert(self); + } + + #[cfg(feature = "parallel")] + pub fn invert(&mut self) { + let ops = ParallelBspOps::new(); + ops.invert(self); + } + + /// Pick the best splitting plane using the default strategy + pub fn pick_best_splitting_plane(&self, polygons: &[Polygon]) -> Plane { + let strategy = BalancedSplittingStrategy::default(); + strategy.pick_best_splitting_plane(polygons) + } + + /// Recursively remove all polygons that are inside this BSP tree + #[cfg(not(feature = "parallel"))] + pub fn clip_polygons(&self, polygons: &[Polygon]) -> Vec> { + let ops = SerialBspOps::new(); + ops.clip_polygons(self, polygons) + } + + #[cfg(feature = "parallel")] + pub fn clip_polygons(&self, polygons: &[Polygon]) -> Vec> { + let ops = ParallelBspOps::new(); + ops.clip_polygons(self, polygons) + } + + /// Remove all polygons in this BSP tree that are inside the other BSP tree + #[cfg(not(feature = "parallel"))] + pub fn clip_to(&mut self, bsp: &Node) { + let ops = SerialBspOps::new(); + ops.clip_to(self, bsp); + } + + #[cfg(feature = "parallel")] + pub fn clip_to(&mut self, bsp: &Node) { + let ops = ParallelBspOps::new(); + ops.clip_to(self, bsp); + } + + /// Return all polygons in this BSP tree + pub fn all_polygons(&self) -> Vec> { + #[cfg(not(feature = "parallel"))] + let ops = SerialBspOps::new(); + #[cfg(feature = "parallel")] + let ops = ParallelBspOps::new(); + + ops.all_polygons(self) + } + + /// Build a BSP tree from the given polygons + #[cfg(not(feature = "parallel"))] + pub fn build(&mut self, polygons: &[Polygon]) { + let ops = SerialBspOps::new(); + ops.build(self, polygons); + } + + #[cfg(feature = "parallel")] + pub fn build(&mut self, polygons: &[Polygon]) { + let ops = ParallelBspOps::new(); + ops.build(self, polygons); + } + + /// Slices this BSP node with the given plane + #[cfg(not(feature = "parallel"))] + pub fn slice(&self, slicing_plane: &Plane) -> (Vec>, Vec<[Vertex; 2]>) { + let ops = SerialBspOps::new(); + ops.slice(self, slicing_plane) + } + + #[cfg(feature = "parallel")] + pub fn slice(&self, slicing_plane: &Plane) -> (Vec>, Vec<[Vertex; 2]>) { + let ops = ParallelBspOps::new(); + ops.slice(self, slicing_plane) + } +} diff --git a/src/mesh/bsp/node.rs b/src/mesh/bsp/node.rs new file mode 100644 index 0000000..79ba179 --- /dev/null +++ b/src/mesh/bsp/node.rs @@ -0,0 +1,41 @@ +//! BSP tree node data structure + +use crate::mesh::plane::Plane; +use crate::mesh::polygon::Polygon; +use std::fmt::Debug; + +/// A BSP tree node, containing polygons plus optional front/back subtrees +#[derive(Debug, Clone)] +pub struct Node { + /// Splitting plane for this node *or* **None** for a leaf that + /// only stores polygons. + pub plane: Option, + + /// Polygons in *front* half‑spaces. + pub front: Option>>, + + /// Polygons in *back* half‑spaces. + pub back: Option>>, + + /// Polygons that lie *exactly* on `plane` + /// (after the node has been built). + pub polygons: Vec>, +} + +impl Default for Node { + fn default() -> Self { + Self::new() + } +} + +impl Node { + /// Create a new empty BSP node + pub const fn new() -> Self { + Self { + plane: None, + front: None, + back: None, + polygons: Vec::new(), + } + } +} diff --git a/src/mesh/bsp/parallel.rs b/src/mesh/bsp/parallel.rs new file mode 100644 index 0000000..0de81d6 --- /dev/null +++ b/src/mesh/bsp/parallel.rs @@ -0,0 +1,290 @@ +//! Parallel implementation of BSP operations + +#[cfg(feature = "parallel")] +use rayon::prelude::*; + +use crate::float_types::EPSILON; +use crate::mesh::bsp::node::Node; +use crate::mesh::bsp::traits::{BalancedSplittingStrategy, BspOps, SplittingPlaneStrategy}; +use crate::mesh::plane::{BACK, COPLANAR, FRONT, Plane, SPANNING}; +use crate::mesh::polygon::Polygon; +use crate::mesh::vertex::Vertex; +use std::fmt::Debug; + +/// Parallel implementation of BSP operations +#[cfg(feature = "parallel")] +pub struct ParallelBspOps< + SP: SplittingPlaneStrategy = BalancedSplittingStrategy, + S: Clone = (), +> { + splitting_strategy: SP, + _phantom: std::marker::PhantomData, +} + +#[cfg(feature = "parallel")] +impl ParallelBspOps { + pub fn new() -> Self { + Self { + splitting_strategy: BalancedSplittingStrategy::default(), + _phantom: std::marker::PhantomData, + } + } +} + +#[cfg(feature = "parallel")] +impl Default for ParallelBspOps { + fn default() -> Self { + Self::new() + } +} + +#[cfg(feature = "parallel")] +impl, S: Clone> ParallelBspOps { + pub const fn with_strategy(strategy: SP) -> Self { + Self { + splitting_strategy: strategy, + _phantom: std::marker::PhantomData, + } + } +} + +#[cfg(feature = "parallel")] +impl + Sync, S: Clone + Send + Sync + Debug> BspOps + for ParallelBspOps +{ + fn invert(&self, node: &mut Node) { + // Use iterative approach with a stack to avoid stack overflow + let mut stack = vec![node]; + + while let Some(current) = stack.pop() { + // Flip all polygons and plane in this node + current.polygons.par_iter_mut().for_each(|p| p.flip()); + if let Some(ref mut plane) = current.plane { + plane.flip(); + } + + // Swap front and back children + std::mem::swap(&mut current.front, &mut current.back); + + // Add children to stack for processing + if let Some(ref mut front) = current.front { + stack.push(front.as_mut()); + } + if let Some(ref mut back) = current.back { + stack.push(back.as_mut()); + } + } + } + + fn clip_polygons(&self, node: &Node, polygons: &[Polygon]) -> Vec> { + // If this node has no plane, just return the original set + if node.plane.is_none() { + return polygons.to_vec(); + } + let plane = node.plane.as_ref().unwrap(); + + // Split each polygon in parallel; gather results + let (coplanar_front, coplanar_back, mut front, mut back) = polygons + .par_iter() + .map(|poly| plane.split_polygon(poly)) + .reduce( + || (Vec::new(), Vec::new(), Vec::new(), Vec::new()), + |mut acc, x| { + acc.0.extend(x.0); + acc.1.extend(x.1); + acc.2.extend(x.2); + acc.3.extend(x.3); + acc + }, + ); + + // Decide where to send the coplanar polygons + coplanar_front + .into_iter() + .chain(coplanar_back) + .for_each(|cp| { + if plane.orient_plane(&cp.plane) == FRONT { + front.push(cp); + } else { + back.push(cp); + } + }); + + // Process front and back using parallel iterators to avoid recursive join + let mut result = if let Some(ref f) = node.front { + self.clip_polygons(f, &front) + } else { + front + }; + + if let Some(ref b) = node.back { + result.extend(self.clip_polygons(b, &back)); + } + + result + } + + fn clip_to(&self, node: &mut Node, bsp: &Node) { + // Use iterative approach with a stack to avoid recursive stack overflow + let mut stack = vec![node]; + + while let Some(current) = stack.pop() { + // Clip polygons at this node + current.polygons = self.clip_polygons(bsp, ¤t.polygons); + + // Add children to stack for processing + if let Some(ref mut front) = current.front { + stack.push(front.as_mut()); + } + if let Some(ref mut back) = current.back { + stack.push(back.as_mut()); + } + } + } + + fn build(&self, node: &mut Node, polygons: &[Polygon]) { + if polygons.is_empty() { + return; + } + + // Choose splitting plane if not already set + if node.plane.is_none() { + node.plane = Some(self.splitting_strategy.pick_best_splitting_plane(polygons)); + } + let plane = node.plane.as_ref().unwrap(); + + // Split polygons in parallel + let (mut coplanar_front, mut coplanar_back, front, back) = + polygons.par_iter().map(|p| plane.split_polygon(p)).reduce( + || (Vec::new(), Vec::new(), Vec::new(), Vec::new()), + |mut acc, x| { + acc.0.extend(x.0); + acc.1.extend(x.1); + acc.2.extend(x.2); + acc.3.extend(x.3); + acc + }, + ); + + // Append coplanar fronts/backs to node.polygons + node.polygons.append(&mut coplanar_front); + node.polygons.append(&mut coplanar_back); + + // Build children sequentially to avoid stack overflow from recursive join + if !front.is_empty() { + let mut front_node = node.front.take().unwrap_or_else(|| Box::new(Node::new())); + self.build(&mut front_node, &front); + node.front = Some(front_node); + } + + if !back.is_empty() { + let mut back_node = node.back.take().unwrap_or_else(|| Box::new(Node::new())); + self.build(&mut back_node, &back); + node.back = Some(back_node); + } + } + + fn all_polygons(&self, node: &Node) -> Vec> { + // Use serial version as parallel collection doesn't provide significant benefit here + let mut result = Vec::new(); + let mut stack = vec![node]; + + while let Some(current) = stack.pop() { + result.extend_from_slice(¤t.polygons); + + // Use iterator to add child nodes more efficiently + stack.extend( + [¤t.front, ¤t.back] + .iter() + .filter_map(|child| child.as_ref().map(|boxed| boxed.as_ref())), + ); + } + result + } + + fn slice( + &self, + node: &Node, + slicing_plane: &Plane, + ) -> (Vec>, Vec<[Vertex; 2]>) { + // Collect all polygons + let all_polys = self.all_polygons(node); + + // Process polygons in parallel + let (coplanar_polygons, intersection_edges) = all_polys + .par_iter() + .map(|poly| { + let vcount = poly.vertices.len(); + if vcount < 2 { + // Degenerate => skip + return (Vec::new(), Vec::new()); + } + + let types: Vec<_> = poly + .vertices + .iter() + .map(|vertex| slicing_plane.orient_point(&vertex.pos)) + .collect(); + + let polygon_type = types.iter().fold(0, |acc, &vertex_type| acc | vertex_type); + + match polygon_type { + COPLANAR => { + // Entire polygon in plane + (vec![poly.clone()], Vec::new()) + }, + FRONT | BACK => { + // Entirely on one side => no intersection + (Vec::new(), Vec::new()) + }, + SPANNING => { + // The polygon crosses the plane => gather intersection edges + let crossing_points: Vec<_> = (0..vcount) + .filter_map(|i| { + let j = (i + 1) % vcount; + let ti = types[i]; + let tj = types[j]; + let vi = &poly.vertices[i]; + let vj = &poly.vertices[j]; + + if (ti | tj) == SPANNING { + // The param intersection at which plane intersects the edge + let denom = slicing_plane.normal().dot(&(vj.pos - vi.pos)); + if denom.abs() > EPSILON { + let intersection = (slicing_plane.offset() + - slicing_plane.normal().dot(&vi.pos.coords)) + / denom; + // Interpolate: + Some(vi.interpolate(vj, intersection)) + } else { + None + } + } else { + None + } + }) + .collect(); + + // Pair up intersection points => edges + let edges: Vec<_> = crossing_points + .chunks_exact(2) + .map(|chunk| [chunk[0].clone(), chunk[1].clone()]) + .collect(); + + (Vec::new(), edges) + }, + _ => (Vec::new(), Vec::new()), + } + }) + .reduce( + || (Vec::new(), Vec::new()), + |mut acc, x| { + acc.0.extend(x.0); + acc.1.extend(x.1); + acc + }, + ); + + (coplanar_polygons, intersection_edges) + } +} diff --git a/src/mesh/bsp/serial.rs b/src/mesh/bsp/serial.rs new file mode 100644 index 0000000..e6b5ee5 --- /dev/null +++ b/src/mesh/bsp/serial.rs @@ -0,0 +1,259 @@ +//! Serial implementation of BSP operations + +use crate::float_types::EPSILON; +use crate::mesh::bsp::node::Node; +use crate::mesh::bsp::traits::{BalancedSplittingStrategy, BspOps, SplittingPlaneStrategy}; +use crate::mesh::plane::{BACK, COPLANAR, FRONT, Plane, SPANNING}; +use crate::mesh::polygon::Polygon; +use crate::mesh::vertex::Vertex; +use std::fmt::Debug; + +/// Serial implementation of BSP operations +pub struct SerialBspOps< + SP: SplittingPlaneStrategy = BalancedSplittingStrategy, + S: Clone = (), +> { + splitting_strategy: SP, + _phantom: std::marker::PhantomData, +} + +impl Default for SerialBspOps { + fn default() -> Self { + Self::new() + } +} + +impl SerialBspOps { + pub fn new() -> Self { + Self { + splitting_strategy: BalancedSplittingStrategy::default(), + _phantom: std::marker::PhantomData, + } + } +} + +impl, S: Clone> SerialBspOps { + pub const fn with_strategy(strategy: SP) -> Self { + Self { + splitting_strategy: strategy, + _phantom: std::marker::PhantomData, + } + } +} + +impl, S: Clone + Send + Sync + Debug> BspOps + for SerialBspOps +{ + fn invert(&self, node: &mut Node) { + // Use iterative approach with a stack + let mut stack = vec![node]; + + while let Some(current) = stack.pop() { + // Flip all polygons and plane in this node + current.polygons.iter_mut().for_each(|p| p.flip()); + if let Some(ref mut plane) = current.plane { + plane.flip(); + } + + // Swap front and back + std::mem::swap(&mut current.front, &mut current.back); + + // Add children to stack + if let Some(ref mut front) = current.front { + stack.push(front.as_mut()); + } + if let Some(ref mut back) = current.back { + stack.push(back.as_mut()); + } + } + } + + fn clip_polygons(&self, node: &Node, polygons: &[Polygon]) -> Vec> { + // If this node has no plane, just return + if node.plane.is_none() { + return polygons.to_vec(); + } + + let plane = node.plane.as_ref().unwrap(); + + // Pre-allocate for better performance + let mut front_polys = Vec::with_capacity(polygons.len()); + let mut back_polys = Vec::with_capacity(polygons.len()); + + // Optimized polygon splitting with iterator patterns + for polygon in polygons { + let (coplanar_front, coplanar_back, mut front_parts, mut back_parts) = + plane.split_polygon(polygon); + + // Efficient coplanar polygon classification using iterator chain + coplanar_front + .into_iter() + .chain(coplanar_back.into_iter()) + .for_each(|coplanar_poly| { + if plane.orient_plane(&coplanar_poly.plane) == FRONT { + front_parts.push(coplanar_poly); + } else { + back_parts.push(coplanar_poly); + } + }); + + front_polys.append(&mut front_parts); + back_polys.append(&mut back_parts); + } + + // Recursively clip with optimized pattern + let mut result = if let Some(front_node) = &node.front { + self.clip_polygons(front_node, &front_polys) + } else { + front_polys + }; + + if let Some(back_node) = &node.back { + result.extend(self.clip_polygons(back_node, &back_polys)); + } + + result + } + + fn clip_to(&self, node: &mut Node, bsp: &Node) { + node.polygons = self.clip_polygons(bsp, &node.polygons); + + if let Some(ref mut front) = node.front { + self.clip_to(front, bsp); + } + + if let Some(ref mut back) = node.back { + self.clip_to(back, bsp); + } + } + + fn all_polygons(&self, node: &Node) -> Vec> { + let mut result = Vec::new(); + let mut stack = vec![node]; + + while let Some(current) = stack.pop() { + result.extend_from_slice(¤t.polygons); + + // Use iterator to add child nodes more efficiently + stack.extend( + [¤t.front, ¤t.back] + .iter() + .filter_map(|child| child.as_ref().map(|boxed| boxed.as_ref())), + ); + } + result + } + + fn build(&self, node: &mut Node, polygons: &[Polygon]) { + if polygons.is_empty() { + return; + } + + // Choose the best splitting plane if not already set + if node.plane.is_none() { + node.plane = Some(self.splitting_strategy.pick_best_splitting_plane(polygons)); + } + let plane = node.plane.as_ref().unwrap(); + + // Pre-allocate with estimated capacity + let mut front = Vec::with_capacity(polygons.len() / 2); + let mut back = Vec::with_capacity(polygons.len() / 2); + + // Optimized polygon classification using iterator pattern + for polygon in polygons { + let (coplanar_front, coplanar_back, mut front_parts, mut back_parts) = + plane.split_polygon(polygon); + + // Extend collections efficiently with iterator chains + node.polygons.extend(coplanar_front); + node.polygons.extend(coplanar_back); + front.append(&mut front_parts); + back.append(&mut back_parts); + } + + // Build child nodes using lazy initialization pattern + if !front.is_empty() { + node.front.get_or_insert_with(|| Box::new(Node::new())); + self.build(node.front.as_mut().unwrap(), &front); + } + + if !back.is_empty() { + node.back.get_or_insert_with(|| Box::new(Node::new())); + self.build(node.back.as_mut().unwrap(), &back); + } + } + + fn slice( + &self, + node: &Node, + slicing_plane: &Plane, + ) -> (Vec>, Vec<[Vertex; 2]>) { + let all_polys = self.all_polygons(node); + + let mut coplanar_polygons = Vec::new(); + let mut intersection_edges = Vec::new(); + + for poly in &all_polys { + let vcount = poly.vertices.len(); + if vcount < 2 { + continue; // degenerate polygon => skip + } + + // Use iterator to compute vertex types + let types: Vec<_> = poly + .vertices + .iter() + .map(|vertex| slicing_plane.orient_point(&vertex.pos)) + .collect(); + + let polygon_type = types.iter().fold(0, |acc, &vertex_type| acc | vertex_type); + + match polygon_type { + COPLANAR => { + coplanar_polygons.push(poly.clone()); + }, + FRONT | BACK => { + // Entirely on one side => no intersection + }, + SPANNING => { + // The polygon crosses the plane + let crossing_points: Vec<_> = (0..vcount) + .filter_map(|i| { + let j = (i + 1) % vcount; + let ti = types[i]; + let tj = types[j]; + let vi = &poly.vertices[i]; + let vj = &poly.vertices[j]; + + if (ti | tj) == SPANNING { + let denom = slicing_plane.normal().dot(&(vj.pos - vi.pos)); + if denom.abs() > EPSILON { + let intersection = (slicing_plane.offset() + - slicing_plane.normal().dot(&vi.pos.coords)) + / denom; + Some(vi.interpolate(vj, intersection)) + } else { + None + } + } else { + None + } + }) + .collect(); + + // Convert crossing points to intersection edges + intersection_edges.extend( + crossing_points + .chunks_exact(2) + .map(|chunk| [chunk[0].clone(), chunk[1].clone()]), + ); + }, + _ => { + // Shouldn't happen in typical classification + }, + } + } + + (coplanar_polygons, intersection_edges) + } +} diff --git a/src/mesh/bsp/traits.rs b/src/mesh/bsp/traits.rs new file mode 100644 index 0000000..97d9cd1 --- /dev/null +++ b/src/mesh/bsp/traits.rs @@ -0,0 +1,86 @@ +//! Traits defining BSP tree operations for dependency inversion + +use crate::float_types::Real; +use crate::mesh::bsp::node::Node; +use crate::mesh::plane::Plane; +use crate::mesh::polygon::Polygon; +use crate::mesh::vertex::Vertex; + +/// Core BSP operations trait - implements algorithms on BSP nodes +pub trait BspOps { + /// Invert all polygons in the BSP tree + fn invert(&self, node: &mut Node); + + /// Recursively remove all polygons that are inside this BSP tree + fn clip_polygons(&self, node: &Node, polygons: &[Polygon]) -> Vec>; + + /// Remove all polygons in this BSP tree that are inside the other BSP tree + fn clip_to(&self, node: &mut Node, other: &Node); + + /// Build a BSP tree from the given polygons + fn build(&self, node: &mut Node, polygons: &[Polygon]); + + /// Return all polygons in this BSP tree + fn all_polygons(&self, node: &Node) -> Vec>; + + /// Slices this BSP node with the given plane + fn slice( + &self, + node: &Node, + slicing_plane: &Plane, + ) -> (Vec>, Vec<[Vertex; 2]>); +} + +/// Trait for picking optimal splitting planes +pub trait SplittingPlaneStrategy { + /// Pick the best splitting plane from a set of polygons + fn pick_best_splitting_plane(&self, polygons: &[Polygon]) -> Plane; +} + +/// Default splitting plane strategy using balanced heuristic +pub struct BalancedSplittingStrategy { + pub span_weight: Real, + pub balance_weight: Real, +} + +impl Default for BalancedSplittingStrategy { + fn default() -> Self { + Self { + span_weight: 8.0, + balance_weight: 1.0, + } + } +} + +impl SplittingPlaneStrategy for BalancedSplittingStrategy { + fn pick_best_splitting_plane(&self, polygons: &[Polygon]) -> Plane { + let mut best_plane = polygons[0].plane.clone(); + let mut best_score = Real::MAX; + + // Take a sample of polygons as candidate planes + let sample_size = polygons.len().min(20); + + polygons.iter().take(sample_size).for_each(|p| { + let plane = &p.plane; + let (num_front, num_back, num_spanning) = polygons + .iter() + .map(|poly| match plane.classify_polygon(poly) { + crate::mesh::plane::COPLANAR => (0, 0, 0), + crate::mesh::plane::FRONT => (1, 0, 0), + crate::mesh::plane::BACK => (0, 1, 0), + _ => (0, 0, 1), + }) + .fold((0, 0, 0), |acc, x| (acc.0 + x.0, acc.1 + x.1, acc.2 + x.2)); + + let score = self.span_weight * num_spanning as Real + + self.balance_weight * ((num_front - num_back) as Real).abs(); + + if score < best_score { + best_score = score; + best_plane = plane.clone(); + } + }); + + best_plane + } +} diff --git a/src/mesh/bsp_parallel.rs b/src/mesh/bsp_parallel.rs deleted file mode 100644 index b1d5220..0000000 --- a/src/mesh/bsp_parallel.rs +++ /dev/null @@ -1,247 +0,0 @@ -//! Parallel versions of [BSP](https://en.wikipedia.org/wiki/Binary_space_partitioning) operations - -use crate::mesh::bsp::Node; -use std::fmt::Debug; - -#[cfg(feature = "parallel")] -use crate::mesh::plane::{BACK, COPLANAR, FRONT, Plane, SPANNING}; - -#[cfg(feature = "parallel")] -use rayon::prelude::*; - -#[cfg(feature = "parallel")] -use crate::mesh::Polygon; - -#[cfg(feature = "parallel")] -use crate::mesh::Vertex; - -#[cfg(feature = "parallel")] -use crate::float_types::EPSILON; - -impl Node { - /// Invert all polygons in the BSP tree using iterative approach to avoid stack overflow - #[cfg(feature = "parallel")] - pub fn invert(&mut self) { - // Use iterative approach with a stack to avoid recursive stack overflow - let mut stack = vec![self]; - - while let Some(node) = stack.pop() { - // Flip all polygons and plane in this node - node.polygons.par_iter_mut().for_each(|p| p.flip()); - if let Some(ref mut plane) = node.plane { - plane.flip(); - } - - // Swap front and back children - std::mem::swap(&mut node.front, &mut node.back); - - // Add children to stack for processing - if let Some(ref mut front) = node.front { - stack.push(front.as_mut()); - } - if let Some(ref mut back) = node.back { - stack.push(back.as_mut()); - } - } - } - - /// Parallel version of clip Polygons - #[cfg(feature = "parallel")] - pub fn clip_polygons(&self, polygons: &[Polygon]) -> Vec> { - // If this node has no plane, just return the original set - if self.plane.is_none() { - return polygons.to_vec(); - } - let plane = self.plane.as_ref().unwrap(); - - // Split each polygon in parallel; gather results - let (coplanar_front, coplanar_back, mut front, mut back) = polygons - .par_iter() - .map(|poly| plane.split_polygon(poly)) // <-- just pass poly - .reduce( - || (Vec::new(), Vec::new(), Vec::new(), Vec::new()), - |mut acc, x| { - acc.0.extend(x.0); - acc.1.extend(x.1); - acc.2.extend(x.2); - acc.3.extend(x.3); - acc - }, - ); - - // Decide where to send the coplanar polygons - for cp in coplanar_front { - if plane.orient_plane(&cp.plane) == FRONT { - front.push(cp); - } else { - back.push(cp); - } - } - for cp in coplanar_back { - if plane.orient_plane(&cp.plane) == FRONT { - front.push(cp); - } else { - back.push(cp); - } - } - - // Process front and back using parallel iterators to avoid recursive join - let mut result = if let Some(ref f) = self.front { - f.clip_polygons(&front) - } else { - front - }; - - if let Some(ref b) = self.back { - result.extend(b.clip_polygons(&back)); - } - // If there's no back node, we simply don't extend (effectively discarding back polygons) - - result - } - - /// Parallel version of `clip_to` using iterative approach to avoid stack overflow - #[cfg(feature = "parallel")] - pub fn clip_to(&mut self, bsp: &Node) { - // Use iterative approach with a stack to avoid recursive stack overflow - let mut stack = vec![self]; - - while let Some(node) = stack.pop() { - // Clip polygons at this node - node.polygons = bsp.clip_polygons(&node.polygons); - - // Add children to stack for processing - if let Some(ref mut front) = node.front { - stack.push(front.as_mut()); - } - if let Some(ref mut back) = node.back { - stack.push(back.as_mut()); - } - } - } - - /// Parallel version of `build`. - #[cfg(feature = "parallel")] - pub fn build(&mut self, polygons: &[Polygon]) { - if polygons.is_empty() { - return; - } - - // Choose splitting plane if not already set - if self.plane.is_none() { - self.plane = Some(self.pick_best_splitting_plane(polygons)); - } - let plane = self.plane.as_ref().unwrap(); - - // Split polygons in parallel - let (mut coplanar_front, mut coplanar_back, front, back) = - polygons.par_iter().map(|p| plane.split_polygon(p)).reduce( - || (Vec::new(), Vec::new(), Vec::new(), Vec::new()), - |mut acc, x| { - acc.0.extend(x.0); - acc.1.extend(x.1); - acc.2.extend(x.2); - acc.3.extend(x.3); - acc - }, - ); - - // Append coplanar fronts/backs to self.polygons - self.polygons.append(&mut coplanar_front); - self.polygons.append(&mut coplanar_back); - - // Build children sequentially to avoid stack overflow from recursive join - // The polygon splitting above already uses parallel iterators for the heavy work - if !front.is_empty() { - let mut front_node = self.front.take().unwrap_or_else(|| Box::new(Node::new())); - front_node.build(&front); - self.front = Some(front_node); - } - - if !back.is_empty() { - let mut back_node = self.back.take().unwrap_or_else(|| Box::new(Node::new())); - back_node.build(&back); - self.back = Some(back_node); - } - } - - // Parallel slice - #[cfg(feature = "parallel")] - pub fn slice(&self, slicing_plane: &Plane) -> (Vec>, Vec<[Vertex; 2]>) { - // Collect all polygons (this can be expensive, but let's do it). - let all_polys = self.all_polygons(); - - // Process polygons in parallel - let (coplanar_polygons, intersection_edges) = all_polys - .par_iter() - .map(|poly| { - let vcount = poly.vertices.len(); - if vcount < 2 { - // Degenerate => skip - return (Vec::new(), Vec::new()); - } - let mut polygon_type = 0; - let mut types = Vec::with_capacity(vcount); - - for vertex in &poly.vertices { - let vertex_type = slicing_plane.orient_point(&vertex.pos); - polygon_type |= vertex_type; - types.push(vertex_type); - } - - match polygon_type { - COPLANAR => { - // Entire polygon in plane - (vec![poly.clone()], Vec::new()) - }, - FRONT | BACK => { - // Entirely on one side => no intersection - (Vec::new(), Vec::new()) - }, - SPANNING => { - // The polygon crosses the plane => gather intersection edges - let mut crossing_points = Vec::new(); - for i in 0..vcount { - let j = (i + 1) % vcount; - let ti = types[i]; - let tj = types[j]; - let vi = &poly.vertices[i]; - let vj = &poly.vertices[j]; - - if (ti | tj) == SPANNING { - // The param intersection at which plane intersects the edge [vi -> vj]. - // Avoid dividing by zero: - let denom = slicing_plane.normal().dot(&(vj.pos - vi.pos)); - if denom.abs() > EPSILON { - let intersection = (slicing_plane.offset() - - slicing_plane.normal().dot(&vi.pos.coords)) - / denom; - // Interpolate: - let intersect_vert = vi.interpolate(vj, intersection); - crossing_points.push(intersect_vert); - } - } - } - - // Pair up intersection points => edges - let mut edges = Vec::new(); - for chunk in crossing_points.chunks_exact(2) { - edges.push([chunk[0].clone(), chunk[1].clone()]); - } - (Vec::new(), edges) - }, - _ => (Vec::new(), Vec::new()), - } - }) - .reduce( - || (Vec::new(), Vec::new()), - |mut acc, x| { - acc.0.extend(x.0); - acc.1.extend(x.1); - acc - }, - ); - - (coplanar_polygons, intersection_edges) - } -} diff --git a/src/mesh/connectivity.rs b/src/mesh/connectivity.rs index 7770557..e0bbb89 100644 --- a/src/mesh/connectivity.rs +++ b/src/mesh/connectivity.rs @@ -46,6 +46,16 @@ impl VertexIndexMap { new_index } + /// Get the index for a vertex position + pub fn get_index(&self, pos: &Point3) -> Option { + for (existing_pos, existing_index) in &self.position_to_index { + if (pos - existing_pos).norm() < self.epsilon { + return Some(*existing_index); + } + } + None + } + /// Get the position for a given index pub fn get_position(&self, index: usize) -> Option> { self.index_to_position.get(&index).copied() diff --git a/src/mesh/metaballs.rs b/src/mesh/metaballs.rs index 1636dd4..ef7c6bc 100644 --- a/src/mesh/metaballs.rs +++ b/src/mesh/metaballs.rs @@ -126,6 +126,7 @@ impl Mesh { } impl fast_surface_nets::ndshape::Shape<3> for GridShape { type Coord = u32; + #[inline] fn as_array(&self) -> [Self::Coord; 3] { [self.nx, self.ny, self.nz] diff --git a/src/mesh/mod.rs b/src/mesh/mod.rs index d657de3..2d71edd 100644 --- a/src/mesh/mod.rs +++ b/src/mesh/mod.rs @@ -12,7 +12,13 @@ use crate::float_types::{ }, {EPSILON, Real}, }; -use crate::mesh::{bsp::Node, plane::Plane, polygon::Polygon, vertex::Vertex}; +use crate::mesh::{ + algorithm::{convex_hull::traits::ConvexHullOps, smoothing::traits::SmoothingOps}, + bsp::Node, + plane::Plane, + polygon::Polygon, + vertex::Vertex, +}; use crate::sketch::Sketch; use crate::traits::CSG; use geo::{CoordsIter, Geometry, Polygon as GeoPolygon}; @@ -24,11 +30,9 @@ use std::{cmp::PartialEq, fmt::Debug, num::NonZeroU32, sync::OnceLock}; #[cfg(feature = "parallel")] use rayon::{iter::IntoParallelRefIterator, prelude::*}; +pub mod algorithm; pub mod bsp; -pub mod bsp_parallel; -#[cfg(feature = "chull")] -pub mod convex_hull; pub mod flatten_slice; #[cfg(feature = "metaballs")] @@ -42,7 +46,6 @@ pub mod quality; #[cfg(feature = "sdf")] pub mod sdf; pub mod shapes; -pub mod smoothing; #[cfg(feature = "sdf")] pub mod tpms; pub mod vertex; @@ -93,7 +96,7 @@ impl Mesh { mesh } - /// Split polygons into (may_touch, cannot_touch) using bounding‑box tests + /// Split polygons into (may_touch, cannot_touch) using bounding-box tests fn partition_polys( polys: &[Polygon], other_bb: &Aabb, @@ -236,7 +239,7 @@ impl Mesh { /// The angle is computed as the angle between the normal vectors of the two polygons. /// /// Returns the angle in radians. - fn dihedral_angle(p1: &Polygon, p2: &Polygon) -> Real { + pub(crate) fn dihedral_angle(p1: &Polygon, p2: &Polygon) -> Real { let n1 = p1.plane.normal(); let n2 = p2.plane.normal(); let dot = n1.dot(&n2).clamp(-1.0, 1.0); @@ -448,6 +451,87 @@ impl Mesh { mesh } + + /// Applies Laplacian smoothing to the mesh. + pub fn laplacian_smooth( + &self, + lambda: Real, + iterations: usize, + preserve_boundaries: bool, + ) -> Mesh { + #[cfg(not(feature = "parallel"))] + let ops = algorithm::smoothing::SerialSmoothingOps::new(); + #[cfg(feature = "parallel")] + let ops = algorithm::smoothing::ParallelSmoothingOps::new(); + + ops.laplacian_smooth(self, lambda, iterations, preserve_boundaries) + } + + /// Applies Taubin smoothing to the mesh. + pub fn taubin_smooth( + &self, + lambda: Real, + mu: Real, + iterations: usize, + preserve_boundaries: bool, + ) -> Mesh { + #[cfg(not(feature = "parallel"))] + let ops = algorithm::smoothing::SerialSmoothingOps::new(); + #[cfg(feature = "parallel")] + let ops = algorithm::smoothing::ParallelSmoothingOps::new(); + + ops.taubin_smooth(self, lambda, mu, iterations, preserve_boundaries) + } + + /// Adaptively refines the mesh based on geometric criteria. + pub fn adaptive_refine( + &self, + quality_threshold: Real, + max_edge_length: Real, + curvature_threshold_deg: Real, + ) -> Mesh { + #[cfg(not(feature = "parallel"))] + let ops = algorithm::smoothing::SerialSmoothingOps::new(); + #[cfg(feature = "parallel")] + let ops = algorithm::smoothing::ParallelSmoothingOps::new(); + + ops.adaptive_refine( + self, + quality_threshold, + max_edge_length, + curvature_threshold_deg, + ) + } + + /// Removes poor-quality triangles from the mesh. + pub fn remove_poor_triangles(&self, min_quality: Real) -> Mesh { + #[cfg(not(feature = "parallel"))] + let ops = algorithm::smoothing::SerialSmoothingOps::new(); + #[cfg(feature = "parallel")] + let ops = algorithm::smoothing::ParallelSmoothingOps::new(); + + ops.remove_poor_triangles(self, min_quality) + } + + /// Computes the convex hull of the mesh. + pub fn convex_hull(&self) -> Mesh { + #[cfg(not(feature = "parallel"))] + let ops = algorithm::convex_hull::SerialConvexHullOps::new(); + #[cfg(feature = "parallel")] + let ops = algorithm::convex_hull::ParallelConvexHullOps::new(); + + ops.convex_hull(self) + } + + /// Computes the Minkowski sum of two meshes. + pub fn minkowski_sum(&self, other: &Mesh) -> Mesh { + #[cfg(not(feature = "parallel"))] + let ops = algorithm::convex_hull::SerialConvexHullOps::new(); + #[cfg(feature = "parallel")] + let ops = algorithm::convex_hull::ParallelConvexHullOps::new(); + + ops.minkowski_sum(self, other) + } } impl CSG for Mesh { @@ -474,7 +558,7 @@ impl CSG for Mesh { /// +-------+ +-------+ /// ``` fn union(&self, other: &Mesh) -> Mesh { - // avoid splitting obvious non‑intersecting faces + // avoid splitting obvious non-intersecting faces let (a_clip, a_passthru) = Self::partition_polys(&self.polygons, &other.bounding_box()); let (b_clip, b_passthru) = @@ -516,7 +600,7 @@ impl CSG for Mesh { /// +-------+ /// ``` fn difference(&self, other: &Mesh) -> Mesh { - // avoid splitting obvious non‑intersecting faces + // avoid splitting obvious non-intersecting faces let (a_clip, a_passthru) = Self::partition_polys(&self.polygons, &other.bounding_box()); let (b_clip, _b_passthru) = diff --git a/src/mesh/sdf.rs b/src/mesh/sdf.rs deleted file mode 100644 index 7d20d8e..0000000 --- a/src/mesh/sdf.rs +++ /dev/null @@ -1,218 +0,0 @@ -//! Create `Mesh`s by meshing signed distance fields ([sdf](https://en.wikipedia.org/wiki/Signed_distance_function)) within a bounding box. - -use crate::float_types::Real; -use crate::mesh::Mesh; -use crate::mesh::polygon::Polygon; -use crate::mesh::vertex::Vertex; -use fast_surface_nets::{SurfaceNetsBuffer, surface_nets}; -use nalgebra::{Point3, Vector3}; -use std::fmt::Debug; - -impl Mesh { - /// Return a Mesh created by meshing a signed distance field within a bounding box - /// - /// ``` - /// # use csgrs::{mesh::Mesh, float_types::Real}; - /// # use nalgebra::Point3; - /// // Example SDF for a sphere of radius 1.5 centered at (0,0,0) - /// let my_sdf = |p: &Point3| p.coords.norm() - 1.5; - /// - /// let resolution = (60, 60, 60); - /// let min_pt = Point3::new(-2.0, -2.0, -2.0); - /// let max_pt = Point3::new( 2.0, 2.0, 2.0); - /// let iso_value = 0.0; // Typically zero for SDF-based surfaces - /// - /// let mesh_shape = Mesh::<()>::sdf(my_sdf, resolution, min_pt, max_pt, iso_value, None); - /// - /// // Now `mesh_shape` is your polygon mesh as a Mesh you can union, subtract, or export: - /// let _ = std::fs::write("stl/sdf_sphere.stl", mesh_shape.to_stl_binary("sdf_sphere").unwrap()); - pub fn sdf( - sdf: F, - resolution: (usize, usize, usize), - min_pt: Point3, - max_pt: Point3, - iso_value: Real, - metadata: Option, - ) -> Mesh - where - // F is a closure or function that takes a 3D point and returns the signed distance. - // Must be `Sync`/`Send` if you want to parallelize the sampling. - F: Fn(&Point3) -> Real + Sync + Send, - { - // Early return if resolution is degenerate - let nx = resolution.0.max(2) as u32; - let ny = resolution.1.max(2) as u32; - let nz = resolution.2.max(2) as u32; - - // Determine grid spacing based on bounding box and resolution - let dx = (max_pt.x - min_pt.x) / (nx as Real - 1.0); - let dy = (max_pt.y - min_pt.y) / (ny as Real - 1.0); - let dz = (max_pt.z - min_pt.z) / (nz as Real - 1.0); - - // Allocate storage for field values: - let array_size = (nx * ny * nz) as usize; - let mut field_values = vec![0.0_f32; array_size]; - - // Optimized finite value checking with iterator patterns - // **Mathematical Foundation**: Ensures all coordinates are finite real numbers - #[inline] - fn point_finite(p: &Point3) -> bool { - p.coords.iter().all(|&c| c.is_finite()) - } - - #[inline] - fn vec_finite(v: &Vector3) -> bool { - v.iter().all(|&c| c.is_finite()) - } - - // Sample the SDF at each grid cell with optimized iteration pattern: - // **Mathematical Foundation**: For SDF f(p), we sample at regular intervals - // and store (f(p) - iso_value) so surface_nets finds zero-crossings at iso_value. - // **Optimization**: Linear memory access pattern with better cache locality. - for i in 0..(nx * ny * nz) { - let iz = i / (nx * ny); - let remainder = i % (nx * ny); - let iy = remainder / nx; - let ix = remainder % nx; - - let xf = min_pt.x + (ix as Real) * dx; - let yf = min_pt.y + (iy as Real) * dy; - let zf = min_pt.z + (iz as Real) * dz; - - let p = Point3::new(xf, yf, zf); - let sdf_val = sdf(&p); - - // Robust finite value handling with mathematical correctness - field_values[i as usize] = if sdf_val.is_finite() { - (sdf_val - iso_value) as f32 - } else { - // For infinite/NaN values, use large positive value to indicate "far outside" - // This preserves the mathematical properties of the distance field - 1e10_f32 - }; - } - - // The shape describing our discrete grid for Surface Nets: - #[derive(Clone, Copy)] - struct GridShape { - nx: u32, - ny: u32, - nz: u32, - } - - impl fast_surface_nets::ndshape::Shape<3> for GridShape { - type Coord = u32; - - #[inline] - fn as_array(&self) -> [Self::Coord; 3] { - [self.nx, self.ny, self.nz] - } - - fn size(&self) -> Self::Coord { - self.nx * self.ny * self.nz - } - - fn usize(&self) -> usize { - (self.nx * self.ny * self.nz) as usize - } - - fn linearize(&self, coords: [Self::Coord; 3]) -> u32 { - let [x, y, z] = coords; - (z * self.ny + y) * self.nx + x - } - - fn delinearize(&self, i: u32) -> [Self::Coord; 3] { - let x = i % self.nx; - let yz = i / self.nx; - let y = yz % self.ny; - let z = yz / self.ny; - [x, y, z] - } - } - - let shape = GridShape { nx, ny, nz }; - - // `SurfaceNetsBuffer` collects the positions, normals, and triangle indices - let mut sn_buffer = SurfaceNetsBuffer::default(); - - // The max valid coordinate in each dimension - let max_x = nx - 1; - let max_y = ny - 1; - let max_z = nz - 1; - - // Run surface nets - surface_nets( - &field_values, - &shape, - [0, 0, 0], - [max_x, max_y, max_z], - &mut sn_buffer, - ); - - // Convert the resulting triangles into Mesh polygons - let mut triangles = Vec::with_capacity(sn_buffer.indices.len() / 3); - - for tri in sn_buffer.indices.chunks_exact(3) { - let i0 = tri[0] as usize; - let i1 = tri[1] as usize; - let i2 = tri[2] as usize; - - let p0i = sn_buffer.positions[i0]; - let p1i = sn_buffer.positions[i1]; - let p2i = sn_buffer.positions[i2]; - - // Convert from [u32; 3] to real coordinates: - let p0 = Point3::new( - min_pt.x + p0i[0] as Real * dx, - min_pt.y + p0i[1] as Real * dy, - min_pt.z + p0i[2] as Real * dz, - ); - let p1 = Point3::new( - min_pt.x + p1i[0] as Real * dx, - min_pt.y + p1i[1] as Real * dy, - min_pt.z + p1i[2] as Real * dz, - ); - let p2 = Point3::new( - min_pt.x + p2i[0] as Real * dx, - min_pt.y + p2i[1] as Real * dy, - min_pt.z + p2i[2] as Real * dz, - ); - - // Retrieve precomputed normal from Surface Nets: - let n0 = sn_buffer.normals[i0]; - let n1 = sn_buffer.normals[i1]; - let n2 = sn_buffer.normals[i2]; - - // Normals come out as [f32;3] – promote to `Real` - let n0v = Vector3::new(n0[0] as Real, n0[1] as Real, n0[2] as Real); - let n1v = Vector3::new(n1[0] as Real, n1[1] as Real, n1[2] as Real); - let n2v = Vector3::new(n2[0] as Real, n2[1] as Real, n2[2] as Real); - - // ── « gate » ──────────────────────────────────────────────── - if !(point_finite(&p0) - && point_finite(&p1) - && point_finite(&p2) - && vec_finite(&n0v) - && vec_finite(&n1v) - && vec_finite(&n2v)) - { - // at least one coordinate was NaN/±∞ – ignore this triangle - continue; - } - - let v0 = - Vertex::new(p0, Vector3::new(n0[0] as Real, n0[1] as Real, n0[2] as Real)); - let v1 = - Vertex::new(p1, Vector3::new(n1[0] as Real, n1[1] as Real, n1[2] as Real)); - let v2 = - Vertex::new(p2, Vector3::new(n2[0] as Real, n2[1] as Real, n2[2] as Real)); - - // Note: reverse v1, v2 if you need to fix winding - let poly = Polygon::new(vec![v0, v1, v2], metadata.clone()); - triangles.push(poly); - } - - // Return as a Mesh - Mesh::from_polygons(&triangles, metadata) - } -} diff --git a/src/mesh/sdf/grid.rs b/src/mesh/sdf/grid.rs new file mode 100644 index 0000000..794938d --- /dev/null +++ b/src/mesh/sdf/grid.rs @@ -0,0 +1,39 @@ +//! Grid shape definition for surface nets + +/// The shape describing our discrete grid for Surface Nets: +#[derive(Clone, Copy)] +pub struct GridShape { + pub nx: u32, + pub ny: u32, + pub nz: u32, +} + +impl fast_surface_nets::ndshape::Shape<3> for GridShape { + type Coord = u32; + + #[inline] + fn as_array(&self) -> [Self::Coord; 3] { + [self.nx, self.ny, self.nz] + } + + fn size(&self) -> Self::Coord { + self.nx * self.ny * self.nz + } + + fn usize(&self) -> usize { + (self.nx * self.ny * self.nz) as usize + } + + fn linearize(&self, coords: [Self::Coord; 3]) -> u32 { + let [x, y, z] = coords; + (z * self.ny + y) * self.nx + x + } + + fn delinearize(&self, i: u32) -> [Self::Coord; 3] { + let x = i % self.nx; + let yz = i / self.nx; + let y = yz % self.ny; + let z = yz / self.ny; + [x, y, z] + } +} diff --git a/src/mesh/sdf/mod.rs b/src/mesh/sdf/mod.rs new file mode 100644 index 0000000..ff7cee1 --- /dev/null +++ b/src/mesh/sdf/mod.rs @@ -0,0 +1,50 @@ +//! Signed Distance Field (SDF) meshing +//! +//! This module provides SDF meshing operations with dependency inversion, +//! allowing for different algorithm implementations (serial/parallel). + +pub mod grid; +pub mod traits; + +#[cfg(not(feature = "parallel"))] +pub mod serial; + +#[cfg(feature = "parallel")] +pub mod parallel; + +// Re-export core types +pub use grid::GridShape; +pub use traits::SdfOps; + +#[cfg(not(feature = "parallel"))] +pub use serial::SerialSdfOps; + +#[cfg(feature = "parallel")] +pub use parallel::ParallelSdfOps; + +use crate::float_types::Real; +use crate::mesh::Mesh; +use nalgebra::Point3; +use std::fmt::Debug; + +impl Mesh { + /// Return a Mesh created by meshing a signed distance field within a bounding box + pub fn sdf( + sdf: F, + resolution: (usize, usize, usize), + min_pt: Point3, + max_pt: Point3, + iso_value: Real, + metadata: Option, + ) -> Mesh + where + F: Fn(&Point3) -> Real + Sync + Send, + { + #[cfg(not(feature = "parallel"))] + let ops = SerialSdfOps::new(); + #[cfg(feature = "parallel")] + let ops = ParallelSdfOps::new(); + + ops.mesh(sdf, resolution, min_pt, max_pt, iso_value, metadata) + } +} diff --git a/src/mesh/sdf/parallel.rs b/src/mesh/sdf/parallel.rs new file mode 100644 index 0000000..94d9f52 --- /dev/null +++ b/src/mesh/sdf/parallel.rs @@ -0,0 +1,161 @@ +//! Parallel implementation of SDF meshing + +use crate::float_types::Real; +use crate::mesh::Mesh; +use crate::mesh::polygon::Polygon; +use crate::mesh::sdf::grid::GridShape; +use crate::mesh::sdf::traits::SdfOps; +use crate::mesh::vertex::Vertex; +use fast_surface_nets::{SurfaceNetsBuffer, surface_nets}; +use nalgebra::{Point3, Vector3}; +use rayon::prelude::*; +use std::fmt::Debug; + +/// Parallel implementation of SDF meshing +pub struct ParallelSdfOps; + +impl ParallelSdfOps { + pub fn new() -> Self { + Self + } +} + +impl Default for ParallelSdfOps { + fn default() -> Self { + Self::new() + } +} + +impl SdfOps for ParallelSdfOps { + fn mesh( + &self, + sdf: F, + resolution: (usize, usize, usize), + min_pt: Point3, + max_pt: Point3, + iso_value: Real, + metadata: Option, + ) -> Mesh + where + F: Fn(&Point3) -> Real + Sync + Send, + { + // Early return if resolution is degenerate + let nx = resolution.0.max(2) as u32; + let ny = resolution.1.max(2) as u32; + let nz = resolution.2.max(2) as u32; + + // Determine grid spacing based on bounding box and resolution + let dx = (max_pt.x - min_pt.x) / (nx as Real - 1.0); + let dy = (max_pt.y - min_pt.y) / (ny as Real - 1.0); + let dz = (max_pt.z - min_pt.z) / (nz as Real - 1.0); + + // Allocate storage for field values: + let array_size = (nx * ny * nz) as usize; + let mut field_values = vec![0.0_f32; array_size]; + + // Sample the SDF at each grid cell in parallel + field_values + .par_iter_mut() + .enumerate() + .for_each(|(i, value)| { + let iz = i / (nx * ny) as usize; + let remainder = i % (nx * ny) as usize; + let iy = remainder / nx as usize; + let ix = remainder % nx as usize; + + let xf = min_pt.x + (ix as Real) * dx; + let yf = min_pt.y + (iy as Real) * dy; + let zf = min_pt.z + (iz as Real) * dz; + + let p = Point3::new(xf, yf, zf); + let sdf_val = sdf(&p); + + *value = if sdf_val.is_finite() { + (sdf_val - iso_value) as f32 + } else { + 1e10_f32 + }; + }); + + // Optimized finite value checking with iterator patterns + #[inline] + fn point_finite(p: &Point3) -> bool { + p.coords.iter().all(|&c| c.is_finite()) + } + + #[inline] + fn vec_finite(v: &Vector3) -> bool { + v.iter().all(|&c| c.is_finite()) + } + + let shape = GridShape { nx, ny, nz }; + let mut sn_buffer = SurfaceNetsBuffer::default(); + let max_x = nx - 1; + let max_y = ny - 1; + let max_z = nz - 1; + + surface_nets( + &field_values, + &shape, + [0, 0, 0], + [max_x, max_y, max_z], + &mut sn_buffer, + ); + + let triangles: Vec> = sn_buffer + .indices + .par_chunks_exact(3) + .filter_map(|tri| { + let i0 = tri[0] as usize; + let i1 = tri[1] as usize; + let i2 = tri[2] as usize; + + let p0i = sn_buffer.positions[i0]; + let p1i = sn_buffer.positions[i1]; + let p2i = sn_buffer.positions[i2]; + + let p0 = Point3::new( + min_pt.x + p0i[0] as Real * dx, + min_pt.y + p0i[1] as Real * dy, + min_pt.z + p0i[2] as Real * dz, + ); + let p1 = Point3::new( + min_pt.x + p1i[0] as Real * dx, + min_pt.y + p1i[1] as Real * dy, + min_pt.z + p1i[2] as Real * dz, + ); + let p2 = Point3::new( + min_pt.x + p2i[0] as Real * dx, + min_pt.y + p2i[1] as Real * dy, + min_pt.z + p2i[2] as Real * dz, + ); + + let n0 = sn_buffer.normals[i0]; + let n1 = sn_buffer.normals[i1]; + let n2 = sn_buffer.normals[i2]; + + let n0v = Vector3::new(n0[0] as Real, n0[1] as Real, n0[2] as Real); + let n1v = Vector3::new(n1[0] as Real, n1[1] as Real, n1[2] as Real); + let n2v = Vector3::new(n2[0] as Real, n2[1] as Real, n2[2] as Real); + + if !(point_finite(&p0) + && point_finite(&p1) + && point_finite(&p2) + && vec_finite(&n0v) + && vec_finite(&n1v) + && vec_finite(&n2v)) + { + return None; + } + + let v0 = Vertex::new(p0, n0v); + let v1 = Vertex::new(p1, n1v); + let v2 = Vertex::new(p2, n2v); + + Some(Polygon::new(vec![v0, v1, v2], metadata.clone())) + }) + .collect(); + + Mesh::from_polygons(&triangles, metadata) + } +} diff --git a/src/mesh/sdf/serial.rs b/src/mesh/sdf/serial.rs new file mode 100644 index 0000000..139960a --- /dev/null +++ b/src/mesh/sdf/serial.rs @@ -0,0 +1,157 @@ +//! Serial implementation of SDF meshing + +use crate::float_types::Real; +use crate::mesh::Mesh; +use crate::mesh::polygon::Polygon; +use crate::mesh::sdf::grid::GridShape; +use crate::mesh::sdf::traits::SdfOps; +use crate::mesh::vertex::Vertex; +use fast_surface_nets::{SurfaceNetsBuffer, surface_nets}; +use nalgebra::{Point3, Vector3}; +use std::fmt::Debug; + +/// Serial implementation of SDF meshing +pub struct SerialSdfOps; + +impl SerialSdfOps { + pub const fn new() -> Self { + Self + } +} + +impl Default for SerialSdfOps { + fn default() -> Self { + Self::new() + } +} + +impl SdfOps for SerialSdfOps { + fn mesh( + &self, + sdf: F, + resolution: (usize, usize, usize), + min_pt: Point3, + max_pt: Point3, + iso_value: Real, + metadata: Option, + ) -> Mesh + where + F: Fn(&Point3) -> Real + Sync + Send, + { + // Early return if resolution is degenerate + let nx = resolution.0.max(2) as u32; + let ny = resolution.1.max(2) as u32; + let nz = resolution.2.max(2) as u32; + + // Determine grid spacing based on bounding box and resolution + let dx = (max_pt.x - min_pt.x) / (nx as Real - 1.0); + let dy = (max_pt.y - min_pt.y) / (ny as Real - 1.0); + let dz = (max_pt.z - min_pt.z) / (nz as Real - 1.0); + + // Allocate storage for field values: + let array_size = (nx * ny * nz) as usize; + let mut field_values = vec![0.0_f32; array_size]; + + // Optimized finite value checking with iterator patterns + #[inline] + fn point_finite(p: &Point3) -> bool { + p.coords.iter().all(|&c| c.is_finite()) + } + + #[inline] + fn vec_finite(v: &Vector3) -> bool { + v.iter().all(|&c| c.is_finite()) + } + + // Sample the SDF at each grid cell + field_values.iter_mut().enumerate().for_each(|(i, value)| { + let iz = i / (nx * ny) as usize; + let remainder = i % (nx * ny) as usize; + let iy = remainder / nx as usize; + let ix = remainder % nx as usize; + + let xf = min_pt.x + (ix as Real) * dx; + let yf = min_pt.y + (iy as Real) * dy; + let zf = min_pt.z + (iz as Real) * dz; + + let p = Point3::new(xf, yf, zf); + let sdf_val = sdf(&p); + + *value = if sdf_val.is_finite() { + (sdf_val - iso_value) as f32 + } else { + 1e10_f32 + }; + }); + + let shape = GridShape { nx, ny, nz }; + let mut sn_buffer = SurfaceNetsBuffer::default(); + let max_x = nx - 1; + let max_y = ny - 1; + let max_z = nz - 1; + + surface_nets( + &field_values, + &shape, + [0, 0, 0], + [max_x, max_y, max_z], + &mut sn_buffer, + ); + + let triangles: Vec> = sn_buffer + .indices + .chunks_exact(3) + .filter_map(|tri| { + let i0 = tri[0] as usize; + let i1 = tri[1] as usize; + let i2 = tri[2] as usize; + + let p0i = sn_buffer.positions[i0]; + let p1i = sn_buffer.positions[i1]; + let p2i = sn_buffer.positions[i2]; + + let p0 = Point3::new( + min_pt.x + p0i[0] as Real * dx, + min_pt.y + p0i[1] as Real * dy, + min_pt.z + p0i[2] as Real * dz, + ); + let p1 = Point3::new( + min_pt.x + p1i[0] as Real * dx, + min_pt.y + p1i[1] as Real * dy, + min_pt.z + p1i[2] as Real * dz, + ); + let p2 = Point3::new( + min_pt.x + p2i[0] as Real * dx, + min_pt.y + p2i[1] as Real * dy, + min_pt.z + p2i[2] as Real * dz, + ); + + let n0 = sn_buffer.normals[i0]; + let n1 = sn_buffer.normals[i1]; + let n2 = sn_buffer.normals[i2]; + + let n0v = Vector3::new(n0[0] as Real, n0[1] as Real, n0[2] as Real); + let n1v = Vector3::new(n1[0] as Real, n1[1] as Real, n1[2] as Real); + let n2v = Vector3::new(n2[0] as Real, n2[1] as Real, n2[2] as Real); + + if !(point_finite(&p0) + && point_finite(&p1) + && point_finite(&p2) + && vec_finite(&n0v) + && vec_finite(&n1v) + && vec_finite(&n2v)) + { + return None; + } + + let v0 = Vertex::new(p0, n0v); + let v1 = Vertex::new(p1, n1v); + let v2 = Vertex::new(p2, n2v); + + Some(Polygon::new(vec![v0, v1, v2], metadata.clone())) + }) + .collect(); + + Mesh::from_polygons(&triangles, metadata) + } +} diff --git a/src/mesh/sdf/traits.rs b/src/mesh/sdf/traits.rs new file mode 100644 index 0000000..8dc1d79 --- /dev/null +++ b/src/mesh/sdf/traits.rs @@ -0,0 +1,22 @@ +//! Traits defining SDF meshing operations for dependency inversion + +use crate::float_types::Real; +use crate::mesh::Mesh; +use nalgebra::Point3; +use std::fmt::Debug; + +/// Core SDF meshing operations trait +pub trait SdfOps { + /// Create a mesh from an SDF + fn mesh( + &self, + sdf: F, + resolution: (usize, usize, usize), + min_pt: Point3, + max_pt: Point3, + iso_value: Real, + metadata: Option, + ) -> Mesh + where + F: Fn(&Point3) -> Real + Sync + Send; +} diff --git a/src/nurbs/mod.rs b/src/nurbs/mod.rs index cd79ddd..4942a7d 100644 --- a/src/nurbs/mod.rs +++ b/src/nurbs/mod.rs @@ -1 +1 @@ -//pub mod nurbs; +// pub mod nurbs; diff --git a/src/sketch/extrudes.rs b/src/sketch/extrudes.rs index 0dc41cc..c30bb29 100644 --- a/src/sketch/extrudes.rs +++ b/src/sketch/extrudes.rs @@ -306,179 +306,177 @@ impl Sketch { Ok(Mesh::from_polygons(&polygons, bottom.metadata.clone())) } - /* - /// Perform a linear extrusion along some axis, with optional twist, center, slices, scale, etc. - /// - /// # Parameters - /// - `direction`: Direction vector for the extrusion. - /// - `twist`: Total twist in degrees around the extrusion axis from bottom to top. - /// - `segments`: Number of intermediate subdivisions. - /// - `scale`: A uniform scale factor to apply at the top slice (bottom is scale=1.0). - /// - /// # Assumptions - /// - This CSG is assumed to represent one or more 2D polygons lying in or near the XY plane. - /// - The resulting shape is extruded *initially* along +Z, then finally rotated if `v != [0,0,1]`. - /// - /// # Returns - /// A new 3D CSG. - /// - /// # Example - /// ``` - /// let shape_2d = CSG::square(2.0, None); // a 2D square in XY - /// let extruded = shape_2d.linear_extrude( - /// direction = Vector3::new(0.0, 0.0, 10.0), - /// twist = 360.0, - /// segments = 32, - /// scale = 1.2, - /// ); - /// ``` - pub fn linear_extrude( - shape: &CCShape, - direction: Vector3, - twist_degs: Real, - segments: usize, - scale_top: Real, - metadata: Option, - ) -> CSG { - let mut polygons_3d = Vec::new(); - if segments < 1 { - return CSG::new(); - } - let height = direction.norm(); - if height < EPSILON { - // no real extrusion - return CSG::new(); - } - - // Step 1) Build a series of “transforms” from bottom=0..top=height, subdivided into `segments`. - // For each i in [0..=segments], compute fraction f and: - // - scale in XY => s_i - // - twist about Z => rot_i - // - translate in Z => z_i - // - // We'll store each “slice” in 3D form as a Vec>>, - // i.e. one 3D polyline for each boundary or hole in the shape. - let mut slices: Vec>>> = Vec::with_capacity(segments + 1); - // The axis to rotate around is the unit of `direction`. We'll do final alignment after constructing them along +Z. - let axis_dir = direction.normalize(); - - for i in 0..=segments { - let f = i as Real / segments as Real; - let s_i = 1.0 + (scale_top - 1.0) * f; // lerp(1, scale_top, f) - let twist_rad = twist_degs.to_radians() * f; - let z_i = height * f; - - // Build transform T = Tz * Rz * Sxy - // - scale in XY - // - twist around Z - // - translate in Z - let mat_scale = Matrix4::new_nonuniform_scaling(&Vector3::new(s_i, s_i, 1.0)); - let mat_rot = Rotation3::from_axis_angle(&Vector3::z_axis(), twist_rad).to_homogeneous(); - let mat_trans = Translation3::new(0.0, 0.0, z_i).to_homogeneous(); - let slice_mat = mat_trans * mat_rot * mat_scale; - - let slice_3d = project_shape_3d(shape, &slice_mat); - slices.push(slice_3d); - } - - // Step 2) “Stitch” consecutive slices to form side polygons. - // For each pair of slices[i], slices[i+1], for each boundary polyline j, - // connect edges. We assume each polyline has the same vertex_count in both slices. - // (If the shape is closed, we do wrap edges [n..0].) - // Then we optionally build bottom & top caps if the polylines are closed. - - // a) bottom + top caps, similar to extrude_vector approach - // For slices[0], build a “bottom” by triangulating in XY, flipping normal. - // For slices[segments], build a “top” by normal up. - // - // But we only do it if each boundary is closed. - // We must group CCW with matching holes. This is the same logic as `extrude_vector`. - - // We'll do a small helper that triangulates shape in 2D, then lifts that triangulation to slice_3d. - // You can re‐use the logic from `extrude_vector`. - - // Build the “bottom” from slices[0] if polylines are all or partially closed - polygons_3d.extend( - build_caps_from_slice(shape, &slices[0], true, metadata.clone()) - ); - // Build the “top” from slices[segments] - polygons_3d.extend( - build_caps_from_slice(shape, &slices[segments], false, metadata.clone()) - ); - - // b) side walls - for i in 0..segments { - let bottom_slice = &slices[i]; - let top_slice = &slices[i + 1]; - - // We know bottom_slice has shape.ccw_plines.len() + shape.cw_plines.len() polylines - // in the same order. Each polyline has the same vertex_count as in top_slice. - // So we can do a direct 1:1 match: bottom_slice[j] <-> top_slice[j]. - for (pline_idx, bot3d) in bottom_slice.iter().enumerate() { - let top3d = &top_slice[pline_idx]; - if bot3d.len() < 2 { - continue; - } - // is it closed? We can check shape’s corresponding polyline - let is_closed = if pline_idx < shape.ccw_plines.len() { - shape.ccw_plines[pline_idx].polyline.is_closed() - } else { - shape.cw_plines[pline_idx - shape.ccw_plines.len()].polyline.is_closed() - }; - let n = bot3d.len(); - let edge_count = if is_closed { n } else { n - 1 }; - - for k in 0..edge_count { - let k_next = (k + 1) % n; - let b_i = bot3d[k]; - let b_j = bot3d[k_next]; - let t_i = top3d[k]; - let t_j = top3d[k_next]; - - let poly_side = Polygon::new( - vec![ - Vertex::new(b_i, Vector3::zeros()), - Vertex::new(b_j, Vector3::zeros()), - Vertex::new(t_j, Vector3::zeros()), - Vertex::new(t_i, Vector3::zeros()), - ], - metadata.clone(), - ); - polygons_3d.push(poly_side); - } - } - } - - // Step 3) If direction is not along +Z, rotate final mesh so +Z aligns with your direction - // (This is optional or can be done up front. Typical OpenSCAD style is to do everything - // along +Z, then rotate the final.) - if (axis_dir - Vector3::z()).norm() > EPSILON { - // rotate from +Z to axis_dir - let rot_axis = Vector3::z().cross(&axis_dir); - let sin_theta = rot_axis.norm(); - if sin_theta > EPSILON { - let cos_theta = Vector3::z().dot(&axis_dir); - let angle = cos_theta.acos(); - let rot = Rotation3::from_axis_angle(&Unit::new_normalize(rot_axis), angle); - let mat = rot.to_homogeneous(); - // transform the polygons - let mut final_polys = Vec::with_capacity(polygons_3d.len()); - for mut poly in polygons_3d { - for v in &mut poly.vertices { - let pos4 = mat * nalgebra::Vector4::new(v.pos.x, v.pos.y, v.pos.z, 1.0); - v.pos = Point3::new(pos4.x / pos4.w, pos4.y / pos4.w, pos4.z / pos4.w); - } - poly.set_new_normal(); - final_polys.push(poly); - } - return CSG::from_polygons(&final_polys); - } - } - - // otherwise, just return as is - CSG::from_polygons(&polygons_3d) - } - */ + // Perform a linear extrusion along some axis, with optional twist, center, slices, scale, etc. + // + // # Parameters + // - `direction`: Direction vector for the extrusion. + // - `twist`: Total twist in degrees around the extrusion axis from bottom to top. + // - `segments`: Number of intermediate subdivisions. + // - `scale`: A uniform scale factor to apply at the top slice (bottom is scale=1.0). + // + // # Assumptions + // - This CSG is assumed to represent one or more 2D polygons lying in or near the XY plane. + // - The resulting shape is extruded *initially* along +Z, then finally rotated if `v != [0,0,1]`. + // + // # Returns + // A new 3D CSG. + // + // # Example + // ``` + // let shape_2d = CSG::square(2.0, None); // a 2D square in XY + // let extruded = shape_2d.linear_extrude( + // direction = Vector3::new(0.0, 0.0, 10.0), + // twist = 360.0, + // segments = 32, + // scale = 1.2, + // ); + // ``` + // pub fn linear_extrude( + // shape: &CCShape, + // direction: Vector3, + // twist_degs: Real, + // segments: usize, + // scale_top: Real, + // metadata: Option, + // ) -> CSG { + // let mut polygons_3d = Vec::new(); + // if segments < 1 { + // return CSG::new(); + // } + // let height = direction.norm(); + // if height < EPSILON { + // no real extrusion + // return CSG::new(); + // } + // + // Step 1) Build a series of “transforms” from bottom=0..top=height, subdivided into `segments`. + // For each i in [0..=segments], compute fraction f and: + // - scale in XY => s_i + // - twist about Z => rot_i + // - translate in Z => z_i + // + // We'll store each “slice” in 3D form as a Vec>>, + // i.e. one 3D polyline for each boundary or hole in the shape. + // let mut slices: Vec>>> = Vec::with_capacity(segments + 1); + // The axis to rotate around is the unit of `direction`. We'll do final alignment after constructing them along +Z. + // let axis_dir = direction.normalize(); + // + // for i in 0..=segments { + // let f = i as Real / segments as Real; + // let s_i = 1.0 + (scale_top - 1.0) * f; // lerp(1, scale_top, f) + // let twist_rad = twist_degs.to_radians() * f; + // let z_i = height * f; + // + // Build transform T = Tz * Rz * Sxy + // - scale in XY + // - twist around Z + // - translate in Z + // let mat_scale = Matrix4::new_nonuniform_scaling(&Vector3::new(s_i, s_i, 1.0)); + // let mat_rot = Rotation3::from_axis_angle(&Vector3::z_axis(), twist_rad).to_homogeneous(); + // let mat_trans = Translation3::new(0.0, 0.0, z_i).to_homogeneous(); + // let slice_mat = mat_trans * mat_rot * mat_scale; + // + // let slice_3d = project_shape_3d(shape, &slice_mat); + // slices.push(slice_3d); + // } + // + // Step 2) “Stitch” consecutive slices to form side polygons. + // For each pair of slices[i], slices[i+1], for each boundary polyline j, + // connect edges. We assume each polyline has the same vertex_count in both slices. + // (If the shape is closed, we do wrap edges [n..0].) + // Then we optionally build bottom & top caps if the polylines are closed. + // + // a) bottom + top caps, similar to extrude_vector approach + // For slices[0], build a “bottom” by triangulating in XY, flipping normal. + // For slices[segments], build a “top” by normal up. + // + // But we only do it if each boundary is closed. + // We must group CCW with matching holes. This is the same logic as `extrude_vector`. + // + // We'll do a small helper that triangulates shape in 2D, then lifts that triangulation to slice_3d. + // You can re‐use the logic from `extrude_vector`. + // + // Build the “bottom” from slices[0] if polylines are all or partially closed + // polygons_3d.extend( + // build_caps_from_slice(shape, &slices[0], true, metadata.clone()) + // ); + // Build the “top” from slices[segments] + // polygons_3d.extend( + // build_caps_from_slice(shape, &slices[segments], false, metadata.clone()) + // ); + // + // b) side walls + // for i in 0..segments { + // let bottom_slice = &slices[i]; + // let top_slice = &slices[i + 1]; + // + // We know bottom_slice has shape.ccw_plines.len() + shape.cw_plines.len() polylines + // in the same order. Each polyline has the same vertex_count as in top_slice. + // So we can do a direct 1:1 match: bottom_slice[j] <-> top_slice[j]. + // for (pline_idx, bot3d) in bottom_slice.iter().enumerate() { + // let top3d = &top_slice[pline_idx]; + // if bot3d.len() < 2 { + // continue; + // } + // is it closed? We can check shape’s corresponding polyline + // let is_closed = if pline_idx < shape.ccw_plines.len() { + // shape.ccw_plines[pline_idx].polyline.is_closed() + // } else { + // shape.cw_plines[pline_idx - shape.ccw_plines.len()].polyline.is_closed() + // }; + // let n = bot3d.len(); + // let edge_count = if is_closed { n } else { n - 1 }; + // + // for k in 0..edge_count { + // let k_next = (k + 1) % n; + // let b_i = bot3d[k]; + // let b_j = bot3d[k_next]; + // let t_i = top3d[k]; + // let t_j = top3d[k_next]; + // + // let poly_side = Polygon::new( + // vec![ + // Vertex::new(b_i, Vector3::zeros()), + // Vertex::new(b_j, Vector3::zeros()), + // Vertex::new(t_j, Vector3::zeros()), + // Vertex::new(t_i, Vector3::zeros()), + // ], + // metadata.clone(), + // ); + // polygons_3d.push(poly_side); + // } + // } + // } + // + // Step 3) If direction is not along +Z, rotate final mesh so +Z aligns with your direction + // (This is optional or can be done up front. Typical OpenSCAD style is to do everything + // along +Z, then rotate the final.) + // if (axis_dir - Vector3::z()).norm() > EPSILON { + // rotate from +Z to axis_dir + // let rot_axis = Vector3::z().cross(&axis_dir); + // let sin_theta = rot_axis.norm(); + // if sin_theta > EPSILON { + // let cos_theta = Vector3::z().dot(&axis_dir); + // let angle = cos_theta.acos(); + // let rot = Rotation3::from_axis_angle(&Unit::new_normalize(rot_axis), angle); + // let mat = rot.to_homogeneous(); + // transform the polygons + // let mut final_polys = Vec::with_capacity(polygons_3d.len()); + // for mut poly in polygons_3d { + // for v in &mut poly.vertices { + // let pos4 = mat * nalgebra::Vector4::new(v.pos.x, v.pos.y, v.pos.z, 1.0); + // v.pos = Point3::new(pos4.x / pos4.w, pos4.y / pos4.w, pos4.z / pos4.w); + // } + // poly.set_new_normal(); + // final_polys.push(poly); + // } + // return CSG::from_polygons(&final_polys); + // } + // } + // + // otherwise, just return as is + // CSG::from_polygons(&polygons_3d) + // } /// **Mathematical Foundation: Surface of Revolution Generation** /// @@ -826,163 +824,161 @@ impl Sketch { // // # Returns // A new 3D `CSG` that is the swept volume. - /* - pub fn sweep(shape_2d: &Polygon, path_2d: &Polygon) -> CSG { - // Gather the path’s vertices in XY - if path_2d.vertices.len() < 2 { - // Degenerate path => no sweep - return CSG::new(); - } - let path_is_closed = !path_2d.open; // If false => open path, if true => closed path - - // Extract path points (x,y,0) from path_2d - let mut path_points = Vec::with_capacity(path_2d.vertices.len()); - for v in &path_2d.vertices { - // We only take X & Y; Z is typically 0 for a 2D path - path_points.push(Point3::new(v.pos.x, v.pos.y, 0.0)); - } - - // Convert the shape_2d into a list of its vertices in local coords (usually in XY). - // We assume shape_2d is a single polygon (can also handle multiple if needed). - let shape_is_closed = !shape_2d.open && shape_2d.vertices.len() >= 3; - let shape_count = shape_2d.vertices.len(); - - // For each path vertex, compute the orientation that aligns +Z to the path tangent. - // Then transform the shape’s 2D vertices into 3D “slice[i]”. - let n_path = path_points.len(); - let mut slices: Vec>> = Vec::with_capacity(n_path); - - for i in 0..n_path { - // The path tangent is p[i+1] - p[i] (or wrap if path is closed) - // If open and i == n_path-1 => we’ll copy the tangent from the last segment - let next_i = if i == n_path - 1 { - if path_is_closed { 0 } else { i - 1 } // if closed, wrap, else reuse the previous - } else { - i + 1 - }; - - let mut dir = path_points[next_i] - path_points[i]; - if dir.norm_squared() < EPSILON { - // Degenerate segment => fallback to the previous direction or just use +Z - dir = Vector3::z(); - } else { - dir.normalize_mut(); - } - - // Build a rotation that maps +Z to `dir`. - // We'll rotate the z-axis (0,0,1) onto `dir`. - let z = Vector3::z(); - let dot = z.dot(&dir); - // If dir is basically the same as z, no rotation needed - if (dot - 1.0).abs() < EPSILON { - return Matrix4::identity(); - } - // If dir is basically opposite z - if (dot + 1.0).abs() < EPSILON { - // 180 deg around X or Y axis - let rot180 = Rotation3::from_axis_angle(&Unit::new_normalize(Vector3::x()), PI); - return rot180.to_homogeneous(); - } - // Otherwise, general axis = z × dir - let axis = z.cross(&dir).normalize(); - let angle = z.dot(&dir).acos(); - let initial_rot = Rotation3::from_axis_angle(&Unit::new_unchecked(axis), angle); - let rot = initial_rot.to_homogeneous() - - // Build a translation that puts shape origin at path_points[i] - let trans = Translation3::from(path_points[i].coords); - - // Combined transform = T * R - let mat = trans.to_homogeneous() * rot; - - // Apply that transform to all shape_2d vertices => slice[i] - let mut slice_i = Vec::with_capacity(shape_count); - for sv in &shape_2d.vertices { - let local_pt = sv.pos; // (x, y, z=0) - let p4 = local_pt.to_homogeneous(); - let p4_trans = mat * p4; - slice_i.push(Point3::from_homogeneous(p4_trans).unwrap()); - } - slices.push(slice_i); - } - - // Build polygons for the new 3D swept solid. - // - (A) “Cap” polygons at start & end if path is open. - // - (B) “Side wall” quads between slice[i] and slice[i+1]. - // - // We’ll gather them all into a Vec>, then make a CSG. - - let mut all_polygons = Vec::new(); - - // Caps if path is open - // We replicate the shape_2d as polygons at slice[0] and slice[n_path-1]. - // We flip the first one so its normal faces outward. The last we keep as is. - if !path_is_closed { - // “Bottom” cap = slice[0], but we flip its winding so outward normal is “down” the path - if shape_is_closed { - let bottom_poly = polygon_from_slice( - &slices[0], - true, // flip - shape_2d.metadata.clone(), - ); - all_polygons.push(bottom_poly); - } - // “Top” cap = slice[n_path-1] (no flip) - if shape_is_closed { - let top_poly = polygon_from_slice( - &slices[n_path - 1], - false, // no flip - shape_2d.metadata.clone(), - ); - all_polygons.push(top_poly); - } - } - - // Side walls: For i in [0..n_path-1], or [0..n_path] if closed - let end_index = if path_is_closed { n_path } else { n_path - 1 }; - - for i in 0..end_index { - let i_next = (i + 1) % n_path; // wraps if closed - let slice_i = &slices[i]; - let slice_next = &slices[i_next]; - - // For each edge in the shape, connect vertices k..k+1 - // shape_2d may be open or closed. If open, we do shape_count-1 edges; if closed, shape_count edges. - let edge_count = if shape_is_closed { - shape_count // because last edge wraps - } else { - shape_count - 1 - }; - - for k in 0..edge_count { - let k_next = (k + 1) % shape_count; - - let v_i_k = slice_i[k]; - let v_i_knext = slice_i[k_next]; - let v_next_k = slice_next[k]; - let v_next_knext = slice_next[k_next]; - - // Build a quad polygon in CCW order for outward normal - // or you might choose a different ordering. Typically: - // [v_i_k, v_i_knext, v_next_knext, v_next_k] - // forms an outward-facing side wall if the shape_2d was originally CCW in XY. - let side_poly = Polygon::new( - vec![ - Vertex::new(v_i_k, Vector3::zeros()), - Vertex::new(v_i_knext, Vector3::zeros()), - Vertex::new(v_next_knext, Vector3::zeros()), - Vertex::new(v_next_k, Vector3::zeros()), - ], - shape_2d.metadata.clone(), - ); - all_polygons.push(side_poly); - } - } - - // Combine into a final CSG - CSG::from_polygons(&all_polygons) - } - */ + // pub fn sweep(shape_2d: &Polygon, path_2d: &Polygon) -> CSG { + // Gather the path’s vertices in XY + // if path_2d.vertices.len() < 2 { + // Degenerate path => no sweep + // return CSG::new(); + // } + // let path_is_closed = !path_2d.open; // If false => open path, if true => closed path + // + // Extract path points (x,y,0) from path_2d + // let mut path_points = Vec::with_capacity(path_2d.vertices.len()); + // for v in &path_2d.vertices { + // We only take X & Y; Z is typically 0 for a 2D path + // path_points.push(Point3::new(v.pos.x, v.pos.y, 0.0)); + // } + // + // Convert the shape_2d into a list of its vertices in local coords (usually in XY). + // We assume shape_2d is a single polygon (can also handle multiple if needed). + // let shape_is_closed = !shape_2d.open && shape_2d.vertices.len() >= 3; + // let shape_count = shape_2d.vertices.len(); + // + // For each path vertex, compute the orientation that aligns +Z to the path tangent. + // Then transform the shape’s 2D vertices into 3D “slice[i]”. + // let n_path = path_points.len(); + // let mut slices: Vec>> = Vec::with_capacity(n_path); + // + // for i in 0..n_path { + // The path tangent is p[i+1] - p[i] (or wrap if path is closed) + // If open and i == n_path-1 => we’ll copy the tangent from the last segment + // let next_i = if i == n_path - 1 { + // if path_is_closed { 0 } else { i - 1 } // if closed, wrap, else reuse the previous + // } else { + // i + 1 + // }; + // + // let mut dir = path_points[next_i] - path_points[i]; + // if dir.norm_squared() < EPSILON { + // Degenerate segment => fallback to the previous direction or just use +Z + // dir = Vector3::z(); + // } else { + // dir.normalize_mut(); + // } + // + // Build a rotation that maps +Z to `dir`. + // We'll rotate the z-axis (0,0,1) onto `dir`. + // let z = Vector3::z(); + // let dot = z.dot(&dir); + // If dir is basically the same as z, no rotation needed + // if (dot - 1.0).abs() < EPSILON { + // return Matrix4::identity(); + // } + // If dir is basically opposite z + // if (dot + 1.0).abs() < EPSILON { + // 180 deg around X or Y axis + // let rot180 = Rotation3::from_axis_angle(&Unit::new_normalize(Vector3::x()), PI); + // return rot180.to_homogeneous(); + // } + // Otherwise, general axis = z × dir + // let axis = z.cross(&dir).normalize(); + // let angle = z.dot(&dir).acos(); + // let initial_rot = Rotation3::from_axis_angle(&Unit::new_unchecked(axis), angle); + // let rot = initial_rot.to_homogeneous() + // + // Build a translation that puts shape origin at path_points[i] + // let trans = Translation3::from(path_points[i].coords); + // + // Combined transform = T * R + // let mat = trans.to_homogeneous() * rot; + // + // Apply that transform to all shape_2d vertices => slice[i] + // let mut slice_i = Vec::with_capacity(shape_count); + // for sv in &shape_2d.vertices { + // let local_pt = sv.pos; // (x, y, z=0) + // let p4 = local_pt.to_homogeneous(); + // let p4_trans = mat * p4; + // slice_i.push(Point3::from_homogeneous(p4_trans).unwrap()); + // } + // slices.push(slice_i); + // } + // + // Build polygons for the new 3D swept solid. + // - (A) “Cap” polygons at start & end if path is open. + // - (B) “Side wall” quads between slice[i] and slice[i+1]. + // + // We’ll gather them all into a Vec>, then make a CSG. + // + // let mut all_polygons = Vec::new(); + // + // Caps if path is open + // We replicate the shape_2d as polygons at slice[0] and slice[n_path-1]. + // We flip the first one so its normal faces outward. The last we keep as is. + // if !path_is_closed { + // “Bottom” cap = slice[0], but we flip its winding so outward normal is “down” the path + // if shape_is_closed { + // let bottom_poly = polygon_from_slice( + // &slices[0], + // true, // flip + // shape_2d.metadata.clone(), + // ); + // all_polygons.push(bottom_poly); + // } + // “Top” cap = slice[n_path-1] (no flip) + // if shape_is_closed { + // let top_poly = polygon_from_slice( + // &slices[n_path - 1], + // false, // no flip + // shape_2d.metadata.clone(), + // ); + // all_polygons.push(top_poly); + // } + // } + // + // Side walls: For i in [0..n_path-1], or [0..n_path] if closed + // let end_index = if path_is_closed { n_path } else { n_path - 1 }; + // + // for i in 0..end_index { + // let i_next = (i + 1) % n_path; // wraps if closed + // let slice_i = &slices[i]; + // let slice_next = &slices[i_next]; + // + // For each edge in the shape, connect vertices k..k+1 + // shape_2d may be open or closed. If open, we do shape_count-1 edges; if closed, shape_count edges. + // let edge_count = if shape_is_closed { + // shape_count // because last edge wraps + // } else { + // shape_count - 1 + // }; + // + // for k in 0..edge_count { + // let k_next = (k + 1) % shape_count; + // + // let v_i_k = slice_i[k]; + // let v_i_knext = slice_i[k_next]; + // let v_next_k = slice_next[k]; + // let v_next_knext = slice_next[k_next]; + // + // Build a quad polygon in CCW order for outward normal + // or you might choose a different ordering. Typically: + // [v_i_k, v_i_knext, v_next_knext, v_next_k] + // forms an outward-facing side wall if the shape_2d was originally CCW in XY. + // let side_poly = Polygon::new( + // vec![ + // Vertex::new(v_i_k, Vector3::zeros()), + // Vertex::new(v_i_knext, Vector3::zeros()), + // Vertex::new(v_next_knext, Vector3::zeros()), + // Vertex::new(v_next_k, Vector3::zeros()), + // ], + // shape_2d.metadata.clone(), + // ); + // all_polygons.push(side_poly); + // } + // } + // + // Combine into a final CSG + // CSG::from_polygons(&all_polygons) + // } } /// Helper to build a single Polygon from a “slice” of 3D points. diff --git a/src/sketch/hershey.rs b/src/sketch/hershey.rs index 564ccfc..452bbfc 100644 --- a/src/sketch/hershey.rs +++ b/src/sketch/hershey.rs @@ -21,7 +21,6 @@ impl Sketch { /// /// # Returns /// A new `Sketch` where each glyph stroke is a `Geometry::LineString` in `geometry`. - /// pub fn from_hershey( text: &str, font: &Font, diff --git a/src/sketch/shapes.rs b/src/sketch/shapes.rs index 0203642..23cc5d7 100644 --- a/src/sketch/shapes.rs +++ b/src/sketch/shapes.rs @@ -21,7 +21,8 @@ impl Sketch { /// /// # Example /// ``` - /// let sq2 = Sketch::rectangle(2.0, 3.0, None); + /// use csgrs::sketch::Sketch; + /// let sq2: Sketch<()> = Sketch::rectangle(2.0, 3.0, None); /// ``` pub fn rectangle(width: Real, length: Real, metadata: Option) -> Self { // In geo, a Polygon is basically (outer: LineString, Vec for holes). diff --git a/src/tests.rs b/src/tests.rs index 36a6947..ae43e76 100644 --- a/src/tests.rs +++ b/src/tests.rs @@ -1401,7 +1401,7 @@ fn test_same_number_of_vertices() { // - 3 side polygons (one for each edge of the triangle) assert_eq!( csg.polygons.len(), - 1 /*bottom*/ + 1 /*top*/ + 3 /*sides*/ + 1 /*bottom*/ + 1 /*top*/ + 3 // sides ); }