rscad_csg_bsp/

csg.rs

use nalgebra::{Affine3, Vector3};

use crate::*;
use core::iter::IntoIterator;

#[derive(Clone)]
pub struct Csg<T: Number> {
    polygons: Vec<Polygon<T, 3>>,
    /// Embedded metadata factory. Polygons reference it via `face_key`.
    factory: CsgFactory<T>,
}

impl<T: Number> Csg<T> {
    pub fn new(polygons: Vec<Polygon<T, 3>>) -> Self {
        Self {
            polygons,
            factory: CsgFactory::new(),
        }
    }

    pub fn from_polygons(polygons: impl IntoIterator<Item = Polygon<T, 3>>) -> Self {
        Self {
            polygons: polygons.into_iter().collect(),
            factory: CsgFactory::new(),
        }
    }

    pub fn factory(&self) -> &CsgFactory<T> {
        &self.factory
    }

    pub fn factory_mut(&mut self) -> &mut CsgFactory<T> {
        &mut self.factory
    }

    /// Replace the embedded factory (e.g. after building via an external [`CsgFactory`]).
    pub fn with_factory(mut self, factory: CsgFactory<T>) -> Self {
        self.factory = factory;
        self
    }

    pub fn union(mut self, other: Csg<T>) -> Self {
        let remap = self.factory.absorb(&other.factory);
        let other_polygons = remap_polygons(other.polygons, &remap);
        let mut a = BspNode::build(self.polygons);
        let mut b = BspNode::build(other_polygons);
        a.clip_to(&b);
        b.clip_to(&a);
        b.invert();
        b.clip_to(&a);
        b.invert();
        a.build_into(b.polygons());
        Self {
            polygons: a.polygons(),
            factory: self.factory,
        }
    }

    pub fn intersection(mut self, other: Csg<T>) -> Self {
        let remap = self.factory.absorb(&other.factory);
        let other_polygons = remap_polygons(other.polygons, &remap);
        let mut a = BspNode::build(self.polygons);
        let mut b = BspNode::build(other_polygons);
        a.invert();
        b.clip_to(&a);
        b.invert();
        a.clip_to(&b);
        b.clip_to(&a);
        a.build_into(b.polygons());
        a.invert();
        Self {
            polygons: a.polygons(),
            factory: self.factory,
        }
    }

    pub fn polygons(&self) -> &[Polygon<T, 3>] {
        self.polygons.as_slice()
    }

    /// Run a cleanup pass: remove degenerate polygons, collapse collinear
    /// vertices, and merge coplanar adjacent polygons. Call once before
    /// mesh generation rather than after every CSG operation.
    pub fn cleanup(mut self) -> Self {
        self.polygons = crate::cleanup::cleanup_polygons(self.polygons);
        self
    }

    /// Apply a function to every polygon in place.
    pub fn mapv_inplace<F>(&mut self, f: F)
    where
        F: Fn(Polygon<T, 3>) -> Polygon<T, 3>,
    {
        let polys = core::mem::take(&mut self.polygons);
        self.polygons = polys.into_iter().map(f).collect();
    }

    /// Apply a function to every polygon in place, in parallel.
    #[cfg(feature = "rayon")]
    pub fn par_mapv_inplace<F>(&mut self, f: F)
    where
        F: Fn(Polygon<T, 3>) -> Polygon<T, 3> + Send + Sync,
    {
        use rayon::prelude::*;
        let polys = core::mem::take(&mut self.polygons);
        self.polygons = polys.into_par_iter().map(f).collect();
    }

    /// Apply a function to every vertex of every polygon in place.
    pub fn mapv_inplace_points<F>(&mut self, f: F)
    where
        F: Fn(nalgebra::Point3<T>) -> nalgebra::Point3<T>,
    {
        self.mapv_inplace(|p| p.map_vertices(&f));
    }

    /// Apply a function to every vertex of every polygon in place, in parallel.
    #[cfg(feature = "rayon")]
    pub fn par_mapv_inplace_points<F>(&mut self, f: F)
    where
        F: Fn(nalgebra::Point3<T>) -> nalgebra::Point3<T> + Send + Sync,
    {
        self.par_mapv_inplace(|p| p.map_vertices(&f));
    }

    /// Translate all polygons by an offset.
    pub fn translate(mut self, offset: nalgebra::Vector3<T>) -> Self {
        self.mapv_inplace(|p| p.translate(offset));
        self
    }

