math_utils/geometry/primitive/

mod.rs

//! Affine, convex, and combinatorial spaces

use rand;

use crate::*;
use super::{intersect, shape};
use super::mesh::VertexEdgeTriangleMesh;

#[cfg(feature = "derive_serdes")]
use serde::{Deserialize, Serialize};

/// Primitives over signed integers
pub mod integer;

mod hull;
mod simplex;

pub use self::hull::*;
pub use self::simplex::{simplex2, simplex3, Simplex2, Simplex3};
pub use simplex2::{
  Segment       as Segment2,
  Triangle      as Triangle2,
  SegmentPoint  as Segment2Point,
  TrianglePoint as Triangle2Point
};
pub use simplex3::{
  Segment          as Segment3,
  Triangle         as Triangle3,
  Tetrahedron      as Tetrahedron3,
  SegmentPoint     as Segment3Point,
  TrianglePoint    as Triangle3Point,
  TetrahedronPoint as Tetrahedron3Point
};

/// A parameterized point on a 2D line
pub type Line2Point <S> = (S, Point2 <S>);
/// A parameterized point on a 3D line
pub type Line3Point <S> = (S, Point3 <S>);
/// A parameterized point on a 2D ray
pub type Ray2Point <S> = (NonNegative <S>, Point2 <S>);
/// A parameterized point on a 3D ray
pub type Ray3Point <S> = (NonNegative <S>, Point3 <S>);

/// Geometric primitives generic over dimension
pub trait Primitive <S, P> : std::fmt::Debug where
  S : Field,
  P : AffineSpace <S>
{
  fn translate (&mut self, displacement : P::Translation);
  fn scale     (&mut self, scale : Positive <S>);
  /// Returns the point with largeest dot product in search direction
  fn support   (&self, direction : <P::Translation as VectorSpace <S>>::NonZero)
    -> (P, S);
}

////////////////////////////////////////////////////////////////////////////////////////
//  structs                                                                           //
////////////////////////////////////////////////////////////////////////////////////////

/// 1D interval. Points are allowed to be equal.
#[cfg_attr(feature = "derive_serdes", derive(Serialize, Deserialize))]
#[derive(Clone, Copy, Debug, Eq, PartialEq)]
pub struct Interval <S> {
  min : S,
  max : S
}

/// 2D axis-aligned bounding box.
///
/// Min/max can be co-linear:
/// ```
/// # use math_utils::geometry::Aabb2;
/// let x = Aabb2::with_minmax (
///   [0.0, 0.0].into(),
///   [0.0, 1.0].into()
/// ).unwrap();
/// ```
/// But not identical:
/// ```should_panic
/// # use math_utils::geometry::Aabb2;
/// // panic!
/// let x = Aabb2::with_minmax (
///   [0.0, 0.0].into(),
///   [0.0, 0.0].into()
/// ).unwrap();
/// ```
#[cfg_attr(feature = "derive_serdes", derive(Serialize, Deserialize))]
#[derive(Clone, Copy, Debug, Eq, PartialEq)]
pub struct Aabb2 <S> {
  min : Point2 <S>,
  max : Point2 <S>
}

/// 3D axis-aligned bounding box.
///
/// See also `shape::Aabb` trait for primitive and shape types that can compute a 3D
/// AABB.
///
/// Min/max can be co-planar or co-linear:
/// ```
/// # use math_utils::geometry::Aabb3;
/// let x = Aabb3::with_minmax (
///   [0.0, 0.0, 0.0].into(),
///   [0.0, 1.0, 1.0].into()
/// ).unwrap();
/// let y = Aabb3::with_minmax (
///   [0.0, 0.0, 0.0].into(),
///   [0.0, 1.0, 0.0].into()
/// ).unwrap();
/// ```
/// But not identical:
/// ```should_panic
/// # use math_utils::geometry::Aabb3;
/// let x = Aabb3::with_minmax (
///   [0.0, 0.0, 0.0].into(),
///   [0.0, 0.0, 0.0].into()
/// ).unwrap();
/// ```
#[cfg_attr(feature = "derive_serdes", derive(Serialize, Deserialize))]
#[derive(Clone, Copy, Debug, Eq, PartialEq)]
pub struct Aabb3 <S> {
  min : Point3 <S>,
  max : Point3 <S>
}

/// 2D oriented bounding box
#[cfg_attr(feature = "derive_serdes", derive(Serialize, Deserialize))]
#[derive(Clone, Copy, Debug, Eq, PartialEq)]
pub struct Obb2 <S> {
  pub center      : Point2 <S>,
  pub extents     : Vector2 <Positive <S>>,
  pub orientation : AngleWrapped <S>
}

/// 3D oriented bounding box
#[cfg_attr(feature = "derive_serdes", derive(Serialize, Deserialize))]
#[derive(Clone, Copy, Debug, Eq, PartialEq)]
pub struct Obb3 <S> {
  pub center      : Point3 <S>,
  pub extents     : Vector3 <Positive <S>>,
  pub orientation : Angles3 <S>
}

/// 3D Z-axis-aligned cylinder
#[cfg_attr(feature = "derive_serdes", derive(Serialize, Deserialize))]
#[derive(Clone, Copy, Debug, Eq, PartialEq)]
pub struct Cylinder3 <S> {
  pub center      : Point3   <S>,
  pub half_height : Positive <S>,
  pub radius      : Positive <S>
}

/// 3D Z-axis-aligned capsule
#[cfg_attr(feature = "derive_serdes", derive(Serialize, Deserialize))]
#[derive(Clone, Copy, Debug, Eq, PartialEq)]
pub struct Capsule3 <S> {
  pub center      : Point3   <S>,
  pub half_height : Positive <S>,
  pub radius      : Positive <S>
}

/// 3D Z-axis-aligned cone
#[cfg_attr(feature = "derive_serdes", derive(Serialize, Deserialize))]
#[derive(Clone, Copy, Debug, Eq, PartialEq)]
pub struct Cone3 <S> {
  pub center      : Point3   <S>,
  pub half_height : Positive <S>,
  pub radius      : Positive <S>
}

/// An infinitely extended line in 2D space defined by a base point and normalized
/// direction
#[cfg_attr(feature = "derive_serdes", derive(Serialize, Deserialize))]
#[derive(Clone, Copy, Debug, Eq, PartialEq)]
pub struct Line2 <S> {
  pub base      : Point2 <S>,
  pub direction : Unit2  <S>
}

/// An infinitely extended ray in 2D space defined by a base point and normalized
/// direction
#[cfg_attr(feature = "derive_serdes", derive(Serialize, Deserialize))]
#[derive(Clone, Copy, Debug, Eq, PartialEq)]
pub struct Ray2 <S> {
  pub base      : Point2 <S>,
  pub direction : Unit2  <S>
}

/// An infinitely extended line in 3D space defined by a base point and normalized
/// direction
// TODO: currently most algorithms for lines are written to work with a line defined by
// a line segment (no need to normalize); can we create another line type to improve
// type checking?
#[cfg_attr(feature = "derive_serdes", derive(Serialize, Deserialize))]
#[derive(Clone, Copy, Debug, Eq, PartialEq)]
pub struct Line3 <S> {
  pub base      : Point3 <S>,
  pub direction : Unit3  <S>
}

/// An infinitely extended ray in 3D space defined by a base point and normalized
/// direction
#[cfg_attr(feature = "derive_serdes", derive(Serialize, Deserialize))]
#[derive(Clone, Copy, Debug, Eq, PartialEq)]
pub struct Ray3 <S> {
  pub base      : Point3 <S>,
  pub direction : Unit3  <S>
}

/// A plane in 3D space defined by a base point and (unit) normal vector
#[cfg_attr(feature = "derive_serdes", derive(Serialize, Deserialize))]
#[derive(Clone, Copy, Debug, Eq, PartialEq)]
pub struct Plane3 <S> {
  pub base   : Point3 <S>,
  pub normal : Unit3  <S>
}

/// Sphere in 2D space (a circle)
#[cfg_attr(feature = "derive_serdes", derive(Serialize, Deserialize))]
#[derive(Clone, Copy, Debug, Eq, PartialEq)]
pub struct Sphere2 <S> {
  pub center : Point2   <S>,
  pub radius : Positive <S>
}

/// Sphere in 3D space
#[cfg_attr(feature = "derive_serdes", derive(Serialize, Deserialize))]
#[derive(Clone, Copy, Debug, Eq, PartialEq)]
pub struct Sphere3 <S> {
  pub center : Point3   <S>,
  pub radius : Positive <S>
}

/// Describes points of a 2D convex hull that support a rectangle and the area of the
/// rectangle
#[derive(Clone, Debug)]
struct Rect <S> {
  pub edge    : Vector2 <S>,
  pub perp    : Vector2 <S>,
  /// bottom, right, top, left
  pub indices : [u32; 4],
  pub area    : S,
  pub edge_len_squared : S
}

////////////////////////////////////////////////////////////////////////////////////////
//  functions                                                                         //
////////////////////////////////////////////////////////////////////////////////////////

/// Returns true when three points lie on the same line in 2D space.
///
/// ```
/// # use math_utils::geometry::colinear_2d;
/// assert!(colinear_2d (
///   [-1.0, -1.0].into(),
///   [ 0.0,  0.0].into(),
///   [ 1.0,  1.0].into())
/// );
/// assert!(!colinear_2d (
///   [-1.0, -1.0].into(),
///   [ 0.0,  1.0].into(),
///   [ 1.0, -1.0].into())
/// );
/// ```
pub fn colinear_2d <S> (a : Point2 <S>, b : Point2 <S>, c : Point2 <S>) -> bool where
  S : OrderedField + approx::AbsDiffEq <Epsilon = S>
{
  use num::Zero;
  let max_point_max_norm = a.0.norm_max().max (b.0.norm_max()).max (c.0.norm_max());
  let min_edge_max_norm  = (b - a).norm_max().min ((c - a).norm_max());
  let ratio = if min_edge_max_norm == NonNegative::zero() {
    NonNegative::zero()
  } else {
    max_point_max_norm / min_edge_max_norm
  };
  let epsilon = S::default_epsilon() * S::eight() * (S::one() + *ratio);
  colinear_2d_epsilon (a, b, c, epsilon)
}

pub fn colinear_2d_epsilon <S> (
  a : Point2 <S>, b : Point2 <S>, c : Point2 <S>, epsilon : S
) -> bool where S : OrderedField {
  // early exit: a pair of points are equal
  if a == b || a == c || b == c {
    return true
  }
  let e1 = b - a;
  let e2 = c - a;
  // pre-scale by max norm: components will lie in [-1, 1]
  let e1_max = e1.norm_max();
  let e2_max = e2.norm_max();
  if *e1_max == S::zero() || *e2_max == S::zero() {
    return true
  }
  let u        = e1 / *e1_max;
  let v        = e2 / *e2_max;
  let u_mag_sq = u.magnitude_squared();
  let v_mag_sq = v.magnitude_squared();
  let cross    = u.exterior_product (v);
  cross * cross <= epsilon * epsilon * u_mag_sq * v_mag_sq
}

/// Returns true when three points lie on the same line in 3D space.
///
/// ```
/// # use math_utils::geometry::colinear_3d;
/// assert!(colinear_3d (
///   [-1.0, -1.0, -1.0].into(),
///   [ 0.0,  0.0,  0.0].into(),
///   [ 1.0,  1.0,  1.0].into())
/// );
/// ```
pub fn colinear_3d <S> (a : Point3 <S>, b : Point3 <S>, c : Point3 <S>) -> bool where
  S : OrderedField + approx::AbsDiffEq <Epsilon = S>
{
  use num::Zero;
  let max_epsilon = S::two().powi (-3);
  let max_point_max_norm = a.0.norm_max().max (b.0.norm_max()).max (c.0.norm_max());
  let min_edge_max_norm  = (b - a).norm_max().min ((c - a).norm_max());
  let ratio = if min_edge_max_norm == NonNegative::zero() {
    NonNegative::zero()
  } else {
    max_point_max_norm / min_edge_max_norm
  };
  let epsilon = (S::default_epsilon() * S::eight() * (S::one() + *ratio))
    .min (max_epsilon);
  colinear_3d_epsilon (a, b, c, epsilon)
}

pub fn colinear_3d_epsilon <S> (
  a : Point3 <S>, b : Point3 <S>, c : Point3 <S>, epsilon : S
) -> bool where
  S : OrderedField
{
  // early exit: a pair of points are equal
  if a == b || a == c || b == c {
    return true
  }
  let e1 = b - a;
  let e2 = c - a;
  // pre-scale by max norm: components will lie in [-1, 1]
  let e1_max = e1.norm_max();
  let e2_max = e2.norm_max();
  if *e1_max == S::zero() || *e2_max == S::zero() {
    return true
  }
  let u        = e1 / *e1_max;
  let v        = e2 / *e2_max;
  let u_mag_sq = u.magnitude_squared();
  let v_mag_sq = v.magnitude_squared();
  let cross_sq = u.cross (v).magnitude_squared();
  cross_sq <= epsilon * epsilon * u_mag_sq * v_mag_sq
}

/// Returns true when four points lie on the same plane in 3D space.
///
/// ```
/// # use math_utils::geometry::coplanar_3d;
/// assert!(coplanar_3d (
///   [-1.0, -1.0, -1.0].into(),
///   [ 1.0,  1.0,  1.0].into(),
///   [-1.0,  1.0,  0.0].into(),
///   [ 1.0, -1.0,  0.0].into()
/// ));
/// ```
pub fn coplanar_3d <S> (a : Point3 <S>, b : Point3 <S>, c : Point3 <S>, d : Point3 <S>)
  -> bool
where S : OrderedField + approx::AbsDiffEq <Epsilon = S> {
  use num::Zero;
  let max_epsilon = S::two().powi (-2);
  let max_point_max_norm =
    a.0.norm_max().max (b.0.norm_max()).max (c.0.norm_max()).max (d.0.norm_max());
  let min_edge_max_norm  =
    (b - a).norm_max().min ((c - a).norm_max()).min ((d - a).norm_max());
  let ratio = if min_edge_max_norm == NonNegative::zero() {
    NonNegative::zero()
  } else {
    max_point_max_norm / min_edge_max_norm
  };
  let epsilon = (S::default_epsilon() * S::eight() * (S::one() + *ratio))
    .min (max_epsilon);
  coplanar_3d_epsilon (a, b, c, d, epsilon)
}

pub fn coplanar_3d_epsilon <S> (
  a : Point3 <S>, b : Point3 <S>, c : Point3 <S>, d : Point3 <S>, epsilon : S
) -> bool where
  S : OrderedField + approx::AbsDiffEq <Epsilon = S>
{
  let ab = b - a;
  let ac = c - a;
  let ad = d - a;
  let ab_max_norm = ab.norm_max();
  let ac_max_norm = ac.norm_max();
  let ad_max_norm = ad.norm_max();
  if *ab_max_norm == S::zero() || *ac_max_norm == S::zero()
    || *ad_max_norm == S::zero()
  {
    // early exit: a pair of points are equal
    return true
  }
  let u = ab / *ab_max_norm;
  let v = ac / *ac_max_norm;
  let w = ad / *ad_max_norm;
  let n = u.cross (v);
  let n_sq = n.magnitude_squared();
  if n_sq == S::zero() {
    // a, b, c are colinear
    return true
  }
  let tp   = n.dot (w);
  let w_sq = w.magnitude_squared();
  tp * tp <= epsilon * epsilon * n_sq * w_sq
}

/// Return a signum indicating which side of the triangle ABC that point D lies on,
/// where +1 indicates that D lies on the side that the normal $(B - A) \times (C - A)$
/// points to, -1 if D lies on the opposite side, and 0 if D lies on the plane ABC,
/// within tolerance
pub fn plane_side_3d <S> (a : Point3 <S>, b : Point3 <S>, c : Point3 <S>, d : Point3 <S>)
  -> S
where S : OrderedField + approx::AbsDiffEq <Epsilon = S> {
  use num::Zero;
  let max_epsilon = S::two().powi (-2);
  let max_point_max_norm =
    a.0.norm_max().max (b.0.norm_max()).max (c.0.norm_max()).max (d.0.norm_max());
  let min_edge_max_norm  =
    (b - a).norm_max().min ((c - a).norm_max()).min ((d - a).norm_max());
  let ratio = if min_edge_max_norm == NonNegative::zero() {
    NonNegative::zero()
  } else {
    max_point_max_norm / min_edge_max_norm
  };
  let epsilon = (S::default_epsilon() * S::eight() * (S::one() + *ratio))
    .min (max_epsilon);
  plane_side_3d_epsilon (a, b, c, d, epsilon)
}

pub fn plane_side_3d_epsilon <S> (
  a : Point3 <S>, b : Point3 <S>, c : Point3 <S>, d : Point3 <S>, epsilon : S
) -> S where
  S : OrderedField + approx::AbsDiffEq <Epsilon = S>
{
  let ab = b - a;
  let ac = c - a;
  let ad = d - a;
  let ab_max_norm = ab.norm_max();
  let ac_max_norm = ac.norm_max();
  let ad_max_norm = ad.norm_max();
  if *ab_max_norm == S::zero() || *ac_max_norm == S::zero()
    || *ad_max_norm == S::zero()
  {
    // a pair of points are equal => degenerate
    return S::zero()
  }
  let u = ab / *ab_max_norm;
  let v = ac / *ac_max_norm;
  let w = ad / *ad_max_norm;
  let n = u.cross (v);
  let n_sq = n.magnitude_squared();
  if n_sq == S::zero() {
    // a, b, c are colinear
    return S::zero()
  }
  let tp   = n.dot (w);
  let w_sq = w.magnitude_squared();
  if tp * tp <= epsilon * epsilon * n_sq * w_sq {
    S::zero()
  } else if tp > S::zero() {
    S::one()
  } else {
    -S::one()
  }
}

