rscad_csg_bsp/

plane.rs

use crate::{Number, Polygon, PolygonClassification};
use core::ops::{AddAssign, MulAssign, SubAssign};
use nalgebra::{Point, Point3, SVector};
use num_traits::real::Real;

pub type Plane3<T> = Plane<T, 3>;

/// A plane in 3D space, defined by a normal vector and a distance from the origin
#[repr(C)]
#[derive(Clone, Copy, Debug, PartialOrd, PartialEq)]
pub struct Plane<T: Number, const N: usize> {
    distance: T,
    normal: SVector<T, N>,
}

#[derive(Copy, Clone, Debug, PartialEq, Eq, PartialOrd, Hash)]
pub enum PlaneSide {
    Front,
    Back,
    On,
}

impl core::ops::Neg for PlaneSide {
    type Output = Self;

    fn neg(self) -> Self::Output {
        match self {
            PlaneSide::Front => PlaneSide::Back,
            PlaneSide::Back => PlaneSide::Front,
            PlaneSide::On => PlaneSide::On,
        }
    }
}

impl<T: Number> Plane<T, 3> {
    pub fn xy_plane() -> Self {
        Self {
            normal: *SVector::ith_axis(2),
            distance: T::zero(),
        }
    }

    pub fn xz_plane() -> Self {
        Self {
            normal: *SVector::ith_axis(1),
            distance: T::zero(),
        }
    }

    pub fn yz_plane() -> Self {
        Self {
            normal: *SVector::ith_axis(0),
            distance: T::zero(),
        }
    }
}

impl<T, const N: usize> Plane<T, N>
where
    T: Number,
{
    pub fn new(normal: SVector<T, N>, distance: T) -> Self {
        let len = normal.magnitude();
        if len > T::default_epsilon() {
            Self {
                normal: normal / len,
                distance: distance / len,
            }
        } else {
            Self { normal, distance }
        }
    }

    #[inline]
    pub fn normal(&self) -> SVector<T, N> {
        self.normal
    }

    pub fn with_normal(mut self, normal: SVector<T, N>) -> Self {
        self.normal = normal;
        self
    }

    /// Returns the signed distance from the origin to the plane along the normal.
    #[inline]
    pub fn distance(&self) -> T {
        self.distance
    }

    /// Approximate equality check using epsilon tolerance.
    pub fn roughly_eq(&self, other: &Self) -> bool {
        let eps = T::default_epsilon();
        (self.normal - other.normal).norm() < eps && Real::abs(self.distance - other.distance) < eps
    }

    /// Returns the point on the plane closest to the origin.
    /// Requires the normal to be a unit vector.
    #[inline]
    pub fn origin_projection(&self) -> Point<T, N> {
        (self.normal * self.distance).into()
    }

    #[inline]
    pub fn classify_point(&self, point: &Point<T, N>) -> PlaneSide {
        // Use a geometry-appropriate epsilon. default_epsilon() (~1.19e-7 for f32) is
        // too tight after multiple BSP split/clip cycles where floating-point drift
        // accumulates, causing misclassification and degenerate slivers.
        let eps = T::from(1e-5).expect("BUG: 1e-5 is representable as T");
        let distance = self.normal.dot(&point.coords) - self.distance;
        if distance > eps {
            PlaneSide::Front
        } else if distance < -eps {
            PlaneSide::Back
        } else {
            PlaneSide::On
        }
    }

    #[inline(always)]
    pub fn has_point(&self, point: &Point<T, N>) -> bool {
        self.classify_point(point) == PlaneSide::On
    }

    /// Cast the plane's numeric type from `T` to `T2`.
    pub fn cast<T2: Number>(self) -> Plane<T2, N>
    where
        T: num_traits::cast::AsPrimitive<T2>,
    {
        Plane {
            normal: self.normal.map(|x| x.as_()),
            distance: self.distance.as_(),
        }
    }

    #[inline]
    #[must_use = "flip takes ownership of the original plane and returns a new one with the flipped normal and distance"]
    pub fn flip(mut self) -> Self {
        self.normal *= -T::one();
        self.distance *= -T::one();
        self
    }

    /// Classifies a polygon
    ///
    /// Computes signed distances for all vertices in a single pass, tracking
    /// min/max to derive the classification without per-vertex branching.
    pub fn classify_polygon(&self, polygon: &Polygon<T, N>) -> PolygonClassification {
        let eps = T::from(1e-5).expect("BUG: 1e-5 is representable as T");
        let neg_eps = -eps;
        let mut min_dist =
            T::from(f64::INFINITY).expect("BUG: f64::INFINITY is representable as T");
        let mut max_dist =
            T::from(f64::NEG_INFINITY).expect("BUG: f64::NEG_INFINITY is representable as T");

        for vertex in polygon.vertices().iter() {
            let d = self.normal.dot(&vertex.coords) - self.distance;
            if d < min_dist {
                min_dist = d;
            }
            if d > max_dist {
                max_dist = d;
            }
        }

        match (max_dist > eps, min_dist < neg_eps) {
            (true, true) => PolygonClassification::Spanning,
            (true, false) => PolygonClassification::Front,
            (false, true) => PolygonClassification::Back,
            (false, false) => PolygonClassification::Coplanar,
        }
    }
}

