nsys-math-utils 1.1.1

pub use self::simplex2::Simplex2;
pub use self::simplex3::Simplex3;

/// 2D simplex types
pub mod simplex2 {
  #[cfg(feature = "derive_serdes")]
  use serde::{Deserialize, Serialize};

  use crate::*;
  use crate::geometry::*;

  /// Parameterized point on a 2D segment
  pub type SegmentPoint <S>  = (Normalized <S>, Point2 <S>);
  /// Parameterized point on a 2D triangle
  pub type TrianglePoint <S> = ([Normalized <S>; 2], Point2 <S>);

  /// A $n<3$-simplex in 2-dimensional space.
  ///
  /// Individual simplex variants should fail to construct in debug builds when points
  /// are degenerate.
  #[cfg_attr(feature = "derive_serdes", derive(Serialize, Deserialize))]
  #[derive(Clone, Copy, Debug, Eq, PartialEq)]
  pub enum Simplex2 <S> {
    Segment  (Segment  <S>),
    Triangle (Triangle <S>)
  }

  /// A 1-simplex or line segment in 2D space.
  ///
  /// Creation methods should fail with a debug assertion if the points are identical.
  #[cfg_attr(feature = "derive_serdes", derive(Serialize, Deserialize))]
  #[derive(Clone, Copy, Debug, Eq, PartialEq)]
  pub struct Segment <S> {
    a : Point2 <S>,
    b : Point2 <S>
  }

  /// A 2-simplex or triangle in 2D space
  ///
  /// Creation methods should fail with a debug assertion if the points are colinear.
  #[cfg_attr(feature = "derive_serdes", derive(Serialize, Deserialize))]
  #[derive(Clone, Copy, Debug, Eq, PartialEq)]
  pub struct Triangle <S> {
    a : Point2 <S>,
    b : Point2 <S>,
    c : Point2 <S>
  }

  impl <S : OrderedRing> Segment <S> {
    /// Returns `None` if the points are identical:
    ///
    /// ```
    /// # use math_utils::geometry::simplex2::Segment;
    /// assert!(Segment::new ([0.0, 0.0].into(), [0.0, 0.0].into()).is_none());
    /// ```
    pub fn new (a : Point2 <S>, b : Point2 <S>) -> Option <Self> {
      if a == b {
        None
      } else {
        Some (Segment { a, b })
      }
    }
    /// Panics if the points are identical:
    ///
    /// ```should_panic
    /// # use math_utils::geometry::simplex2::Segment;
    /// let s = Segment::noisy ([0.0, 0.0].into(), [0.0, 0.0].into());
    /// ```
    pub fn noisy (a : Point2 <S>, b : Point2 <S>) -> Self {
      assert_ne!(a, b);
      Segment { a, b }
    }
    /// Debug panic if the points are identical:
    ///
    /// ```should_panic
    /// # use math_utils::geometry::simplex2::Segment;
    /// let s = Segment::unchecked ([0.0, 0.0].into(), [0.0, 0.0].into());
    /// ```
    pub fn unchecked (a : Point2 <S>, b : Point2 <S>) -> Self {
      debug_assert_ne!(a, b);
      Segment { a, b }
    }
    #[inline]
    pub fn from_array ([a, b] : [Point2 <S>; 2]) -> Option <Self> {
      Self::new (a, b)
    }
    #[inline]
    pub fn numcast <T> (self) -> Option <Segment <T>> where
      S : num::NumCast,
      T : num::NumCast
    {
      Some (Segment {
        a: self.a.numcast()?,
        b: self.b.numcast()?
      })
    }
    #[inline]
    pub const fn point_a (self) -> Point2 <S> {
      self.a
    }
    #[inline]
    pub const fn point_b (self) -> Point2 <S> {
      self.b
    }
    #[inline]
    pub const fn points (self) -> [Point2 <S>; 2] {
      [self.a, self.b]
    }
    /// Return a parameterized point along the segment
    pub fn point (self, param : Normalized <S>) -> Point2 <S> {
      self.a + *self.vector() * *param
    }
    #[inline]
    pub fn length (self) -> NonNegative <S> where S : Field + Sqrt {
      self.vector().norm()
    }
    #[inline]
    pub fn length_squared (self) -> NonNegative <S> {
      self.vector().self_dot()
    }
    /// Returns `point_b() - point_a()`
    #[inline]
    pub fn vector (self) -> NonZero2 <S> {
      NonZero2::unchecked (self.b - self.a)
    }
    #[inline]
    pub fn perpendicular (self) -> NonZero2 <S> {
      self.vector().map_unchecked (Vector2::yx)
    }
    #[inline]
    pub fn translate (&mut self, displacement : Vector2 <S>) {
      self.a += displacement;
      self.b += displacement;
    }
    #[inline]
    pub fn midpoint (self) -> Point2 <S> where S : Field {
      ((self.a.0 + self.b.0) * S::half()).into()
    }
    #[inline]
    pub fn aabb2 (self) -> Aabb2 <S> {
      Aabb2::from_points_unchecked (self.a, self.b)
    }
    #[inline]
    pub fn line (self) -> Line2 <S> where S : Real {
      Line2::new (self.a, self.vector().normalize())
    }
    pub fn affine_line (self) -> frame::Line2 <S> {
      frame::Line2 {
        origin: self.a,
        basis:  self.vector()
      }
    }
    #[inline]
    pub fn nearest_point (self, point : Point2 <S>) -> SegmentPoint <S> where
      S : Field
    {
      distance::nearest_segment2_point2 (self, point)
    }
    #[inline]
    pub fn intersect_aabb (self, aabb : Aabb2 <S>)
      -> Option <(SegmentPoint <S>, SegmentPoint <S>)>
    where S : Real {
      intersect::continuous_segment2_aabb2 (self, aabb)
    }
    #[inline]
    pub fn intersect_sphere (self, sphere : Sphere2 <S>)
      -> Option <(SegmentPoint <S>, SegmentPoint <S>)>
    where S : Real {
      intersect::continuous_segment2_sphere2 (self, sphere)
    }
  }
  impl <S : Ring> Default for Segment <S> {
    /// A default simplex is arbitrarily chosen to be the simplex from -1.0 to 1.0 lying
    /// on the X axis:
    ///
    /// ```
    /// # use math_utils::geometry::simplex2::Segment;
    /// assert_eq!(
    ///   Segment::default(),
    ///   Segment::noisy ([-1.0, 0.0].into(), [1.0, 0.0].into())
    /// );
    /// ```
    fn default() -> Self {
      Segment {
        a: [-S::one(), S::zero()].into(),
        b: [ S::one(), S::zero()].into()
      }
    }
  }
  impl <S> std::ops::Index <usize> for Segment <S> where S : Copy {
    type Output = Point2 <S>;
    fn index (&self, index : usize) -> &Point2 <S> {
      [&self.a, &self.b][index]
    }
  }
  impl <S> TryFrom <[Point2 <S>; 2]> for Segment <S> where S : OrderedRing {
    type Error = [Point2 <S>; 2];
    fn try_from ([a, b] : [Point2 <S>; 2]) -> Result <Self, Self::Error> {
      Self::new (a, b).ok_or([a, b])
    }
  }