    /// Scale all polygon vertices by per-axis factors.
    ///
    /// Planes are recomputed from the scaled vertices: a non-uniform scale
    /// changes plane normals (inverse-transpose transform), so the cached
    /// plane can't just be carried over.
    pub fn scale(mut self, factors: nalgebra::Vector3<T>) -> Self {
        self.mapv_inplace(|mut p| {
            p.vertices
                .mapv_inplace(|v| v.coords.component_mul(&factors).into());
            p.recompute_plane();
            p
        });
        self
    }

    /// Rotate all polygons by Euler angles (degrees), applied as Z × Y × X.
    pub fn rotate_euler(mut self, angles_deg: nalgebra::Vector3<T>) -> Self {
        let rotation = crate::euler_rotation_matrix(angles_deg);
        self.mapv_inplace(|p| p.rotate(rotation));
        self
    }

    /// Mirror across axes where the corresponding component is non-zero.
    ///
    /// Reverses polygon winding to maintain correct normals.
    pub fn mirror(mut self, axes: nalgebra::Vector3<T>) -> Self {
        let signs = axes.map(|a| if a != T::zero() { -T::one() } else { T::one() });
        self.mapv_inplace(|mut p| {
            p.vertices
                .mapv_inplace(|v| v.coords.component_mul(&signs).into());
            p.vertices.invert_axis(ndarray::Axis(0));
            let normal = p.plane.normal().component_mul(&signs);
            let d = normal.dot(&p.vertices[0].coords);
            p.plane = Plane::new(normal, d);
            p
        });
        self
    }

    /// Symmetric difference: `(a - b) ∪ (b - a)`.
    pub fn xor(self, other: Csg<T>) -> Self {
        let a2 = self.clone();
        let b2 = other.clone();
        self.difference(b2).union(other.difference(a2))
    }

    /// Different (a - (a ^ b))
    pub fn difference(mut self, other: Csg<T>) -> Self {
        let remap = self.factory.absorb(&other.factory);
        let other_polygons = remap_polygons(other.polygons, &remap);
        let mut a = BspNode::build(self.polygons);
        let mut b = BspNode::build(other_polygons);
        a.invert();
        a.clip_to(&b);
        b.clip_to(&a);
        b.invert();
        b.clip_to(&a);
        b.invert();
        a.build_into(b.polygons());
        a.invert();
        Self {
            polygons: a.polygons(),
            factory: self.factory,
        }
    }

    /// Extrudes a planar shape along a vector.
    ///
    /// The extrusion is subdivided into `slices` levels. At each level the
    /// rotation angle and scale factor are linearly interpolated from identity
    /// (at the bottom) to the full `transform` (at the top).
    ///
    /// The `transform` is decomposed into a rotation (axis-angle) and uniform
    /// scale. Non-uniform scale is not supported.
    ///
    /// Side walls are triangulated to guarantee planarity when twist is used.
    pub fn extrude(
        self,
        along: Vector3<T>,
        transform: Affine3<T>,
        slices: usize,
    ) -> Result<Self, Error> {
        use nalgebra::Rotation3;
        use num_traits::Float;

        if !self.coplanar() {
            return Err(Error::NonPlanarExtrude);
        }

        let slices = slices.max(1);
        let slices_t = T::from(slices).expect("BUG: slices must convert to T");

        // Decompose transform into rotation (axis-angle) + uniform scale.
        let m = transform.matrix().fixed_view::<3, 3>(0, 0).into_owned();
        let scale_factor = Float::cbrt(m.determinant());
        let rot_matrix = m / scale_factor;
        let rotation = Rotation3::from_matrix_unchecked(rot_matrix);
        let axis_angle = rotation.scaled_axis();

        let vertex_at = |v: nalgebra::Point3<T>, t: T| {
            let s = T::one() + (scale_factor - T::one()) * t;
            let rot = Rotation3::from_scaled_axis(axis_angle * t);
            let offset = along * t;
            nalgebra::Point3::from((rot * (v.coords * s)) + offset)
        };

        use itertools::Itertools;

        let polygons: Vec<Polygon<T, 3>> =
            self.polygons
                .iter()
                .flat_map(|base_poly| {
                    let aligned = base_poly.plane().normal().dot(&along) > T::zero();

                    let base_verts: Vec<nalgebra::Point3<T>> =
                        base_poly.vertices().iter().copied().collect();

                    let levels: Vec<Vec<nalgebra::Point3<T>>> = (0..=slices)
                        .map(|i| {
                            let t = T::from(i).expect("BUG: slice index") / slices_t;
                            base_verts.iter().map(|&v| vertex_at(v, t)).collect()
                        })
                        .collect();

                    let sides = levels
                        .iter()
                        .map(|level| {
                            level
                                .iter()
                                .copied()
                                .circular_tuple_windows::<(_, _)>()
                                .collect::<Vec<_>>()
                        })
                        .tuple_windows()
                        .flat_map(move |(bot_edges, top_edges)| {
                            bot_edges.into_iter().zip(top_edges).flat_map(
                                move |((a, b), (d, c))| {
                                    if aligned {
                                        [
                                            Polygon::from_points([a, b, c]).ok(),
                                            Polygon::from_points([a, c, d]).ok(),
                                        ]
                                    } else {
                                        [
                                            Polygon::from_points([a, c, b]).ok(),
                                            Polygon::from_points([a, d, c]).ok(),
                                        ]
                                    }
                                    .into_iter()
                                    .flatten()
                                },
                            )
                        });

                    let first_level = levels.first().expect("BUG: levels is non-empty");
                    let last_level = levels.last().expect("BUG: levels is non-empty");
                    let (bottom_cap, top_cap) = if aligned {
                        (
                            Polygon::from_points(first_level.iter().rev().copied()),
                            Polygon::from_points(last_level.iter().copied()),
                        )
                    } else {
                        (
                            Polygon::from_points(first_level.iter().copied()),
                            Polygon::from_points(last_level.iter().rev().copied()),
                        )
                    };

                    sides.chain(bottom_cap).chain(top_cap).collect::<Vec<_>>()
                })
                .collect();

        Ok(Csg::from_polygons(polygons))
    }