impl<T: Number, const N: usize> Plane<T, N> {
    #[allow(clippy::type_complexity)]
    pub fn split_edge(
        &self,
        a: &Point<T, N>,
        b: &Point<T, N>,
    ) -> Result<(Point<T, N>, Point<T, N>, Point<T, N>), PlaneSide> {
        // Classify both points
        let a_class = self.classify_point(a);
        let b_class = self.classify_point(b);
        if a_class == b_class {
            Err(a_class)?;
        }
        let t = (self.distance() - self.normal().dot(&a.coords)) / self.normal().dot(&(b - a));
        let p = a + (b - a) * t;
        Ok((*a, p, *b))
    }

    /// Split a spanning polygon into a front fragment and a back fragment.
    ///
    /// Algorithm:
    /// 1. Classify each vertex with `classify_point`.
    /// 2. For each consecutive edge (a, b):
    ///    - If `a` is `Front`, add `a` to the front list.
    ///    - If `a` is `Back`, add `a` to the back list.
    ///    - If `a` is `On`, add `a` to both lists.
    ///    - If edge (a, b) crosses the plane (a and b have opposite Front/Back classification):
    ///      compute the intersection point `p = a + t*(b-a)` where `t = (d - n·a) / (n·(b-a))`,
    ///      then add `p` to both lists.
    /// 3. Build two new `Polygon`s from the front and back vertex lists, keeping the original plane.
    ///
    /// Metadata is preserved: both fragments inherit the source polygon's metadata.
    pub fn split_polygon(&self, polygon: Polygon<T, N>) -> (Polygon<T, N>, Polygon<T, N>) {
        let verts = polygon.vertices();
        let n = verts.len();
        let cap = n + 1;
        let mut front = Vec::with_capacity(cap);
        let mut back = Vec::with_capacity(cap);

        // Pre-compute per-vertex classifications to avoid redundant classify_point calls.
        // Each vertex is classified once instead of twice (as this + as next).
        let classes: Vec<PlaneSide> = verts.iter().map(|v| self.classify_point(v)).collect();

        for i in 0..n {
            let j = if i + 1 < n { i + 1 } else { 0 };
            let this = &verts[i];
            let this_class = classes[i];
            let next_class = classes[j];

            match this_class {
                PlaneSide::Front => front.push(*this),
                PlaneSide::Back => back.push(*this),
                PlaneSide::On => {
                    front.push(*this);
                    back.push(*this);
                }
            }
            if this_class != PlaneSide::On
                && next_class != PlaneSide::On
                && this_class != next_class
            {
                let next = &verts[j];
                let t = (self.distance() - self.normal().dot(&this.coords))
                    / self.normal().dot(&(next - this));
                let p = this + (next - this) * t;
                front.push(p);
                back.push(p);
            }
        }

        let a = ndarray::Array1::from_vec(front);
        let b = ndarray::Array1::from_vec(back);
        (polygon.clone_with_points(a), polygon.clone_with_points(b))
    }
}