  impl <S : OrderedField> Triangle <S> {
    /// Returns `None` if the points are colinear:
    ///
    /// ```
    /// # use math_utils::geometry::simplex2::Triangle;
    /// assert!(Triangle::new (
    ///   [-1.0, -1.0].into(),
    ///   [ 0.0,  0.0].into(),
    ///   [ 1.0,  1.0].into()
    /// ).is_none());
    /// ```
    pub fn new (a : Point2 <S>, b : Point2 <S>, c : Point2 <S>) -> Option <Self> where
      S : approx::AbsDiffEq <Epsilon = S>
    {
      if colinear_2d (a, b, c) {
        None
      } else {
        Some (Triangle { a, b, c })
      }
    }
    /// Panics if the points are colinear:
    ///
    /// ```should_panic
    /// # use math_utils::geometry::simplex2::Triangle;
    /// let s = Triangle::noisy (
    ///   [-1.0, -1.0].into(),
    ///   [ 0.0,  0.0].into(),
    ///   [ 1.0,  1.0].into());
    /// ```
    pub fn noisy (a : Point2 <S>, b : Point2 <S>, c : Point2 <S>) -> Self where
      S : approx::AbsDiffEq <Epsilon = S>
    {
      debug_assert_ne!(a, b);
      debug_assert_ne!(a, c);
      debug_assert_ne!(b, c);
      assert!(!colinear_2d (a, b, c));
      Triangle { a, b, c }
    }
    /// Debug panic if the points are colinear:
    ///
    /// ```should_panic
    /// # use math_utils::geometry::simplex2::Triangle;
    /// let s = Triangle::unchecked (
    ///   [-1.0, -1.0].into(),
    ///   [ 0.0,  0.0].into(),
    ///   [ 1.0,  1.0].into());
    /// ```
    pub fn unchecked (a : Point2 <S>, b : Point2 <S>, c : Point2 <S>) -> Self where
      S : approx::AbsDiffEq <Epsilon = S>
    {
      debug_assert_ne!(a, b);
      debug_assert_ne!(a, c);
      debug_assert_ne!(b, c);
      debug_assert!(!colinear_2d (a, b, c));
      Triangle { a, b, c }
    }
    #[inline]
    pub fn from_array ([a, b, c] : [Point2 <S>; 3]) -> Option <Self> where
      S : approx::AbsDiffEq <Epsilon = S>
    {
      Self::new (a, b, c)
    }
    #[inline]
    pub const fn point_a (self) -> Point2 <S> {
      self.a
    }
    #[inline]
    pub const fn point_b (self) -> Point2 <S> {
      self.b
    }
    #[inline]
    pub const fn point_c (self) -> Point2 <S> {
      self.c
    }
    #[inline]
    pub const fn points (self) -> [Point2 <S>; 3] {
      [self.a, self.b, self.c]
    }
    /// Return a parameterized point on the triangle.
    ///
    /// Returns `None` if `s + t > 1.0`.
    pub fn point (self, s : Normalized <S>, t : Normalized <S>) -> Option <Point2 <S>> {
      if *s + *t > S::one() {
        None
      } else {
        Some (self.a + (self.b - self.a) * *s + (self.c - self.a) * *t)
      }
    }
    #[inline]
    pub const fn edge_ab (self) -> Segment2 <S> {
      Segment { a: self.a, b: self.b }
    }
    #[inline]
    pub const fn edge_bc (self) -> Segment2 <S> {
      Segment { a: self.b, b: self.c }
    }
    #[inline]
    pub const fn edge_ca (self) -> Segment2 <S> {
      Segment { a: self.c, b: self.a }
    }
    /// Returns edges `[AB, BC, CA]`
    #[inline]
    pub const fn edges (self) -> [Segment2 <S>; 3] {
      [ self.edge_ab(), self.edge_bc(), self.edge_ca() ]
    }
    /// Change triangle winding by swapping points B and C.
    ///
    /// Equivalent to `triangle.acb()`.
    #[inline]
    pub const fn rewind (self) -> Self {
      self.acb()
    }
    /// Change triangle winding by swapping points B and C
    #[inline]
    pub const fn acb (self) -> Self {
      Triangle { a: self.a, b: self.c, c: self.b }
    }
    #[inline]
    pub fn translate (&mut self, displacement : Vector2 <S>) {
      self.a += displacement;
      self.b += displacement;
      self.c += displacement;
    }
    /// Returns perpendicular vector to edge `AB`.
    ///
    /// ```
    /// # use math_utils::geometry::simplex2::Triangle;
    /// let s = Triangle::unchecked (
    ///   [0.0, 0.0].into(),
    ///   [0.0, 2.0].into(),
    ///   [2.0, 1.0].into());
    /// assert_eq!(*s.perpendicular_ab(), [-2.0, 0.0].into());
    /// let s = Triangle::unchecked (
    ///   [ 0.0, 0.0].into(),
    ///   [ 0.0, 2.0].into(),
    ///   [-2.0, 1.0].into());
    /// assert_eq!(*s.perpendicular_ab(), [2.0, 0.0].into());
    /// ```
    pub fn perpendicular_ab (self) -> NonZero2 <S> {
      let ac = self.c - self.a;
      let mut perp = self.edge_ab().perpendicular();
      if perp.dot (ac) > S::zero() {
        perp = -perp
      }
      perp
    }
    /// Returns perpendicular vector to edge `BC`.
    ///
    /// ```
    /// # use math_utils::geometry::simplex2::Triangle;
    /// let s = Triangle::unchecked (
    ///   [2.0, 1.0].into(),
    ///   [0.0, 0.0].into(),
    ///   [0.0, 2.0].into());
    /// assert_eq!(*s.perpendicular_bc(), [-2.0, 0.0].into());
    /// let s = Triangle::unchecked (
    ///   [-2.0, 1.0].into(),
    ///   [ 0.0, 0.0].into(),
    ///   [ 0.0, 2.0].into());
    /// assert_eq!(*s.perpendicular_bc(), [2.0, 0.0].into());
    /// ```
    pub fn perpendicular_bc (self) -> NonZero2 <S> {
      let bc = self.c - self.b;
      let ba = self.a - self.b;
      let mut perp = bc.yx();
      if perp.dot (ba) > S::zero() {
        perp = -perp
      }
      NonZero2::unchecked (perp)
    }
    /// Returns perpendicular vector to edge `CA`.
    ///
    /// ```
    /// # use math_utils::geometry::simplex2::Triangle;
    /// let s = Triangle::unchecked (
    ///   [0.0, 2.0].into(),
    ///   [2.0, 1.0].into(),
    ///   [0.0, 0.0].into());
    /// assert_eq!(*s.perpendicular_ca(), [-2.0, 0.0].into());
    /// let s = Triangle::unchecked (
    ///   [ 0.0, 2.0].into(),
    ///   [-2.0, 1.0].into(),
    ///   [ 0.0, 0.0].into());
    /// assert_eq!(*s.perpendicular_ca(), [2.0, 0.0].into());
    /// ```
    pub fn perpendicular_ca (self) -> NonZero2 <S> {
      let ca = self.a - self.c;
      let cb = self.b - self.c;
      let mut perp = ca.yx();
      if perp.dot (cb) > S::zero() {
        perp = -perp
      }
      NonZero2::unchecked (perp)
    }
    #[inline]
    pub fn area (self) -> Positive <S> where S : num::real::Real {
      let ab = self.b - self.a;
      let ac = self.c - self.a;
      Positive::unchecked (ab.exterior_product (ac).abs())
    }
    pub fn barycentric (self, coords : Vector3 <S>) -> Point2 <S> {
      (self.a.0 * coords.x + self.b.0 * coords.y + self.c.0 * coords.z).into()
    }
    pub fn trilinear (self, coords : Vector3 <S>) -> Point2 <S> where
      S : Real + num::real::Real
    {
      let (ab, bc, ca) = (self.b - self.a, self.c - self.b, self.a - self.c);
      let len_ab = ab.magnitude();
      let len_bc = bc.magnitude();
      let len_ca = ca.magnitude();
      let ax = len_bc * coords.x;
      let by = len_ca * coords.y;
      let cz = len_ab * coords.z;
      let denom = S::one() / (ax + by + cz);
      (self.a.0 * ax * denom + self.b.0 * by * denom + self.c.0 * cz * denom).into()
    }
    #[inline]
    pub fn affine_plane (self) -> frame::Plane2 <S> {
      frame::Plane2 {
        origin: self.a,
        basis:  frame::PlanarBasis2::new (Matrix2::from_col_arrays ([
          (self.b - self.a).into_array(),
          (self.c - self.a).into_array()
        ])).unwrap()
      }
    }
    #[inline]
    pub fn nearest_point (self, point : Point2 <S>) -> TrianglePoint <S> {
      distance::nearest_triangle2_point2 (self, point)
    }
  }
  impl <S : Field + Sqrt> Default for Triangle <S> {
    /// A default triangle is arbitrarily chosen to be the equilateral triangle with
    /// point at 1.0 on the Y axis:
    ///
    /// ```
    /// # use math_utils::approx;
    /// # use math_utils::geometry::simplex2::Triangle;
    /// use std::f64::consts::FRAC_1_SQRT_2;
    /// let s = Triangle::default();
    /// let t = Triangle::noisy (
    ///   [           0.0,            1.0].into(),
    ///   [-FRAC_1_SQRT_2, -FRAC_1_SQRT_2].into(),
    ///   [ FRAC_1_SQRT_2, -FRAC_1_SQRT_2].into()
    /// );
    /// approx::assert_relative_eq!(s.point_a(), t.point_a());
    /// approx::assert_relative_eq!(s.point_b(), t.point_b());
    /// approx::assert_relative_eq!(s.point_c(), t.point_c());
    /// ```
    fn default() -> Self {
      let frac_1_sqrt_2 = S::one() / (S::one() + S::one()).sqrt();
      Triangle {
        a: [     S::zero(),       S::one()].into(),
        b: [-frac_1_sqrt_2, -frac_1_sqrt_2].into(),
        c: [ frac_1_sqrt_2, -frac_1_sqrt_2].into()
      }
    }
  }
  impl <S> std::ops::Index <usize> for Triangle <S> where S : Copy {
    type Output = Point2 <S>;
    fn index (&self, index : usize) -> &Point2 <S> {
      [&self.a, &self.b, &self.c][index]
    }
  }
  impl <S> TryFrom <[Point2 <S>; 3]> for Triangle <S>
    where S : OrderedField + approx::AbsDiffEq <Epsilon = S>
  {
    type Error = [Point2 <S>; 3];
    fn try_from ([a, b, c] : [Point2 <S>; 3]) -> Result <Self, Self::Error> {
      Self::new (a, b, c).ok_or([a, b, c])
    }
  }
}

