nsys-math-utils 0.1.0

//! Affine, convex, and combinatorial spaces

use approx;
use num_traits as num;
use rand;

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

#[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;
pub use simplex3::Segment     as Segment3;
pub use simplex2::Triangle    as Triangle2;
pub use simplex3::Triangle    as Triangle3;
pub use simplex3::Tetrahedron as Tetrahedron3;

pub trait Primitive <S, V> where
  S : Field,
  V : VectorSpace <S>
{
  fn translate (self, displacement : V) -> Self;
  fn scale     (self, scale : NonZero <S>) -> Self;
}

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

/// 1D interval
#[cfg_attr(feature = "derive_serdes", derive(Serialize, Deserialize))]
#[derive(Clone, Copy, Debug, 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()
/// );
/// ```
/// 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()
/// );
/// ```
#[cfg_attr(feature = "derive_serdes", derive(Serialize, Deserialize))]
#[derive(Clone, Copy, Debug, 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()
/// );
/// let y = Aabb3::with_minmax (
///   [0.0, 0.0, 0.0].into(),
///   [0.0, 1.0, 0.0].into()
/// );
/// ```
/// 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()
/// );
/// ```
#[cfg_attr(feature = "derive_serdes", derive(Serialize, Deserialize))]
#[derive(Clone, Copy, Debug, 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, PartialEq)]
pub struct Obb2 <S> {
  center        : Point2 <S>,
  half_extent_x : Positive <S>,
  half_extent_y : Positive <S>,
  orientation   : AngleWrapped <S>
}

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

/// 3D Z-axis-aligned cylinder
#[cfg_attr(feature = "derive_serdes", derive(Serialize, Deserialize))]
#[derive(Clone, Copy, Debug, 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, PartialEq)]
pub struct Capsule3 <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, PartialEq)]
pub struct Line2 <S> {
  pub base      : Point2 <S>,
  pub direction : Unit2  <S>
}

/// An infinitely extended line in 3D space defined by a base point and
/// normalized direction
#[cfg_attr(feature = "derive_serdes", derive(Serialize, Deserialize))]
#[derive(Clone, Copy, Debug, PartialEq)]
pub struct Line3 <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, 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, 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, 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 : Field + approx::RelativeEq
{
  // early exit: a pair of points are equal
  if a == b || a == c || b == c {
    true
  } else {
    let a = Vector3::from_point_2d (a.0).into();
    let b = Vector3::from_point_2d (b.0).into();
    let c = Vector3::from_point_2d (c.0).into();
    let determinant = determinant_3d (&a, &b, &c);
    // a zero determinant indicates that the points are colinear
    approx::relative_eq!(S::zero(), determinant)
  }
}

/// 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 : Field + Sqrt + approx::RelativeEq
{
  // early exit: a pair of points are equal
  if a == b || a == c || b == c {
    true
  } else {
    let determinant = determinant_3d (a, b, c);
    // a zero determinant indicates that the points are coplanar
    if approx::relative_eq!(S::zero(), determinant) {
      approx::relative_eq!(S::zero(), triangle_3d_area2 (a, b, c))
    } else {
      false
    }
  }
}

/// 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 : Ring + approx::RelativeEq
{
  // early exit: a pair of points are equal
  if a == b || a == c || a == d || b == c || b == d || c == d {
    true
  } else {
    let ab_cross_ac_dot_ad = (*b-a).cross (*c-a).dot (*d-a);
    approx::relative_eq!(S::zero(), ab_cross_ac_dot_ad)
  }
}

/// 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_area2;
/// assert_relative_eq!(
///   3.0/4.0,
///   triangle_3d_area2 (
///     &[-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_area2;
/// assert_eq!(
///   0.0,
///   triangle_3d_area2 (
///     &[-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_area2 <S> (a : &Point3 <S>, b : &Point3 <S>, c : &Point3 <S>)
  -> S
where
  S : Field + Sqrt
{
  if a == b || a == c || b == c {
    return S::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 max >= mid >= min
  let max = S::max (S::max (ab_mag, ac_mag), bc_mag);
  let min = S::min (S::min (ab_mag, ac_mag), bc_mag);
  let mid = if min <= ab_mag && ab_mag <= max {
    ab_mag
  } else if min <= ac_mag && ac_mag <= max {
    ac_mag
  } else {
    bc_mag
  };
  let two = S::one() + S::one();
  let sixteen = two * two * two * two;
  // 1.0/16.0
  let frac = S::one() / sixteen;
  frac
    * (max + (mid + min))
    * (min - (max - mid))
    * (min + (max - mid))
    * (max + (mid - min))
}

/// 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()
}