impl<T> Plane<T, 3>
where
    T: Number + AddAssign + SubAssign + MulAssign + nalgebra::RealField,
{
    /// Construct a plane from three non-collinear points.
    /// The normal is oriented by the right-hand rule: (p2 - p1) × (p3 - p1).
    pub fn from_points(p1: Point3<T>, p2: Point3<T>, p3: Point3<T>) -> Self {
        let v1 = p2 - p1;
        let v2 = p3 - p1;
        let normal = v1.cross(&v2).normalize();
        let distance = normal.dot(&p1.coords);
        Self { normal, distance }
    }

    /// Compute the 3×3 rotation matrix **R** such that `R * self.normal() == other.normal()`.
    ///
    /// In other words, given a point `p` that lies in the plane described by `self`, the product
    /// `R * p` will lie in the plane described by `other`.  This is the unique shortest-arc rotation
    /// (minimum rotation angle) that carries one normal onto the other, computed via the
    /// **Rodrigues rotation formula** in its trigonometry-free form.
    ///
    /// # Mathematical derivation
    ///
    /// Let **a** = `self.normal()` and **b** = `other.normal()` (both assumed to be unit vectors).
    ///
    /// The angle θ between them satisfies:
    /// ```text
    ///   cos θ  = a · b           (dot product)
    ///   sin θ  = ‖a × b‖         (magnitude of cross product)
    /// ```
    ///
    /// The rotation axis is `v = a × b` (unnormalized; its magnitude is sin θ).
    ///
    /// The skew-symmetric cross-product matrix `[v]×` encodes the "cross with v" operation:
    /// ```text
    ///         ⎡  0   −v₃   v₂ ⎤
    /// [v]× =  ⎢  v₃   0   −v₁ ⎥
    ///         ⎣ −v₂   v₁   0  ⎦
    /// ```
    ///
    /// Rodrigues' formula then gives:
    /// ```text
    /// R = I  +  [v]×  +  [v]ײ  ·  1/(1 + cos θ)
    /// ```
    ///
    /// This avoids calling `acos`/`sin`/`cos` and is numerically stable for all θ except
    /// θ = 180° (anti-parallel normals), which is handled separately below.
    ///
    /// # Special case: anti-parallel normals (θ ≈ 180°)
    ///
    /// When `a · b ≈ −1` the denominator `(1 + cos θ)` approaches zero, so the formula above
    /// is undefined.  Geometrically, any 180° rotation whose axis is perpendicular to **a** is a
    /// valid solution; we break the degeneracy by choosing the world axis that is *least* aligned
    /// with **a** (found with `iamin()` — the index of the smallest absolute component) as a
    /// helper to build a perpendicular rotation axis:
    /// ```text
    /// rot_axis = (a × world_axis).normalize()
    /// R = −I + 2 · rot_axis · rot_axisᵀ
    /// ```
    /// (This is Rodrigues with θ = π: cos π = −1, sin π = 0.)
    pub fn rotation_to(&self, other: Plane<T, 3>) -> nalgebra::Matrix3<T>
    where
        T: num_traits::Signed,
    {
        // Unit normals of the source and destination planes.
        let a = self.normal().normalize();
        let b = other.normal().normalize();

        // cos θ = a · b  and  v = a × b  (unnormalized rotation axis, ‖v‖ = sin θ).
        let cos_theta = a.dot(&b);
        let v = a.cross(&b);

        let one = T::one();
        let two = one + one;

        // -----------------------------------------------------------------------
        // Anti-parallel special case: a · b ≈ −1  (θ ≈ 180°).
        //
        // The general Rodrigues formula below scales `[v]ײ` by 1/(1 + cos θ). As the
        // normals approach anti-parallel, (1 + cos θ) → 0, so that factor explodes and
        // amplifies the rounding noise in `v = a × b` (which is itself ≈ 0 here). Writing
        // δ for the angular deviation from exactly 180°, 1 + cos θ ≈ δ²/2, so the matrix
        // error grows like ε/δ². The old cutoff only kicked in for δ ≲ 3e-8 rad, leaving
        // a wide band of nearly-anti-parallel inputs to fall through to the unstable
        // formula and produce garbage matrices (non-orthogonal, det ≠ 1).
        //
        // We therefore branch to the exact 180° construction long before the denominator
        // becomes dangerous. The `1e-6` threshold corresponds to δ ≈ 1.4e-3 rad, where the
        // Rodrigues term still carries only ~1e-10 error — comfortably below the crossover
        // while covering every input the formula cannot handle cleanly.
        //
        // Any 180° rotation whose axis is perpendicular to `a` sends a → −a ≈ b, so we
        // build one explicitly. `iamin()` returns the index of the component of `a` with
        // the smallest absolute value, i.e. the world axis (X=0, Y=1, Z=2) that is *least*
        // parallel to `a`; crossing `a` with it yields a well-conditioned perpendicular.
        // -----------------------------------------------------------------------
        let anti_parallel_threshold = T::from(1e-6).expect("BUG: 1e-6 is representable as T");
        if cos_theta < -one + anti_parallel_threshold {
            // Pick the world axis least aligned with `a` to avoid a near-zero cross product.
            let least_aligned_axis = a.iamin();
            let world_axis = SVector::<T, 3>::ith_axis(least_aligned_axis);

            // rot_axis is perpendicular to `a` and unit-length.
            let rot_axis = a.cross(&world_axis.into_inner()).normalize();

            // 180° Rodrigues: R = −I + 2 * rot_axis * rot_axisᵀ
            // (cos π = −1 cancels the identity, sin π = 0 removes the [v]× term)
            return -nalgebra::Matrix3::identity() + rot_axis * rot_axis.transpose() * two;
        }

        // -----------------------------------------------------------------------
        // General case via the trigonometry-free Rodrigues formula:
        //
        //   R = I  +  [v]×  +  [v]ײ  ·  1 / (1 + cos θ)
        //
        // where [v]× is the skew-symmetric matrix of v = a × b.
        // -----------------------------------------------------------------------

        // Build the skew-symmetric cross-product matrix [v]×.
        let vx = nalgebra::Matrix3::new(
            T::zero(),
            -v[2],
            v[1],
            v[2],
            T::zero(),
            -v[0],
            -v[1],
            v[0],
            T::zero(),
        );

        // R = I + [v]× + [v]ײ / (1 + cos θ)
        nalgebra::Matrix3::identity() + vx + vx * vx * (one / (one + cos_theta))
    }
}