/// 3D simplex types
pub mod simplex3 {
  #[cfg(feature = "derive_serdes")]
  use serde::{Deserialize, Serialize};

  use crate::*;
  use crate::geometry::*;

  /// Parameterized point on a 3D segment
  pub type SegmentPoint <S>     = (Normalized <S>, Point3 <S>);
  /// Parameterized point on a 3D triangle
  pub type TrianglePoint <S>    = ([Normalized <S>; 2], Point3 <S>);
  /// Parameterized point on a 3D tetrahedron
  pub type TetrahedronPoint <S> = ([Normalized <S>; 3], Point3 <S>);

  /// A $n<4$-simplex in 3-dimensional space.
  ///
  /// Individual simplex variants should fail to construct in debug builds when points
  /// are degenerate.
  #[cfg_attr(feature = "derive_serdes", derive(Serialize, Deserialize))]
  #[derive(Clone, Copy, Debug, Eq, PartialEq)]
  pub enum Simplex3 <S> {
    Segment     (Segment     <S>),
    Triangle    (Triangle    <S>),
    Tetrahedron (Tetrahedron <S>)
  }

  /// A 1-simplex or line segment in 3D space.
  ///
  /// Creation methods should fail with a debug assertion if the points are identical.
  #[cfg_attr(feature = "derive_serdes", derive(Serialize, Deserialize))]
  #[derive(Clone, Copy, Debug, Eq, PartialEq)]
  pub struct Segment <S> {
    a : Point3 <S>,
    b : Point3 <S>
  }

  /// A 2-simplex or triangle in 3D space
  ///
  /// Creation methods should fail with a debug assertion if the points are colinear.
  #[cfg_attr(feature = "derive_serdes", derive(Serialize, Deserialize))]
  #[derive(Clone, Copy, Debug, Eq, PartialEq)]
  pub struct Triangle <S> {
    a : Point3 <S>,
    b : Point3 <S>,
    c : Point3 <S>
  }

  /// A 3-simplex or tetrahedron in 3D space
  ///
  /// Creation methods will fail with a debug assertion if the points are coplanar.
  #[cfg_attr(feature = "derive_serdes", derive(Serialize, Deserialize))]
  #[derive(Clone, Copy, Debug, Eq, PartialEq)]
  pub struct Tetrahedron <S> {
    a : Point3 <S>,
    b : Point3 <S>,
    c : Point3 <S>,
    d : Point3 <S>
  }