/// Computes the determinant of a matrix formed by the three points as columns
#[inline]
pub fn determinant_3d <S : Ring> (a : Point3 <S>, b : Point3 <S>, c : Point3 <S>) -> S {
  Matrix3::from_col_arrays ([a.0.into_array(), b.0.into_array(), c.0.into_array()])
    .determinant()
}

/// Create an orthonormal basis given the first basis component
pub fn orthonormal_basis <S> (e1 : Unit3 <S>) -> [Unit3 <S>; 3] where S : Real {
  use num::Zero;
  let u = Unit3::axis_x();
  let e2 = {
    let mut e2 = (*e1).cross (*u);
    e2.z += if e2.is_zero() {
      S::one()
    } else {
      S::zero()
    };
    Unit3::new (e2).unwrap()
  };
  let e3 = Unit3::new ((*e1).cross (*e2)).unwrap();
  [e1, e2, e3]
}

/// Coordinate-wise min
#[inline]
pub fn point2_min <S : PartialOrd> (a : Point2 <S>, b : Point2 <S>) -> Point2 <S> {
  Vector2::partial_min (a.0, b.0).into()
}

/// Coordinate-wise max
#[inline]
pub fn point2_max <S : PartialOrd> (a : Point2 <S>, b : Point2 <S>) -> Point2 <S> {
  Vector2::partial_max (a.0, b.0).into()
}

/// Coordinate-wise min
#[inline]
pub fn point3_min <S : PartialOrd> (a : Point3 <S>, b : Point3 <S>) -> Point3 <S> {
  Vector3::partial_min (a.0, b.0).into()
}

/// Coordinate-wise max
#[inline]
pub fn point3_max <S : PartialOrd> (a : Point3 <S>, b : Point3 <S>) -> Point3 <S> {
  Vector3::partial_max (a.0, b.0).into()
}

/// Given a 2D point and a 2D line, returns the nearest point on the line to the given
/// point
#[inline]
pub fn project_2d_point_on_line <S : OrderedRing> (point : Point2 <S>, line : Line2 <S>)
  -> Line2Point <S>
{
  let dot_dir = line.direction.dot (point - line.base);
  (dot_dir, line.point (dot_dir))
}

/// Given a 3D point and a 3D line, returns the nearest point on the line to the given
/// point
#[inline]
pub fn project_3d_point_on_line <S : OrderedRing> (point : Point3 <S>, line : Line3 <S>)
  -> Line3Point <S>
{
  let dot_dir = line.direction.dot (point - line.base);
  (dot_dir, line.point (dot_dir))
}

/// Given a 3D point and a 3D plane, returns the nearest point on the plane to the given
/// point
#[inline]
pub fn project_3d_point_on_plane <S : Real> (point : Point3 <S>, plane : Plane3 <S>)
  -> Point3 <S>
{
  let v = point - plane.base;
  let distance = (*plane.normal).dot (v);
  point - *plane.normal * distance
}

/// Returns point furthest along given direction together with dot product with the
/// direction vector
#[inline]
pub fn support2 <S> (points : &[Point2 <S>], direction : NonZero2 <S>) -> (Point2 <S>, S)
  where S : Ring + PartialOrd
{
  debug_assert!(!points.is_empty());
  points.iter()
    .map (|point| (*point, direction.dot (point.0)))
    .max_by (|(_, dot_a), (_, dot_b)| dot_a.partial_cmp (dot_b).unwrap()).unwrap()
}

/// Returns point furthest along given direction together with dot product with the
/// direction vector
#[inline]
pub fn support3 <S> (points : &[Point3 <S>], direction : NonZero3 <S>) -> (Point3 <S>, S)
  where S : Ring + PartialOrd
{
  debug_assert!(!points.is_empty());
  points.iter()
    .map (|point| (*point, direction.dot (point.0)))
    .max_by (|(_, dot_a), (_, dot_b)| dot_a.partial_cmp (dot_b).unwrap()).unwrap()
}

/// Square area of three points in 3D space.
///
/// Uses a numerically stable Heron's formula:
/// <https://en.wikipedia.org/wiki/Heron%27s_formula#Numerical_stability>
///
/// ```
/// # use math_utils::approx::assert_relative_eq;
/// # use math_utils::geometry::triangle_3d_area_squared;
/// assert_relative_eq!(
///   3.0/4.0,
///   *triangle_3d_area_squared (
///     [-1.0,  0.0,  0.0].into(),
///     [ 0.0,  0.0,  1.0].into(),
///     [ 0.0,  1.0,  0.0].into())
/// );
/// ```
///
/// If the area squared is zero then the points are colinear:
///
/// ```
/// # use math_utils::geometry::triangle_3d_area_squared;
/// assert_eq!(
///   0.0,
///   *triangle_3d_area_squared (
///     [-1.0, -1.0, -1.0].into(),
///     [ 0.0,  0.0,  0.0].into(),
///     [ 1.0,  1.0,  1.0].into())
/// );
/// ```
pub fn triangle_3d_area_squared <S> (a : Point3 <S>, b : Point3 <S>, c : Point3 <S>)
  -> NonNegative <S>
where
  S : OrderedField + Sqrt + approx::AbsDiffEq <Epsilon = S>
{
  use num::Zero;
  if a == b || a == c || b == c {
    return NonNegative::zero()
  }
  // compute the length of each side
  let ab_mag = *(b-a).norm();
  let ac_mag = *(c-a).norm();
  let bc_mag = *(c-b).norm();
  // order as min <= mid <= max
  let min_ab_ac = S::min (ab_mag, ac_mag);
  let max_ab_ac = S::max (ab_mag, ac_mag);
  let (min, mid, max) = if bc_mag <= min_ab_ac {
    (bc_mag, min_ab_ac, max_ab_ac)
  } else if bc_mag >= max_ab_ac {
    (min_ab_ac, max_ab_ac, bc_mag)
  } else {
    (min_ab_ac, bc_mag, max_ab_ac)
  };
  // 1.0/16.0
  let frac = S::two().powi (-4);
  let mut area_squared = frac
    * (max + (mid + min))
    * (min - (max - mid))
    * (min + (max - mid))
    * (max + (mid - min));
  // numerical accuracy may result in a negative value; from tests this is always a
  // relatively small value, in pathological cases (e.g. colinear points at several
  // thousand units) can result in a value on the order of -0.00001
  if area_squared < S::zero() {
    if cfg!(debug_assertions) {
      let lengths_squared = (b - a).magnitude_squared() + (c - b).magnitude_squared()
        + (a - c).magnitude_squared();
      approx::assert_abs_diff_eq!(area_squared, S::zero(),
        epsilon = S::default_epsilon() * S::two().powi (44) * lengths_squared)
    }
    area_squared = S::zero();
  }
  NonNegative::unchecked (area_squared)
}

// private

/// Smallest rectangle for given edge indices
fn smallest_rect <S : Real> (i0 : u32, i1 : u32, hull : &Hull2 <S>) -> Rect <S> {
  // NOTE: the following algorithm assumes that the points in the hull are
  // counter-clockwise sorted and that the hull does not contain any colinear points.
  // The construction of the hull should ensure these conditions.
  let points = hull.points();
  let edge             = points[i1 as usize] - points[i0 as usize];
  let perp             = vector2 (-edge.y, edge.x);
  let edge_len_squared = edge.magnitude_squared();
  let origin           = points[i1 as usize];
  let mut support      = [Point2::<S>::origin(); 4];
  let mut rect_index   = [i1, i1, i1, i1];
  #[expect(clippy::cast_possible_truncation)]
  for (i, point) in points.iter().enumerate() {
    let diff = point - origin;
    let v    = point2 (edge.dot (diff), perp.dot (diff));
    if v.x() > support[1].x() || v.x() == support[1].x() && v.y() > support[1].y() {
      rect_index[1] = i as u32;
      support[1]    = v;
    }
    if v.y() > support[2].y() || v.y() == support[2].y() && v.x() < support[2].x() {
      rect_index[2] = i as u32;
      support[2]    = v;
    }
    if v.x() < support[3].x() || v.x() == support[3].x() && v.y() < support[3].y() {
      rect_index[3] = i as u32;
      support[3]    = v;
    }
  }
  let scaled_width  = support[1].x() - support[3].x();
  let scaled_height = support[2].y();
  let area = scaled_width * scaled_height / edge_len_squared;

  Rect {
    edge, perp, area, edge_len_squared,
    indices: [rect_index[0], rect_index[1], rect_index[2], rect_index[3]]
  }
}

////////////////////////////////////////////////////////////////////////////////
//  impls                                                                     //
////////////////////////////////////////////////////////////////////////////////

impl <S> Interval <S> where S : Copy {
  #[inline]
  pub const fn min (self) -> S {
    self.min
  }
  #[inline]
  pub const fn max (self) -> S {
    self.max
  }
}
impl <S> Interval <S> where
  S : MinMax + Copy + PartialOrd + std::fmt::Debug
{
  #[inline]
  pub fn from_points (a : S, b : S) -> Self {
    let min = S::min (a, b);
    let max = S::max (a, b);
    Interval { min, max }
  }
  /// Returns `None` if points are not min/max
  #[inline]
  pub fn with_minmax (min : S, max : S) -> Option <Self> {
    if min > max {
      None
    } else {
      Some (Interval { min, max })
    }
  }
  /// Construct a new AABB with given min and max points.
  ///
  /// Debug panic if points are not min/max:
  ///
  /// ```should_panic
  /// # use math_utils::geometry::*;
  /// let b = Interval::with_minmax_unchecked (1.0, 0.0);  // panic!
  /// ```
  #[inline]
  pub fn with_minmax_unchecked (min : S, max : S) -> Self {
    debug_assert_eq!(min, S::min (min, max));
    debug_assert_eq!(max, S::max (min, max));
    Interval { min, max }
  }
  /// Construct the minimum AABB containing the given set of points.
  ///
  /// Returns `None` fewer than 2 points are given or all points are identical:
  ///
  /// ```
  /// # use math_utils::geometry::*;
  /// assert!(Interval::<f32>::containing (&[]).is_none());
  /// ```
  #[inline]
  pub fn containing (points : &[S]) -> Option <Self> {
    if points.len() < 2 {
      None
    } else {
      debug_assert!(points.len() >= 2);
      let mut min = S::min (points[0], points[1]);
      let mut max = S::max (points[0], points[1]);
      for point in points.iter().skip (2) {
        min = S::min (min, *point);
        max = S::max (max, *point);
      }
      Interval::with_minmax (min, max)
    }
  }
  #[inline]
  pub fn numcast <T> (self) -> Option <Interval <T>> where
    S : num::ToPrimitive,
    T : num::NumCast
  {
    Some (Interval {
      min: T::from (self.min)?,
      max: T::from (self.max)?
    })
  }
  /// Create a new AABB that is the union of the two input AABBs
  #[inline]
  pub fn union (a : Interval <S>, b : Interval <S>) -> Self {
    Interval::with_minmax_unchecked (
      S::min (a.min(), b.min()), S::max (a.max(), b.max()))
  }
  #[inline]
  pub fn length (self) -> NonNegative <S> where
    S : num::Zero + std::ops::Sub <Output=S>
  {
    self.width()
  }
  #[inline]
  pub fn width (self) -> NonNegative <S> where
    S : num::Zero + std::ops::Sub <Output=S>
  {
    NonNegative::unchecked (self.max - self.min)
  }
  #[inline]
  pub fn contains (self, point : S) -> bool {
    self.min < point && point < self.max
  }
  /// Clamp a given point to the AABB interval.
  ///
  /// ```
  /// # use math_utils::geometry::*;
  /// let b = Interval::with_minmax_unchecked (-1.0, 1.0);
  /// assert_eq!(b.clamp (-2.0), -1.0);
  /// assert_eq!(b.clamp ( 2.0),  1.0);
  /// assert_eq!(b.clamp ( 0.0),  0.0);
  /// ```
  #[inline]
  pub fn clamp (self, point : S) -> S {
    point.clamp (self.min, self.max)
  }
  /// Generate a random point contained in the AABB
  ///
  /// ```
  /// # use rand;
  /// # use rand_xorshift;
  /// # use math_utils::geometry::*;
  /// # fn main () {
  /// use rand_xorshift;
  /// use rand::SeedableRng;
  /// // random sequence will be the same each time this is run
  /// let mut rng = rand_xorshift::XorShiftRng::seed_from_u64 (0);
  /// let aabb = Interval::<f32>::with_minmax_unchecked (-10.0, 10.0);
  /// let point = aabb.rand_point (&mut rng);
  /// assert!(aabb.contains (point));
  /// let point = aabb.rand_point (&mut rng);
  /// assert!(aabb.contains (point));
  /// let point = aabb.rand_point (&mut rng);
  /// assert!(aabb.contains (point));
  /// let point = aabb.rand_point (&mut rng);
  /// assert!(aabb.contains (point));
  /// let point = aabb.rand_point (&mut rng);
  /// assert!(aabb.contains (point));
  /// # }
  /// ```
  #[inline]
  pub fn rand_point <R> (self, rng : &mut R) -> S where
    S : rand::distr::uniform::SampleUniform,
    R : rand::RngExt
  {
    rng.random_range (self.min..self.max)
  }
  #[inline]
  pub fn intersects (self, other : Interval <S>) -> bool {
    intersect::discrete_interval (self, other)
  }
  #[inline]
  pub fn intersection (self, other : Interval <S>) -> Option <Interval <S>> {
    intersect::continuous_interval (self, other)
  }

  /// Increase or decrease each endpoint by the given amount.
  ///
  /// Returns `None` if the interval width is less than twice the given amount:
  ///
  /// ```
  /// # use math_utils::geometry::*;
  /// let x = Interval::<f32>::with_minmax_unchecked (-1.0, 1.0);
  /// assert!(x.dilate (-2.0).is_none());
  /// ```
  pub fn dilate (mut self, amount : S) -> Option <Self> where
    S : std::ops::AddAssign + std::ops::SubAssign
  {
    self.min -= amount;
    self.max += amount;
    if self.min > self.max {
      None
    } else {
      Some (self)
    }
  }
}

impl <S> Primitive <S, S> for Interval <S> where
  S : OrderedField + AffineSpace <S, Translation=S>
    + VectorSpace <S, NonZero=NonZero <S>>
{
  fn translate (&mut self, displacement : S) {
    self.min += displacement;
    self.max += displacement;
  }
  fn scale (&mut self, scale : Positive <S>) {
    let old_width = *self.width();
    let new_width = old_width * *scale;
    let half_difference = (new_width - old_width) / S::two();
    self.min -= half_difference;
    self.max += half_difference;
  }
  fn support (&self, direction : NonZero <S>) -> (S, S) {
    if *direction > S::zero() {
      (self.max, self.max * *direction)
    } else {
      (self.min, self.min * *direction)
    }
  }
}

impl <S> std::ops::Add <S> for Interval <S> where
  S    : Field + AffineSpace <S, Translation=S> + VectorSpace <S>,
  Self : Primitive <S, S>
{
  type Output = Self;
  #[expect(clippy::renamed_function_params)]
  fn add (mut self, displacement : S) -> Self {
    self.translate (displacement);
    self
  }
}

impl <S> std::ops::Sub <S> for Interval <S> where
  S    : Field + AffineSpace <S, Translation=S> + VectorSpace <S>,
  Self : Primitive <S, S>
{
  type Output = Self;
  #[expect(clippy::renamed_function_params)]
  fn sub (mut self, displacement : S) -> Self {
    self.translate (-displacement);
    self
  }
}