    /// Revolves a planar XZ-plane profile about the +Y axis (lathe).
    ///
    /// Profile vertices `(x, 0, z)` sweep to `(x·cosθ, z, x·sinθ)`; profile-x
    /// is the radius (x ≥ 0 expected), profile-z the height along the axis.
    /// Full revolutions weld the seam; partial sweeps are closed with start
    /// and end caps. Side walls are triangulated (swept quads are non-planar),
    /// and triangles degenerated by on-axis vertices are dropped.
    pub fn lathe(self, angle_degrees: T, segments: usize) -> Result<Self, Error> {
        use itertools::Itertools;
        use num_traits::Float;

        if !self.coplanar() {
            return Err(Error::NonPlanarExtrude);
        }

        let segments = segments.max(3);
        let full_turn = T::from(360.0).expect("BUG: const");
        let full_circle = angle_degrees >= full_turn;
        let sweep = if full_circle {
            full_turn
        } else {
            angle_degrees
        };
        let step =
            Float::to_radians(sweep) / T::from(segments).expect("BUG: segments must convert to T");

        let vertex_at = |v: nalgebra::Point3<T>, i: usize| {
            let theta = step * T::from(i).expect("BUG: segment index");
            nalgebra::Point3::new(v.x * Float::cos(theta), v.z, v.x * Float::sin(theta))
        };

        // Full circles reuse level 0 as the final level so the seam welds exactly.
        let level_count = if full_circle { segments } else { segments + 1 };

        let polygons: Vec<Polygon<T, 3>> = self
            .polygons
            .iter()
            .flat_map(|base_poly| {
                // Profile polygons wind with normal −Y in the standard case
                // (see `rect`); a +Y normal flips every winding below.
                let aligned = base_poly.plane().normal().y < T::zero();

                let base_verts: Vec<nalgebra::Point3<T>> =
                    base_poly.vertices().iter().copied().collect();

                let levels: Vec<Vec<nalgebra::Point3<T>>> = (0..level_count)
                    .map(|i| base_verts.iter().map(|&v| vertex_at(v, i)).collect())
                    .collect();

                let edges_of = |level: &Vec<nalgebra::Point3<T>>| {
                    level
                        .iter()
                        .copied()
                        .circular_tuple_windows::<(_, _)>()
                        .collect::<Vec<_>>()
                };

                let sides =
                    (0..segments)
                        .map(|i| (i, (i + 1) % level_count))
                        .flat_map(|(bot, top)| {
                            edges_of(&levels[bot])
                                .into_iter()
                                .zip(edges_of(&levels[top]))
                                .flat_map(move |((a, b), (d, c))| {
                                    if aligned {
                                        [
                                            Polygon::from_points([a, b, c]).ok(),
                                            Polygon::from_points([a, c, d]).ok(),
                                        ]
                                    } else {
                                        [
                                            Polygon::from_points([a, c, b]).ok(),
                                            Polygon::from_points([a, d, c]).ok(),
                                        ]
                                    }
                                    .into_iter()
                                    .flatten()
                                })
                                .collect::<Vec<_>>()
                        });

                let caps = if full_circle {
                    [None, None]
                } else {
                    let first_level = levels.first().expect("BUG: levels is non-empty");
                    let last_level = levels.last().expect("BUG: levels is non-empty");
                    if aligned {
                        [
                            Polygon::from_points(first_level.iter().rev().copied()).ok(),
                            Polygon::from_points(last_level.iter().copied()).ok(),
                        ]
                    } else {
                        [
                            Polygon::from_points(first_level.iter().copied()).ok(),
                            Polygon::from_points(last_level.iter().rev().copied()).ok(),
                        ]
                    }
                };

                sides.chain(caps.into_iter().flatten()).collect::<Vec<_>>()
            })
            .collect();

        Ok(Csg::from_polygons(polygons))
    }