  impl <S : OrderedRing> Segment <S> {
    /// Returns `None` if the points are identical:
    ///
    /// ```
    /// # use math_utils::geometry::simplex3::Segment;
    /// assert!(Segment::new ([0.0, 0.0, 1.0].into(), [0.0, 0.0, 1.0].into()).is_none());
    /// ```
    pub fn new (a : Point3 <S>, b : Point3 <S>) -> Option <Self> {
      if a == b {
        None
      } else {
        Some (Segment { a, b })
      }
    }
    /// Panics if the points are identical:
    ///
    /// ```should_panic
    /// # use math_utils::geometry::simplex3::Segment;
    /// let s = Segment::noisy ([0.0, 0.0, 1.0].into(), [0.0, 0.0, 1.0].into());
    /// ```
    pub fn noisy (a : Point3 <S>, b : Point3 <S>) -> Self {
      assert_ne!(a, b);
      Segment { a, b }
    }
    /// Debug panic if the points are identical:
    ///
    /// ```should_panic
    /// # use math_utils::geometry::simplex3::Segment;
    /// let s = Segment::unchecked ([0.0, 0.0, 1.0].into(), [0.0, 0.0, 1.0].into());
    /// ```
    pub fn unchecked (a : Point3 <S>, b : Point3 <S>) -> Self {
      debug_assert_ne!(a, b);
      Segment { a, b }
    }
    #[inline]
    pub fn from_array ([a, b] : [Point3 <S>; 2]) -> Option <Self> {
      Self::new (a, b)
    }
    #[inline]
    pub fn numcast <T> (self) -> Option <Segment <T>> where
      S : num::NumCast,
      T : num::NumCast
    {
      Some (Segment {
        a: self.a.numcast()?,
        b: self.b.numcast()?
      })
    }
    #[inline]
    pub const fn point_a (self) -> Point3 <S> {
      self.a
    }
    #[inline]
    pub const fn point_b (self) -> Point3 <S> {
      self.b
    }
    #[inline]
    pub const fn points (self) -> [Point3 <S>; 2] {
      [self.a, self.b]
    }
    /// Return a parameterized point along the segment
    pub fn point (self, param : Normalized <S>) -> Point3 <S> {
      self.a + *self.vector() * *param
    }
    #[inline]
    pub fn length (self) -> NonNegative <S> where S : Field + Sqrt {
      self.vector().norm()
    }
    #[inline]
    pub fn length_squared (self) -> NonNegative <S> {
      self.vector().self_dot()
    }
    /// Returns `point_b() - point_a()`
    #[inline]
    pub fn vector (self) -> NonZero3 <S> {
      NonZero3::unchecked (self.b - self.a)
    }
    /// Equivalent to `self.vector().orthogonal()`
    #[inline]
    pub fn perpendicular (self) -> NonZero3 <S> where
      S : Field + approx::AbsDiffEq <Epsilon = S>
    {
      NonZero3::unchecked (self.vector().orthogonal())
    }
    #[inline]
    pub fn translate (&mut self, displacement : Vector3 <S>) {
      self.a += displacement;
      self.b += displacement;
    }
    #[inline]
    pub fn midpoint (self) -> Point3 <S> where S : Field {
      ((self.a.0 + self.b.0) * S::half()).into()
    }
    #[inline]
    pub fn aabb3 (self) -> Aabb3 <S> {
      Aabb3::from_points_unchecked (self.a, self.b)
    }
    #[inline]
    pub fn line (self) -> Line3 <S> where S : Field + Sqrt {
      Line3::new (self.a, Unit3::normalize (*self.vector()))
    }
    pub fn affine_line (self) -> frame::Line3 <S> {
      frame::Line3 {
        origin: self.a,
        basis:  self.vector()
      }
    }
    #[inline]
    pub fn nearest_point (self, point : Point3 <S>) -> SegmentPoint <S> where
      S : Field
    {
      distance::nearest_segment3_point3 (self, point)
    }
    #[inline]
    pub fn intersect_aabb (self, aabb : Aabb3 <S>)
      -> Option <(SegmentPoint <S>, SegmentPoint <S>)>
    where S : Real + num::Float + approx::RelativeEq <Epsilon=S> {
      intersect::continuous_segment3_aabb3 (self, aabb)
    }
    #[inline]
    pub fn intersect_triangle (self, triangle : Triangle3 <S>)
      -> Option <SegmentPoint <S>>
    where S : Real + num::real::Real + approx::AbsDiffEq <Epsilon=S> {
      intersect::continuous_triangle3_segment3 (triangle, self)
    }
    #[inline]
    pub fn intersect_sphere (self, sphere : Sphere3 <S>)
      -> Option <(SegmentPoint <S>, SegmentPoint <S>)>
    where S : Field + Sqrt {
      intersect::continuous_segment3_sphere3 (self, sphere)
    }
    #[inline]
    pub fn intersect_cylinder (self, sphere : Cylinder3 <S>)
      -> Option <(SegmentPoint <S>, SegmentPoint <S>)>
    where S : Real {
      intersect::continuous_segment3_cylinder3 (self, sphere)
    }
    #[inline]
    pub fn intersect_capsule (self, capsule : Capsule3 <S>)
      -> Option <(SegmentPoint <S>, SegmentPoint <S>)>
    where S : Real {
      intersect::continuous_segment3_capsule3 (self, capsule)
    }
    #[inline]
    pub fn intersect_hull (self, hull : &Hull3 <S>)
      -> Option <(SegmentPoint <S>, SegmentPoint <S>)>
    {
      intersect::continuous_segment3_hull3 (self, hull)
    }
  }
  impl <S : Ring> Default for Segment <S> {
    /// A default simplex is arbitrarily chosen to be the simplex from -1.0 to 1.0 lying
    /// on the X axis:
    ///
    /// ```
    /// # use math_utils::geometry::simplex3::Segment;
    /// assert_eq!(
    ///   Segment::default(),
    ///   Segment::noisy ([-1.0, 0.0, 0.0].into(), [1.0, 0.0, 0.0].into())
    /// );
    /// ```
    fn default() -> Self {
      Segment {
        a: [-S::one(), S::zero(), S::zero()].into(),
        b: [ S::one(), S::zero(), S::zero()].into()
      }
    }
  }
  impl <S> std::ops::Index <usize> for Segment <S> where S : Copy {
    type Output = Point3 <S>;
    fn index (&self, index : usize) -> &Point3 <S> {
      [&self.a, &self.b][index]
    }
  }
  impl <S> TryFrom <[Point3 <S>; 2]> for Segment <S> where S : OrderedRing {
    type Error = [Point3 <S>; 2];
    fn try_from ([a, b] : [Point3 <S>; 2]) -> Result <Self, Self::Error> {
      Self::new (a, b).ok_or([a, b])
    }
  }