impl_numcast_fields!(Aabb2, min, max);
impl <S> Aabb2 <S> where S : Copy {
  #[inline]
  pub const fn min (self) -> Point2 <S> {
    self.min
  }
  #[inline]
  pub const fn max (self) -> Point2 <S> {
    self.max
  }
}
impl <S> Aabb2 <S> where
  S : Copy + PartialOrd + std::fmt::Debug
{
  /// Returns `None` if points are not min/max or points are identical
  #[inline]
  pub fn with_minmax (min : Point2 <S>, max : Point2 <S>) -> Option <Self> {
    if min == max || min != point2_min (min, max) || max != point2_max (min, max) {
      None
    } else {
      Some (Aabb2 { min, max })
    }
  }
  /// Returns `None` if points are identical
  #[inline]
  pub fn from_points (a : Point2 <S>, b : Point2 <S>) -> Option <Self> {
    if a == b {
      None
    } else {
      let min = point2_min (a, b);
      let max = point2_max (a, b);
      Some (Aabb2 { min, max })
    }
  }
  /// Construct a new AABB with given min and max points.
  ///
  /// Debug panic if points are not min/max:
  ///
  /// ```should_panic
  /// # use math_utils::geometry::*;
  /// let b = Aabb2::with_minmax_unchecked (
  ///   [1.0, 1.0].into(),
  ///   [0.0, 0.0].into());  // panic!
  /// ```
  ///
  /// Debug panic if points are identical:
  ///
  /// ```should_panic
  /// # use math_utils::geometry::*;
  /// let b = Aabb2::with_minmax_unchecked (
  ///   [0.0, 0.0].into(),
  ///   [0.0, 0.0].into());  // panic!
  /// ```
  #[inline]
  pub fn with_minmax_unchecked (min : Point2 <S>, max : Point2 <S>) -> Self {
    debug_assert_ne!(min, max);
    debug_assert_eq!(min, point2_min (min, max));
    debug_assert_eq!(max, point2_max (min, max));
    Aabb2 { min, max }
  }
  /// Construct a new AABB using the two given points to determine min/max.
  ///
  /// Debug panic if points are identical:
  ///
  /// ```should_panic
  /// # use math_utils::geometry::*;
  /// let b = Aabb2::from_points_unchecked (
  ///   [0.0, 0.0].into(),
  ///   [0.0, 0.0].into());  // panic!
  /// ```
  #[inline]
  pub fn from_points_unchecked (a : Point2 <S>, b : Point2 <S>) -> Self {
    debug_assert_ne!(a, b);
    let min = point2_min (a, b);
    let max = point2_max (a, b);
    Aabb2 { min, max }
  }
  /// Construct the minimum AABB containing the given set of points.
  ///
  /// Returns `None` if fewer than 2 points are given or all points are identical:
  ///
  /// ```
  /// # use math_utils::geometry::*;
  /// assert!(Aabb2::<f32>::containing (&[]).is_none());
  /// ```
  ///
  /// ```
  /// # use math_utils::geometry::*;
  /// assert!(Aabb2::containing (&[[0.0, 0.0].into()]).is_none());
  /// ```
  #[inline]
  pub fn containing (points : &[Point2 <S>]) -> Option <Self> {
    if points.len() < 2 {
      None
    } else {
      debug_assert!(points.len() >= 2);
      let mut min = point2_min (points[0], points[1]);
      let mut max = point2_max (points[0], points[1]);
      for point in points.iter().skip (2) {
        min = point2_min (min, *point);
        max = point2_max (max, *point);
      }
      Aabb2::with_minmax (min, max)
    }
  }
  /// Create a new AABB that is the union of the two input AABBs
  #[inline]
  pub fn union (a : Aabb2 <S>, b : Aabb2 <S>) -> Self {
    Aabb2::with_minmax_unchecked (
      point2_min (a.min(), b.min()), point2_max (a.max(), b.max()))
  }
  #[inline]
  pub fn center (self) -> Point2 <S> where S : Ring + std::ops::Div <Output=S> {
    Point2 ((self.min.0 + self.max.0) / S::two())
  }
  #[inline]
  pub fn dimensions (self) -> Vector2 <NonNegative <S>> where
    S : num::Zero + std::ops::Sub <Output=S>
  {
    (self.max.0 - self.min.0).map (NonNegative::unchecked)
  }
  #[inline]
  pub fn extents (self) -> Vector2 <NonNegative <S>> where S : Field {
    self.dimensions().map (|d| d * NonNegative::unchecked (S::half()))
  }
  /// X dimension
  #[inline]
  pub fn width (self) -> NonNegative <S> where S : num::Zero + std::ops::Sub <Output=S> {
    NonNegative::unchecked (self.max.x() - self.min.x())
  }
  /// Y dimension
  #[inline]
  pub fn height (self) -> NonNegative <S> where
    S : num::Zero + std::ops::Sub <Output=S>
  {
    NonNegative::unchecked (self.max.y() - self.min.y())
  }
  #[inline]
  pub fn area (self) -> NonNegative <S> where
    S : num::Zero + std::ops::Sub <Output=S> + std::ops::Mul <Output=S>
  {
    self.width() * self.height()
  }
  #[inline]
  pub fn contains (self, point : Point2 <S>) -> bool {
    self.min.x() < point.x() && point.x() < self.max.x() &&
    self.min.y() < point.y() && point.y() < self.max.y()
  }
  /// Clamp a given point to the AABB.
  ///
  /// ```
  /// # use math_utils::geometry::*;
  /// let b = Aabb2::with_minmax ([-1.0, -1.0].into(), [1.0, 1.0].into()).unwrap();
  /// assert_eq!(b.clamp ([-2.0, 0.0].into()), [-1.0, 0.0].into());
  /// assert_eq!(b.clamp ([ 2.0, 2.0].into()), [ 1.0, 1.0].into());
  /// assert_eq!(b.clamp ([-0.5, 0.5].into()), [-0.5, 0.5].into());
  /// ```
  pub fn clamp (self, point : Point2 <S>) -> Point2 <S> where S : MinMax {
    [ point.x().clamp (self.min.x(), self.max.x()),
      point.y().clamp (self.min.y(), self.max.y())
    ].into()
  }
  /// Generate a random point contained in the AABB
  ///
  /// ```
  /// # use rand;
  /// # use rand_xorshift;
  /// # use math_utils::geometry::*;
  /// # fn main () {
  /// use rand_xorshift;
  /// use rand::SeedableRng;
  /// // random sequence will be the same each time this is run
  /// let mut rng = rand_xorshift::XorShiftRng::seed_from_u64 (0);
  /// let aabb = Aabb2::<f32>::with_minmax (
  ///   [-10.0, -10.0].into(),
  ///   [ 10.0,  10.0].into()
  /// ).unwrap();
  /// let point = aabb.rand_point (&mut rng);
  /// assert!(aabb.contains (point));
  /// let point = aabb.rand_point (&mut rng);
  /// assert!(aabb.contains (point));
  /// let point = aabb.rand_point (&mut rng);
  /// assert!(aabb.contains (point));
  /// let point = aabb.rand_point (&mut rng);
  /// assert!(aabb.contains (point));
  /// let point = aabb.rand_point (&mut rng);
  /// assert!(aabb.contains (point));
  /// # }
  /// ```
  #[inline]
  pub fn rand_point <R> (self, rng : &mut R) -> Point2 <S> where
    S : rand::distr::uniform::SampleUniform,
    R : rand::RngExt
  {
    let x_range = self.min.x()..self.max.x();
    let y_range = self.min.y()..self.max.y();
    let x = if x_range.is_empty() {
      self.min.x()
    } else {
      rng.random_range (x_range)
    };
    let y = if y_range.is_empty() {
      self.min.y()
    } else {
      rng.random_range (y_range)
    };
    [x, y].into()
  }
  #[inline]
  pub fn intersects (self, other : Aabb2 <S>) -> bool {
    intersect::discrete_aabb2_aabb2 (self, other)
  }
  #[inline]
  pub fn intersection (self, other : Aabb2 <S>) -> Option <Aabb2 <S>> {
    intersect::continuous_aabb2_aabb2 (self, other)
  }

  pub fn corner (self, quadrant : Quadrant) -> Point2 <S> {
    match quadrant {
      // upper-right
      Quadrant::PosPos => self.max,
      // lower-right
      Quadrant::PosNeg => [self.max.x(), self.min.y()].into(),
      // upper-left
      Quadrant::NegPos => [self.min.x(), self.max.y()].into(),
      // lower-left
      Quadrant::NegNeg => self.min
    }
  }

  /// Corner points clockwise from min
  pub fn corners (self) -> [Point2 <S>; 4] {
    [ self.min, [self.min.x(), self.max.y()].into(),
      self.max, [self.max.x(), self.min.y()].into()
    ]
  }

  pub fn edge (self, direction : SignedAxis2) -> Self {
    let (min, max) = match direction {
      SignedAxis2::PosX => (
        self.corner (Quadrant::PosNeg),
        self.corner (Quadrant::PosPos) ),
      SignedAxis2::NegX => (
        self.corner (Quadrant::NegNeg),
        self.corner (Quadrant::NegPos) ),
      SignedAxis2::PosY => (
        self.corner (Quadrant::NegPos),
        self.corner (Quadrant::PosPos) ),
      SignedAxis2::NegY => (
        self.corner (Quadrant::NegNeg),
        self.corner (Quadrant::PosNeg) )
    };
    Aabb2::with_minmax_unchecked (min, max)
  }

  /// Extrude (or inset) a new Aabb from a face:
  ///
  /// ```
  /// # use math_utils::*;
  /// # use math_utils::geometry::*;
  /// let x = Aabb2::with_minmax ([-1.0, -1.0].into(), [1.0, 1.0].into())
  ///   .unwrap();
  /// assert_eq!(x.extrude (SignedAxis2::PosY, 0.5),
  ///   Aabb2::with_minmax ([-1.0, 1.0].into(), [1.0, 1.5].into()));
  /// assert_eq!(x.extrude (SignedAxis2::PosY, -0.5),
  ///   Aabb2::with_minmax ([-1.0, 0.5].into(), [1.0, 1.0].into()));
  /// ```
  #[must_use]
  pub fn extrude (self, axis : SignedAxis2, amount : S) -> Option <Self> where S : Ring {
    let (p1, p2) = match axis {
      SignedAxis2::PosX => (
        self.min + Vector2::zero().with_x (*self.width()),
        self.max + Vector2::zero().with_x (amount) ),
      SignedAxis2::NegX => (
        self.min - Vector2::zero().with_x (amount),
        self.max - Vector2::zero().with_x (*self.height()) ),
      SignedAxis2::PosY => (
        self.min + Vector2::zero().with_y (*self.height()),
        self.max + Vector2::zero().with_y (amount) ),
      SignedAxis2::NegY => (
        self.min - Vector2::zero().with_y (amount),
        self.max - Vector2::zero().with_y (*self.height()) )
    };
    Self::from_points (p1, p2)
  }

  /// Extend (or contract) one of the extents
  #[must_use]
  pub fn extend (mut self, axis : SignedAxis2, amount : S) -> Option <Self>
    where S : std::ops::AddAssign + std::ops::SubAssign
  {
    match axis {
      SignedAxis2::PosX => self.max.0.x += amount,
      SignedAxis2::NegX => self.min.0.x -= amount,
      SignedAxis2::PosY => self.max.0.y += amount,
      SignedAxis2::NegY => self.min.0.y -= amount
    }
    Self::with_minmax (self.min, self.max)
  }

  pub fn with_z (self, z : Interval <S>) -> Aabb3 <S> {
    Aabb3::with_minmax_unchecked (
      self.min.0.with_z (z.min()).into(),
      self.max.0.with_z (z.max()).into()
    )
  }

  /// Increase or decrease each extent by the given amount.
  ///
  /// Returns `None` if any extent would become negative or if all dimensions become
  /// zero:
  /// ```
  /// # use math_utils::geometry::*;
  /// let x = Aabb2::with_minmax ([0.0, 0.0].into(), [1.0, 2.0].into()).unwrap();
  /// assert!(x.dilate (-1.0).is_none());
  /// ```
  #[must_use]
  pub fn dilate (mut self, amount : S) -> Option <Self> where S : Ring {
    self.min -= Vector2::broadcast (amount);
    self.max += Vector2::broadcast (amount);
    if self.min == self.max || self.min != point2_min (self.min, self.max)
      || self.max != point2_max (self.min, self.max)
    {
      None
    } else {
      Some (self)
    }
  }

  /// Project *along* the given axis.
  ///
  /// For example, projecting *along* the X-axis projects *onto* the Y-axis.
  pub fn project (self, axis : Axis2) -> Interval <S> where S : MinMax {
    let (min, max) = match axis {
      Axis2::X => (self.min.y(), self.max.y()),
      Axis2::Y => (self.min.x(), self.max.x())
    };
    Interval::with_minmax_unchecked (min, max)
  }
}

impl <S> Primitive <S, Point2 <S>> for Aabb2 <S> where S : OrderedField {
  fn translate (&mut self, displacement : Vector2 <S>) {
    self.min += displacement;
    self.max += displacement;
  }
  fn scale (&mut self, scale : Positive <S>) {
    let center = self.center();
    self.translate (-center.0);
    self.min.0 *= *scale;
    self.max.0 *= *scale;
    self.translate (center.0)
  }
  fn support (&self, direction : NonZero2 <S>) -> (Point2 <S>, S) {
    let out = if let Some (quadrant) = Quadrant::from_vec_strict (*direction) {
      self.corner (quadrant)
    } else if direction.x > S::zero() {
      debug_assert_eq!(direction.y, S::zero());
      point2 (self.max.x(), (self.min.y() + self.max.y()) / S::two())
    } else if direction.x < S::zero() {
      debug_assert_eq!(direction.y, S::zero());
      point2 (self.min.x(), (self.min.y() + self.max.y()) / S::two())
    } else {
      debug_assert_eq!(direction.x, S::zero());
      if direction.y > S::zero() {
        point2 ((self.min.x() + self.max.x()) / S::two(), self.max.y())
      } else {
        debug_assert!(direction.y < S::zero());
        point2 ((self.min.x() + self.max.x()) / S::two(), self.min.y())
      }
    };
    (out, out.0.dot (*direction))
  }
}

impl <S> std::ops::Add <Vector2 <S>> for Aabb2 <S> where
  S    : Field,
  Self : Primitive <S, Point2 <S>>
{
  type Output = Self;
  #[expect(clippy::renamed_function_params)]
  fn add (mut self, displacement : Vector2 <S>) -> Self {
    self.translate (displacement);
    self
  }
}

impl <S> std::ops::Sub <Vector2 <S>> for Aabb2 <S> where
  S    : Field,
  Self : Primitive <S, Point2 <S>>
{
  type Output = Self;
  #[expect(clippy::renamed_function_params)]
  fn sub (mut self, displacement : Vector2 <S>) -> Self {
    self.translate (-displacement);
    self
  }
}