/// Coordinate-wise min
pub fn point2_min <S : Ring> (a : &Point2 <S>, b : &Point2 <S>) -> Point2 <S> {
  Vector2::partial_min (a.0, b.0).into()
}

/// Coordinate-wise max
pub fn point2_max <S : Ring> (a : &Point2 <S>, b : &Point2 <S>) -> Point2 <S> {
  Vector2::partial_max (a.0, b.0).into()
}

/// Coordinate-wise min
pub fn point3_min <S : Ring> (a : &Point3 <S>, b : &Point3 <S>) -> Point3 <S> {
  Vector3::partial_min (a.0, b.0).into()
}

/// Coordinate-wise max
pub fn point3_max <S : Ring> (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
pub fn project_2d_point_on_line <S : Real>
  (point : &Point2 <S>, line : &Line2 <S>) -> Point2 <S>
{
  let dot_dir = line.direction.dot (*point - line.base);
  line.point (dot_dir)
}

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

// 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.norm_squared();
  let origin           = points[i1 as usize];
  let mut support      = [Point2::<S>::origin(); 4];
  let mut rect_index   = [i1, i1, i1, i1];
  for (i, point) in points.iter().enumerate() {
    let diff = point - origin;
    let v    = point2 (edge.dot (diff), perp.dot (diff));
    if v.0.x > support[1].0.x ||
      v.0.x == support[1].0.x && v.0.y > support[1].0.y
    {
      rect_index[1] = i as u32;
      support[1]    = v;
    }
    if v.0.y > support[2].0.y ||
      v.0.y == support[2].0.y && v.0.x < support[2].0.x
    {
      rect_index[2] = i as u32;
      support[2]    = v;
    }
    if v.0.x < support[3].0.x ||
      v.0.x == support[3].0.x && v.0.y < support[3].0.y
    {
      rect_index[3] = i as u32;
      support[3]    = v;
    }
  }
  let scaled_width  = support[1].0.x - support[3].0.x;
  let scaled_height = support[2].0.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 : Ring> Interval <S> {
  /// 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 (1.0, 0.0);  // panic!
  /// ```
  ///
  /// Debug panic if points are identical:
  ///
  /// ```should_panic
  /// # use math_utils::geometry::*;
  /// let b = Interval::with_minmax (0.0, 0.0);  // panic!
  /// ```
  #[inline]
  pub fn with_minmax (min : S, max : S) -> Self where S : std::fmt::Debug {
    debug_assert_ne!(min, max);
    debug_assert_eq!(min, S::min (min, max));
    debug_assert_eq!(max, S::max (min, max));
    Interval { 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 = Interval::from_points (0.0, 0.0);  // panic!
  /// ```
  #[inline]
  pub fn from_points (a : S, b : S) -> Self where S : std::fmt::Debug {
    debug_assert_ne!(a, b);
    let min = S::min (a, b);
    let max = S::max (a, b);
    Interval { min, max }
  }
  /// Construct the minimum AABB containing the given set of points.
  ///
  /// Debug panic if fewer than 2 points are given:
  ///
  /// ```should_panic
  /// # use math_utils::geometry::*;
  /// let b = Interval::<f32>::containing (&[]);  // panic!
  /// ```
  ///
  /// ```should_panic
  /// # use math_utils::geometry::*;
  /// let b = Interval::containing (&[0.0]);  // panic!
  /// ```
  #[inline]
  pub fn containing (points : &[S]) -> Self where S : std::fmt::Debug {
    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 where
    S : std::fmt::Debug
  {
    Interval::with_minmax (
      S::min (a.min(), b.min()),
      S::max (a.max(), b.max()))
  }
  #[inline]
  pub fn min (&self) -> S {
    self.min
  }
  #[inline]
  pub fn max (&self) -> S {
    self.max
  }
  #[inline]
  pub fn length (&self) -> NonNegative <S> {
    self.width()
  }
  #[inline]
  pub fn width (&self) -> NonNegative <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::from_points (-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);
  /// ```
  pub fn clamp (&self, point : &S) -> S {
    S::max (S::min (self.max, *point), self.min)
  }
  /// 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 (-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::Rng
  {
    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>>
    where S : std::fmt::Debug
  {
    intersect::continuous_interval (&self, &other)
  }

  /// Increase or decrease each endpoint by the given amount.
  ///
  /// Debug panic if the interval width is less than or equal to twice the given
  /// amount:
  /// ```should_panic
  /// # use math_utils::geometry::*;
  /// let x = Interval::<f32>::with_minmax (-1.0, 1.0);
  /// let y = x.dilate (-1.0); // panic!
  /// ```
  pub fn dilate (mut self, amount : S) -> Self {
    self.min -= amount;
    self.max += amount;
    debug_assert!(self.min < self.max);
    self
  }
}

impl <S> Primitive <S, S> for Interval <S> where
  S : Field + VectorSpace <S>
{
  fn translate (mut self, displacement : S) -> Self {
    self.min += displacement;
    self.max += displacement;
    self
  }
  fn scale (mut self, scale : NonZero <S>) -> Self {
    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;
    self
  }
}

impl <S> std::ops::Add <S> for Interval <S> where
  S    : Field + VectorSpace <S>,
  Self : Primitive <S, S>
{
  type Output = Self;
  fn add (self, displacement : S) -> Self {
    self.translate (displacement)
  }
}

impl <S> std::ops::Sub <S> for Interval <S> where
  S    : Field + VectorSpace <S>,
  Self : Primitive <S, S>
{
  type Output = Self;
  fn sub (self, displacement : S) -> Self {
    self.translate (-displacement)
  }
}

impl_numcast_fields!(Aabb2, min, max);
impl <S : Ring> Aabb2 <S> {
  /// 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 ([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 ([0.0, 0.0].into(), [0.0, 0.0].into());  // panic!
  /// ```
  #[inline]
  pub fn with_minmax (min : Point2 <S>, max : Point2 <S>) -> Self where
    S : std::fmt::Debug
  {
    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 ([0.0, 0.0].into(), [0.0, 0.0].into());  // panic!
  /// ```
  #[inline]
  pub fn from_points (a : Point2 <S>, b : Point2 <S>) -> Self where
    S : std::fmt::Debug
  {
    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.
  ///
  /// Debug panic if fewer than 2 points are given:
  ///
  /// ```should_panic
  /// # use math_utils::geometry::*;
  /// let b = Aabb2::<f32>::containing (&[]);  // panic!
  /// ```
  ///
  /// ```should_panic
  /// # use math_utils::geometry::*;
  /// let b = Aabb2::containing (&[[0.0, 0.0].into()]);  // panic!
  /// ```
  #[inline]
  pub fn containing (points : &[Point2 <S>]) -> Self where S : std::fmt::Debug {
    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 where
    S : std::fmt::Debug
  {
    Aabb2::with_minmax (
      point2_min (a.min(), b.min()),
      point2_max (a.max(), b.max())
    )
  }
  #[inline]
  pub fn min (&self) -> &Point2 <S> {
    &self.min
  }
  #[inline]
  pub fn max (&self) -> &Point2 <S> {
    &self.max
  }
  #[inline]
  pub fn center (&self) -> Point2 <S> where S : Field {
    Point2 ((self.min.0 + self.max.0) / S::two())
  }
  #[inline]
  pub fn dimensions (&self) -> Vector2 <S> {
    self.max.0 - self.min.0
  }
  #[inline]
  pub fn extents (&self) -> Vector2 <S> where S : Field {
    self.dimensions() * S::half()
  }
  /// X dimension
  #[inline]
  pub fn width (&self) -> NonNegative <S> {
    NonNegative::unchecked (self.max.0.x - self.min.0.x)
  }
  /// Y dimension
  #[inline]
  pub fn height (&self) -> NonNegative <S> {
    NonNegative::unchecked (self.max.0.y - self.min.0.y)
  }
  #[inline]
  pub fn area (&self) -> NonNegative <S> {
    self.width() * self.height()
  }
  #[inline]
  pub fn contains (&self, point : &Point2 <S>) -> bool {
    self.min.0.x < point.0.x && point.0.x < self.max.0.x &&
    self.min.0.y < point.0.y && point.0.y < self.max.0.y
  }
  /// Clamp a given point to the AABB.
  ///
  /// ```
  /// # use math_utils::geometry::*;
  /// let b = Aabb2::from_points ([-1.0, -1.0].into(), [1.0, 1.0].into());
  /// 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> {
    [ S::max (S::min (self.max.0.x, point.0.x), self.min.0.x),
      S::max (S::min (self.max.0.y, point.0.y), self.min.0.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());
  /// 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::Rng
  {
    let x_range = self.min.0.x..self.max.0.x;
    let y_range = self.min.0.y..self.max.0.y;
    let x = if x_range.is_empty() {
      self.min.0.x
    } else {
      rng.random_range (x_range)
    };
    let y = if y_range.is_empty() {
      self.min.0.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>> where
    S : std::fmt::Debug
  {
    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.0.x, self.min.0.y].into(),
      // upper-left
      Quadrant::NegPos => [self.min.0.x, self.max.0.y].into(),
      // lower-left
      Quadrant::NegNeg => self.min
    }
  }

  pub fn edge (&self, direction : SignedAxis2) -> Self where
    S : std::fmt::Debug
  {
    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 (min, max)
  }

  pub fn extrude (&self, axis : SignedAxis2, amount : Positive <S>) -> Self {
    let (min, max) = 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())
      )
    };
    Aabb2 { min, max }
  }

  pub fn extend (&mut self, axis : SignedAxis2, amount : Positive <S>) {
    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
    }
  }

  pub fn with_z (self, z : Interval <S>) -> Aabb3 <S> where S : std::fmt::Debug {
    Aabb3::with_minmax (
      self.min.0.with_z (z.min()).into(),
      self.max.0.with_z (z.max()).into()
    )
  }

  /// Increase or decrease each dimension by the given amount.
  ///
  /// Debug panic if any dimension would become negative:
  /// ```should_panic
  /// # use math_utils::geometry::*;
  /// let x = Aabb2::with_minmax ([0.0, 0.0].into(), [1.0, 2.0].into());
  /// let y = x.dilate (-1.0); // panic!
  /// ```
  pub fn dilate (mut self, amount : S) -> Self where S : std::fmt::Debug {
    self.min -= Vector2::broadcast (amount);
    self.max += Vector2::broadcast (amount);
    debug_assert_ne!(self.min, self.max);
    debug_assert_eq!(self.min, point2_min (&self.min, &self.max));
    debug_assert_eq!(self.max, point2_max (&self.min, &self.max));
    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 : std::fmt::Debug {
    let (min, max) = match axis {
      Axis2::X => (self.min.0.y, self.max.0.y),
      Axis2::Y => (self.min.0.x, self.max.0.x)
    };
    Interval::with_minmax (min, max)
  }
}

impl <S : Field> Primitive <S, Vector2 <S>> for Aabb2 <S> {
  fn translate (mut self, displacement : Vector2 <S>) -> Self {
    self.min.0 += displacement;
    self.max.0 += displacement;
    self
  }
  fn scale (mut self, scale : NonZero <S>) -> Self {
    let center = self.center();
    self = self.translate (-center.0);
    self.min.0 *= *scale;
    self.max.0 *= *scale;
    self.translate (center.0)
  }
}

impl <S> std::ops::Add <Vector2 <S>> for Aabb2 <S> where
  S    : Field,
  Self : Primitive <S, Vector2 <S>>
{
  type Output = Self;
  fn add (self, displacement : Vector2 <S>) -> Self {
    self.translate (displacement)
  }
}

impl <S> std::ops::Sub <Vector2 <S>> for Aabb2 <S> where
  S    : Field,
  Self : Primitive <S, Vector2 <S>>
{
  type Output = Self;
  fn sub (self, displacement : Vector2 <S>) -> Self {
    self.translate (-displacement)
  }
}

impl_numcast_fields!(Aabb3, min, max);
impl <S : Ring> Aabb3 <S> {
  /// 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 ([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 ([0.0, 0.0, 0.0].into(), [0.0, 0.0, 0.0].into());  // panic!
  /// ```
  #[inline]
  pub fn with_minmax (min : Point3 <S>, max : Point3 <S>) -> Self where
    S : std::fmt::Debug
  {
    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 (
  ///   [0.0, 0.0, 0.0].into(),
  ///   [0.0, 0.0, 0.0].into());  // panic!
  /// ```
  #[inline]
  pub fn from_points (a : Point3 <S>, b : Point3 <S>) -> Self where
    S : std::fmt::Debug
  {
    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.
  ///
  /// Debug panic if fewer than 2 points are given:
  ///
  /// ```should_panic
  /// # use math_utils::geometry::*;
  /// let b = Aabb3::<f32>::containing (&[]);  // panic!
  /// ```
  ///
  /// ```should_panic
  /// # use math_utils::geometry::*;
  /// let b = Aabb3::containing (&[[0.0, 0.0, 0.0].into()]);  // panic!
  /// ```
  #[inline]
  pub fn containing (points : &[Point3 <S>]) -> Self where S : std::fmt::Debug {
    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 where
    S : std::fmt::Debug
  {
    Aabb3::with_minmax (
      point3_min (a.min(), b.min()),
      point3_max (a.max(), b.max())
    )
  }
  #[inline]
  pub fn min (&self) -> &Point3 <S> {
    &self.min
  }
  #[inline]
  pub fn max (&self) -> &Point3 <S> {
    &self.max
  }
  #[inline]
  pub fn center (&self) -> Point3 <S> where S : Field {
    Point3 ((self.min.0 + self.max.0) / S::two())
  }
  #[inline]
  pub fn dimensions (&self) -> Vector3 <S> {
    self.max.0 - self.min.0
  }
  #[inline]
  pub fn extents (&self) -> Vector3 <S> where S : Field {
    self.dimensions() * S::half()
  }
  /// X dimension
  #[inline]
  pub fn width (&self) -> NonNegative <S> {
    NonNegative::unchecked (self.max.0.x - self.min.0.x)
  }
  /// Y dimension
  #[inline]
  pub fn depth (&self) -> NonNegative <S> {
    NonNegative::unchecked (self.max.0.y - self.min.0.y)
  }
  /// Z dimension
  #[inline]
  pub fn height (&self) -> NonNegative <S> {
    NonNegative::unchecked (self.max.0.z - self.min.0.z)
  }
  #[inline]
  pub fn volume (&self) -> NonNegative <S> {
    self.width() * self.depth() * self.height()
  }
  #[inline]
  pub fn contains (&self, point : &Point3 <S>) -> bool {
    self.min.0.x < point.0.x && point.0.x < self.max.0.x &&
    self.min.0.y < point.0.y && point.0.y < self.max.0.y &&
    self.min.0.z < point.0.z && point.0.z < self.max.0.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());
  /// 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> {
    [ S::max (S::min (self.max.0.x, point.0.x), self.min.0.x),
      S::max (S::min (self.max.0.y, point.0.y), self.min.0.y),
      S::max (S::min (self.max.0.z, point.0.z), self.min.0.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());
  /// 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::Rng
  {
    let x_range = self.min.0.x..self.max.0.x;
    let y_range = self.min.0.y..self.max.0.y;
    let z_range = self.min.0.z..self.max.0.z;
    let x = if x_range.is_empty() {
      self.min.0.x
    } else {
      rng.random_range (x_range)
    };
    let y = if y_range.is_empty() {
      self.min.0.y
    } else {
      rng.random_range (y_range)
    };
    let z = if z_range.is_empty() {
      self.min.0.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>> where
    S : std::fmt::Debug
  {
    intersect::continuous_aabb3_aabb3 (&self, &other)
  }

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

  pub fn face (&self, direction : SignedAxis3) -> Self where
    S : std::fmt::Debug
  {
    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 (min, max)
  }

  pub fn extrude (&self, axis : SignedAxis3, amount : Positive <S>) -> Self {
    let (min, max) = 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 { min, max }
  }

  pub fn extend (&mut self, axis : SignedAxis3, amount : Positive <S>) {
    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
    }
  }

  /// Increase or decrease each dimension by the given amount.
  ///
  /// Debug panic if any dimension would become negative:
  /// ```should_panic
  /// # use math_utils::geometry::*;
  /// let x = Aabb3::with_minmax ([0.0, 0.0, 0.0].into(), [1.0, 2.0, 3.0].into());
  /// let y = x.dilate (-1.0); // panic!
  /// ```
  pub fn dilate (mut self, amount : S) -> Self where S : std::fmt::Debug {
    self.min -= Vector3::broadcast (amount);
    self.max += Vector3::broadcast (amount);
    debug_assert_ne!(self.min, self.max);
    debug_assert_eq!(self.min, point3_min (&self.min, &self.max));
    debug_assert_eq!(self.max, point3_max (&self.min, &self.max));
    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> where S : std::fmt::Debug {
    let (min, max) = match axis {
      Axis3::X => ([self.min.0.y, self.min.0.z], [self.max.0.y, self.max.0.z]),
      Axis3::Y => ([self.min.0.x, self.min.0.z], [self.max.0.x, self.max.0.z]),
      Axis3::Z => ([self.min.0.x, self.min.0.y], [self.max.0.x, self.max.0.y]),
    };
    Aabb2::with_minmax (min.into(), max.into())
  }
}

impl <S : Field> Primitive <S, Vector3 <S>> for Aabb3 <S> {
  fn translate (mut self, displacement : Vector3 <S>) -> Self {
    self.min.0 += displacement;
    self.max.0 += displacement;
    self
  }
  fn scale (mut self, scale : NonZero <S>) -> Self {
    let center = self.center();
    self = self.translate (-center.0);
    self.min.0 *= *scale;
    self.max.0 *= *scale;
    self.translate (center.0)
  }
}

impl <S> std::ops::Add <Vector3 <S>> for Aabb3 <S> where
  S    : Field,
  Self : Primitive <S, Vector3 <S>>
{
  type Output = Self;
  fn add (self, displacement : Vector3 <S>) -> Self {
    self.translate (displacement)
  }
}

impl <S> std::ops::Sub <Vector3 <S>> for Aabb3 <S> where
  S    : Field,
  Self : Primitive <S, Vector3 <S>>
{
  type Output = Self;
  fn sub (self, displacement : Vector3 <S>) -> Self {
    self.translate (-displacement)
  }
}

impl_numcast_fields!(Obb2, center, half_extent_x, half_extent_y, orientation);
impl <S : Ring> Obb2 <S> {
  /// Construct a new OBB
  #[inline]
  pub fn new (
    center        : Point2 <S>,
    half_extent_x : Positive <S>,
    half_extent_y : Positive <S>,
    orientation   : AngleWrapped <S>
  ) -> Self {
    Obb2 { center, half_extent_x, half_extent_y, 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
  {
    if points.len() < 3 {
      return None
    }
    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 if points span an empty set:
  ///
  /// ```
  /// # use math_utils::*;
  /// # use math_utils::geometry::*;
  /// 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()];
    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 (&mut 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 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 half_extent_x = Positive::unchecked (S::sqrt (
      (S::half() * min_rect.edge.dot (diff[0])).squared()
      / min_rect.edge_len_squared));
    let half_extent_y = Positive::unchecked (S::sqrt (
      (S::half() * min_rect.perp.dot (diff[1])).squared()
      / min_rect.edge_len_squared));
    let orientation =
      AngleWrapped::wrap (Rad (min_rect.edge.y.atan2 (min_rect.edge.x)));
    let obb = Obb2 { center, half_extent_x, half_extent_y, 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;
          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.norm_squared();
          let sin_theta_squared = (dp * dp) / e_lensq;
          angles[*nangles as usize] = (sin_theta_squared, k0 as u32);
          *nangles += 1;
        }
        k0 = k1;
        k1 += 1;
      }
      *nangles > 0
    }
    fn sort_angles <S : PartialOrd> (angles : &mut [(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;
          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.norm_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 center (&self) -> Point2 <S> {
    self.center
  }
  #[inline]
  pub fn half_extent_x (&self) -> Positive <S> {
    self.half_extent_x
  }
  #[inline]
  pub fn half_extent_y (&self) -> Positive <S> {
    self.half_extent_y
  }
  #[inline]
  pub fn orientation (&self) -> AngleWrapped <S> {
    self.orientation
  }
  pub fn dimensions (&self) -> Vector2 <S> where S : Field {
    self.extents() * S::two()
  }
  pub fn extents (&self) -> Vector2 <S> {
    [*self.half_extent_x, *self.half_extent_y].into()
  }
  /// X dimension
  #[inline]
  pub fn width (&self) -> Positive <S> where S : Field {
    self.half_extent_x * Positive::unchecked (S::two())
  }
  /// Y dimension
  #[inline]
  pub fn height (&self) -> Positive <S> where S : Field {
    self.half_extent_y * Positive::unchecked (S::two())
  }
  #[inline]
  pub fn area (&self) -> Positive <S> where S : Field {
    self.width() * self.height()
  }
  #[inline]
  pub fn corners (&self) -> [Point2 <S>; 4] where
    S : Real + num::real::Real + MaybeSerDes
  {
    let extents = self.extents();
    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] = [extents.x * sign_x, extents.y * sign_y].into();
      }
    }
    let rotation = Rotation2::from_angle (self.orientation.angle());
    points.map (|point| rotation.rotate (point) + self.center.0)
  }

  /// 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(),
  ///   Positive::unchecked (2.0),
  ///   Positive::unchecked (0.75),
  ///   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.0.x.abs() < *self.half_extent_x && p.0.y.abs() < *self.half_extent_y
  }

  /// Increase or decrease each dimension by the given amount.
  ///
  /// Panics if any dimension would become zero or negative:
  ///
  /// ```should_panic
  /// # use math_utils::*;
  /// # use math_utils::geometry::*;
  /// # use math_utils::num_traits::Zero;
  /// let x = Obb2::new (
  ///   Point2::origin(),
  ///   Positive::unchecked (0.5),
  ///   Positive::unchecked (0.5),
  ///   AngleWrapped::zero()
  /// );
  /// let y = x.dilate (-1.0); // panic!
  /// ```
  pub fn dilate (mut self, amount : S) -> Self where S : Field {
    self.half_extent_x = Positive::noisy (*self.half_extent_x + amount);
    self.half_extent_y = Positive::noisy (*self.half_extent_y + 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_traits::Zero;
  /// let x = Aabb2::with_minmax (
  ///   [0.0, 0.0].into(),
  ///   [0.0, 1.0].into()
  /// );
  /// let y = Obb2::from (x); // panic!
  /// ```
  fn from (aabb : Aabb2 <S>) -> Self {
    use num::Zero;
    let center        = aabb.center();
    let half_extent_x = Positive::noisy (*aabb.width()  / S::two());
    let half_extent_y = Positive::noisy (*aabb.height() / S::two());
    let orientation   = AngleWrapped::zero();
    Obb2::new (center, half_extent_x, half_extent_y, orientation)
  }
}

impl_numcast_fields!(Obb3, center, half_extent_x, half_extent_y, half_extent_z,
  orientation);
impl <S : Ring> Obb3 <S> {
  /// Construct a new OBB
  #[inline]
  pub fn new (
    center        : Point3 <S>,
    half_extent_x : Positive <S>,
    half_extent_y : Positive <S>,
    half_extent_z : Positive <S>,
    orientation   : Angles3 <S>
  ) -> Self {
    Obb3 { center, half_extent_x, half_extent_y, half_extent_z, orientation }
  }
  /// Construct the minimum OBB containing the convex hull of a given set of
  /// points
  pub fn containing_points (_points : &[Point3 <S>]) -> Option <Self> where
    S : Real
  {
    unimplemented!("TODO")
  }
  /// Construct the minimum OBB containing the given convex hull of points
  pub fn containing (_hull : &Hull3 <S>) -> Option <Self> where S : Real {
    unimplemented!("TODO")
  }
  #[inline]
  pub fn center (&self) -> Point3 <S> {
    self.center
  }
  #[inline]
  pub fn half_extent_x (&self) -> Positive <S> {
    self.half_extent_x
  }
  #[inline]
  pub fn half_extent_y (&self) -> Positive <S> {
    self.half_extent_y
  }
  #[inline]
  pub fn half_extent_z (&self) -> Positive <S> {
    self.half_extent_z
  }
  #[inline]
  pub fn orientation (&self) -> Angles3 <S> {
    self.orientation
  }
  pub fn dimensions (&self) -> Vector3 <S> where S : Field {
    self.extents() * S::two()
  }
  pub fn extents (&self) -> Vector3 <S> where S : Field {
    [*self.half_extent_x, *self.half_extent_y, *self.half_extent_z].into()
  }
  /// X dimension
  #[inline]
  pub fn width (&self) -> Positive <S> where S : Field {
    self.half_extent_x * Positive::unchecked (S::two())
  }
  /// Y dimension
  #[inline]
  pub fn depth (&self) -> Positive <S> where S : Field {
    self.half_extent_y * Positive::unchecked (S::two())
  }
  /// Z dimension
  #[inline]
  pub fn height (&self) -> Positive <S> where S : Field {
    self.half_extent_z * Positive::unchecked (S::two())
  }
  #[inline]
  pub fn volume (&self) -> Positive <S> where S : Field {
    self.width() * self.depth() * self.height()
  }
  #[inline]
  pub fn corners (&self) -> [Point3 <S>; 8] where
    S : Real + num::real::Real + MaybeSerDes
  {
    let extents = self.extents();
    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] =
            [extents.x * sign_x, extents.y * sign_y, extents.z * sign_z].into();
        }
      }
    }
    let rotation = Rotation3::from (self.orientation);
    points.map (|point| rotation.rotate (point) + self.center.0)
  }

  /// Check if point is contained by the bounding box
  pub fn contains (&self, _point : &Point2 <S>) -> bool where
    S : Real + num::real::Real + MaybeSerDes
  {
    unimplemented!("TODO")
  }

  /// Increase or decrease each dimension by the given amount.
  ///
  /// Panics if any dimension would become zero or negative:
  ///
  /// ```should_panic
  /// # use math_utils::*;
  /// # use math_utils::geometry::*;
  /// let x = Obb3::new (
  ///   Point3::origin(),
  ///   Positive::unchecked (0.5),
  ///   Positive::unchecked (0.5),
  ///   Positive::unchecked (0.5),
  ///   Angles3::zero()
  /// );
  /// let y = x.dilate (-1.0); // panic!
  /// ```
  pub fn dilate (mut self, amount : S) -> Self where S : Field {
    self.half_extent_x = Positive::noisy (*self.half_extent_x + amount);
    self.half_extent_y = Positive::noisy (*self.half_extent_y + amount);
    self.half_extent_z = Positive::noisy (*self.half_extent_z + 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_traits::Zero;
  /// let x = Aabb3::with_minmax (
  ///   [0.0, 0.0, 0.0].into(),
  ///   [0.0, 1.0, 1.0].into()
  /// );
  /// let y = Obb3::from (x); // panic!
  /// ```
  fn from (aabb : Aabb3 <S>) -> Self {
    let center = aabb.center();
    let [half_extent_x, half_extent_y, half_extent_z] =
      aabb.extents().into_array().map (Positive::noisy);
    let orientation = Angles3::zero();
    Obb3::new (center, half_extent_x, half_extent_y, half_extent_z, orientation)
  }
}

impl_numcast_fields!(Capsule3, center, half_height, radius);
impl <S : Ring> Capsule3 <S> {
  pub fn aabb3 (&self) -> Aabb3 <S> where S : Real + std::fmt::Debug {
    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 (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)
  }
}

// TODO: impl Primitive for Capsule3

impl_numcast_fields!(Cylinder3, center, half_height, radius);
impl <S : Ring> Cylinder3 <S> {
  pub fn aabb3 (&self) -> Aabb3 <S> where S : Real + std::fmt::Debug {
    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 (min, max)
  }
}

// TODO: impl Primitive for Cylinder3

impl_numcast_fields!(Line2, base, direction);
impl <S : Real> Line2 <S> {
  /// Construct a new 2D line
  #[inline]
  pub 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]
  #[allow(clippy::type_complexity)]
  pub fn intersect_aabb (&self, aabb : &Aabb2 <S>)
    -> Option <((S, Point2 <S>), (S, Point2 <S>))>
  where S : std::fmt::Debug {
    intersect::continuous_line2_aabb2 (self, aabb)
  }
  #[inline]
  #[allow(clippy::type_complexity)]
  pub fn intersect_sphere (&self, sphere : &Sphere2 <S>)
    -> Option <((S, Point2 <S>), (S, Point2 <S>))>
  {
    intersect::continuous_line2_sphere2 (self, sphere)
  }
}
impl <S : Real> Default for Line2 <S> {
  fn default() -> Self {
    Line2 {
      base:      Point2::origin(),
      direction: Unit2::axis_y()
    }
  }
}

// TODO: impl Primitive for Line2

impl_numcast_fields!(Line3, base, direction);
impl <S : Real> Line3 <S> {
  /// Consructs a new 3D line with given base point and direction vector
  #[inline]
  pub 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 intersect_plane (&self, plane : &Plane3 <S>)
    -> Option <(S, Point3 <S>)>
  where
    S : approx::RelativeEq
  {
    intersect::continuous_line3_plane3 (self, plane)
  }
  #[inline]
  #[allow(clippy::type_complexity)]
  pub fn intersect_aabb (&self, aabb : &Aabb3 <S>)
    -> Option <((S, Point3 <S>), (S, Point3 <S>))>
  where S : num::Float + approx::RelativeEq <Epsilon=S> + std::fmt::Debug {
    intersect::continuous_line3_aabb3 (self, aabb)
  }
  #[inline]
  #[allow(clippy::type_complexity)]
  pub fn intersect_sphere (&self, sphere : &Sphere3 <S>)
    -> Option <((S, Point3 <S>), (S, Point3 <S>))>
  {
    intersect::continuous_line3_sphere3 (self, sphere)
  }
}
impl <S : Real> Default for Line3 <S> {
  fn default() -> Self {
    Line3 {
      base:      Point3::origin(),
      direction: Unit3::axis_z()
    }
  }
}

// TODO: impl Primitive for Line3

impl_numcast_fields!(Plane3, base, normal);
impl <S : Real> Plane3 <S> {
  /// Consructs a new 3D plane with given base point and normal vector
  #[inline]
  pub 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()
    }
  }
}

// TODO: impl Primitive for Plane3

impl_numcast_fields!(Sphere2, center, radius);
impl <S : Ring> Sphere2 <S> {
  /// Unit circle
  #[inline]
  pub fn unit() -> Self where S : Field {
    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)
  }
}

// TODO: impl Primitive for Sphere2

impl_numcast_fields!(Sphere3, center, radius);
impl <S : Ring> 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)
  }
}

// TODO: impl Primitive for Sphere3

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

#[cfg(test)]
mod tests {
  use approx::{assert_relative_eq, assert_ulps_eq};
  use super::*;

  #[test]
  fn test_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), [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), [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.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),
      project_2d_point_on_line (&point, &line_b));
    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),
      project_2d_point_on_line (&point, &line_b)
    );
    assert_relative_eq!(
      project_2d_point_on_line (&point, &line_a),
      [2.0, 2.0].into());
  }

  #[test]
  fn test_project_3d_point_on_line() {
    use Unit3;
    // 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), [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), [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.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),
      project_3d_point_on_line (&point, &line_b));
    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),
      project_3d_point_on_line (&point, &line_b));
    assert_relative_eq!(
      project_3d_point_on_line (&point, &line_a),
      [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),
      project_3d_point_on_line (&point, &line_b),
      epsilon = 0.00000000000001
    );
    assert_relative_eq!(
      project_3d_point_on_line (&point, &line_a),
      [2.0/3.0, 2.0/3.0, 2.0/3.0].into(),
      epsilon = 0.00000000000001
    );
  }

  #[test]
  fn test_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 test_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);
  }
} // end tests