Refactor/mesh algorithm dependency inversion

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::<usize>() 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 {

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,

src/io/stl.rs

@@ -406,7 +406,6 @@ impl<S: Clone + Debug + Send + Sync> Sketch<S> {
             }
         }
 
-        //
         // (C) Encode into a binary STL buffer
         //
         let mut cursor = Cursor::new(Vec::new());

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

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"),
     //);

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;

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<S: Clone + Debug + Send + Sync> ConvexHullOps<S> for ParallelConvexHullOps {
+    fn convex_hull(&self, mesh: &Mesh<S>) -> Mesh<S> {
+        // Gather all (x, y, z) coordinates from the polygons in parallel
+        let points: Vec<Point3<Real>> = mesh
+            .polygons
+            .par_iter()
+            .flat_map(|poly| poly.vertices.par_iter().map(|v| v.pos))
+            .collect();
+
+        let points_for_hull: Vec<Vec<Real>> =
+            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<S>, other: &Mesh<S>) -> Mesh<S> {
+        // Collect all vertices (x, y, z) from self
+        let verts_a: Vec<Point3<Real>> = 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<Point3<Real>> = 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<Polygon<S>> = 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())
+    }
+}

src/mesh/convex_hull.rs → 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<S: Clone + Debug + Send + Sync> Mesh<S> {
-    /// Compute the convex hull of all vertices in this Mesh.
-    pub fn convex_hull(&self) -> Mesh<S> {
+/// 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<S: Clone + Debug + Send + Sync> ConvexHullOps<S> for SerialConvexHullOps {
+    fn convex_hull(&self, mesh: &Mesh<S>) -> Mesh<S> {
         // Gather all (x, y, z) coordinates from the polygons
-        let points: Vec<Point3<Real>> = self
+        let points: Vec<Point3<Real>> = mesh
             .polygons
             .iter()
             .flat_map(|poly| poly.vertices.iter().map(|v| v.pos))
@@ -54,19 +58,12 @@ impl<S: Clone + Debug + Send + Sync> Mesh<S> {
             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<S>) -> Mesh<S> {
+    fn minkowski_sum(&self, mesh: &Mesh<S>, other: &Mesh<S>) -> Mesh<S> {
         // Collect all vertices (x, y, z) from self
-        let verts_a: Vec<Point3<Real>> = self
+        let verts_a: Vec<Point3<Real>> = mesh
             .polygons
             .iter()
             .flat_map(|poly| poly.vertices.iter().map(|v| v.pos))
@@ -132,6 +129,6 @@ impl<S: Clone + Debug + Send + Sync> Mesh<S> {
             })
             .collect();
 
-        Mesh::from_polygons(&polygons, self.metadata.clone())
+        Mesh::from_polygons(&polygons, mesh.metadata.clone())
     }
 }

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<S: Clone + Debug + Send + Sync> {
+    /// Computes the convex hull of the mesh.
+    fn convex_hull(&self, mesh: &Mesh<S>) -> Mesh<S>;
+
+    /// Computes the Minkowski sum of two meshes.
+    fn minkowski_sum(&self, mesh: &Mesh<S>, other: &Mesh<S>) -> Mesh<S>;
+}