impl_numcast_fields!(Aabb3, min, max);
impl <S> Aabb3 <S> where S : Copy {
  #[inline]
  pub const fn min (self) -> Point3 <S> {
    self.min
  }
  #[inline]
  pub const fn max (self) -> Point3 <S> {
    self.max
  }
}
impl <S> Aabb3 <S> where S : Copy + PartialOrd + std::fmt::Debug {
  /// Returns `None` if points are not min/max or points are identical
  #[inline]
  pub fn with_minmax (min : Point3 <S>, max : Point3 <S>) -> Option <Self> {
    if min == max || min != point3_min (min, max) || max != point3_max (min, max) {
      None
    } else {
      Some (Aabb3 { min, max })
    }
  }
  /// Returns `None` if points are identical
  #[inline]
  pub fn from_points (a : Point3 <S>, b : Point3 <S>) -> Option <Self> {
    if a == b {
      None
    } else {
      let min = point3_min (a, b);
      let max = point3_max (a, b);
      Some (Aabb3 { min, max })
    }
  }
  /// Construct a new AABB with given min and max points.
  ///
  /// Debug panic if points are not min/max:
  ///
  /// ```should_panic
  /// # use math_utils::geometry::*;
  /// let b = Aabb3::with_minmax_unchecked (
  ///   [1.0, 1.0, 1.0].into(),
  ///   [0.0, 0.0, 0.0].into());  // panic!
  /// ```
  ///
  /// Debug panic if points are identical:
  ///
  /// ```should_panic
  /// # use math_utils::geometry::*;
  /// let b = Aabb3::with_minmax_unchecked (
  ///   [0.0, 0.0, 0.0].into(),
  ///   [0.0, 0.0, 0.0].into());  // panic!
  /// ```
  #[inline]
  pub fn with_minmax_unchecked (min : Point3 <S>, max : Point3 <S>) -> Self {
    debug_assert_ne!(min, max);
    debug_assert_eq!(min, point3_min (min, max));
    debug_assert_eq!(max, point3_max (min, max));
    Aabb3 { min, max }
  }
  /// Construct a new AABB using the two given points to determine min/max.
  ///
  /// Panic if points are identical:
  ///
  /// ```should_panic
  /// # use math_utils::geometry::*;
  /// let b = Aabb3::from_points_unchecked (
  ///   [0.0, 0.0, 0.0].into(),
  ///   [0.0, 0.0, 0.0].into());  // panic!
  /// ```
  #[inline]
  pub fn from_points_unchecked (a : Point3 <S>, b : Point3 <S>) -> Self {
    debug_assert_ne!(a, b);
    let min = point3_min (a, b);
    let max = point3_max (a, b);
    Aabb3 { min, max }
  }
  /// Construct the minimum AABB containing the given set of points.
  ///
  /// Returns `None` if fewer than 2 points are given or if all points are identical:
  ///
  /// ```should_panic
  /// # use math_utils::geometry::*;
  /// assert!(Aabb3::<f32>::containing (&[]).is_none());
  /// ```
  ///
  /// ```should_panic
  /// # use math_utils::geometry::*;
  /// assert!(Aabb3::containing (&[[0.0, 0.0, 0.0].into()]).is_none());
  /// ```
  #[inline]
  pub fn containing (points : &[Point3 <S>]) -> Option <Self> {
    debug_assert!(points.len() >= 2);
    let mut min = point3_min (points[0], points[1]);
    let mut max = point3_max (points[0], points[1]);
    for point in points.iter().skip (2) {
      min = point3_min (min, *point);
      max = point3_max (max, *point);
    }
    Aabb3::with_minmax (min, max)
  }
  /// Create a new AABB that is the union of the two input AABBs
  #[inline]
  pub fn union (a : Aabb3 <S>, b : Aabb3 <S>) -> Self {
    Aabb3::with_minmax_unchecked (
      point3_min (a.min(), b.min()), point3_max (a.max(), b.max()))
  }
  #[inline]
  pub fn center (self) -> Point3 <S> where S : Ring + std::ops::Div <Output=S> {
    Point3 ((self.min.0 + self.max.0) / S::two())
  }
  #[inline]
  pub fn dimensions (self) -> Vector3 <NonNegative <S>> where
    S : num::Zero + std::ops::Sub <Output=S>
  {
    (self.max.0 - self.min.0).map (NonNegative::unchecked)
  }
  #[inline]
  pub fn extents (self) -> Vector3 <NonNegative <S>> where S : Field {
    self.dimensions().map (|d| d * NonNegative::unchecked (S::half()))
  }
  /// X dimension
  #[inline]
  pub fn width (self) -> NonNegative <S> where S : num::Zero + std::ops::Sub <Output=S> {
    NonNegative::unchecked (self.max.x() - self.min.x())
  }
  /// Y dimension
  #[inline]
  pub fn depth (self) -> NonNegative <S> where S : num::Zero + std::ops::Sub <Output=S> {
    NonNegative::unchecked (self.max.y() - self.min.y())
  }
  /// Z dimension
  #[inline]
  pub fn height (self) -> NonNegative <S> where
    S : num::Zero + std::ops::Sub <Output=S>
  {
    NonNegative::unchecked (self.max.z() - self.min.z())
  }
  #[inline]
  pub fn volume (self) -> NonNegative <S> where
    S : num::Zero + std::ops::Sub <Output=S> + std::ops::Mul <Output=S>
  {
    self.width() * self.depth() * self.height()
  }
  #[inline]
  pub fn contains (self, point : Point3 <S>) -> bool {
    self.min.x() < point.x() && point.x() < self.max.x() &&
    self.min.y() < point.y() && point.y() < self.max.y() &&
    self.min.z() < point.z() && point.z() < self.max.z()
  }
  /// Clamp a given point to the AABB.
  ///
  /// ```
  /// # use math_utils::geometry::*;
  /// let b = Aabb3::with_minmax ([-1.0, -1.0, -1.0].into(), [1.0, 1.0, 1.0].into())
  ///   .unwrap();
  /// assert_eq!(b.clamp ([-2.0, 0.0, 0.0].into()), [-1.0, 0.0, 0.0].into());
  /// assert_eq!(b.clamp ([ 2.0, 2.0, 0.0].into()), [ 1.0, 1.0, 0.0].into());
  /// assert_eq!(b.clamp ([-1.0, 2.0, 3.0].into()), [-1.0, 1.0, 1.0].into());
  /// assert_eq!(b.clamp ([-0.5, 0.5, 0.0].into()), [-0.5, 0.5, 0.0].into());
  /// ```
  pub fn clamp (self, point : Point3 <S>) -> Point3 <S> where S : MinMax {
    [ point.x().clamp (self.min.x(), self.max.x()),
      point.y().clamp (self.min.y(), self.max.y()),
      point.z().clamp (self.min.z(), self.max.z())
    ].into()
  }
  /// Generate a random point contained in the AABB
  ///
  /// ```
  /// # use rand;
  /// # use rand_xorshift;
  /// # use math_utils::geometry::*;
  /// # fn main () {
  /// use rand_xorshift;
  /// use rand::SeedableRng;
  /// // random sequence will be the same each time this is run
  /// let mut rng = rand_xorshift::XorShiftRng::seed_from_u64 (0);
  /// let aabb = Aabb3::<f32>::with_minmax (
  ///   [-10.0, -10.0, -10.0].into(),
  ///   [ 10.0,  10.0,  10.0].into()
  /// ).unwrap();
  /// let point = aabb.rand_point (&mut rng);
  /// assert!(aabb.contains (point));
  /// let point = aabb.rand_point (&mut rng);
  /// assert!(aabb.contains (point));
  /// let point = aabb.rand_point (&mut rng);
  /// assert!(aabb.contains (point));
  /// let point = aabb.rand_point (&mut rng);
  /// assert!(aabb.contains (point));
  /// let point = aabb.rand_point (&mut rng);
  /// assert!(aabb.contains (point));
  /// # }
  /// ```
  #[inline]
  pub fn rand_point <R> (self, rng : &mut R) -> Point3 <S> where
    S : rand::distr::uniform::SampleUniform,
    R : rand::RngExt
  {
    let x_range = self.min.x()..self.max.x();
    let y_range = self.min.y()..self.max.y();
    let z_range = self.min.z()..self.max.z();
    let x = if x_range.is_empty() {
      self.min.x()
    } else {
      rng.random_range (x_range)
    };
    let y = if y_range.is_empty() {
      self.min.y()
    } else {
      rng.random_range (y_range)
    };
    let z = if z_range.is_empty() {
      self.min.z()
    } else {
      rng.random_range (z_range)
    };
    [x, y, z].into()
  }
  #[inline]
  pub fn intersects (self, other : Aabb3 <S>) -> bool {
    intersect::discrete_aabb3_aabb3 (self, other)
  }
  #[inline]
  pub fn intersection (self, other : Aabb3 <S>) -> Option <Aabb3 <S>> {
    intersect::continuous_aabb3_aabb3 (self, other)
  }

  pub fn corner (self, octant : Octant) -> Point3 <S> {
    match octant {
      Octant::PosPosPos => self.max,
      Octant::NegPosPos => [self.min.x(), self.max.y(), self.max.z()].into(),
      Octant::PosNegPos => [self.max.x(), self.min.y(), self.max.z()].into(),
      Octant::NegNegPos => [self.min.x(), self.min.y(), self.max.z()].into(),
      Octant::PosPosNeg => [self.max.x(), self.max.y(), self.min.z()].into(),
      Octant::NegPosNeg => [self.min.x(), self.max.y(), self.min.z()].into(),
      Octant::PosNegNeg => [self.max.x(), self.min.y(), self.min.z()].into(),
      Octant::NegNegNeg => self.min
    }
  }

  pub fn corners (self) -> [Point3 <S>; 8] {
    [ self.min,
      [self.min.x(), self.max.y(), self.max.z()].into(),
      [self.max.x(), self.min.y(), self.max.z()].into(),
      [self.min.x(), self.min.y(), self.max.z()].into(),
      [self.max.x(), self.max.y(), self.min.z()].into(),
      [self.min.x(), self.max.y(), self.min.z()].into(),
      [self.max.x(), self.min.y(), self.min.z()].into(),
      self.max
    ]
  }

  pub fn face (self, direction : SignedAxis3) -> Self {
    let (min, max) = match direction {
      SignedAxis3::PosX => (
        self.corner (Octant::PosNegNeg),
        self.corner (Octant::PosPosPos) ),
      SignedAxis3::NegX => (
        self.corner (Octant::NegNegNeg),
        self.corner (Octant::NegPosPos) ),
      SignedAxis3::PosY => (
        self.corner (Octant::NegPosNeg),
        self.corner (Octant::PosPosPos) ),
      SignedAxis3::NegY => (
        self.corner (Octant::NegNegNeg),
        self.corner (Octant::PosNegPos) ),
      SignedAxis3::PosZ => (
        self.corner (Octant::NegNegPos),
        self.corner (Octant::PosPosPos) ),
      SignedAxis3::NegZ => (
        self.corner (Octant::NegNegNeg),
        self.corner (Octant::PosPosNeg) )
    };
    Aabb3::with_minmax_unchecked (min, max)
  }

  /// Extrude (or inset) a new Aabb from a face:
  ///
  /// ```
  /// # use math_utils::*;
  /// # use math_utils::geometry::*;
  /// let x = Aabb3::with_minmax ([-1.0, -1.0, -1.0].into(), [1.0, 1.0, 1.0].into())
  ///   .unwrap();
  /// assert_eq!(x.extrude (SignedAxis3::PosZ, 0.5),
  ///   Aabb3::with_minmax ([-1.0, -1.0, 1.0].into(), [1.0, 1.0, 1.5].into()));
  /// assert_eq!(x.extrude (SignedAxis3::PosZ, -0.5),
  ///   Aabb3::with_minmax ([-1.0, -1.0, 0.5].into(), [1.0, 1.0, 1.0].into()));
  /// ```
  #[must_use]
  pub fn extrude (self, axis : SignedAxis3, amount : S) -> Option <Self> where S : Ring {
    let (p1, p2) = match axis {
      SignedAxis3::PosX => (
        self.min + Vector3::zero().with_x (*self.width()),
        self.max + Vector3::zero().with_x (amount)
      ),
      SignedAxis3::NegX => (
        self.min - Vector3::zero().with_x (amount),
        self.max - Vector3::zero().with_x (*self.width())
      ),
      SignedAxis3::PosY => (
        self.min + Vector3::zero().with_y (*self.depth()),
        self.max + Vector3::zero().with_y (amount)
      ),
      SignedAxis3::NegY => (
        self.min - Vector3::zero().with_y (amount),
        self.max - Vector3::zero().with_y (*self.depth())
      ),
      SignedAxis3::PosZ => (
        self.min + Vector3::zero().with_z (*self.height()),
        self.max + Vector3::zero().with_z (amount)
      ),
      SignedAxis3::NegZ => (
        self.min - Vector3::zero().with_z (amount),
        self.max - Vector3::zero().with_z (*self.height())
      )
    };
    Aabb3::from_points (p1, p2)
  }

  /// Extend (or contract) one of the extents
  #[must_use]
  pub fn extend (mut self, axis : SignedAxis3, amount : S) -> Option <Self>
    where S : std::ops::AddAssign + std::ops::SubAssign
  {
    match axis {
      SignedAxis3::PosX => self.max.0.x += amount,
      SignedAxis3::NegX => self.min.0.x -= amount,
      SignedAxis3::PosY => self.max.0.y += amount,
      SignedAxis3::NegY => self.min.0.y -= amount,
      SignedAxis3::PosZ => self.max.0.z += amount,
      SignedAxis3::NegZ => self.min.0.z -= amount
    }
    Self::with_minmax (self.min, self.max)
  }

  /// Increase or decrease each extent by the given amount.
  ///
  /// Returns `None` if any extent would become negative:
  /// ```
  /// # use math_utils::geometry::*;
  /// let x = Aabb3::with_minmax ([0.0, 0.0, 0.0].into(), [1.0, 2.0, 3.0].into())
  ///   .unwrap();
  /// assert!(x.dilate (-1.0).is_none());
  /// ```
  #[must_use]
  pub fn dilate (mut self, amount : S) -> Option <Self> where S : Ring {
    self.min -= Vector3::broadcast (amount);
    self.max += Vector3::broadcast (amount);
    if self.min == self.max || self.min != point3_min (self.min, self.max)
      || self.max != point3_max (self.min, self.max)
    {
      None
    } else {
      Some (self)
    }
  }

  /// Project *along* the given axis.
  ///
  /// For example, projecting *along* the X-axis projects *onto* the YZ-plane.
  pub fn project (self, axis : Axis3) -> Aabb2 <S> {
    let (min, max) = match axis {
      Axis3::X => ([self.min.y(), self.min.z()], [self.max.y(), self.max.z()]),
      Axis3::Y => ([self.min.x(), self.min.z()], [self.max.x(), self.max.z()]),
      Axis3::Z => ([self.min.x(), self.min.y()], [self.max.x(), self.max.y()]),
    };
    Aabb2::with_minmax_unchecked (min.into(), max.into())
  }
}

impl <S> Primitive <S, Point3 <S>> for Aabb3 <S> where S : OrderedField {
  fn translate (&mut self, displacement : Vector3 <S>) {
    self.min += displacement;
    self.max += displacement;
  }
  fn scale (&mut self, scale : Positive <S>) {
    let center = self.center();
    self.translate (-center.0);
    self.min.0 *= *scale;
    self.max.0 *= *scale;
    self.translate (center.0);
  }
  fn support (&self, direction : NonZero3 <S>) -> (Point3 <S>, S) {
    let out = if let Some (octant) = Octant::from_vec_strict (*direction) {
      self.corner (octant)
    } else {
      use Sign::*;
      let center = self.center();
      match direction.map (S::sign).into_array() {
        // faces
        [Positive, Zero, Zero] => center.0.with_x (self.max.x()).into(),
        [Negative, Zero, Zero] => center.0.with_x (self.min.x()).into(),
        [Zero, Positive, Zero] => center.0.with_y (self.max.y()).into(),
        [Zero, Negative, Zero] => center.0.with_y (self.min.y()).into(),
        [Zero, Zero, Positive] => center.0.with_z (self.max.z()).into(),
        [Zero, Zero, Negative] => center.0.with_z (self.min.z()).into(),
        // edges
        [Positive, Positive, Zero] =>
          center.0.with_x (self.max.x()).with_y (self.max.y()).into(),
        [Positive, Negative, Zero] =>
          center.0.with_x (self.max.x()).with_y (self.min.y()).into(),
        [Negative, Positive, Zero] =>
          center.0.with_x (self.min.x()).with_y (self.max.y()).into(),
        [Negative, Negative, Zero] =>
          center.0.with_x (self.min.x()).with_y (self.min.y()).into(),

        [Positive, Zero, Positive] =>
          center.0.with_x (self.max.x()).with_z (self.max.z()).into(),
        [Positive, Zero, Negative] =>
          center.0.with_x (self.max.x()).with_z (self.min.z()).into(),
        [Negative, Zero, Positive] =>
          center.0.with_x (self.min.x()).with_z (self.max.z()).into(),
        [Negative, Zero, Negative] =>
          center.0.with_x (self.min.x()).with_z (self.min.z()).into(),

        [Zero, Positive, Positive] =>
          center.0.with_y (self.max.y()).with_z (self.max.z()).into(),
        [Zero, Positive, Negative] =>
          center.0.with_y (self.max.y()).with_z (self.min.z()).into(),
        [Zero, Negative, Positive] =>
          center.0.with_y (self.min.y()).with_z (self.max.z()).into(),
        [Zero, Negative, Negative] =>
          center.0.with_y (self.min.y()).with_z (self.min.z()).into(),
        [_, _, _] => unreachable!()
      }
    };
    (out, out.0.dot (*direction))
  }
}

impl <S> std::ops::Add <Vector3 <S>> for Aabb3 <S> where
  S    : Field,
  Self : Primitive <S, Point3 <S>>
{
  type Output = Self;
  #[expect(clippy::renamed_function_params)]
  fn add (mut self, displacement : Vector3 <S>) -> Self {
    self.translate (displacement);
    self
  }
}

impl <S> std::ops::Sub <Vector3 <S>> for Aabb3 <S> where
  S    : Field,
  Self : Primitive <S, Point3 <S>>
{
  type Output = Self;
  #[expect(clippy::renamed_function_params)]
  fn sub (mut self, displacement : Vector3 <S>) -> Self {
    self.translate (-displacement);
    self
  }
}