#[cfg(test)]
mod tests {
    use super::*;

    #[test]
    fn classify_point_front() {
        let plane = Plane::xy_plane();
        let point = Point::from([0f32, 0f32, 1f32]);
        assert_eq!(plane.classify_point(&point), PlaneSide::Front);
    }

    #[test]
    fn classify_point_back() {
        let plane = Plane::xy_plane();
        let point = Point::from([0f32, 0f32, -1f32]);
        assert_eq!(plane.classify_point(&point), PlaneSide::Back);
    }

    #[test]
    fn classify_point_on_plane() {
        let plane = Plane::xy_plane();
        let point = Point::from([5f32, -3f32, 0f32]);
        assert_eq!(plane.classify_point(&point), PlaneSide::On);
    }

    #[test]
    fn classify_point_offset_plane_front() {
        let normal = SVector::from([0f32, 0f32, 1f32]);
        let plane = Plane::new(normal, 5f32);
        let point = Point::from([10f32, 20f32, 6f32]);
        assert_eq!(plane.classify_point(&point), PlaneSide::Front);
    }

    #[test]
    fn classify_point_offset_plane_back() {
        let normal = SVector::from([0f32, 0f32, 1f32]);
        let plane = Plane::new(normal, 5f32);
        let point = Point::from([10f32, 20f32, 4f32]);
        assert_eq!(plane.classify_point(&point), PlaneSide::Back);
    }