    pub fn coplanar(&self) -> bool {
        use itertools::Itertools;
        self.polygons()
            .iter()
            .tuple_windows()
            .all(|(last, next)| last.plane().roughly_eq(&next.plane()))
    }

    /// Cast the CSG solid's numeric type from `T` to `T2`, e.g. `f64` → `f32`.
    ///
    /// Face keys are preserved unchanged; only vertex coordinates and plane
    /// equations are converted via [`num_traits::cast::AsPrimitive`].
    pub fn cast<T2: Number>(self) -> Csg<T2>
    where
        T: num_traits::cast::AsPrimitive<T2>,
    {
        let (new_factory, face_remap) = self.factory.cast::<T2>();
        let polygons = self
            .polygons
            .into_iter()
            .map(|p| {
                let mut p2 = p.cast::<T2>();
                if let Some(fk) = p2.face_key()
                    && let Some(&new_fk) = face_remap.0.get(&fk)
                {
                    p2.set_face_key(Some(new_fk));
                }
                p2
            })
            .collect();
        Csg {
            polygons,
            factory: new_factory,
        }
    }
}

/// Update polygon `face_key`s using a remap from a factory merge.
fn remap_polygons<T: Number>(
    mut polygons: Vec<Polygon<T, 3>>,
    remap: &FaceKeyRemap,
) -> Vec<Polygon<T, 3>> {
    if remap.0.is_empty() {
        return polygons;
    }
    for p in &mut polygons {
        if let Some(fk) = p.face_key()
            && let Some(&new_fk) = remap.0.get(&fk)
        {
            p.set_face_key(Some(new_fk));
        }
    }
    polygons
}

#[cfg(feature = "bevy")]
impl Csg<f32> {
    /// Look up the color for a polygon via face_key → FaceMeta → PrimitiveMeta → color.
    fn polygon_color(&self, polygon: &Polygon<f32, 3>) -> Option<[f32; 4]> {
        let fk = polygon.face_key()?;
        let fm = self.factory.face(fk)?;
        let pm = self.factory.primitive(fm.primitive_key)?;
        pm.color
    }

    pub fn to_bevy_mesh(&self) -> bevy_mesh::Mesh {
        use bevy_asset::RenderAssetUsages;
        use bevy_mesh::{Indices, Mesh};
        use wgpu::PrimitiveTopology;

        let has_colors = self
            .polygons
            .iter()
            .any(|p| self.polygon_color(p).is_some());
        let color_fn = |p: &Polygon<f32, 3>| self.polygon_color(p);

        let (vertices, normals, indices, colors) =
            build_mesh_buffers(&self.polygons, has_colors, &color_fn);

        let mut mesh = Mesh::new(
            PrimitiveTopology::TriangleList,
            RenderAssetUsages::default(),
        );
        mesh.insert_attribute(Mesh::ATTRIBUTE_POSITION, vertices);
        mesh.insert_attribute(Mesh::ATTRIBUTE_NORMAL, normals);
        if let Some(colors) = colors {
            mesh.insert_attribute(Mesh::ATTRIBUTE_COLOR, colors);
        }
        mesh.insert_indices(Indices::U32(indices));
        mesh
    }
}