impl_numcast_fields!(Obb2, center, extents, orientation);
impl <S : OrderedRing> Obb2 <S> {
  /// Construct a new OBB
  #[inline]
  pub const fn new (
    center      : Point2 <S>,
    extents     : Vector2 <Positive <S>>,
    orientation : AngleWrapped <S>
  ) -> Self {
    Obb2 { center, extents, orientation }
  }
  /// Compute the OBB of a given orientation.
  ///
  /// Returns None if the points span an empty set.
  pub fn containing_points_with_orientation (
    points : &[Point2 <S>], orientation : AngleWrapped <S>
  ) -> Option <Self> where
    S : Real + num::real::Real + MaybeSerDes
  {
    let [x, y] = Rotation2::from_angle (orientation.angle()).cols
      .map (NonZero2::unchecked).into_array();
    let min_dot_x = -support2 (points, -x).1;
    let max_dot_x =  support2 (points,  x).1;
    let min_dot_y = -support2 (points, -y).1;
    let max_dot_y =  support2 (points,  y).1;
    let center = Point2 (
      (*x * min_dot_x + *x * max_dot_x) * S::half() +
      (*y * min_dot_y + *y * max_dot_y) * S::half());
    let extents = vector2 (
      Positive::new (S::half() * (max_dot_x - min_dot_x))?,
      Positive::new (S::half() * (max_dot_y - min_dot_y))?);
    Some (Obb2 { center, extents, orientation })
  }
  /// Construct the minimum OBB containing the convex hull of a given set of points
  /// using rotating calipers algorithm. To construct the OBB containing a convex hull
  /// directoy use the `Obb2::containing` method.
  ///
  /// Returns None if fewer than 3 points are given or if points span an empty set:
  ///
  /// ```
  /// # use math_utils::*;
  /// # use math_utils::geometry::*;
  /// assert!(Obb2::containing_points (
  ///   &[[0.0, 1.0], [1.0, 0.0]].map (Point2::from)
  /// ).is_none());
  ///
  /// assert!(Obb2::containing_points (
  ///   &[[1.0, 1.0], [0.0, 0.0], [-1.0, -1.0]].map (Point2::from)
  /// ).is_none());
  /// ```
  pub fn containing_points (points : &[Point2 <S>]) -> Option <Self> where S : Real {
    let hull = Hull2::from_points (points)?;
    Self::containing (&hull)
  }
  /// Construct the minimum OBB containing the given convex hull of points using
  /// rotating calipers algorithm.
  ///
  /// Returns None if fewer than 3 points are given or if points span an empty set:
  ///
  /// ```
  /// # use math_utils::*;
  /// # use math_utils::geometry::*;
  /// let hull = Hull2::from_points (
  ///   &[[0.0, 1.0], [1.0, 0.0]].map (Point2::from)).unwrap();
  /// assert!(Obb2::containing (&hull).is_none());
  ///
  /// let hull = Hull2::from_points (
  ///   &[[1.0, 1.0], [0.0, 0.0], [-1.0, -1.0]].map (Point2::from)).unwrap();
  /// assert!(Obb2::containing (&hull).is_none());
  /// ```
  // code follows GeometricTools:
  // <https://github.com/davideberly/GeometricTools/blob/f78dd0b65bc3a0a4723586de6991dd2c339b08ad/GTE/Mathematics/MinimumAreaBox2.h>
  pub fn containing (hull : &Hull2 <S>) -> Option <Self> where S : Real {
    let points = hull.points();
    if points.len() < 3 {
      return None
    }
    let mut visited  = vec![false; points.len()];
    #[expect(clippy::cast_possible_truncation)]
    let mut min_rect = smallest_rect (points.len() as u32 - 1, 0, hull);
    visited[min_rect.indices[0] as usize] = true;
    // rotating calipers
    let mut rect = min_rect.clone();
    for _ in 0..points.len() {
      let mut angles  = [(S::zero(), u32::MAX); 4];
      let mut nangles = u32::MAX;
      if !compute_angles (hull, &rect, &mut angles, &mut nangles) {
        break
      }
      let sorted = sort_angles (&angles, nangles);
      if !update_support (&angles, nangles, &sorted, hull, &mut visited, &mut rect) {
        break
      }
      if rect.area < min_rect.area {
        min_rect = rect.clone();
      }
    }
    // construct obb from min rect
    let sum = [
      points[min_rect.indices[1] as usize].0 + points[min_rect.indices[3] as usize].0,
      points[min_rect.indices[2] as usize].0 + points[min_rect.indices[0] as usize].0
    ];
    let diff = [
      points[min_rect.indices[1] as usize].0 - points[min_rect.indices[3] as usize].0,
      points[min_rect.indices[2] as usize].0 - points[min_rect.indices[0] as usize].0
    ];
    let obb = {
      let center = Point2::from ((min_rect.edge * min_rect.edge.dot (sum[0]) +
        min_rect.perp * min_rect.perp.dot (sum[1]))
        * S::half() / min_rect.edge_len_squared);
      let extents = {
        let extent_x = Positive::unchecked (S::sqrt (
          (S::half() * min_rect.edge.dot (diff[0])).squared()
          / min_rect.edge_len_squared));
        let extent_y = Positive::unchecked (S::sqrt (
          (S::half() * min_rect.perp.dot (diff[1])).squared()
          / min_rect.edge_len_squared));
        vector2 (extent_x, extent_y)
      };
      let orientation =
        AngleWrapped::wrap (Rad (min_rect.edge.y.atan2 (min_rect.edge.x)));
      Obb2 { center, extents, orientation }
    };
    fn compute_angles <S : Real> (
      hull    : &Hull2 <S>,
      rect    : &Rect <S>,
      angles  : &mut [(S, u32); 4],
      nangles : &mut u32
    ) -> bool {
      *nangles = 0;
      let points = hull.points();
      let mut k0 = 3;
      let mut k1 = 0;
      while k1 < 4 {
        if rect.indices[k0] != rect.indices[k1] {
          let d = {
            let mut d = if k0 & 1 != 0 {
              rect.perp
            } else {
              rect.edge
            };
            if k0 & 2 != 0 {
              d = -d;
            }
            d
          };
          let j0 = rect.indices[k0];
          let mut j1 = j0 + 1;
          #[expect(clippy::cast_possible_truncation)]
          if j1 == points.len() as u32 {
            j1 = 0;
          }
          let e = points[j1 as usize] - points[j0 as usize];
          let dp = d.dot ([e.y, -e.x].into());
          let e_lensq = e.magnitude_squared();
          let sin_theta_squared = (dp * dp) / e_lensq;
          angles[*nangles as usize] =
            #[expect(clippy::cast_possible_truncation)]
            (sin_theta_squared, k0 as u32);
          *nangles += 1;
        }
        k0 = k1;
        k1 += 1;
      }
      *nangles > 0
    }
    fn sort_angles <S : PartialOrd> (angles : &[(S, u32); 4], nangles : u32)
      -> [u32; 4]
    {
      let mut sorted = [0u32, 1, 2, 3];
      match nangles {
        0 | 1 => {}
        2 => if angles[sorted[0] as usize].0 > angles[sorted[1] as usize].0 {
          sorted.swap (0, 1)
        }
        3 => {
          if angles[sorted[0] as usize].0 > angles[sorted[1] as usize].0 {
            sorted.swap (0, 1)
          }
          if angles[sorted[0] as usize].0 > angles[sorted[2] as usize].0 {
            sorted.swap (0, 2)
          }
          if angles[sorted[1] as usize].0 > angles[sorted[2] as usize].0 {
            sorted.swap (1, 2)
          }
        }
        4 => {
          if angles[sorted[0] as usize].0 > angles[sorted[1] as usize].0 {
            sorted.swap (0, 1)
          }
          if angles[sorted[2] as usize].0 > angles[sorted[3] as usize].0 {
            sorted.swap (2, 3)
          }
          if angles[sorted[0] as usize].0 > angles[sorted[2] as usize].0 {
            sorted.swap (0, 2)
          }
          if angles[sorted[1] as usize].0 > angles[sorted[3] as usize].0 {
            sorted.swap (1, 3)
          }
          if angles[sorted[1] as usize].0 > angles[sorted[2] as usize].0 {
            sorted.swap (1, 2)
          }
        }
        _ => unreachable!()
      }
      sorted
    }
    fn update_support <S : Real> (
      angles  : &[(S, u32); 4],
      nangles : u32,
      sorted  : &[u32; 4],
      hull    : &Hull2 <S>,
      visited : &mut [bool],
      rect    : &mut Rect <S>
    ) -> bool {
      let points = hull.points();
      let min_angle = angles[sorted[0] as usize];
      for k in 0..nangles as usize {
        let (angle, index) = angles[sorted[k] as usize];
        if angle == min_angle.0 {
          rect.indices[index as usize] += 1;
          #[expect(clippy::cast_possible_truncation)]
          if rect.indices[index as usize] == points.len() as u32 {
            rect.indices[index as usize] = 0;
          }
        }
      }
      let bottom = rect.indices[min_angle.1 as usize];
      if visited[bottom as usize] {
        return false
      }
      visited[bottom as usize] = true;
      let mut next_index = [u32::MAX; 4];
      for k in 0u32..4 {
        next_index[k as usize] = rect.indices[((min_angle.1 + k) % 4) as usize];
      }
      rect.indices = next_index;
      let j1 = rect.indices[0] as usize;
      let j0 = if j1 == 0 {
        points.len() - 1
      } else {
        j1 - 1
      };
      rect.edge = points[j1] - points[j0];
      rect.perp = vector2 (-rect.edge.y, rect.edge.x);
      let edge_len_squared = rect.edge.magnitude_squared();
      let diff = [
        points[rect.indices[1] as usize] - points[rect.indices[3] as usize],
        points[rect.indices[2] as usize] - points[rect.indices[0] as usize]
      ];
      rect.area = rect.edge.dot (diff[0]) * rect.perp.dot (diff[1]) / edge_len_squared;
      true
    }
    Some (obb)
  }
  #[inline]
  pub fn dimensions (self) -> Vector2 <Positive <S>> {
    self.extents * Positive::unchecked (S::two())
  }
  /// X dimension
  #[inline]
  pub fn width (self) -> Positive <S> {
    self.extents.x * Positive::unchecked (S::two())
  }
  /// Y dimension
  #[inline]
  pub fn height (self) -> Positive <S> {
    self.extents.y * Positive::unchecked (S::two())
  }
  #[inline]
  pub fn area (self) -> Positive <S> {
    self.width() * self.height()
  }
  #[inline]
  pub fn corners (self) -> [Point2 <S>; 4] where
    S : Real + num::real::Real + MaybeSerDes
  {
    let mut points = [Point2::origin(); 4];
    for i in 0..2 {
      for j in 0..2 {
        let sign_x = if i % 2 == 0 { -S::one() } else { S::one() };
        let sign_y = if j % 2 == 0 { -S::one() } else { S::one() };
        points[i * 2 + j] = [*self.extents.x * sign_x, *self.extents.y * sign_y].into();
      }
    }
    let rotation = Rotation2::from_angle (self.orientation.angle());
    points.map (|point| rotation.rotate (point) + self.center.0)
  }

  pub fn aabb (self) -> Aabb2 <S> where S : Real + num::real::Real + MaybeSerDes {
    Aabb2::containing (&self.corners()).unwrap()
  }

  /// Check if point is contained by the bounding box.
  ///
  /// ```
  /// # use math_utils::*;
  /// # use math_utils::geometry::*;
  /// let x = Obb2::new (
  ///   [3.0, 6.0].into(),
  ///   [2.0, 0.75].map (Positive::noisy).into(),
  ///   AngleWrapped::wrap (Rad (std::f32::consts::FRAC_PI_4)));
  /// assert!(x.contains ([2.0f32, 5.0].into()));
  /// assert!(!x.contains ([3.0f32, 2.0].into()));
  /// ```
  pub fn contains (self, point : Point2 <S>) -> bool where
    S : Real + num::real::Real + MaybeSerDes
  {
    // re-center
    let p = point - self.center.0;
    // reverse rotation
    let p = Rotation2::from_angle (-self.orientation.angle()).rotate (p);
    p.x().abs() < *self.extents.x && p.y().abs() < *self.extents.y
  }

  /// Increase or decrease each extent by the given amount.
  ///
  /// ```
  /// # use math_utils::*;
  /// # use math_utils::geometry::*;
  /// # use math_utils::num::Zero;
  /// let x = Obb2::new (
  ///   Point2::origin(),
  ///   [0.5, 0.5].map (Positive::noisy).into(),
  ///   AngleWrapped::zero()
  /// );
  /// let y = x.dilate (1.0);
  /// assert_eq!(*y.dimensions().x, *x.dimensions().x + 2.0);
  /// assert_eq!(*y.dimensions().y, *x.dimensions().y + 2.0);
  /// ```
  ///
  /// Panics if any extent would become zero or negative:
  ///
  /// ```should_panic
  /// # use math_utils::*;
  /// # use math_utils::geometry::*;
  /// # use math_utils::num::Zero;
  /// let x = Obb2::new (
  ///   Point2::origin(),
  ///   [0.5, 0.5].map (Positive::noisy).into(),
  ///   AngleWrapped::zero()
  /// );
  /// let y = x.dilate (-1.0); // panic!
  /// ```
  pub fn dilate (mut self, amount : S) -> Self {
    self.extents = self.extents.map (|e| e.map_noisy (|e| e + amount));
    self
  }
}
impl <S : Real> From <Aabb2 <S>> for Obb2 <S> {
  /// Panics if AABB has a zero dimension:
  ///
  /// ```should_panic
  /// # use math_utils::*;
  /// # use math_utils::geometry::*;
  /// # use math_utils::num::Zero;
  /// let x = Aabb2::with_minmax (
  ///   [0.0, 0.0].into(),
  ///   [0.0, 1.0].into()
  /// ).unwrap();
  /// let y = Obb2::from (x); // panic!
  /// ```
  fn from (aabb : Aabb2 <S>) -> Self {
    use num::Zero;
    let center        = aabb.center();
    let extents = vector2 (
      Positive::noisy (*aabb.width()  / S::two()),
      Positive::noisy (*aabb.height() / S::two()));
    let orientation   = AngleWrapped::zero();
    Obb2::new (center, extents, orientation)
  }
}
impl <S> Primitive <S, Point2 <S>> for Obb2 <S> where
  S : Real + num::real::Real + MaybeSerDes
{
  fn translate (&mut self, displacement : Vector2 <S>) {
    self.center += displacement;
  }
  fn scale (&mut self, scale : Positive <S>) {
    self.extents *= scale;
  }
  fn support (&self, direction : NonZero2 <S>) -> (Point2 <S>, S) {
    use num::Inv;
    let rotation = Rotation2::from_angle (self.orientation.angle());
    let aabb = Aabb2::with_minmax_unchecked (
      self.center - self.extents.map (|e| *e),
      self.center + self.extents.map (|e| *e));
    let aabb_direction = rotation.inv() * direction;
    let (aabb_support, _) = aabb.support (aabb_direction);
    let out = rotation.rotate_around (aabb_support, self.center);
    (out, out.0.dot (*direction))
  }
}

impl_numcast_fields!(Obb3, center, extents, orientation);
impl <S : OrderedRing> Obb3 <S> {
  /// Construct a new OBB
  #[inline]
  pub const fn new (
    center      : Point3 <S>,
    extents     : Vector3 <Positive <S>>,
    orientation : Angles3 <S>
  ) -> Self {
    Obb3 { center, extents, orientation }
  }
  /// Compute the OBB of a given orientation.
  ///
  /// Returns None if the points span an empty set.
  pub fn containing_points_with_orientation (
    points : &[Point3 <S>], orientation : Angles3 <S>
  ) -> Option <Self> where
    S : Real + num::real::Real + MaybeSerDes
  {
    let [x, y, z] = Rotation3::from (orientation).cols.map (NonZero3::unchecked)
      .into_array();
    let min_dot_x = -support3 (points, -x).1;
    let max_dot_x =  support3 (points,  x).1;
    let min_dot_y = -support3 (points, -y).1;
    let max_dot_y =  support3 (points,  y).1;
    let min_dot_z = -support3 (points, -z).1;
    let max_dot_z =  support3 (points,  z).1;
    let center = Point3 (
      (*x * min_dot_x + *x * max_dot_x) * S::half() +
      (*y * min_dot_y + *y * max_dot_y) * S::half() +
      (*z * min_dot_z + *z * max_dot_z) * S::half());
    let extents = vector3 (
      Positive::new (S::half() * (max_dot_x - min_dot_x))?,
      Positive::new (S::half() * (max_dot_y - min_dot_y))?,
      Positive::new (S::half() * (max_dot_z - min_dot_z))?);
    Some (Obb3 { center, extents, orientation })
  }
  /// Construct an approximate minimum OBB containing the convex hull of a given set of
  /// points.
  ///
  /// Returns None if fewer than 4 points are given or if points span an empty set:
  /// ```
  /// # use math_utils::*;
  /// # use math_utils::geometry::*;
  /// assert!(Obb3::containing_points_approx (
  ///   &[[0.0, 1.0, 0.0], [1.0, 0.0, 0.0], [0.0, 0.0, 1.0]].map (Point3::from)
  /// ).is_none());
  ///
  /// assert!(Obb3::containing_points_approx (
  ///   &[[1.0, 0.0, 1.0], [0.0, 1.0, 0.0], [0.0, -1.0, 0.0], [-1.0,  0.0, -1.0]]
  ///     .map (Point3::from)
  /// ).is_none());
  /// ```
  pub fn containing_points_approx (points : &[Point3 <S>]) -> Option <Self> where
    S : Real + num::real::Real + approx::RelativeEq <Epsilon=S> + MaybeSerDes
  {
    if points.len() < 4 {
      return None
    }
    let (hull, mesh) = Hull3::from_points_with_mesh (points)?;
    Self::containing_approx (&hull, &mesh)
  }
  /// Construct an approximate minimum OBB containing the given convex hull and
  /// corresponding mesh.
  ///
  /// Returns None if fewer than 4 points are given or if points span an empty set:
  /// ```
  /// # use math_utils::*;
  /// # use math_utils::geometry::*;
  /// let (hull, mesh) = Hull3::from_points_with_mesh (
  ///   &[[0.0, 1.0, 0.0], [1.0, 0.0, 0.0], [0.0, 0.0, 1.0]].map (Point3::from)
  /// ).unwrap();
  /// assert!(Obb3::containing_approx (&hull, &mesh).is_none());
  ///
  /// let (hull, mesh) = Hull3::from_points_with_mesh (
  ///   &[[1.0, 0.0, 1.0], [0.0, 1.0, 0.0], [0.0, -1.0, 0.0], [-1.0,  0.0, -1.0]]
  ///     .map (Point3::from)
  /// ).unwrap();
  /// assert!(Obb3::containing_approx (&hull, &mesh).is_none());
  /// ```
  pub fn containing_approx (hull : &Hull3 <S>, mesh : &VertexEdgeTriangleMesh)
    -> Option <Self>
  where S : Real + num::NumCast + approx::RelativeEq <Epsilon=S> + MaybeSerDes {
    // algorithm from:
    // <https://github.com/isl-org/Open3D/blob/fb7088ceebef38d54c575b47935e568979b50954/cpp/open3d/t/geometry/kernel/MinimumOBB.cpp>
    use num::Inv;
    if hull.points().len() < 4 {
      return None
    }
    let mut min_obb : Obb3 <S> = Aabb3::containing (hull.points()).unwrap().into();
    let mut min_vol = min_obb.volume();
    for triangle in mesh.triangles().values() {
      let [a, b, c] = triangle.vertices().map (|i| hull.points()[i as usize]);
      let mut u = b - a;
      let mut v = c - a;
      let mut w = u.cross (v);
      v = w.cross (u);
      u = *u.normalize();
      v = *v.normalize();
      w = *w.normalize();
      let m = Matrix3::from_row_arrays (std::array::from_fn (|i| [u[i], v[i], w[i]]));
      let m_rot       = Rotation3::new_approx (m).unwrap();
      let rot_inv     = m_rot.inv();
      let center      = a;
      let rotated     = hull.points().iter().copied()
        .map (|p| rot_inv.rotate_around (p, center))
        .collect::<Vec <_>>();
      let aabb        = Aabb3::containing (rotated.as_slice()).unwrap();
      let aabb_volume = aabb.volume();
      if approx::abs_diff_eq!(*aabb_volume, S::zero()) {
        // hull contains no volume
        return None
      }
      let volume = Positive::new (*aabb_volume).unwrap();
      if volume < min_vol {
        min_vol = volume;
        min_obb = aabb.into();
        min_obb.orientation = m_rot.into();
        min_obb.center = m_rot.rotate_around (min_obb.center, center);
      }
    }
    Some (min_obb)
  }
  /// Construct an approximate minimum OBB containing the given convex hull of points
  /// oriented according to principal component analysis bases.
  ///
  /// Returns None if fewer than 4 points are given or if points span an empty set:
  /// ```
  /// # use math_utils::*;
  /// # use math_utils::geometry::*;
  /// let hull = Hull3::from_points (
  ///   &[[0.0, 1.0, 0.0], [1.0, 0.0, 0.0], [0.0, 0.0, 1.0]].map (Point3::from)
  /// ).unwrap();
  /// assert!(Obb3::containing_approx_pca (&hull).is_none());
  ///
  /// let hull = Hull3::from_points (
  ///   &[[1.0, 0.0, 1.0], [0.0, 1.0, 0.0], [0.0, -1.0, 0.0], [-1.0,  0.0, -1.0]]
  ///     .map (Point3::from)
  /// ).unwrap();
  /// assert!(Obb3::containing_approx_pca (&hull).is_none());
  /// ```
  pub fn containing_approx_pca (hull : &Hull3 <S>) -> Option <Self> where
    S : Real + num::real::Real + num::NumCast + approx::RelativeEq <Epsilon=S>
      + MaybeSerDes
  {
    if hull.points().len() < 4 {
      return None
    }
    // coplanar check
    // TODO: should we store dimensionality information in the hull to avoid doing this
    // calculation ?
    let [p1, p2, p3] = &hull.points()[..3] else { unreachable!() };
    let mut coplanar = true;
    for p in &hull.points()[3..] {
      if !coplanar_3d (*p1, *p2, *p3, *p) {
        coplanar = false;
        break
      }
    }
    if coplanar {
      return None
    }
    // do PCA analysis of points
    let orientation = data::pca_3 (hull.points()).into();
    Some (Self::containing_points_with_orientation (hull.points(), orientation).unwrap())
  }
  #[inline]
  pub fn dimensions (self) -> Vector3 <Positive <S>> {
    self.extents * Positive::unchecked (S::two())
  }
  /// X dimension
  #[inline]
  pub fn width (self) -> Positive <S> {
    self.extents.x * Positive::unchecked (S::two())
  }
  /// Y dimension
  #[inline]
  pub fn depth (self) -> Positive <S> {
    self.extents.y * Positive::unchecked (S::two())
  }
  /// Z dimension
  #[inline]
  pub fn height (self) -> Positive <S> {
    self.extents.z * Positive::unchecked (S::two())
  }
  #[inline]
  pub fn volume (self) -> Positive <S> {
    self.width() * self.depth() * self.height()
  }
  #[inline]
  pub fn corners (self) -> [Point3 <S>; 8] where
    S : Real + num::real::Real + MaybeSerDes
  {
    let mut points = [Point3::origin(); 8];
    for i in 0..2 {
      for j in 0..2 {
        for k in 0..2 {
          let sign_x = if i % 2 == 0 { -S::one() } else { S::one() };
          let sign_y = if j % 2 == 0 { -S::one() } else { S::one() };
          let sign_z = if k % 2 == 0 { -S::one() } else { S::one() };
          points[i * 4 + j * 2 + k] = [
            *self.extents.x * sign_x,
            *self.extents.y * sign_y,
            *self.extents.z * sign_z
          ].into();
        }
      }
    }
    let rotation = Rotation3::from (self.orientation);
    points.map (|point| rotation.rotate (point) + self.center.0)
  }