    #[test]
    fn classify_point_offset_plane_on() {
        let normal = SVector::from([0f32, 0f32, 1f32]);
        let plane = Plane::new(normal, 5f32);
        let point = Point::from([-10f32, 30f32, 5f32]);
        assert_eq!(plane.classify_point(&point), PlaneSide::On);
    }

    #[test]
    fn edge_both_front() {
        let plane = Plane::xy_plane();
        let a = Point::from([0f32, 0f32, 5f32]);
        let b = Point::from([10f32, 10f32, 3f32]);
        assert_eq!(plane.split_edge(&a, &b), Err(PlaneSide::Front));
    }

    #[test]
    fn edge_both_back() {
        let plane = Plane::xy_plane();
        let a = Point::from([0f32, 0f32, -5f32]);
        let b = Point::from([10f32, 10f32, -3f32]);
        assert_eq!(plane.split_edge(&a, &b), Err(PlaneSide::Back));
    }

    #[test]
    fn edge_both_on_plane() {
        let plane = Plane::xy_plane();
        let a = Point::from([10f32, 10f32, 0f32]);
        let b = Point::from([-10f32, -10f32, 0f32]);
        assert_eq!(plane.split_edge(&a, &b), Err(PlaneSide::On));
    }

    #[test]
    fn edge_crossing_front_to_back() {
        let plane = Plane::xy_plane();
        let a = Point::from([0f32, 0f32, 2f32]);
        let b = Point::from([0f32, 0f32, -2f32]);
        let result = plane.split_edge(&a, &b);
        assert!(result.is_ok());
        let (pa, intersection, pb) = result.unwrap();
        assert_eq!(pa, a);
        assert_eq!(pb, b);
        assert_eq!(intersection, Point::from([0f32, 0f32, 0f32]));
    }

    #[test]
    fn edge_crossing_back_to_front() {
        let plane = Plane::xy_plane();
        let a = Point::from([0f32, 0f32, -2f32]);
        let b = Point::from([0f32, 0f32, 2f32]);
        let result = plane.split_edge(&a, &b);
        assert!(result.is_ok());
        let (pa, intersection, pb) = result.unwrap();
        assert_eq!(pa, a);
        assert_eq!(pb, b);
        assert_eq!(intersection, Point::from([0f32, 0f32, 0f32]));
    }

    #[test]
    fn edge_crossing_offset_plane() {
        let normal = SVector::from([0f32, 0f32, 1f32]);
        let plane = Plane::new(normal, 3f32);
        let a = Point::from([1f32, 2f32, 1f32]);
        let b = Point::from([1f32, 2f32, 5f32]);
        let result = plane.split_edge(&a, &b);
        assert!(result.is_ok());
        let (_, intersection, _) = result.unwrap();
        assert_eq!(intersection, Point::from([1f32, 2f32, 3f32]));
    }

    #[test]
    fn edge_crossing_offset_plane_negative_distance() {
        let normal = SVector::from([0f32, 0f32, 1f32]);
        let plane = Plane::new(normal, -2f32);
        let a = Point::from([5f32, 7f32, -5f32]);
        let b = Point::from([5f32, 7f32, 1f32]);
        let result = plane.split_edge(&a, &b);
        assert!(result.is_ok());
        let (_, intersection, _) = result.unwrap();
        assert_eq!(intersection, Point::from([5f32, 7f32, -2f32]));
    }