#[cfg(all(feature = "bevy", feature = "rayon"))]
#[allow(clippy::type_complexity)]
fn build_mesh_buffers(
    polygons: &[Polygon<f32, 3>],
    emit_colors: bool,
    color_fn: &dyn Fn(&Polygon<f32, 3>) -> Option<[f32; 4]>,
) -> (
    Vec<[f32; 3]>,
    Vec<[f32; 3]>,
    Vec<u32>,
    Option<Vec<[f32; 4]>>,
) {
    use rayon::prelude::*;

    #[allow(clippy::type_complexity)]
    let chunks: Vec<(Vec<[f32; 3]>, [f32; 3], Vec<u32>, Option<[f32; 4]>)> = polygons
        .par_iter()
        .map(|polygon| {
            let normal: [f32; 3] = bytemuck::cast(polygon.plane().normal());
            let verts: Vec<[f32; 3]> = polygon
                .vertices()
                .iter()
                .map(|v| bytemuck::cast(*v))
                .collect();
            let tris: Vec<u32> = if polygon.vertices().len() == 3 {
                vec![0, 1, 2]
            } else {
                polygon
                    .tessellate_earcut::<u32>()
                    .triangles()
                    .iter()
                    .copied()
                    .collect()
            };
            // color_fn can't be called from par_iter (not Sync), so we use face_key as a proxy
            // and look up in the sequential pass below
            (verts, normal, tris, None)
        })
        .collect();

    let mut offset = 0u32;
    let offsets: Vec<u32> = chunks
        .iter()
        .map(|(verts, _, _, _)| {
            let o = offset;
            offset += verts.len() as u32;
            o
        })
        .collect();

    let total_verts: usize = chunks.iter().map(|(v, _, _, _)| v.len()).sum();
    let total_tris: usize = chunks.iter().map(|(_, _, t, _)| t.len()).sum();

    let mut vertices = Vec::with_capacity(total_verts);
    let mut normals = Vec::with_capacity(total_verts);
    let mut indices = Vec::with_capacity(total_tris);
    let mut colors = if emit_colors {
        Some(Vec::with_capacity(total_verts))
    } else {
        None
    };

    for (((verts, normal, tris, _), offset), polygon) in
        chunks.into_iter().zip(offsets).zip(polygons)
    {
        let n = verts.len();
        normals.extend(std::iter::repeat_n(normal, n));
        vertices.extend(verts);
        indices.extend(tris.into_iter().map(|i| i + offset));
        if let Some(ref mut color_buf) = colors {
            let c = color_fn(polygon).unwrap_or([1.0, 1.0, 1.0, 1.0]);
            color_buf.extend(std::iter::repeat_n(c, n));
        }
    }

    (vertices, normals, indices, colors)
}

#[cfg(all(feature = "bevy", not(feature = "rayon")))]
#[allow(clippy::type_complexity)]
fn build_mesh_buffers(
    polygons: &[Polygon<f32, 3>],
    emit_colors: bool,
    color_fn: &dyn Fn(&Polygon<f32, 3>) -> Option<[f32; 4]>,
) -> (
    Vec<[f32; 3]>,
    Vec<[f32; 3]>,
    Vec<u32>,
    Option<Vec<[f32; 4]>>,
) {
    let mut vertices: Vec<[f32; 3]> = Vec::new();
    let mut normals: Vec<[f32; 3]> = Vec::new();
    let mut indices: Vec<u32> = Vec::new();
    let mut colors: Option<Vec<[f32; 4]>> = if emit_colors { Some(Vec::new()) } else { None };
    let mut vertex_offset: u32 = 0;

    for polygon in polygons {
        let normal: [f32; 3] = bytemuck::cast(polygon.plane().normal());
        let n_verts = polygon.vertices().len();

        for v in polygon.vertices().iter() {
            vertices.push(bytemuck::cast(*v));
            normals.push(normal);
        }
        if let Some(ref mut color_buf) = colors {
            let c = color_fn(polygon).unwrap_or([1.0, 1.0, 1.0, 1.0]);
            color_buf.extend(std::iter::repeat_n(c, n_verts));
        }
        if n_verts == 3 {
            indices.extend_from_slice(&[vertex_offset, vertex_offset + 1, vertex_offset + 2]);
        } else {
            let tessellated = polygon.tessellate_earcut::<u32>();
            for idx in tessellated.triangles().iter() {
                indices.push(*idx + vertex_offset);
            }
        }
        vertex_offset += n_verts as u32;
    }

    (vertices, normals, indices, colors)
}

#[derive(thiserror::Error, Debug)]
pub enum Error {
    #[error("The shape you tried to extrude is non planar")]
    NonPlanarExtrude,
}

#[cfg(test)]
mod tests {
    /// Non-uniform scale changes plane normals (inverse-transpose transform).
    /// A stale cached plane leaves vertices off their own plane, which makes
    /// BSP classification inconsistent and cascades into unbounded polygon
    /// splitting (stack overflow on export).
    #[test]
    fn non_uniform_scale_keeps_polygon_planes_consistent() {
        let sphere = crate::CsgFactory::<f64>::new().sphere(0.5, 32, 32);
        let scaled = sphere.scale(nalgebra::Vector3::new(3.6, 3.98, 1.0));
        for polygon in scaled.polygons() {
            let plane = polygon.plane();
            for v in polygon.vertices().iter() {
                assert!(
                    plane.has_point(v),
                    "vertex {v:?} off its polygon plane {plane:?}"
                );
            }
        }
    }
}