  impl <S : OrderedField> Triangle <S> {
    /// Returns `None` if points are colinear:
    ///
    /// ```
    /// # use math_utils::geometry::simplex3::Triangle;
    /// assert!(Triangle::new (
    ///   [-1.0, -1.0, -1.0].into(),
    ///   [ 0.0,  0.0,  0.0].into(),
    ///   [ 1.0,  1.0,  1.0].into()
    /// ).is_none());
    /// ```
    pub fn new (a : Point3 <S>, b : Point3 <S>, c : Point3 <S>) -> Option <Self> where
      S : Sqrt + approx::RelativeEq <Epsilon = S>
    {
      if colinear_3d (a, b, c) {
        None
      } else {
        Some (Triangle { a, b, c })
      }
    }
    /// Panics if the points are colinear:
    ///
    /// ```should_panic
    /// # use math_utils::geometry::simplex3::Triangle;
    /// let s = Triangle::noisy (
    ///   [-1.0, -1.0, -1.0].into(),
    ///   [ 0.0,  0.0,  0.0].into(),
    ///   [ 1.0,  1.0,  1.0].into());
    /// ```
    pub fn noisy (a : Point3 <S>, b : Point3 <S>, c : Point3 <S>) -> Self where
      S : Sqrt + approx::RelativeEq <Epsilon = S>
    {
      debug_assert_ne!(a, b);
      debug_assert_ne!(a, c);
      debug_assert_ne!(b, c);
      assert!(!colinear_3d (a, b, c));
      Triangle { a, b, c }
    }
    /// Debug panic if the points are colinear:
    ///
    /// ```should_panic
    /// # use math_utils::geometry::simplex3::Triangle;
    /// let s = Triangle::unchecked (
    ///   [-1.0, -1.0, -1.0].into(),
    ///   [ 0.0,  0.0,  0.0].into(),
    ///   [ 1.0,  1.0,  1.0].into());
    /// ```
    pub fn unchecked (a : Point3 <S>, b : Point3 <S>, c : Point3 <S>) -> Self where
      S : Sqrt + approx::RelativeEq <Epsilon = S>
    {
      debug_assert_ne!(a, b);
      debug_assert_ne!(a, c);
      debug_assert_ne!(b, c);
      debug_assert!(!colinear_3d (a, b, c));
      Triangle { a, b, c }
    }
    #[inline]
    pub fn from_array ([a, b, c] : [Point3 <S>; 3]) -> Option <Self> where
      S : Sqrt + approx::RelativeEq <Epsilon = S>
    {
      Self::new (a, b, c)
    }
    #[inline]
    pub const fn point_a (self) -> Point3 <S> {
      self.a
    }
    #[inline]
    pub const fn point_b (self) -> Point3 <S> {
      self.b
    }
    #[inline]
    pub const fn point_c (self) -> Point3 <S> {
      self.c
    }
    #[inline]
    pub const fn points (self) -> [Point3 <S>; 3] {
      [self.a, self.b, self.c]
    }
    /// Return a parameterized point on the triangle.
    ///
    /// Returns `None` if `s + t > 1.0`.
    pub fn point (self, s : Normalized <S>, t : Normalized <S>) -> Option <Point3 <S>> {
      if *s + *t > S::one() {
        None
      } else {
        Some (self.a + (self.b - self.a) * *s + (self.c - self.a) * *t)
      }
    }
    #[inline]
    pub const fn edge_ab (self) -> Segment3 <S> {
      Segment { a: self.a, b: self.b }
    }
    #[inline]
    pub const fn edge_bc (self) -> Segment3 <S> {
      Segment { a: self.b, b: self.c }
    }
    #[inline]
    pub const fn edge_ca (self) -> Segment3 <S> {
      Segment { a: self.c, b: self.a }
    }
    #[inline]
    pub const fn edges (self) -> [Segment3 <S>; 3] {
      [ Segment { a: self.a, b: self.b },
        Segment { a: self.b, b: self.c },
        Segment { a: self.c, b: self.a }
      ]
    }
    /// Change triangle winding by swapping points B and C.
    ///
    /// Equivalent to `triangle.acb()`.
    #[inline]
    pub const fn rewind (self) -> Self {
      self.acb()
    }
    /// Change triangle winding by swapping points B and C
    #[inline]
    pub const fn acb (self) -> Self {
      Triangle { a: self.a, b: self.c, c: self.b }
    }
    pub fn translate (&mut self, displacement : Vector3 <S>) where
      S : Copy + AdditiveGroup
    {
      self.a += displacement;
      self.b += displacement;
      self.c += displacement;
    }
    /// Cross product of edges `AB x AC`
    #[inline]
    pub fn normal (self) -> NonZero3 <S> {
      let ab = self.b - self.a;
      let ac = self.c - self.b;
      NonZero3::unchecked (ab.cross (ac))
    }
    /// Normalized cross product of edges `AB x AC`
    #[inline]
    pub fn unit_normal (self) -> Unit3 <S> where S : Sqrt {
      self.normal().normalize()
    }
    /// Returns perpendicular vector to edge `AB`.
    ///
    /// ```
    /// # use math_utils::*;
    /// # use math_utils::geometry::simplex3::Triangle;
    /// let s = Triangle::unchecked (
    ///   [0.0, 0.0, 0.0].into(),
    ///   [0.0, 2.0, 0.0].into(),
    ///   [2.0, 1.0, 0.0].into());
    /// assert_eq!(*s.perpendicular_ab(), vector![-8.0, 0.0, 0.0]);
    /// let s = Triangle::unchecked (
    ///   [ 0.0, 0.0, 0.0].into(),
    ///   [ 0.0, 2.0, 0.0].into(),
    ///   [-2.0, 1.0, 0.0].into());
    /// assert_eq!(*s.perpendicular_ab(), vector![8.0, 0.0, 0.0]);
    /// ```
    pub fn perpendicular_ab (self) -> NonZero3 <S> {
      let ab = self.b - self.a;
      let ac = self.c - self.a;
      let perp = ab.cross (ab.cross (ac));
      debug_assert!(perp.dot (ac) < S::zero());
      NonZero3::unchecked (perp)
    }
    /// Returns perpendicular vector to edge `BC`.
    ///
    /// ```
    /// # use math_utils::*;
    /// # use math_utils::geometry::simplex3::Triangle;
    /// let s = Triangle::unchecked (
    ///   [2.0, 1.0, 0.0].into(),
    ///   [0.0, 0.0, 0.0].into(),
    ///   [0.0, 2.0, 0.0].into());
    /// assert_eq!(*s.perpendicular_bc(), vector![-8.0, 0.0, 0.0]);
    /// let s = Triangle::unchecked (
    ///   [-2.0, 1.0, 0.0].into(),
    ///   [ 0.0, 0.0, 0.0].into(),
    ///   [ 0.0, 2.0, 0.0].into());
    /// assert_eq!(*s.perpendicular_bc(), vector![8.0, 0.0, 0.0]);
    /// ```
    pub fn perpendicular_bc (self) -> NonZero3 <S> {
      let bc = self.c - self.b;
      let ba = self.a - self.b;
      let perp = bc.cross (bc.cross (ba));
      debug_assert!(perp.dot (ba) < S::zero());
      NonZero3::unchecked (perp)
    }
    /// Returns perpendicular vector to edge `CA`.
    ///
    /// ```
    /// # use math_utils::*;
    /// # use math_utils::geometry::simplex3::Triangle;
    /// let s = Triangle::unchecked (
    ///   [0.0, 2.0, 0.0].into(),
    ///   [2.0, 1.0, 0.0].into(),
    ///   [0.0, 0.0, 0.0].into());
    /// assert_eq!(*s.perpendicular_ca(), vector![-8.0, 0.0, 0.0]);
    /// let s = Triangle::unchecked (
    ///   [ 0.0, 2.0, 0.0].into(),
    ///   [-2.0, 1.0, 0.0].into(),
    ///   [ 0.0, 0.0, 0.0].into());
    /// assert_eq!(*s.perpendicular_ca(), vector![8.0, 0.0, 0.0]);
    /// ```
    pub fn perpendicular_ca (self) -> NonZero3 <S> {
      let ca = self.a - self.c;
      let cb = self.b - self.c;
      let perp = ca.cross (ca.cross (cb));
      debug_assert!(perp.dot (cb) < S::zero());
      NonZero3::unchecked (perp)
    }
    #[inline]
    pub fn area (self) -> Positive <S> where S : Sqrt {
      Positive::unchecked (*self.normal().norm() * S::half())
    }
    pub fn barycentric (self, coords : Vector3 <S>) -> Point3 <S> {
      (self.a.0 * coords.x + self.b.0 * coords.y + self.c.0 * coords.z).into()
    }
    pub fn trilinear (self, coords : Vector3 <S>) -> Point3 <S> where
      S : Real + num::real::Real
    {
      let (ab, bc, ca) = (self.b - self.a, self.c - self.b, self.a - self.c);
      let len_ab = ab.magnitude();
      let len_bc = bc.magnitude();
      let len_ca = ca.magnitude();
      let ax = len_bc * coords.x;
      let by = len_ca * coords.y;
      let cz = len_ab * coords.z;
      let denom = S::one() / (ax + by + cz);
      (self.a.0 * ax * denom + self.b.0 * by * denom + self.c.0 * cz * denom).into()
    }
    #[inline]
    pub fn affine_plane (self) -> frame::Plane3 <S> {
      frame::Plane3 {
        origin: self.a,
        basis:  frame::PlanarBasis3::new ([
          (self.b - self.a),
          (self.c - self.a)
        ].map (NonZero3::unchecked)).unwrap()
      }
    }
    #[inline]
    pub fn nearest_point (self, point : Point3 <S>) -> TrianglePoint <S> {
      distance::nearest_triangle3_point3 (self, point)
    }
  }
  impl <S : Field + Sqrt> Default for Triangle <S> {
    /// A default simplex is arbitrarily chosen to be the simplex on the XY plane formed
    /// by an equilateral triangle with point at 1.0 on the Y axis:
    ///
    /// ```
    /// # use math_utils::approx;
    /// # use math_utils::geometry::simplex3::Triangle;
    /// # use math_utils::*;
    /// use std::f64::consts::FRAC_1_SQRT_2;
    /// let s = Triangle::default();
    /// let t = Triangle::noisy (
    ///   [           0.0,            1.0, 0.0].into(),
    ///   [-FRAC_1_SQRT_2, -FRAC_1_SQRT_2, 0.0].into(),
    ///   [ FRAC_1_SQRT_2, -FRAC_1_SQRT_2, 0.0].into()
    /// );
    /// approx::assert_relative_eq!(
    ///   Matrix3::from_col_arrays ([
    ///     s.point_a().0.into_array(),
    ///     s.point_b().0.into_array(),
    ///     s.point_c().0.into_array()
    ///   ]),
    ///   Matrix3::from_col_arrays ([
    ///     t.point_a().0.into_array(),
    ///     t.point_b().0.into_array(),
    ///     t.point_c().0.into_array()
    ///   ]));
    /// ```
    fn default() -> Self {
      let frac_1_sqrt_2 = S::one() / (S::one() + S::one()).sqrt();
      Triangle {
        a: [     S::zero(),       S::one(), S::zero()].into(),
        b: [-frac_1_sqrt_2, -frac_1_sqrt_2, S::zero()].into(),
        c: [ frac_1_sqrt_2, -frac_1_sqrt_2, S::zero()].into()
      }
    }
  }
  impl <S> std::ops::Index <usize> for Triangle <S> where S : Copy {
    type Output = Point3 <S>;
    fn index (&self, index : usize) -> &Point3 <S> {
      [&self.a, &self.b, &self.c][index]
    }
  }
  impl <S> TryFrom <[Point3 <S>; 3]> for Triangle <S> where
    S : OrderedField + Sqrt + approx::RelativeEq <Epsilon = S>
  {
    type Error = [Point3 <S>; 3];
    fn try_from ([a, b, c] : [Point3 <S>; 3]) -> Result <Self, Self::Error> {
      Self::new (a, b, c).ok_or([a, b, c])
    }
  }