  pub fn aabb (self) -> Aabb3 <S> where S : Real + num::real::Real + MaybeSerDes {
    Aabb3::containing (&self.corners()).unwrap()
  }

  /// Check if point is contained by the bounding box
  pub fn contains (self, point : Point3 <S>) -> bool where
    S : Real + num::real::Real + MaybeSerDes
  {
    use num::Inv;
    // re-center
    let p = point - self.center.0;
    // reverse rotation
    let p = Rotation3::from (self.orientation).inv().rotate (p);
    p.x().abs() < *self.extents.x &&
    p.y().abs() < *self.extents.y &&
    p.z().abs() < *self.extents.z
  }

  /// Increase or decrease each extent by the given amount.
  ///
  /// Panics if any extent would become zero or negative:
  ///
  /// ```should_panic
  /// # use math_utils::*;
  /// # use math_utils::geometry::*;
  /// let x = Obb3::new (
  ///   Point3::origin(),
  ///   [0.5, 0.5, 0.5].map (Positive::noisy).into(),
  ///   Angles3::zero()
  /// );
  /// let y = x.dilate (-1.0); // panic!
  /// ```
  pub fn dilate (mut self, amount : S) -> Self {
    self.extents = self.extents.map (|e| e.map_noisy (|e| e + amount));
    self
  }
}
impl <S : Real> From <Aabb3 <S>> for Obb3 <S> {
  /// Panics if AABB has a zero dimension:
  ///
  /// ```should_panic
  /// # use math_utils::*;
  /// # use math_utils::geometry::*;
  /// # use math_utils::num::Zero;
  /// let x = Aabb3::with_minmax (
  ///   [0.0, 0.0, 0.0].into(),
  ///   [0.0, 1.0, 1.0].into()
  /// ).unwrap();
  /// let y = Obb3::from (x); // panic!
  /// ```
  fn from (aabb : Aabb3 <S>) -> Self {
    let center = aabb.center();
    let extents = aabb.extents().map (|d| Positive::noisy (*d));
    let orientation = Angles3::zero();
    Obb3::new (center, extents, orientation)
  }
}
impl <S> Primitive <S, Point3 <S>> for Obb3 <S> where
  S : Real + num::real::Real + MaybeSerDes
{
  fn translate (&mut self, displacement : Vector3 <S>) {
    self.center += displacement;
  }
  fn scale (&mut self, scale : Positive <S>) {
    self.extents *= scale;
  }
  fn support (&self, direction : NonZero3 <S>) -> (Point3 <S>, S) {
    use num::Inv;
    let rotation = Rotation3::from (self.orientation);
    let aabb = Aabb3::with_minmax_unchecked (
      self.center - self.extents.map (|e| *e),
      self.center + self.extents.map (|e| *e));
    let aabb_direction = rotation.inv() * direction;
    let (aabb_support, _) = aabb.support (aabb_direction);
    let out = rotation.rotate_around (aabb_support, self.center);
    (out, out.0.dot (*direction))
  }
}

impl_numcast_fields!(Capsule3, center, half_height, radius);
impl <S : OrderedRing> Capsule3 <S> {
  pub fn aabb3 (self) -> Aabb3 <S> where S : Real {
    use shape::Aabb;
    let shape_aabb = shape::Capsule {
      radius:      self.radius,
      half_height: self.half_height
    }.aabb();
    let center_vec = self.center.0;
    let min        = shape_aabb.min() + center_vec;
    let max        = shape_aabb.max() + center_vec;
    Aabb3::with_minmax_unchecked (min, max)
  }
  /// Return the (upper sphere, cylinder, lower sphere) making up this capsule
  pub fn decompose (self) -> (Sphere3 <S>, Option <Cylinder3 <S>>, Sphere3 <S>) {
    let cylinder = if *self.half_height > S::zero() {
      Some (Cylinder3 {
        center:      self.center,
        radius:      self.radius,
        half_height: Positive::unchecked (*self.half_height)
      })
    } else {
      None
    };
    let upper_sphere = Sphere3 {
      center: self.center + (Vector3::unit_z() * *self.half_height),
      radius: self.radius
    };
    let lower_sphere = Sphere3 {
      center: self.center - (Vector3::unit_z() * *self.half_height),
      radius: self.radius
    };
    (upper_sphere, cylinder, lower_sphere)
  }
}
impl <S> Primitive <S, Point3 <S>> for Capsule3 <S> where
  S : Real + num::real::Real + MaybeSerDes
{
  fn translate (&mut self, displacement : Vector3 <S>) {
    self.center += displacement;
  }
  fn scale (&mut self, scale : Positive <S>) {
    self.half_height *= scale;
    self.radius *= scale;
  }
  fn support (&self, direction : NonZero3 <S>) -> (Point3 <S>, S) {
    let dir_unit = direction.normalized();
    let out = self.center + if direction.z > S::zero() {
      Vector3::unit_z() * *self.half_height + dir_unit * *self.radius
    } else if direction.z < S::zero() {
      -Vector3::unit_z() * *self.half_height + dir_unit * *self.radius
    } else {
      dir_unit * *self.radius
    };
    (out, out.0.dot (*direction))
  }
}

impl_numcast_fields!(Cylinder3, center, half_height, radius);
impl <S : OrderedRing> Cylinder3 <S> {
  pub fn aabb3 (self) -> Aabb3 <S> where S : Real {
    use shape::Aabb;
    let shape_aabb = shape::Cylinder::unchecked (*self.radius, *self.half_height).aabb();
    let center_vec = self.center.0;
    let min        = shape_aabb.min() + center_vec;
    let max        = shape_aabb.max() + center_vec;
    Aabb3::with_minmax_unchecked (min, max)
  }
}
impl <S> Primitive <S, Point3 <S>> for Cylinder3 <S> where
  S : Real + num::real::Real + MaybeSerDes
{
  fn translate (&mut self, displacement : Vector3 <S>) {
    self.center += displacement;
  }
  fn scale (&mut self, scale : Positive <S>) {
    self.half_height *= scale;
    self.radius *= scale;
  }
  fn support (&self, direction : NonZero3 <S>) -> (Point3 <S>, S) {
    let out = self.center + if *direction == Vector3::unit_z() {
      Vector3::unit_z() * *self.half_height
    } else if *direction == -Vector3::unit_z() {
      -Vector3::unit_z() * *self.half_height
    } else if direction.z == S::zero() {
      direction.normalized() * *self.radius
    } else if direction.z > S::zero() {
      Vector3::unit_z() * *self.half_height
        + vector3 (direction.x, direction.y, S::zero()).normalized() * *self.radius
    } else {
      debug_assert!(direction.z < S::zero());
      -Vector3::unit_z() * *self.half_height
        + vector3 (direction.x, direction.y, S::zero()).normalized() * *self.radius
    };
    (out, out.0.dot (*direction))
  }
}

impl_numcast_fields!(Cone3, center, half_height, radius);
impl <S : OrderedRing> Cone3 <S> {
  pub fn height (self) -> Positive <S> {
    self.half_height * Positive::unchecked (S::two())
  }
  pub fn aabb3 (self) -> Aabb3 <S> where S : Real {
    use shape::Aabb;
    let shape_aabb = shape::Cone {
      radius:      self.radius,
      half_height: self.half_height
    }.aabb();
    let center_vec = self.center.0;
    let min        = shape_aabb.min() + center_vec;
    let max        = shape_aabb.max() + center_vec;
    Aabb3::with_minmax_unchecked (min, max)
  }
}
impl <S> Primitive <S, Point3 <S>> for Cone3 <S> where
  S : Real + num::real::Real + approx::RelativeEq + MaybeSerDes
{
  fn translate (&mut self, displacement : Vector3 <S>) {
    self.center += displacement;
  }
  fn scale (&mut self, scale : Positive <S>) {
    self.half_height *= scale;
    self.radius *= scale;
  }
  fn support (&self, direction : NonZero3 <S>) -> (Point3 <S>, S) {
    let top = || Vector3::unit_z() * *self.half_height;
    let bottom = || -top();
    let edge = || bottom()
      + vector3 (direction.x, direction.y, S::zero()).normalized() * *self.radius;
    let out = self.center + if *direction == Vector3::unit_z() {
      top()
    } else if *direction == -Vector3::unit_z() {
      bottom()
    } else if direction.z <= S::zero() {
      edge()
    } else {
      // TODO: possibly cheaper with dot products ?
      let angle_dir =
        Trig::atan2 (direction.z, Sqrt::sqrt (S::one() - direction.z.squared()));
      debug_assert!(angle_dir > S::zero());
      let angle_cone = Trig::atan2 (-*self.height(), *self.radius);
      debug_assert!(angle_cone < S::zero());
      let pi_2 = S::pi() / S::two();
      let angle_total = angle_dir - angle_cone;
      if approx::relative_eq!(angle_total, pi_2) {
        (top() + edge()) * S::half()
      } else if angle_total > pi_2 {
        top()
      } else {
        debug_assert!(angle_total < pi_2);
        edge()
      }
    };
    (out, out.0.dot (*direction))
  }
}

impl_numcast_fields!(Line2, base, direction);
impl <S : OrderedRing> Line2 <S> {
  /// Construct a new 2D line
  #[inline]
  pub const fn new (base : Point2 <S>, direction : Unit2 <S>) -> Self {
    Line2 { base, direction }
  }
  #[inline]
  pub fn point (self, t : S) -> Point2 <S> {
    self.base + *self.direction * t
  }
  #[inline]
  pub fn affine_line (self) -> frame::Line2 <S> {
    frame::Line2 {
      origin: self.base,
      basis:  self.direction.into()
    }
  }
  #[inline]
  pub fn nearest_point (self, point : Point2 <S>) -> Line2Point <S> {
    project_2d_point_on_line (point, self)
  }
  #[inline]
  pub fn intersect_aabb (self, aabb : Aabb2 <S>)
    -> Option <(Line2Point <S>, Line2Point <S>)>
  {
    intersect::continuous_line2_aabb2 (self, aabb)
  }
  #[inline]
  pub fn intersect_sphere (self, sphere : Sphere2 <S>)
    -> Option <(Line2Point <S>, Line2Point <S>)>
  where S : Field + Sqrt {
    intersect::continuous_line2_sphere2 (self, sphere)
  }
}
impl <S : Real> Default for Line2 <S> {
  fn default() -> Self {
    Line2 {
      base:      Point2::origin(),
      direction: Unit2::axis_y()
    }
  }
}
impl <S> Primitive <S, Point2 <S>> for Line2 <S> where
  S : OrderedField + num::float::FloatCore + approx::AbsDiffEq
{
  fn translate (&mut self, displacement : Vector2 <S>) {
    self.base += displacement;
  }
  #[expect(unused_variables)]
  fn scale (&mut self, scale : Positive <S>) {
    /* no-op */
  }
  fn support (&self, direction : NonZero2 <S>) -> (Point2 <S>, S) {
    let dot = self.direction.dot (*direction);
    let out = if approx::abs_diff_eq!(dot, S::zero()) {
      self.base
    } else if dot > S::zero() {
      Point2 (self.direction.sigvec().map (|x| if x == S::zero() {
        S::zero()
      } else {
        x * S::infinity()
      }))
    } else {
      debug_assert!(dot < S::zero());
      Point2 ((-self.direction).sigvec().map (|x| if x == S::zero() {
        S::zero()
      } else {
        x * S::infinity()
      }))
    };
    (out, out.0.dot (*direction))
  }
}
impl <S> From <Ray2 <S>> for Line2 <S> {
  fn from (Ray2 { base, direction } : Ray2 <S>) -> Line2 <S> {
    Line2 { base, direction }
  }
}

impl_numcast_fields!(Ray2, base, direction);
impl <S : OrderedRing> Ray2 <S> {
  /// Construct a new 2D ray
  #[inline]
  pub const fn new (base : Point2 <S>, direction : Unit2 <S>) -> Self {
    Ray2 { base, direction }
  }
  #[inline]
  pub fn point (self, t : S) -> Point2 <S> {
    self.base + *self.direction * t
  }
  #[inline]
  pub fn affine_line (self) -> frame::Line2 <S> {
    frame::Line2 {
      origin: self.base,
      basis:  self.direction.into()
    }
  }
  pub fn nearest_point (self, point : Point2 <S>) -> Ray2Point <S> {
    use num::Zero;
    let (t, point) = Line2::from (self).nearest_point (point);
    if t < S::zero() {
      (NonNegative::zero(), self.base)
    } else {
      (NonNegative::unchecked (t), point)
    }
  }
  #[inline]
  pub fn intersect_aabb (self, aabb : Aabb2 <S>)
    -> Option <(Ray2Point <S>, Ray2Point <S>)>
  {
    intersect::continuous_ray2_aabb2 (self, aabb)
  }
  #[inline]
  pub fn intersect_sphere (self, sphere : Sphere2 <S>)
    -> Option <(Ray2Point <S>, Ray2Point <S>)>
  {
    intersect::continuous_ray2_sphere2 (self, sphere)
  }
}
impl <S : Real> Default for Ray2 <S> {
  fn default() -> Self {
    Ray2 {
      base:      Point2::origin(),
      direction: Unit2::axis_y()
    }
  }
}
impl <S> Primitive <S, Point2 <S>> for Ray2 <S> where
  S : OrderedField + num::float::FloatCore + approx::AbsDiffEq <Epsilon = S>
{
  fn translate (&mut self, displacement : Vector2 <S>) {
    self.base += displacement;
  }
  #[expect(unused_variables)]
  fn scale (&mut self, scale : Positive <S>) {
    /* no-op */
  }
  fn support (&self, direction : NonZero2 <S>) -> (Point2 <S>, S) {
    let dot = self.direction.dot (*direction);
    let out = if dot < S::default_epsilon() {
      self.base
    } else {
      Point2 (self.direction.sigvec().map (|x| if x == S::zero() {
        S::zero()
      } else {
        x * S::infinity()
      }))
    };
    (out, out.0.dot (*direction))
  }
}
impl <S> From <Line2 <S>> for Ray2 <S> {
  fn from (Line2 { base, direction } : Line2 <S>) -> Ray2 <S> {
    Ray2 { base, direction }
  }
}