    #[test]
    fn edge_crossing_positive_x_plane() {
        let normal = SVector::from([1f32, 0f32, 0f32]);
        let plane = Plane::new(normal, 4f32);
        let a = Point::from([2f32, 3f32, 5f32]);
        let b = Point::from([8f32, 3f32, 5f32]);
        let result = plane.split_edge(&a, &b);
        assert!(result.is_ok());
        let (_, intersection, _) = result.unwrap();
        assert_eq!(intersection, Point::from([4f32, 3f32, 5f32]));
    }

    #[test]
    fn edge_crossing_positive_y_plane() {
        let normal = SVector::from([0f32, 1f32, 0f32]);
        let plane = Plane::new(normal, -6f32);
        let a = Point::from([1f32, -10f32, 2f32]);
        let b = Point::from([1f32, 0f32, 2f32]);
        let result = plane.split_edge(&a, &b);
        assert!(result.is_ok());
        let (_, intersection, _) = result.unwrap();
        assert_eq!(intersection, Point::from([1f32, -6f32, 2f32]));
    }

    #[test]
    fn edge_crossing_diagonal_plane() {
        let normal = SVector::from([1f32, 0f32, 0f32]);
        let plane = Plane::new(normal, 0f32);
        let a = Point::from([-2f32, 3f32, 7f32]);
        let b = Point::from([4f32, 3f32, 7f32]);
        let result = plane.split_edge(&a, &b);
        assert!(result.is_ok());
        let (_, intersection, _) = result.unwrap();
        assert_eq!(intersection, Point::from([0f32, 3f32, 7f32]));
    }

    #[test]
    fn edge_crossing_arbitrary_points() {
        let normal = SVector::from([0f32, 0f32, 1f32]);
        let plane = Plane::new(normal, 10f32);
        let a = Point::from([100f32, 200f32, 5f32]);
        let b = Point::from([100f32, 200f32, 20f32]);
        let result = plane.split_edge(&a, &b);
        assert!(result.is_ok());
        let (_, intersection, _) = result.unwrap();
        assert_eq!(intersection, Point::from([100f32, 200f32, 10f32]));
    }

    #[test]
    fn edge_crossing_slanted_edge() {
        let normal = SVector::from([0f32, 0f32, 1f32]);
        let plane = Plane::new(normal, 1f32);
        let a = Point::from([0f32, 0f32, 0f32]);
        let b = Point::from([10f32, 20f32, 2f32]);
        let result = plane.split_edge(&a, &b);
        assert!(result.is_ok());
        let (_, intersection, _) = result.unwrap();
        assert_eq!(intersection.x, 5f32);
        assert_eq!(intersection.y, 10f32);
        assert_eq!(intersection.z, 1f32);
    }

    #[test]
    fn edge_crossing_asymmetric() {
        let normal = SVector::from([0f32, 0f32, 1f32]);
        let plane = Plane::new(normal, 5f32);
        let a = Point::from([2f32, -3f32, 2f32]);
        let b = Point::from([2f32, -3f32, 8f32]);
        let result = plane.split_edge(&a, &b);
        assert!(result.is_ok());
        let (_, intersection, _) = result.unwrap();
        assert_eq!(intersection, Point::from([2f32, -3f32, 5f32]));
    }

    #[test]
    fn edge_crossing_diagonal() {
        let plane = Plane::xy_plane();
        let a = Point::from([0f32, 0f32, 1f32]);
        let b = Point::from([0f32, 0f32, -1f32]);
        let result = plane.split_edge(&a, &b);
        assert!(result.is_ok());
        let (_, intersection, _) = result.unwrap();
        assert_eq!(intersection.z, 0f32);
    }

    #[test]
    fn edge_one_on_plane_with_front() {
        let plane = Plane::xy_plane();
        let a = Point::from([0f32, 0f32, 0f32]);
        let b = Point::from([0f32, 0f32, 5f32]);
        let result = plane.split_edge(&a, &b);
        assert!(result.is_ok());
        let (_, intersection, _) = result.unwrap();
        assert_eq!(intersection, Point::from([0f32, 0f32, 0f32]));
    }