  impl <S : OrderedField> Tetrahedron <S> {
    /// Returns `None` if the points are coplanar:
    ///
    /// ```
    /// # use math_utils::geometry::simplex3::Tetrahedron;
    /// assert!(Tetrahedron::new (
    ///   [-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()
    /// ).is_none());
    /// ```
    pub fn new (a : Point3 <S>, b : Point3 <S>, c : Point3 <S>, d : Point3 <S>)
      -> Option <Self>
    where S : approx::AbsDiffEq <Epsilon = S> {
      if coplanar_3d (a, b, c, d) {
        None
      } else {
        Some (Tetrahedron { a, b, c, d })
      }
    }
    /// Panics if the points are coplanar:
    ///
    /// ```should_panic
    /// # use math_utils::geometry::simplex3::Tetrahedron;
    /// let s = Tetrahedron::noisy (
    ///   [-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 noisy (a : Point3 <S>, b : Point3 <S>, c : Point3 <S>, d : Point3 <S>)
      -> Self
    where S : approx::AbsDiffEq <Epsilon = S> {
      assert!(!coplanar_3d (a, b, c, d));
      Tetrahedron { a, b, c, d }
    }
    /// Debug panic if the points are coplanar:
    ///
    /// ```should_panic
    /// # use math_utils::geometry::simplex3::Tetrahedron;
    /// let s = Tetrahedron::unchecked (
    ///   [-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 unchecked (a : Point3 <S>, b : Point3 <S>, c : Point3 <S>, d : Point3 <S>)
      -> Self
    where S : approx::AbsDiffEq <Epsilon = S> {
      debug_assert!(!coplanar_3d (a, b, c, d));
      Tetrahedron { a, b, c, d }
    }
    #[inline]
    pub fn from_array ([a, b, c, d] : [Point3 <S>; 4]) -> Option <Self> where
      S : approx::AbsDiffEq <Epsilon = S>
    {
      Self::new (a, b, c, d)
    }
    #[inline]
    pub const fn point_a (self) -> Point3 <S> {
      self.a
    }
    #[inline]
    pub const fn point_b (self) -> Point3 <S> {
      self.b
    }
    #[inline]
    pub const fn point_c (self) -> Point3 <S> {
      self.c
    }
    #[inline]
    pub const fn point_d (self) -> Point3 <S> {
      self.d
    }
    #[inline]
    pub const fn points (self) -> [Point3 <S>; 4] {
      [self.a, self.b, self.c, self.d]
    }
    #[inline]
    pub const fn edge_ab (self) -> Segment3 <S> {
      Segment { a: self.a, b: self.b }
    }
    #[inline]
    pub const fn edge_ac (self) -> Segment3 <S> {
      Segment { a: self.a, b: self.c }
    }
    #[inline]
    pub const fn edge_ad (self) -> Segment3 <S> {
      Segment { a: self.a, b: self.d }
    }
    #[inline]
    pub const fn edge_bc (self) -> Segment3 <S> {
      Segment { a: self.b, b: self.c }
    }
    #[inline]
    pub const fn edge_bd (self) -> Segment3 <S> {
      Segment { a: self.b, b: self.d }
    }
    #[inline]
    pub const fn edge_cd (self) -> Segment3 <S> {
      Segment { a: self.c, b: self.d }
    }
    #[inline]
    pub const fn edges (self) -> [Segment3 <S>; 6] {
      [ self.edge_ab(),
        self.edge_ac(),
        self.edge_ad(),
        self.edge_bc(),
        self.edge_bd(),
        self.edge_cd()
      ]
    }
    /// Note that the points in the returned triangle will be ordered so that the
    /// `triangle.normal()` will face away from the tetrahedron.
    #[inline]
    pub fn face_abc (self) -> Triangle <S> {
      let triangle = Triangle { a: self.a, b: self.b, c: self.c };
      if triangle.normal().dot (self.d - self.a) > S::zero() {
        triangle.rewind()
      } else {
        triangle
      }
    }
    /// Note that the points in the returned triangle will be ordered so that the
    /// `triangle.normal()` will face away from the tetrahedron.
    #[inline]
    pub fn face_abd (self) -> Triangle <S> {
      let triangle = Triangle { a: self.a, b: self.b, c: self.d };
      if triangle.normal().dot (self.c - self.a) > S::zero() {
        triangle.rewind()
      } else {
        triangle
      }
    }
    /// Note that the points in the returned triangle will be ordered so that the
    /// `triangle.normal()` will face away from the tetrahedron.
    #[inline]
    pub fn face_acd (self) -> Triangle <S> {
      let triangle = Triangle { a: self.a, b: self.c, c: self.d };
      if triangle.normal().dot (self.b - self.a) > S::zero() {
        triangle.rewind()
      } else {
        triangle
      }
    }
    /// Note that the points in the returned triangle will be ordered so that the
    /// `triangle.normal()` will face away from the tetrahedron.
    #[inline]
    pub fn face_bcd (self) -> Triangle <S> {
      let triangle = Triangle { a: self.b, b: self.c, c: self.d };
      if triangle.normal().dot (self.a - self.b) > S::zero() {
        triangle.rewind()
      } else {
        triangle
      }
    }
    #[inline]
    pub fn faces (self) -> [Triangle <S>; 4] {
      [ self.face_abc(),
        self.face_abd(),
        self.face_acd(),
        self.face_bcd()
      ]
    }
    pub fn normal_abc (self) -> NonZero3 <S> {
      let mut abc_normal = self.face_abc().normal();
      let ad = self.d - self.a;
      if abc_normal.dot (ad) > S::zero() {
        abc_normal = -abc_normal;
      }
      abc_normal
    }
    pub fn normal_abd (self) -> NonZero3 <S> {
      let mut abd_normal = self.face_abd().normal();
      let ac = self.c - self.a;
      if abd_normal.dot (ac) > S::zero() {
        abd_normal = -abd_normal;
      }
      abd_normal
    }
    pub fn normal_acd (self) -> NonZero3 <S> {
      let mut acd_normal = self.face_acd().normal();
      let ab = self.b - self.a;
      if acd_normal.dot (ab) > S::zero() {
        acd_normal = -acd_normal;
      }
      acd_normal
    }
    pub fn normal_bcd (self) -> NonZero3 <S> {
      let mut bcd_normal = self.face_bcd().normal();
      let ba = self.a - self.b;
      if bcd_normal.dot (ba) > S::zero() {
        bcd_normal = -bcd_normal;
      }
      bcd_normal
    }
    pub fn volume (self) -> Positive <S> {
      Positive::unchecked (
        (self.b - self.a).cross (self.c - self.a).dot (self.d - self.a).abs()
          * (S::one() / S::six()))
    }
    #[inline]
    pub fn translate (&mut self, displacement : Vector3 <S>) {
      self.a += displacement;
      self.b += displacement;
      self.c += displacement;
      self.d += displacement;
    }
    /// Checks if a point is contained in the tetrahedron:
    ///
    /// ```
    /// # use math_utils::*;
    /// # use math_utils::geometry::simplex3::Tetrahedron;
    /// let s = Tetrahedron::noisy (
    ///   [ 0.0,  0.0,  1.0].into(),
    ///   [ 0.0,  1.0, -1.0].into(),
    ///   [-1.0, -1.0, -1.0].into(),
    ///   [ 1.0, -1.0, -1.0].into()
    /// );
    /// assert!(s.contains (Point3::origin()));
    /// assert!(!s.contains ([0.0, 0.0, 2.0].into()));
    /// ```
    #[inline]
    pub fn contains (self, point : Point3 <S>) -> bool {
      let ap = point - self.a;
      let bp = point - self.b;
      self.normal_abc().dot (ap) < S::zero() &&
      self.normal_abd().dot (ap) < S::zero() &&
      self.normal_acd().dot (ap) < S::zero() &&
      self.normal_bcd().dot (bp) < S::zero()
    }
  }
  impl <S : Field + Sqrt> Default for Tetrahedron <S> {
    /// A default simplex is arbitrarily chosen to be the simplex with vertices at unit
    /// distance from the origin with A at `[0.0, 0.0, 1.0]` and B at
    /// `[0.0, sqrt(8.0/9.0), -1.0/3.0]`, and points C and D at
    /// `[ sqrt(2.0/3.0), -sqrt(2.0/9.0), -1.0/3.0]` and
    /// `[-sqrt(2.0/3.0), -sqrt(2.0/9.0), -1.0/3.0]`.
    ///
    /// ```
    /// # use math_utils::approx;
    /// # use math_utils::geometry::simplex3::Tetrahedron;
    /// # use math_utils::*;
    /// let s = Tetrahedron::default();
    /// let t = Tetrahedron::noisy (
    ///   [                0.0,                 0.0,      1.0].into(),
    ///   [                0.0,  f64::sqrt(8.0/9.0), -1.0/3.0].into(),
    ///   [ f64::sqrt(2.0/3.0), -f64::sqrt(2.0/9.0), -1.0/3.0].into(),
    ///   [-f64::sqrt(2.0/3.0), -f64::sqrt(2.0/9.0), -1.0/3.0].into());
    /// approx::assert_relative_eq!(
    ///   Matrix4::from_col_arrays ([
    ///     s.point_a().0.with_w (0.0).into_array(),
    ///     s.point_b().0.with_w (0.0).into_array(),
    ///     s.point_c().0.with_w (0.0).into_array(),
    ///     s.point_d().0.with_w (0.0).into_array()
    ///   ]),
    ///   Matrix4::from_col_arrays ([
    ///     t.point_a().0.with_w (0.0).into_array(),
    ///     t.point_b().0.with_w (0.0).into_array(),
    ///     t.point_c().0.with_w (0.0).into_array(),
    ///     t.point_d().0.with_w (0.0).into_array()
    ///   ]));
    /// ```
    fn default() -> Self {
      let frac_1_3 = S::one()    / S::three();
      let sqrt_2_3 = (S::two()   / S::three()).sqrt();
      let sqrt_8_9 = (S::eight() / S::nine()).sqrt();
      let sqrt_2_9 = (S::two()   / S::nine()).sqrt();
      Tetrahedron {
        a: [S::zero(), S::zero(),  S::one()].into(),
        b: [S::zero(),  sqrt_8_9, -frac_1_3].into(),
        c: [ sqrt_2_3, -sqrt_2_9, -frac_1_3].into(),
        d: [-sqrt_2_3, -sqrt_2_9, -frac_1_3].into()
      }
    }
  }
  impl <S> std::ops::Index <usize> for Tetrahedron <S> where S : Copy {
    type Output = Point3 <S>;
    fn index (&self, index : usize) -> &Point3 <S> {
      [&self.a, &self.b, &self.c, &self.d][index]
    }
  }