impl_numcast_fields!(Line3, base, direction);
impl <S : OrderedRing> Line3 <S> {
  /// Consructs a new 3D line with given base point and direction vector
  #[inline]
  pub const fn new (base : Point3 <S>, direction : Unit3 <S>) -> Self {
    Line3 { base, direction }
  }
  #[inline]
  pub fn point (self, t : S) -> Point3 <S> {
    self.base + *self.direction * t
  }
  #[inline]
  pub fn affine_line (self) -> frame::Line3 <S> {
    frame::Line3 {
      origin: self.base,
      basis:  self.direction.into()
    }
  }
  #[inline]
  pub fn nearest_point (self, point : Point3 <S>) -> Line3Point <S> {
    project_3d_point_on_line (point, self)
  }
  #[inline]
  pub fn intersect_plane (self, plane : Plane3 <S>) -> Option <Line3Point <S>> where
    S : Field + approx::RelativeEq
  {
    intersect::continuous_line3_plane3 (self, plane)
  }
  #[inline]
  pub fn intersect_aabb (self, aabb : Aabb3 <S>)
    -> Option <(Line3Point <S>, Line3Point <S>)>
  where S : num::Float + approx::RelativeEq <Epsilon=S> {
    intersect::continuous_line3_aabb3 (self, aabb)
  }
  #[inline]
  pub fn intersect_sphere (self, sphere : Sphere3 <S>)
    -> Option <(Line3Point <S>, Line3Point <S>)>
  where S : Field + Sqrt {
    intersect::continuous_line3_sphere3 (self, sphere)
  }
  #[inline]
  pub fn intersect_triangle (self, triangle : Triangle3 <S>)
    -> Option <Line3Point <S>>
  where S : OrderedField + num::real::Real + approx::AbsDiffEq <Epsilon = S> {
    intersect::continuous_line3_triangle3 (self, triangle)
  }
}
impl <S : Real> Default for Line3 <S> {
  fn default() -> Self {
    Line3 {
      base:      Point3::origin(),
      direction: Unit3::axis_z()
    }
  }
}
impl <S> Primitive <S, Point3 <S>> for Line3 <S> where
  S : OrderedField + num::float::FloatCore + approx::AbsDiffEq
{
  fn translate (&mut self, displacement : Vector3 <S>) {
    self.base += displacement;
  }
  #[expect(unused_variables)]
  fn scale (&mut self, scale : Positive <S>) {
    /* no-op */
  }
  fn support (&self, direction : NonZero3 <S>) -> (Point3 <S>, S) {
    let dot = self.direction.dot (*direction);
    let out = if approx::abs_diff_eq!(dot, S::zero()) {
      self.base
    } else if dot > S::zero() {
      Point3 (self.direction.sigvec().map (|x| if x == S::zero() {
        S::zero()
      } else {
        x * S::infinity()
      }))
    } else {
      debug_assert!(dot < S::zero());
      Point3 (-self.direction.sigvec().map (|x| if x == S::zero() {
        S::zero()
      } else {
        x * S::infinity()
      }))
    };
    (out, out.0.dot (*direction))
  }
}
impl <S> From <Ray3 <S>> for Line3 <S> {
  fn from (Ray3 { base, direction } : Ray3 <S>) -> Line3 <S> {
    Line3 { base, direction }
  }
}

impl_numcast_fields!(Ray3, base, direction);
impl <S : OrderedRing> Ray3 <S> {
  /// Consructs a new 3D ray with given base point and direction vector
  #[inline]
  pub const fn new (base : Point3 <S>, direction : Unit3 <S>) -> Self {
    Ray3 { base, direction }
  }
  #[inline]
  pub fn point (self, t : S) -> Point3 <S> {
    self.base + *self.direction * t
  }
  #[inline]
  pub fn affine_line (self) -> frame::Line3 <S> {
    frame::Line3 {
      origin: self.base,
      basis:  self.direction.into()
    }
  }
  pub fn nearest_point (self, point : Point3 <S>) -> Ray3Point <S> {
    use num::Zero;
    let (t, point) = Line3::from (self).nearest_point (point);
    if t < S::zero() {
      (NonNegative::zero(), self.base)
    } else {
      (NonNegative::unchecked (t), point)
    }
  }
  #[inline]
  pub fn intersect_plane (self, plane : Plane3 <S>) -> Option <Ray3Point <S>> where
    S : approx::RelativeEq
  {
    intersect::continuous_ray3_plane3 (self, plane)
  }
  #[inline]
  pub fn intersect_aabb (self, aabb : Aabb3 <S>)
    -> Option <(Ray3Point <S>, Ray3Point <S>)>
  where S : num::Float + approx::RelativeEq <Epsilon=S> {
    intersect::continuous_ray3_aabb3 (self, aabb)
  }
  #[inline]
  pub fn intersect_sphere (self, sphere : Sphere3 <S>)
    -> Option <(Ray3Point <S>, Ray3Point <S>)>
  {
    intersect::continuous_ray3_sphere3 (self, sphere)
  }
  #[inline]
  pub fn intersect_triangle (self, triangle : Triangle3 <S>)
    -> Option <Ray3Point <S>>
  where S : Real + num::real::Real + approx::AbsDiffEq <Epsilon = S> {
    intersect::continuous_ray3_triangle3 (self, triangle)
  }
}
impl <S : Real> Default for Ray3 <S> {
  fn default() -> Self {
    Ray3 {
      base:      Point3::origin(),
      direction: Unit3::axis_z()
    }
  }
}
impl <S> Primitive <S, Point3 <S>> for Ray3 <S> where
  S : OrderedField + num::float::FloatCore + approx::AbsDiffEq <Epsilon = S>
{
  fn translate (&mut self, displacement : Vector3 <S>) {
    self.base += displacement;
  }
  #[expect(unused_variables)]
  fn scale (&mut self, scale : Positive <S>) {
    /* no-op */
  }
  fn support (&self, direction : NonZero3 <S>) -> (Point3 <S>, S) {
    let dot = self.direction.dot (*direction);
    let out = if dot < S::default_epsilon() {
      self.base
    } else {
      Point3 (self.direction.sigvec().map (|x| if x == S::zero() {
        S::zero()
      } else {
        x * S::infinity()
      }))
    };
    (out, out.0.dot (*direction))
  }
}
impl <S> From <Line3 <S>> for Ray3 <S> {
  fn from (Line3 { base, direction } : Line3 <S>) -> Ray3 <S> {
    Ray3 { base, direction }
  }
}

impl_numcast_fields!(Plane3, base, normal);
impl <S> Plane3 <S> {
  /// Consructs a new 3D plane with given base point and normal vector
  #[inline]
  pub const fn new (base : Point3 <S>, normal : Unit3 <S>) -> Self {
    Plane3 { base, normal }
  }
}
impl <S : Real> Default for Plane3 <S> {
  fn default() -> Self {
    Plane3 {
      base:   Point3::origin(),
      normal: Unit3::axis_z()
    }
  }
}
impl <S> Primitive <S, Point3 <S>> for Plane3 <S> where
  S : Real + num::float::FloatCore + approx::RelativeEq
{
  fn translate (&mut self, displacement : Vector3 <S>) {
    self.base += displacement;
  }
  #[expect(unused_variables)]
  fn scale (&mut self, scale : Positive <S>) {
    /* no-op */
  }
  fn support (&self, direction : NonZero3 <S>) -> (Point3 <S>, S) {
    let dot = self.normal.dot (*direction);
    let out = if approx::relative_eq!(dot, S::one()) {
      self.base
    } else {
      let projected =
        project_3d_point_on_plane (Point3 (self.base.0 + *direction), *self).0;
      Point3 (projected.sigvec().map (|x| if x == S::zero() {
        S::zero()
      } else {
        x * S::infinity()
      }))
    };
    (out, out.0.dot (*direction))
  }
}

impl_numcast_fields!(Sphere2, center, radius);
impl <S : OrderedRing> Sphere2 <S> {
  /// Unit circle
  #[inline]
  pub fn unit() -> Self {
    Sphere2 {
      center: Point2::origin(),
      radius: num::One::one()
    }
  }
  /// Discrete intersection test with another sphere
  #[inline]
  pub fn intersects (self, other : Self) -> bool {
    intersect::discrete_sphere2_sphere2 (self, other)
  }
}
impl <S> Primitive <S, Point2 <S>> for Sphere2 <S> where
  S : Real + num::real::Real + MaybeSerDes
{
  fn translate (&mut self, displacement : Vector2 <S>) {
    self.center += displacement;
  }
  fn scale (&mut self, scale : Positive <S>) {
    self.radius *= scale;
  }
  fn support (&self, direction : NonZero2 <S>) -> (Point2 <S>, S) {
    let out = self.center + direction.normalized() * *self.radius;
    (out, out.0.dot (*direction))
  }
}

impl_numcast_fields!(Sphere3, center, radius);
impl <S : OrderedRing> Sphere3 <S> {
  /// Unit sphere
  #[inline]
  pub fn unit() -> Self where S : Field {
    Sphere3 {
      center: Point3::origin(),
      radius: num::One::one()
    }
  }
  /// Discrete intersection test with another sphere
  #[inline]
  pub fn intersects (self, other : Self) -> bool {
    intersect::discrete_sphere3_sphere3 (self, other)
  }
}
impl <S> Primitive <S, Point3 <S>> for Sphere3 <S> where
  S : Real + num::real::Real + std::fmt::Debug + MaybeSerDes
{
  fn translate (&mut self, displacement : Vector3 <S>) {
    self.center += displacement;
  }
  fn scale (&mut self, scale : Positive <S>) {
    self.radius *= scale;
  }
  fn support (&self, direction : NonZero3 <S>) -> (Point3 <S>, S) {
    let out = self.center + direction.normalized() * *self.radius;
    (out, out.0.dot (*direction))
  }
}

////////////////////////////////////////////////////////////////////////////////////////
//  tests                                                                             //
////////////////////////////////////////////////////////////////////////////////////////

#[cfg(test)]
mod tests {
  use approx::{assert_relative_eq, assert_ulps_eq};
  use rand;
  use rand_distr;
  use rand_xorshift;
  use sorted_vec::partial::SortedSet;
  use super::*;

  #[test]
  fn colinear_2() {
    use rand::{RngExt, SeedableRng};
    use rand_xorshift::XorShiftRng;
    use rand_distr::Distribution;
    macro test ($t:ty) {
      let mut rng = XorShiftRng::seed_from_u64 (0);
      let normal  = rand_distr::StandardNormal;
      let cauchy  = rand_distr::Cauchy::new (0.0, 1.0).unwrap();
      let randn   = |rng : &mut XorShiftRng| normal.sample (rng);
      let randf   = |rng : &mut XorShiftRng| cauchy.sample (rng);
      // colinear
      for _ in 0..100000 {
        let dir  = Unit2::<$t>::normalize ([randn (&mut rng), randn (&mut rng)].into());
        let base = point2 (randf (&mut rng), randf (&mut rng));
        let line = Line2::new (base, dir);
        let p1   = line.point (randf (&mut rng));
        let p2   = line.point (randf (&mut rng));
        let p3   = line.point (randf (&mut rng));
        assert!(colinear_2d (p1, p2, p3));
      }
      // not colinear
      for _ in 0..100000 {
        let p1        = point2 (randf (&mut rng), randf (&mut rng));
        let v         = p1.0;
        let len       = 0.1 + randf (&mut rng).abs();
        let angle     = Deg (rng.random_range (-180.0..180.0)).into();
        let rotation  = Rotation2::from_angle (angle);
        let p2        = p1 + rotation * (v.normalized() * len);
        let v         = p2 - p1;
        let len       = 0.1 + randf (&mut rng).abs();
        let mut angle = rng.random_range (5.0..175.0);
        if rng.random_bool (0.5) {
          angle = -angle;
        }
        let rotation = Rotation2::from_angle (Deg (angle).into());
        let p3       = p2 + rotation * v * len;
        assert!(!colinear_2d (p1, p2, p3));
      }
    }
    test!(f32);
    test!(f64);
  }

  #[test]
  fn colinear_3() {
    use rand::{RngExt, SeedableRng};
    use rand_xorshift::XorShiftRng;
    use rand_distr::Distribution;
    macro test ($t:ty) {
      let mut rng = XorShiftRng::seed_from_u64 (0);
      let normal  = rand_distr::StandardNormal;
      let cauchy  = rand_distr::Cauchy::new (0.0, 1.0).unwrap();
      let randf   = |rng : &mut XorShiftRng| cauchy.sample (rng);
      let randn   = |rng : &mut XorShiftRng| normal.sample (rng);
      // colinear
      for _ in 0..50000 {
        let dir = Unit3::<$t>::normalize (
          [randn (&mut rng), randn (&mut rng), randn (&mut rng)].into());
        let base = point3 (randf (&mut rng), randf (&mut rng), randf (&mut rng));
        let line = Line3::new (base, dir);
        let p1   = line.point (randf (&mut rng));
        let p2   = line.point (randf (&mut rng));
        let p3   = line.point (randf (&mut rng));
        assert!(colinear_3d (p1, p2, p3));
      }
      // not colinear
      for _ in 0..50000 {
        let p1       = point3 (randf (&mut rng), randf (&mut rng), randf (&mut rng));
        let v        = p1.0;
        let len      = 0.1 + randf (&mut rng).abs();
        let axis     = vector3 (randn (&mut rng), randn (&mut rng), randn (&mut rng));
        let angle    = Deg (rng.random_range (-180.0..180.0)).into();
        let rotation = Rotation3::from_axis_angle (axis, angle);
        let p2       = p1 + rotation * (v.normalized() * len);
        let v        = p2 - p1;
        let len      = 0.1 + randf (&mut rng).abs();
        let axis     = loop {
          let axis = vector3 (randn (&mut rng), randn (&mut rng), randn (&mut rng));
          if axis.normalized().cross (v.normalized()).magnitude() > 0.1 {
            break axis
          }
        };
        let mut angle = rng.random_range (5.0..175.0);
        if rng.random_bool (0.5) {
          angle = -angle;
        }
        let rotation = Rotation3::from_axis_angle (axis, Deg (angle).into());
        let p3       = p2 + rotation * v * len;
        assert!(!colinear_3d (p1, p2, p3));
      }
    }
    test!(f32);
    test!(f64);
  }

  #[test]
  fn coplanar_3() {
    use rand::SeedableRng;
    use rand_xorshift::XorShiftRng;
    use rand_distr::Distribution;
    macro test ($t:ty) {
      let mut rng = XorShiftRng::seed_from_u64 (0);
      let cauchy  = rand_distr::Cauchy::new (0.0, 1.0).unwrap();
      let randf   = |rng : &mut XorShiftRng| cauchy.sample (rng);
      // coplanar
      for _ in 0..50000 {
        let p1   = point3 (randf (&mut rng), randf (&mut rng), randf (&mut rng));
        let p2   = point3 (randf (&mut rng), randf (&mut rng), randf (&mut rng));
        let p3   = point3 (randf (&mut rng), randf (&mut rng), randf (&mut rng));
        let s    = randf (&mut rng);
        let t    = randf (&mut rng);
        let p4   = p1 + (p2 - p1) * s + (p3 - p1) * t;
        assert!(coplanar_3d::<$t> (p1, p2, p3, p4));
      }
      // not coplanar
      for _ in 0..50000 {
        let p1    = point3 (randf (&mut rng), randf (&mut rng), randf (&mut rng));
        let p2    = point3 (randf (&mut rng), randf (&mut rng), randf (&mut rng));
        let p3    = point3 (randf (&mut rng), randf (&mut rng), randf (&mut rng));
        let s     = randf (&mut rng);
        let t     = randf (&mut rng);
        let e1    = p2 - p1;
        let e2    = p3 - p1;
        let n     = e1.cross (e2).normalized();
        let rn    = randf (&mut rng);
        let pnew  = p1 + e1 * s + e2 * t;
        let n_len = (0.3 * (pnew - p1).magnitude()).mul_add (rn.signum(), rn);
        //let n_len = rn + 0.3 * (pnew - p1).magnitude() * rn.signum();
        let p4    = pnew + n * n_len;
        assert!(!coplanar_3d::<$t> (p1, p2, p3, p4));
      }
    }
    test!(f32);
    test!(f64);
  }

  #[test]
  fn triangle_area_squared() {
    assert_relative_eq!(
      81.0,
      *triangle_3d_area_squared (
        [-2.0, -1.0,  0.0].into(),
        [ 1.0, -1.0,  0.0].into(),
        [ 1.0,  5.0,  0.0].into())
    );
    assert_relative_eq!(
      20.25,
      *triangle_3d_area_squared (
        [-2.0, -1.0,  0.0].into(),
        [ 1.0, -1.0,  0.0].into(),
        [ 1.0,  2.0,  0.0].into())
    );
    assert_relative_eq!(
      20.25,
      *triangle_3d_area_squared (
        [ 1.0, -1.0,  0.0].into(),
        [-2.0, -1.0,  0.0].into(),
        [ 1.0,  2.0,  0.0].into())
    );
    assert_relative_eq!(
      20.25,
      *triangle_3d_area_squared (
        [ 1.0, -1.0,  0.0].into(),
        [ 1.0,  2.0,  0.0].into(),
        [-2.0, -1.0,  0.0].into())
    );
  }