    #[test]
    fn edge_one_on_plane_with_front_reversed() {
        let plane = Plane::xy_plane();
        let a = Point::from([0f32, 0f32, 5f32]);
        let b = Point::from([0f32, 0f32, 0f32]);
        let result = plane.split_edge(&a, &b);
        assert!(result.is_ok());
        let (_, intersection, _) = result.unwrap();
        assert_eq!(intersection, Point::from([0f32, 0f32, 0f32]));
    }

    #[test]
    fn edge_one_on_plane_with_back() {
        let plane = Plane::xy_plane();
        let a = Point::from([0f32, 0f32, 0f32]);
        let b = Point::from([0f32, 0f32, -5f32]);
        let result = plane.split_edge(&a, &b);
        assert!(result.is_ok());
        let (_, intersection, _) = result.unwrap();
        assert_eq!(intersection, Point::from([0f32, 0f32, 0f32]));
    }

    #[test]
    fn edge_one_on_plane_with_back_reversed() {
        let plane = Plane::xy_plane();
        let a = Point::from([0f32, 0f32, -5f32]);
        let b = Point::from([0f32, 0f32, 0f32]);
        let result = plane.split_edge(&a, &b);
        assert!(result.is_ok());
        let (_, intersection, _) = result.unwrap();
        assert_eq!(intersection, Point::from([0f32, 0f32, 0f32]));
    }

    #[test]
    fn split_polygon_triangle_spanning_xy_plane() {
        let plane = Plane::xy_plane();
        let polygon = Polygon::from_points(vec![
            Point3::new(0f32, 0f32, 1f32),
            Point3::new(1f32, 0f32, 1f32),
            Point3::new(0.5f32, 1f32, -1f32),
        ])
        .unwrap();
        assert_eq!(
            plane.classify_polygon(&polygon),
            PolygonClassification::Spanning
        );
        let (front, back) = plane.split_polygon(polygon);
        assert!(front.vertices().len() >= 3);
        assert!(back.vertices().len() >= 3);
        for v in front.vertices().iter() {
            assert!(plane.classify_point(v) != PlaneSide::Back);
        }
        for v in back.vertices().iter() {
            assert!(plane.classify_point(v) != PlaneSide::Front);
        }
    }

    #[test]
    fn split_polygon_square_spanning_xy_plane() {
        let plane = Plane::xy_plane();
        let polygon: Polygon<f32, 3> = Polygon::from_points(vec![
            Point3::new(-1f32, -1f32, 1f32),
            Point3::new(1f32, -1f32, 1f32),
            Point3::new(1f32, 1f32, -1f32),
            Point3::new(-1f32, 1f32, -1f32),
        ])
        .unwrap();
        assert_eq!(
            plane.classify_polygon(&polygon),
            PolygonClassification::Spanning
        );
        let (front, back) = plane.split_polygon(polygon);
        assert!(front.vertices().len() >= 3);
        assert!(back.vertices().len() >= 3);
        for v in front.vertices().iter() {
            assert!(plane.classify_point(v) != PlaneSide::Back);
        }
        for v in back.vertices().iter() {
            assert!(plane.classify_point(v) != PlaneSide::Front);
        }
    }

    #[test]
    fn split_polygon_vertex_counts_sum() {
        let plane = Plane::xy_plane();
        let polygon: Polygon<f32, 3> = Polygon::from_points(vec![
            Point3::new(0f32, 0f32, 2f32),
            Point3::new(2f32, 0f32, 2f32),
            Point3::new(2f32, 2f32, -2f32),
            Point3::new(0f32, 2f32, -2f32),
        ])
        .unwrap();
        let (front, back) = plane.split_polygon(polygon);
        assert_eq!(front.vertices().len() + back.vertices().len(), 8);
    }