  impl <S> TryFrom <[Point3 <S>; 4]> for Tetrahedron <S> where
    S : OrderedField + approx::RelativeEq <Epsilon = S>
  {
    type Error = [Point3 <S>; 4];
    fn try_from ([a, b, c, d] : [Point3 <S>; 4]) -> Result <Self, Self::Error> {
      Self::new (a, b, c, d).ok_or([a, b, c, d])
    }
  }
}

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

  #[test]
  fn barycentric_2() {
    let triangle = simplex2::Triangle::noisy (
      [-2.0, -2.0].into(), [2.0, -2.0].into(), [0.0, 2.0].into());
    assert_eq!(triangle.barycentric ([1.0, 0.0, 0.0].into()), [-2.0, -2.0].into());
    assert_eq!(triangle.barycentric ([0.0, 1.0, 0.0].into()), [ 2.0, -2.0].into());
    assert_eq!(triangle.barycentric ([0.0, 0.0, 1.0].into()), [ 0.0,  2.0].into());
    assert_eq!(triangle.barycentric ([0.5, 0.5, 0.0].into()), [ 0.0, -2.0].into());
    assert_eq!(triangle.barycentric ([0.5, 0.0, 0.5].into()), [-1.0,  0.0].into());
    assert_eq!(triangle.barycentric ([0.0, 0.5, 0.5].into()), [ 1.0,  0.0].into());
  }

  #[test]
  fn barycentric_3() {
    let triangle = simplex3::Triangle::noisy (
      [-2.0, -2.0, 0.0].into(), [2.0, -2.0, 0.0].into(), [0.0, 2.0, 0.0].into());
    assert_eq!(triangle.barycentric ([1.0, 0.0, 0.0].into()), [-2.0, -2.0, 0.0].into());
    assert_eq!(triangle.barycentric ([0.0, 1.0, 0.0].into()), [ 2.0, -2.0, 0.0].into());
    assert_eq!(triangle.barycentric ([0.0, 0.0, 1.0].into()), [ 0.0,  2.0, 0.0].into());
    assert_eq!(triangle.barycentric ([0.5, 0.5, 0.0].into()), [ 0.0, -2.0, 0.0].into());
    assert_eq!(triangle.barycentric ([0.5, 0.0, 0.5].into()), [-1.0,  0.0, 0.0].into());
    assert_eq!(triangle.barycentric ([0.0, 0.5, 0.5].into()), [ 1.0,  0.0, 0.0].into());
  }

  #[test]
  fn trilinear_2() {
    let triangle = simplex2::Triangle::noisy (
      [-2.0, -2.0].into(), [2.0, -2.0].into(), [0.0, 2.0].into());
    assert_eq!(triangle.trilinear ([1.0, 0.0, 0.0].into()), [-2.0, -2.0].into());
    assert_eq!(triangle.trilinear ([0.0, 1.0, 0.0].into()), [ 2.0, -2.0].into());
    assert_eq!(triangle.trilinear ([0.0, 0.0, 1.0].into()), [ 0.0,  2.0].into());
    // centroid
    assert_eq!(triangle.trilinear (
      [1.0 / 20.0f32.sqrt(), 1.0 / 20.0f32.sqrt(), 0.25].into()),
      [0.0,  -2.0/3.0].into());
    // incenter
    approx::assert_relative_eq!(triangle.trilinear (
      [1.0, 1.0, 1.0].into()),
      [ 0.0, 5.0f32.sqrt() - 3.0].into());
  }

  #[test]
  fn trilinear_3() {
    let triangle = simplex3::Triangle::noisy (
      [-2.0, -2.0, 0.0].into(), [2.0, -2.0, 0.0].into(), [0.0, 2.0, 0.0].into());
    assert_eq!(triangle.trilinear ([1.0, 0.0, 0.0].into()), [-2.0, -2.0, 0.0].into());
    assert_eq!(triangle.trilinear ([0.0, 1.0, 0.0].into()), [ 2.0, -2.0, 0.0].into());
    assert_eq!(triangle.trilinear ([0.0, 0.0, 1.0].into()), [ 0.0,  2.0, 0.0].into());
    // centroid
    assert_eq!(triangle.trilinear (
      [1.0 / 20.0f32.sqrt(), 1.0 / 20.0f32.sqrt(), 0.25].into()),
      [0.0,  -2.0/3.0, 0.0].into());
    // incenter
    approx::assert_relative_eq!(triangle.trilinear (
      [1.0, 1.0, 1.0].into()),
      [ 0.0, 5.0f32.sqrt() - 3.0, 0.0].into());
  }

  #[test]
  fn tetrahedron_face_normals() {
    // face normals should face out from tetrahedron
    use rand;
    use rand_xorshift::XorShiftRng;
    use rand::SeedableRng;
    let aabb = Aabb3::<f32>::with_minmax_unchecked (
      [-10.0, -10.0, -10.0].into(),
      [ 10.0,  10.0,  10.0].into());
    let mut rng = XorShiftRng::seed_from_u64 (0);
    for _ in 0..10000 {
      let t = Tetrahedron3::noisy (
        aabb.rand_point (&mut rng),
        aabb.rand_point (&mut rng),
        aabb.rand_point (&mut rng),
        aabb.rand_point (&mut rng));
      for face in t.faces() {
        let other_p = {
          let mut point = None;
          for p in t.points() {
            if !face.points().contains (&p) {
              point = Some (p);
              break
            }
          }
          point.unwrap()
        };
        assert!(face.normal().dot (other_p - face.point_a()) < 0.0);
      }
    }
  }
}