  #[test]
  fn project_2d_point_on_line() {
    use super::project_2d_point_on_line;
    use Unit2;
    let point : Point2 <f64> = [2.0, 2.0].into();
    let line = Line2::<f64>::new ([0.0, 0.0].into(), Unit2::axis_x());
    assert_eq!(project_2d_point_on_line (point, line).1, [2.0, 0.0].into());
    let line = Line2::<f64>::new ([0.0, 0.0].into(), Unit2::axis_y());
    assert_eq!(project_2d_point_on_line (point, line).1, [0.0, 2.0].into());
    let point : Point2 <f64> = [0.0, 1.0].into();
    let line = Line2::<f64>::new (
      [0.0, -1.0].into(), Unit2::normalize ([1.0, 1.0].into()));
    assert_relative_eq!(
      project_2d_point_on_line (point, line).1, [1.0, 0.0].into());
    // the answer should be the same for lines with equivalent definitions
    let point : Point2 <f64> = [1.0, 3.0].into();
    let line_a = Line2::<f64>::new (
      [0.0, -1.0].into(), Unit2::normalize ([2.0, 1.0].into()));
    let line_b = Line2::<f64>::new (
      [2.0, 0.0].into(),  Unit2::normalize ([2.0, 1.0].into()));
    assert_relative_eq!(
      project_2d_point_on_line (point, line_a).1,
      project_2d_point_on_line (point, line_b).1);
    let line_a = Line2::<f64>::new (
      [0.0, 0.0].into(),   Unit2::normalize ([1.0, 1.0].into()));
    let line_b = Line2::<f64>::new (
      [-2.0, -2.0].into(), Unit2::normalize ([1.0, 1.0].into()));
    assert_ulps_eq!(
      project_2d_point_on_line (point, line_a).1,
      project_2d_point_on_line (point, line_b).1
    );
    assert_relative_eq!(
      project_2d_point_on_line (point, line_a).1,
      [2.0, 2.0].into());
  }

  #[test]
  fn project_3d_point_on_line() {
    use Unit3;
    use super::project_3d_point_on_line;
    // all the tests from 2d projection with 0.0 for the Z component
    let point : Point3 <f64> = [2.0, 2.0, 0.0].into();
    let line = Line3::<f64>::new ([0.0, 0.0, 0.0].into(), Unit3::axis_x());
    assert_eq!(project_3d_point_on_line (point, line).1, [2.0, 0.0, 0.0].into());
    let line = Line3::<f64>::new ([0.0, 0.0, 0.0].into(), Unit3::axis_y());
    assert_eq!(project_3d_point_on_line (point, line).1, [0.0, 2.0, 0.0].into());
    let point : Point3 <f64> = [0.0, 1.0, 0.0].into();
    let line = Line3::<f64>::new (
      [0.0, -1.0, 0.0].into(), Unit3::normalize ([1.0, 1.0, 0.0].into()));
    assert_relative_eq!(
      project_3d_point_on_line (point, line).1, [1.0, 0.0, 0.0].into());
    // the answer should be the same for lines with equivalent definitions
    let point : Point3 <f64> = [1.0, 3.0, 0.0].into();
    let line_a = Line3::<f64>::new (
      [0.0, -1.0, 0.0].into(), Unit3::normalize ([2.0, 1.0, 0.0].into()));
    let line_b = Line3::<f64>::new (
      [2.0, 0.0, 0.0].into(),  Unit3::normalize ([2.0, 1.0, 0.0].into()));
    assert_relative_eq!(
      project_3d_point_on_line (point, line_a).1,
      project_3d_point_on_line (point, line_b).1);
    let line_a = Line3::<f64>::new (
      [0.0, 0.0, 0.0].into(), Unit3::normalize ([1.0, 1.0, 0.0].into()));
    let line_b = Line3::<f64>::new (
      [2.0, 2.0, 0.0].into(), Unit3::normalize ([1.0, 1.0, 0.0].into()));
    assert_relative_eq!(
      project_3d_point_on_line (point, line_a).1,
      project_3d_point_on_line (point, line_b).1);
    assert_relative_eq!(
      project_3d_point_on_line (point, line_a).1,
      [2.0, 2.0, 0.0].into());
    // more tests
    let point : Point3 <f64> = [0.0, 0.0, 2.0].into();
    let line_a = Line3::<f64>::new (
      [-4.0, -4.0, -4.0].into(), Unit3::normalize ([1.0, 1.0, 1.0].into()));
    let line_b = Line3::<f64>::new (
      [4.0, 4.0, 4.0].into(),    Unit3::normalize ([1.0, 1.0, 1.0].into()));
    assert_relative_eq!(
      project_3d_point_on_line (point, line_a).1,
      project_3d_point_on_line (point, line_b).1,
      epsilon = 0.000_000_000_000_01
    );
    assert_relative_eq!(
      project_3d_point_on_line (point, line_a).1,
      [2.0/3.0, 2.0/3.0, 2.0/3.0].into(),
      epsilon = 0.000_000_000_000_01
    );
  }

  #[test]
  fn smallest_rect() {
    use super::smallest_rect;
    // points are in counter-clockwise order
    let points : Vec <Point2 <f32>> = [
      [-3.0, -3.0],
      [ 3.0, -3.0],
      [ 3.0,  3.0],
      [ 0.0,  6.0],
      [-3.0,  3.0]
    ].map (Point2::from).into_iter().collect();
    let hull = Hull2::from_points (points.as_slice()).unwrap();
    assert_eq!(hull.points(), points);
    let rect = smallest_rect (0, 1, &hull);
    assert_eq!(rect.area, 54.0);
    // points are in counter-clockwise order
    let points : Vec <Point2 <f32>> = [
      [-1.0, -4.0],
      [ 2.0,  2.0],
      [-4.0, -1.0]
    ].map (Point2::from).into_iter().collect();
    let hull = Hull2::from_points (points.as_slice()).unwrap();
    assert_eq!(hull.points(), points);
    let rect = smallest_rect (0, 1, &hull);
    assert_eq!(rect.area, 27.0);
  }

  #[test]
  fn capsule3() {
    use num::One;
    let capsule : Capsule3 <f32> = Capsule3 {
      center:      Point3::origin(),
      radius:      Positive::noisy (2.0),
      half_height: Positive::one()
    };
    let support = capsule.support (Unit3::axis_z().into()).0;
    assert_eq!(support, point3 (0.0, 0.0, 3.0));
    let support = capsule.support (Unit3::axis_x().into()).0;
    assert_eq!(support, point3 (2.0, 0.0, 0.0));
    let support = capsule.support (NonZero3::noisy (vector3 (1.0, 0.0, 1.0))).0;
    assert_eq!(support,
      point3 (0.0, 0.0, 1.0) + *Unit3::normalize (vector3 (1.0, 0.0, 1.0)) * 2.0);
  }

  #[test]
  fn cylinder3() {
    use num::One;
    let cylinder : Cylinder3 <f32> = Cylinder3 {
      center:      Point3::origin(),
      radius:      Positive::noisy (2.0),
      half_height: Positive::one()
    };
    let support = cylinder.support (Unit3::axis_z().into()).0;
    assert_eq!(support, point3 (0.0, 0.0, 1.0));
    let support = cylinder.support (Unit3::axis_x().into()).0;
    assert_eq!(support, point3 (2.0, 0.0, 0.0));
    let support = cylinder.support (NonZero3::noisy (vector3 (1.0, 0.0, 1.0))).0;
    assert_eq!(support, point3 (2.0, 0.0, 1.0));
  }

  #[test]
  fn obb2() {
    let points : Vec <Point2 <f32>> = [
      [-1.0, -4.0],
      [ 2.0,  2.0],
      [-4.0, -1.0]
    ].map (Point2::from).into_iter().collect();
    let obb = Obb2::containing_points (&points).unwrap();
    let corners = obb.corners();
    assert_eq!(obb.center, [-0.25, -0.25].into());
    approx::assert_relative_eq!(point2 (-4.0, -1.0), corners[0], max_relative = 0.00001);
    approx::assert_relative_eq!(point2 ( 0.5,  3.5), corners[1], max_relative = 0.00001);
    approx::assert_relative_eq!(point2 (-1.0, -4.0), corners[2], max_relative = 0.00001);
    approx::assert_relative_eq!(point2 ( 3.5,  0.5), corners[3], max_relative = 0.00001);
    // test support
    let points : Vec <Point2 <f32>> = [
      [-2.0,  0.0],
      [ 0.0,  2.0],
      [ 2.0,  0.0],
      [ 0.0, -2.0],
    ].map (Point2::from).into_iter().collect();
    let obb = Obb2::containing_points (&points).unwrap();
    let support = obb.support (Unit2::axis_y().into()).0;
    approx::assert_relative_eq!(support, point2 (0.0, 2.0), epsilon = 0.000_001);
  }

  #[test]
  #[expect(clippy::suboptimal_flops)]
  fn obb3() {
    use std::f32::consts::{FRAC_PI_2, FRAC_PI_4, PI};
    use approx::AbsDiffEq;

    let points = [
      [ 1.0, -1.0, -1.0],
      [ 1.0, -1.0,  1.0],
      [ 1.0,  1.0, -1.0],
      [ 1.0,  1.0,  1.0],
      [-1.0, -1.0, -1.0],
      [-1.0, -1.0,  1.0],
      [-1.0,  1.0, -1.0],
      [-1.0,  1.0,  1.0]
    ].map (Point3::<f32>::from);
    let obb1 = Obb3::containing_points_with_orientation (&points, Angles3::zero())
      .unwrap();
    assert_eq!(
      SortedSet::from_unsorted (obb1.corners().to_vec()),
      SortedSet::from_unsorted (
        Aabb3::containing (&points).unwrap().corners().to_vec()));
    let obb2 = Obb3::containing_points_approx (&points).unwrap();
    assert_eq!(obb1, obb2);

    let rotation = Rotation3::from_angle_y (FRAC_PI_4.into());
    let points = points.map (|p| rotation.rotate (p));
    let obb1 = Obb3::containing_points_with_orientation (&points,
      // because the cube is symmetrical, it gets oriented by one of the faces
      Angles3::wrap (
        (3.0 * FRAC_PI_2).into(),
        (2.0 * PI - FRAC_PI_4).into(),
        FRAC_PI_2.into())
    ).unwrap();
    let obb2 = Obb3::containing_points_approx (&points).unwrap();
    approx::assert_abs_diff_eq!(obb1.center, obb2.center);
    assert_eq!(obb1.orientation, obb2.orientation);
    approx::assert_abs_diff_eq!(*obb1.extents.x, *obb2.extents.x,
      epsilon=4.0 * f32::default_epsilon());
    approx::assert_abs_diff_eq!(*obb1.extents.y, *obb2.extents.y,
      epsilon=4.0 * f32::default_epsilon());
    approx::assert_abs_diff_eq!(*obb1.extents.z, *obb2.extents.z,
      epsilon=2.0 * f32::default_epsilon());
    // because of rounding errors we can't sort the points to compare
    let corners = obb2.corners();
    approx::assert_abs_diff_eq!(corners[0], points[7],
      epsilon=8.0 * f32::default_epsilon());
    approx::assert_abs_diff_eq!(corners[1], points[5],
      epsilon=8.0 * f32::default_epsilon());
    approx::assert_abs_diff_eq!(corners[2], points[3],
      epsilon=8.0 * f32::default_epsilon());
    approx::assert_abs_diff_eq!(corners[3], points[1],
      epsilon=8.0 * f32::default_epsilon());
    approx::assert_abs_diff_eq!(corners[4], points[6],
      epsilon=8.0 * f32::default_epsilon());
    approx::assert_abs_diff_eq!(corners[5], points[4],
      epsilon=8.0 * f32::default_epsilon());
    approx::assert_abs_diff_eq!(corners[6], points[2],
      epsilon=8.0 * f32::default_epsilon());
    approx::assert_abs_diff_eq!(corners[7], points[0],
      epsilon=8.0 * f32::default_epsilon());

    // test support
    let points : Vec <Point3 <f32>> = [
      [-2.0,  0.0, -2.0],
      [ 0.0,  2.0, -2.0],
      [ 2.0,  0.0, -2.0],
      [ 0.0, -2.0, -2.0],
      [-2.0,  0.0,  2.0],
      [ 0.0,  2.0,  2.0],
      [ 2.0,  0.0,  2.0],
      [ 0.0, -2.0,  2.0]
    ].map (Point3::from).into_iter().collect();
    let obb = Obb3::containing_points_approx (&points).unwrap();
    let support = obb.support (Unit3::axis_y().into()).0;
    approx::assert_relative_eq!(support, point3 (0.0, 2.0, 0.0), epsilon = 0.000_001);
    let support = obb.support (NonZero3::noisy (vector3 (0.0, 1.0, 1.0))).0;
    approx::assert_relative_eq!(support, point3 (0.0, 2.0, 2.0), epsilon = 0.000_001);
  }

  #[test]
  fn support_aabb3() {
    let aabb = Aabb3::with_minmax (
      [-1.0, -2.0, -3.0].into(),
      [ 1.0,  2.0,  3.0].into()).unwrap();
    assert_eq!(
      [-1.0, -2.0, -3.0],
      aabb.support (NonZero3::noisy ([-1.0, -1.0, -1.0].into())).0.0.into_array());
    assert_eq!(
      [-1.0, -2.0,  3.0],
      aabb.support (NonZero3::noisy ([-1.0, -1.0,  1.0].into())).0.0.into_array());
    assert_eq!(
      [-1.0,  2.0, -3.0],
      aabb.support (NonZero3::noisy ([-1.0,  1.0, -1.0].into())).0.0.into_array());
    assert_eq!(
      [-1.0,  2.0,  3.0],
      aabb.support (NonZero3::noisy ([-1.0,  1.0,  1.0].into())).0.0.into_array());
    assert_eq!(
      [ 1.0, -2.0, -3.0],
      aabb.support (NonZero3::noisy ([ 1.0, -1.0, -1.0].into())).0.0.into_array());
    assert_eq!(
      [ 1.0, -2.0,  3.0],
      aabb.support (NonZero3::noisy ([ 1.0, -1.0,  1.0].into())).0.0.into_array());
    assert_eq!(
      [ 1.0,  2.0, -3.0],
      aabb.support (NonZero3::noisy ([ 1.0,  1.0, -1.0].into())).0.0.into_array());
    assert_eq!(
      [ 1.0,  2.0,  3.0],
      aabb.support (NonZero3::noisy ([ 1.0,  1.0,  1.0].into())).0.0.into_array());

    assert_eq!(
      [-1.0, -2.0,  0.0],
      aabb.support (NonZero3::noisy ([-1.0, -1.0,  0.0].into())).0.0.into_array());
    assert_eq!(
      [-1.0,  2.0,  0.0],
      aabb.support (NonZero3::noisy ([-1.0,  1.0,  0.0].into())).0.0.into_array());
    assert_eq!(
      [ 1.0, -2.0,  0.0],
      aabb.support (NonZero3::noisy ([ 1.0, -1.0,  0.0].into())).0.0.into_array());
    assert_eq!(
      [ 1.0,  2.0,  0.0],
      aabb.support (NonZero3::noisy ([ 1.0,  1.0,  0.0].into())).0.0.into_array());

    assert_eq!(
      [-1.0,  0.0, -3.0],
      aabb.support (NonZero3::noisy ([-1.0,  0.0, -1.0].into())).0.0.into_array());
    assert_eq!(
      [-1.0,  0.0,  3.0],
      aabb.support (NonZero3::noisy ([-1.0,  0.0,  1.0].into())).0.0.into_array());
    assert_eq!(
      [ 1.0,  0.0, -3.0],
      aabb.support (NonZero3::noisy ([ 1.0,  0.0, -1.0].into())).0.0.into_array());
    assert_eq!(
      [ 1.0,  0.0,  3.0],
      aabb.support (NonZero3::noisy ([ 1.0,  0.0,  1.0].into())).0.0.into_array());

    assert_eq!(
      [ 0.0, -2.0, -3.0],
      aabb.support (NonZero3::noisy ([ 0.0, -1.0, -1.0].into())).0.0.into_array());
    assert_eq!(
      [ 0.0, -2.0,  3.0],
      aabb.support (NonZero3::noisy ([ 0.0, -1.0,  1.0].into())).0.0.into_array());
    assert_eq!(
      [ 0.0,  2.0, -3.0],
      aabb.support (NonZero3::noisy ([ 0.0,  1.0, -1.0].into())).0.0.into_array());
    assert_eq!(
      [ 0.0,  2.0,  3.0],
      aabb.support (NonZero3::noisy ([ 0.0,  1.0,  1.0].into())).0.0.into_array());

    assert_eq!(
      [-1.0,  0.0,  0.0],
      aabb.support (NonZero3::noisy ([-1.0,  0.0,  0.0].into())).0.0.into_array());
    assert_eq!(
      [ 1.0,  0.0,  0.0],
      aabb.support (NonZero3::noisy ([ 1.0,  0.0,  0.0].into())).0.0.into_array());
    assert_eq!(
      [ 0.0, -2.0,  0.0],
      aabb.support (NonZero3::noisy ([ 0.0, -1.0,  0.0].into())).0.0.into_array());
    assert_eq!(
      [ 0.0,  2.0,  0.0],
      aabb.support (NonZero3::noisy ([ 0.0,  1.0,  0.0].into())).0.0.into_array());
    assert_eq!(
      [ 0.0,  0.0, -3.0],
      aabb.support (NonZero3::noisy ([ 0.0,  0.0, -1.0].into())).0.0.into_array());
    assert_eq!(
      [ 0.0,  0.0,  3.0],
      aabb.support (NonZero3::noisy ([ 0.0,  0.0,  1.0].into())).0.0.into_array());
  }
} // end tests