    #[test]
    fn split_polygon_offset_plane() {
        let plane = Plane::new(SVector::from([0f32, 0f32, 1f32]), 1f32);
        let polygon: Polygon<f32, 3> = Polygon::from_points(vec![
            Point3::new(0f32, 0f32, 0f32),
            Point3::new(2f32, 0f32, 0f32),
            Point3::new(2f32, 2f32, 2f32),
            Point3::new(0f32, 2f32, 2f32),
        ])
        .unwrap();
        let (front, back) = plane.split_polygon(polygon);
        for v in front.vertices().iter() {
            assert!(plane.classify_point(v) != PlaneSide::Back);
        }
        for v in back.vertices().iter() {
            assert!(plane.classify_point(v) != PlaneSide::Front);
        }
    }

    /// Assert that `rotation_to` produced a proper rotation matrix `R` (`RᵀR = I`,
    /// `det R = 1`) that carries the source normal `a` onto the target normal `b`
    /// within `map_tol`.
    fn assert_valid_rotation(a: SVector<f64, 3>, b: SVector<f64, 3>, map_tol: f64) {
        let src = Plane::new(a, 0.0);
        let dst = Plane::new(b, 0.0);
        let r = src.rotation_to(dst);

        let ortho_err = (r.transpose() * r - nalgebra::Matrix3::<f64>::identity()).norm();
        assert!(
            ortho_err < 1e-9,
            "‖RᵀR−I‖ = {ortho_err} (a = {a:?}, b = {b:?})"
        );

        let det = r.determinant();
        assert!(
            (det - 1.0).abs() < 1e-9,
            "det(R) = {det} (a = {a:?}, b = {b:?})"
        );

        let map_err = (r * src.normal() - dst.normal()).norm();
        assert!(
            map_err < map_tol,
            "‖R·a − b‖ = {map_err} exceeds {map_tol} (a = {a:?}, b = {b:?})"
        );
    }

    #[test]
    fn rotation_to_generic_pair() {
        // A plain 90° rotation must map the source normal onto the target exactly.
        let a = SVector::from([1.0, 0.0, 0.0]);
        let b = SVector::from([0.0, 1.0, 0.0]);
        assert_valid_rotation(a, b, 1e-12);
    }

    #[test]
    fn rotation_to_near_parallel_is_identity() {
        // cos θ ≈ +1: the general Rodrigues path must still yield a valid rotation
        // (≈ identity) and not be disturbed by the anti-parallel widening.
        let a = SVector::from([0.3, -0.5, 0.8]).normalize();
        let perp = a.cross(&SVector::from([1.0, 0.0, 0.0])).normalize();
        let b = (a * 1e-7_f64.cos() + perp * 1e-7_f64.sin()).normalize();
        assert_valid_rotation(a, b, 1e-6);
    }

    #[test]
    fn rotation_to_near_and_exact_anti_parallel() {
        // Base normal (off-axis so no component is trivially zero) plus a unit vector
        // perpendicular to it, used to tilt `-a` by a tiny angle δ away from exactly
        // anti-parallel. Before the fix, every δ below reached the unstable
        // 1/(1 + cos θ) branch and produced a non-orthogonal matrix (det ≠ 1).
        let a = SVector::from([0.3, -0.5, 0.8]).normalize();
        let perp = a.cross(&SVector::from([1.0, 0.0, 0.0])).normalize();

        // δ = 0 (exactly anti-parallel) via cos 0 = 1, sin 0 = 0, plus the pathological
        // deltas flagged in issue #100.
        for &delta in &[0.0_f64, 5e-8, 1e-7, 1e-6] {
            let b = (-a * delta.cos() + perp * delta.sin()).normalize();
            // R maps a to exactly −a, which differs from b by ~δ, so the mapping
            // tolerance absorbs δ; orthogonality and det stay pinned at 1e-9.
            assert_valid_rotation(a, b, 1e-5);
        }
    }
}