rscad_csg_bsp/
polygon.rs

1//! A polygon module, containing the `Polygon` struct and related types and methods.
2//! A polygon is a simple 2d shape in an N dimensional space, defined by its vertices and holes, and the plane it lies on.  It can be tessellated into triangles using the `tessellate_earcut` method, which uses the `earcut` crate, or the `tessellate_delaunay` method, which uses the `spade` crate for Delaunay triangulation.  The tessellated polygon is represented by the `TessellatedPolygon` struct, which contains a reference to the original polygon and a list of triangle indices.
3use crate::metadata::FaceKey;
4use crate::{Index, Number, Plane};
5pub mod tessellate;
6use earcut::Earcut;
7use nalgebra::{Point, Point3, SVector, TAffine, Transform};
8use ndarray::{Array1, ArrayView2, Ix2};
9use spade::{DelaunayTriangulation, HasPosition, Point2 as SpadePoint2, SpadeNum, Triangulation};
10
11#[repr(C)]
12#[derive(Clone, Debug, PartialEq)]
13pub struct TessellatedPolygon<'p, T: Number, const N: usize, S: Index = usize> {
14    polygon: &'p Polygon<T, N>,
15    triangles: Vec<S>,
16}
17
18impl<T, const N: usize, S> TessellatedPolygon<'_, T, N, S>
19where
20    T: Number,
21    S: Index,
22{
23    pub fn polygon(&self) -> &Polygon<T, N> {
24        self.polygon
25    }
26
27    pub fn into_triangles(self) -> Vec<S> {
28        self.triangles
29    }
30    pub(crate) fn as_triangles(&self) -> ArrayView2<'_, S> {
31        let triangles: &[S] = self.triangles.as_slice();
32        ArrayView2::from_shape([triangles.len() / 3, 3], triangles)
33            .expect("BUG: triangles length should be divisible by 3 possible bug in earcut")
34    }
35
36    // pub fn triangles(&self) -> CowArray<'_, S, Ix2> {
37    pub fn triangles(&self) -> ndarray::ArrayView<'_, S, Ix2> {
38        let axis = self.polygon.plane.normal().iamax();
39        let normal_component = self.polygon.plane.normal()[axis];
40        // When the dominant axis is Y (axis==1), the 2D projection is XZ where
41        // X×Z = -Y, so CCW in 2D corresponds to the *negative* Y direction.
42        // For X (Y×Z=+X) and Z (X×Y=+Z), CCW corresponds to the positive direction.
43        let inverse = (axis == 1) != (normal_component < T::zero());
44
45        let mut triangles = self.as_triangles();
46        if inverse {
47            triangles.invert_axis(ndarray::Axis(1));
48            triangles
49        } else {
50            triangles
51        }
52    }
53}
54
55#[cfg(feature = "bevy")]
56impl<'p> TessellatedPolygon<'p, f32, 3, u32> {
57    pub fn to_bevy_mesh(&self) -> bevy_mesh::Mesh {
58        use bevy_asset::RenderAssetUsages;
59        use bevy_mesh::{Indices, Mesh};
60        use wgpu::PrimitiveTopology;
61        let mut mesh = Mesh::new(
62            PrimitiveTopology::TriangleList,
63            RenderAssetUsages::default(),
64        );
65
66        let vertices: Vec<[f32; 3]> = self
67            .polygon
68            .vertices
69            .iter()
70            .copied()
71            .map(|v| [v.x, v.y, v.z])
72            .collect();
73
74        mesh.insert_attribute(bevy_mesh::Mesh::ATTRIBUTE_POSITION, vertices);
75        // mesh.insert_attribute(
76        //     bevy_mesh::Mesh::ATTRIBUTE_NORMAL,
77        //     vec![[0.0, 0.0, 1.0]; self.polygon.vertices.len()],
78        // );
79        mesh.insert_indices(Indices::U32(self.triangles.clone()));
80        mesh
81    }
82}
83
84/// A simple 2d polygon in an N dimensional space, defined by its vertices and holes, and the plane it lies on
85#[repr(C)]
86#[derive(Clone, Debug, PartialEq)]
87pub struct Polygon<T: Number, const N: usize> {
88    pub(crate) vertices: Array1<Point<T, N>>,
89    holes: Vec<Array1<usize>>,
90    pub(crate) plane: Plane<T, N>,
91    face_key: Option<FaceKey>,
92}
93
94impl<T, const N: usize> core::ops::Mul<Polygon<T, N>> for nalgebra::SMatrix<T, N, N>
95where
96    T: Number,
97{
98    type Output = Polygon<T, N>;
99
100    fn mul(self, rhs: Polygon<T, N>) -> Self::Output {
101        // NOTE: per-vertex rayon parallelism benchmarked slower than sequential here
102        // due to dispatch overhead on cheap ops (matrix * point). See benches/rayon_overhead.rs.
103        rhs.map_vertices(|v| self * v)
104    }
105}
106
107impl<T, const N: usize> Polygon<T, N>
108where
109    T: Number,
110{
111    pub fn from_points_and_plane(
112        vertices: impl Into<Vec<Point<T, N>>>,
113        plane: Plane<T, N>,
114    ) -> Result<Self, String> {
115        let vertices = vertices.into();
116        // NOTE: per-vertex rayon parallelism benchmarked slower than sequential here
117        // due to dispatch overhead on cheap ops (point coplanarity check). See benches/rayon_overhead.rs.
118        if vertices.iter().any(|v| !plane.has_point(v)) {
119            return Err("All points must be coplanar with the plane".to_string());
120        }
121        Ok(Self {
122            vertices: Array1::from_vec(vertices),
123            holes: Vec::new(),
124            plane,
125            face_key: None,
126        })
127    }
128
129    pub(crate) fn clone_with_points(&self, vertices: Array1<Point<T, N>>) -> Self {
130        Self {
131            vertices,
132            holes: Vec::new(),
133            plane: self.plane,
134            face_key: self.face_key,
135        }
136    }
137
138    pub fn face_key(&self) -> Option<FaceKey> {
139        self.face_key
140    }
141
142    pub fn set_face_key(&mut self, key: Option<FaceKey>) {
143        self.face_key = key;
144    }
145
146    pub fn holes(&self) -> &[ndarray::Array1<usize>] {
147        self.holes.as_slice()
148    }
149
150    pub fn vertices(&self) -> ndarray::ArrayView1<'_, Point<T, N>> {
151        self.vertices.view()
152    }
153
154    pub fn translate(mut self, offset: impl Into<SVector<T, N>>) -> Self {
155        let offset = offset.into();
156        // NOTE: per-vertex rayon parallelism benchmarked slower than sequential here
157        // due to dispatch overhead on cheap ops (point + vector). See benches/rayon_overhead.rs.
158        self.vertices.mapv_inplace(|v| v + offset);
159        self.plane = Plane::new(
160            self.plane.normal(),
161            self.plane.distance() + self.plane.normal().dot(&offset),
162        );
163        self
164    }
165
166    pub fn plane(&self) -> Plane<T, N> {
167        self.plane
168    }
169
170    /// Takes an affine transformation matrix and applies it to a n dimensional point polygon
171    pub fn transform(mut self, transform: Transform<T, TAffine, N>) -> Self
172    where
173        nalgebra::Const<N>: nalgebra::ToTypenum,
174        <nalgebra::Const<N> as nalgebra::ToTypenum>::Typenum:
175            core::ops::Add<typenum::uint::UInt<typenum::uint::UTerm, typenum::bit::B1>>,
176        nalgebra::Const<N>: nalgebra::DimNameAdd<nalgebra::U1>,
177        nalgebra::DefaultAllocator: nalgebra::allocator::Allocator<
178                nalgebra::DimNameSum<nalgebra::Const<N>, nalgebra::U1>,
179                nalgebra::DimNameSum<nalgebra::Const<N>, nalgebra::U1>,
180            > + nalgebra::allocator::Allocator<nalgebra::DimNameSum<nalgebra::Const<N>, nalgebra::U1>>,
181    {
182        // NOTE: per-vertex rayon parallelism benchmarked slower than sequential here
183        // due to dispatch overhead on cheap ops (affine transform). See benches/rayon_overhead.rs.
184        self.vertices
185            .mapv_inplace(|v| transform.transform_point(&v));
186        self
187    }
188
189    pub fn bounding_box(&self) -> bounding_box::AxisAlignedBoundingBox<T, N> {
190        let mut min = self.vertices[0];
191        let mut max = min;
192        self.vertices.iter().for_each(|v| {
193            min = min.inf(v);
194            max = max.sup(v);
195        });
196        bounding_box::AxisAlignedBoundingBox::new(min, max)
197    }
198
199    #[inline(always)]
200    #[must_use = "map_vertices takes ownership of the original polygon and returns a new one with the mapped vertices"]
201    pub fn map_vertices<F>(mut self, f: F) -> Self
202    where
203        F: Fn(Point<T, N>) -> Point<T, N>,
204    {
205        self.vertices.mapv_inplace(&f);
206        let d = self
207            .plane
208            .normal()
209            .dot(&f(self.plane.origin_projection()).coords);
210        self.plane = Plane::new(self.plane.normal(), d);
211        self
212    }
213
214    /// Cast the polygon's numeric type from `T` to `T2`.
215    pub fn cast<T2: Number>(self) -> Polygon<T2, N>
216    where
217        T: num_traits::cast::AsPrimitive<T2>,
218    {
219        Polygon {
220            vertices: self.vertices.mapv(|v| v.map(|x| x.as_())),
221            holes: self.holes,
222            plane: self.plane.cast(),
223            face_key: self.face_key,
224        }
225    }
226}
227
228impl<T> Polygon<T, 3>
229where
230    T: Number,
231{
232    /// Construct a polygon from a list of coplanar 3D points.
233    /// The plane is derived from the first three vertices.
234    /// Returns an error if fewer than 3 vertices are provided, or if any
235    /// vertex is not coplanar with the plane defined by the first three.
236    pub fn from_points(
237        vertices: impl IntoIterator<Item = Point3<T>>,
238    ) -> Result<Self, PolygonError> {
239        let vertices: Vec<_> = vertices.into_iter().collect();
240        if vertices.len() < 3 {
241            return Err(PolygonError::NotEnoughPoints(vertices.len()));
242        }
243        let plane = Plane::from_points(vertices[0], vertices[1], vertices[2]);
244        if let Some((_i, v)) = vertices
245            .iter()
246            .enumerate()
247            .find(|(_, v)| !plane.has_point(*v))
248        {
249            let dist = plane.normal().dot(&v.coords) - plane.distance();
250            return Err(PolygonError::NonCoplanarPoint(
251                v.map(|c| num_traits::cast(c).unwrap_or(0.0)),
252                num_traits::cast(dist).unwrap_or(0.0),
253            ));
254        }
255        Ok(Self {
256            vertices: Array1::from_vec(vertices),
257            holes: Vec::new(),
258            plane,
259            face_key: None,
260        })
261    }
262
263    pub fn circle(radius: T, segments: usize, plane: Plane<T, 3>) -> Self {
264        let one = T::one();
265        let two = one + one;
266        let two_pi = two * T::PI();
267        let transformation_matrix = Plane::xy_plane().rotation_to(plane);
268        let vertices = (0..segments)
269            .map(|i| {
270                let angle = T::from(i).expect("BUG: loop index is representable as T") * two_pi
271                    / T::from(segments).expect("BUG: segment count is representable as T");
272                let x = radius * num_traits::real::Real::cos(angle);
273                let y = radius * num_traits::real::Real::sin(angle);
274                let v = transformation_matrix * SVector::from([x, y, T::zero()]);
275                v.into()
276            })
277            .collect();
278        Self {
279            vertices: Array1::from_vec(vertices),
280            holes: Vec::new(),
281            plane,
282            face_key: None,
283        }
284    }
285
286    pub fn square(size: T, plane: Plane<T, 3>) -> Self {
287        let half_size = size / (T::one() + T::one());
288        let transformation_matrix = Plane::xy_plane().rotation_to(plane);
289        let vertices = [
290            SVector::from([-half_size, -half_size, T::zero()]),
291            SVector::from([half_size, -half_size, T::zero()]),
292            SVector::from([half_size, half_size, T::zero()]),
293            SVector::from([-half_size, half_size, T::zero()]),
294        ]
295        .into_iter()
296        .map(|v| transformation_matrix * v)
297        .map(|v| v.into())
298        .collect();
299        Self {
300            vertices: Array1::from_vec(vertices),
301            holes: Vec::new(),
302            plane,
303            face_key: None,
304        }
305    }
306
307    /// Rotate all vertices and the plane by a 3×3 rotation matrix.
308    ///
309    /// The rotation matrix must be orthogonal (e.g., from [`euler_rotation_matrix`]).
310    /// For an orthogonal matrix, the plane distance is preserved.
311    pub fn rotate(mut self, rotation: nalgebra::Matrix3<T>) -> Self {
312        // NOTE: per-vertex rayon parallelism benchmarked slower than sequential here
313        // due to dispatch overhead on cheap ops (matrix * vector). See benches/rayon_overhead.rs.
314        self.vertices.mapv_inplace(|v| {
315            let rv = rotation * v.coords;
316            Point3::from(rv)
317        });
318        let new_normal = rotation * self.plane.normal();
319        self.plane = Plane::new(new_normal, self.plane.distance());
320        self
321    }
322
323    /// Rotate by Euler angles (degrees), applied as Z × Y × X.
324    pub fn rotate_euler(self, angles_deg: nalgebra::Vector3<T>) -> Self {
325        self.rotate(crate::euler_rotation_matrix(angles_deg))
326    }
327
328    pub fn flip(mut self) -> Self {
329        self.vertices.invert_axis(ndarray::Axis(0));
330        self.plane = self.plane.flip();
331        self
332    }
333
334    pub fn flip_in_place(&mut self) {
335        self.vertices.invert_axis(ndarray::Axis(0));
336        self.plane = self.plane.flip();
337    }
338
339    /// Recompute the cached plane from the current vertices via Newell's
340    /// method (normal oriented by winding, matching [`Plane::from_points`],
341    /// but robust to near-collinear leading vertices).
342    ///
343    /// Required after any vertex map that doesn't preserve the plane normal
344    /// (e.g. non-uniform scale): a stale plane makes BSP classification
345    /// inconsistent with the vertices, which cascades into unbounded
346    /// polygon splitting.
347    pub fn recompute_plane(&mut self) {
348        let n = self.vertices.len();
349        if n < 3 {
350            return;
351        }
352        let mut normal = SVector::<T, 3>::zeros();
353        let mut centroid = SVector::<T, 3>::zeros();
354        for i in 0..n {
355            let a = self.vertices[i].coords;
356            let b = self.vertices[(i + 1) % n].coords;
357            normal.x += (a.y - b.y) * (a.z + b.z);
358            normal.y += (a.z - b.z) * (a.x + b.x);
359            normal.z += (a.x - b.x) * (a.y + b.y);
360            centroid += a;
361        }
362        let normal = normal.normalize();
363        let distance = normal.dot(&centroid) / T::from(n).expect("BUG: vertex count converts");
364        self.plane = Plane::new(normal, distance);
365    }
366}
367
368#[derive(Debug, Clone, Copy, PartialEq, Eq)]
369pub enum PolygonClassification {
370    Front,
371    Back,
372    Coplanar,
373    Spanning,
374}
375
376#[test]
377fn test_ndarray_invet_axis() {
378    let mut arr = Array1::from_vec(vec![1, 2, 3]);
379    arr.invert_axis(ndarray::Axis(0));
380    assert_eq!(arr, Array1::from_vec(vec![3, 2, 1]));
381}
382
383#[derive(thiserror::Error, Debug)]
384pub enum PolygonError {
385    #[error("Not enough points to construct a polygon from expected > 3 got ({0})")]
386    NotEnoughPoints(usize),
387    #[error("Vertex is not coplanar wit the rest")]
388    NonCoplanarPoint(Point3<f64>, f64),
389}
rscad_csg_bsp/polygon.rs

rscad_csg_bsp/
polygon.rs
