linear-sim 0.7.0

use crate::{collision, component, constraint, object};
use crate::geometry::shape;
use crate::math::*;

mod gjk;

#[cfg(test)]
mod bench;

/// Relation configuration of a pair of bounded objects.
///
/// - 'Disjoint' corresponds to a positive distance `> CONTACT_DISTANCE`
/// - 'Contact' corresponds to a positive distance `<= CONTACT_DISTANCE`
/// - 'Inersect' corresponds to a strictly negative contact distance `< 0.0`
#[derive(Clone, Copy, Debug, Eq, PartialEq)]
pub enum Relation {
  Disjoint,
  Contact,
  Intersect
}

/// Result of a proximity query between a pair of bounded objects.
///
/// Note that the *axis vector* points from object A to object B, while the *normal
/// vector* points from object B to object A, but they should be parallel when the axis
/// vector is non-zero.
#[derive(Clone, Debug, PartialEq)]
pub struct Proximity {
  /// Signed distance.
  ///
  /// A positive distance indicates the separating distance for a disjoint result, and a
  /// negative distance indicates a penetration depth of an intersection result.
  pub distance  : f64,
  /// Non-normalized vector pointing along the axis of separation or penetration from
  /// object A to object B.
  ///
  /// For a disjoint configuration, adding this vector to the midpoint gives the nearest
  /// point on object B and subtracting this vector from the midpoint gives the nearest
  /// point on object A. Translating object A by this vector and translating object B by
  /// its inverse will bring the objects into contact.
  ///
  /// For an intersecting configuration, translating object A by this vector and
  /// translating object B by its inverse will resolve the inter-penetration.
  ///
  /// Note this vector can be zero if the objects are in contact.
  pub half_axis : Vector3 <f64>,
  /// Indicates the midpoint of the proximity query
  pub midpoint  : Point3  <f64>,
  /// The unit normal to the separating plane directed from object B to object A
  pub normal    : Unit3 <f64>
}

impl Proximity {
  pub fn nearest_a (&self) -> Point3 <f64> {
    self.midpoint - self.half_axis
  }

  pub fn nearest_b (&self) -> Point3 <f64> {
    self.midpoint + self.half_axis
  }

  pub fn relation (&self) -> Relation {
    if self.distance < 0.0 {
      Relation::Intersect
    } else if self.distance <= collision::CONTACT_DISTANCE {
      Relation::Contact
    } else {
      debug_assert!(self.distance > collision::CONTACT_DISTANCE);
      Relation::Disjoint
    }
  }

  #[inline]
  pub const fn constraint_planar (&self) -> constraint::Planar {
    constraint::Planar::new (self.midpoint, self.normal)
  }

  /// Flip the axis and normal so that the relation of objects A and B are reversed
  pub fn flip (mut self) -> Self {
    self.half_axis = -self.half_axis;
    self.normal = -self.normal;
    self
  }

  /// Proximity query.
  ///
  /// # Examples
  ///
  /// **Capsule v. capsule**
  ///
  /// By convention if the axes of capsules A and B intersect, the +X axis will be
  /// chosen for the axis vector and normal:
  ///
  /// ```
  /// # extern crate linear_sim;
  /// # use linear_sim::*;
  /// # use collision::*;
  /// # use geometry::shape;
  /// # use math;
  /// # fn main() {
  /// let a = object::Static {
  ///   position:   component::Position ([0.0, 0.0, 1.0].into()),
  ///   bound:      component::Bound::from (shape::Capsule::noisy (2.0, 3.0)),
  ///   material:   component::MATERIAL_STONE,
  ///   collidable: true
  /// };
  /// let b = object::Static {
  ///   position:   component::Position ([0.0, 0.0, -1.0].into()),
  ///   bound:      component::Bound::from (shape::Capsule::noisy (1.0, 2.0)),
  ///   material:   component::MATERIAL_STONE,
  ///   collidable: true
  /// };
  /// assert_eq!(
  ///   Proximity::query (&a, &b),
  ///   Proximity {
  ///     distance:  -3.0,
  ///     half_axis: [ 1.5, 0.0,  0.0].into(),
  ///     normal:    math::Unit3::axis_x(),
  ///     midpoint:  [-0.5, 0.0, -0.5].into()
  ///   }
  /// );
  /// # }
  /// ```
  pub fn query <A, B> (object_a : &A, object_b : &B) -> Self where
    A : object::Bounded, B : object::Bounded
  {
    let component::Position (position_a) = object_a.position();
    let component::Position (position_b) = object_b.position();
    let component::Bound    (shape_a)    = object_a.bound();
    let component::Bound    (shape_b)    = object_b.bound();

    match (shape_a, shape_b) {
      //
      // sphere vs. sphere
      //
      ( shape::Variant::Bounded (shape::Bounded::Sphere (sphere_a)),
        shape::Variant::Bounded (shape::Bounded::Sphere (sphere_b))
      ) => Self::query_sphere_sphere (position_a, sphere_a, position_b, sphere_b),
      //
      // capsule vs. capsule
      //
      ( shape::Variant::Bounded (shape::Bounded::Capsule (capsule_a)),
        shape::Variant::Bounded (shape::Bounded::Capsule (capsule_b))
      ) => Self::query_capsule_capsule (position_a, capsule_a, position_b, capsule_b),
      //
      // cuboid vs. cuboid
      //
      ( shape::Variant::Bounded (shape::Bounded::Cuboid (cuboid_a)),
        shape::Variant::Bounded (shape::Bounded::Cuboid (cuboid_b))
      ) => Self::query_cuboid_cuboid (position_a, cuboid_a, position_b, cuboid_b),
      //
      // sphere vs. capsule
      //
      ( shape::Variant::Bounded (shape::Bounded::Sphere (sphere_a)),
        shape::Variant::Bounded (shape::Bounded::Capsule (capsule_b))
      ) => Self::query_sphere_capsule (position_a, sphere_a, position_b, capsule_b),
      ( shape::Variant::Bounded (shape::Bounded::Capsule (capsule_a)),
        shape::Variant::Bounded (shape::Bounded::Sphere (sphere_b))
      ) => Self::query_sphere_capsule (position_b, sphere_b, position_a, capsule_a)
        .flip(),
      //
      // sphere vs. cuboid
      //
      ( shape::Variant::Bounded (shape::Bounded::Sphere (sphere_a)),
        shape::Variant::Bounded (shape::Bounded::Cuboid (cuboid_b))
      ) => Self::query_sphere_cuboid (
        position_a, sphere_a, position_b, cuboid_b),
      ( shape::Variant::Bounded (shape::Bounded::Cuboid (cuboid_a)),
        shape::Variant::Bounded (shape::Bounded::Sphere (sphere_b))
      ) => Self::query_sphere_cuboid (position_b, sphere_b, position_a, cuboid_a).flip(),
      //
      // sphere vs. orthant
      //
      ( shape::Variant::Bounded   (shape::Bounded::Sphere    (sphere_a)),
        shape::Variant::Unbounded (shape::Unbounded::Orthant (orthant_b))
      ) => Self::query_sphere_orthant (position_a, sphere_a, position_b, orthant_b),
      ( shape::Variant::Unbounded (shape::Unbounded::Orthant (orthant_a)),
        shape::Variant::Bounded   (shape::Bounded::Sphere    (sphere_b))
      ) => Self::query_sphere_orthant (position_b, sphere_b, position_a, orthant_a)
        .flip(),
      //
      // capsule vs. cuboid
      //
      ( shape::Variant::Bounded (shape::Bounded::Capsule (capsule_a)),
        shape::Variant::Bounded (shape::Bounded::Cuboid  (cuboid_b))
      ) => Self::query_capsule_cuboid (position_a, capsule_a, position_b, cuboid_b),
      ( shape::Variant::Bounded (shape::Bounded::Cuboid  (cuboid_a)),
        shape::Variant::Bounded (shape::Bounded::Capsule (capsule_b))
      ) => Self::query_capsule_cuboid (position_b, capsule_b, position_a, cuboid_a)
        .flip(),
      //
      // capsule vs. orthant
      //
      ( shape::Variant::Bounded   (shape::Bounded::Capsule   (capsule_a)),
        shape::Variant::Unbounded (shape::Unbounded::Orthant (orthant_b))
      ) => Self::query_capsule_orthant (position_a, capsule_a, position_b, orthant_b),
      ( shape::Variant::Unbounded (shape::Unbounded::Orthant (orthant_a)),
        shape::Variant::Bounded   (shape::Bounded::Capsule   (capsule_b))
      ) => Self::query_capsule_orthant (position_b, capsule_b, position_a, orthant_a)
        .flip(),
      //
      // cuboid vs. orthant
      //
      ( shape::Variant::Bounded   (shape::Bounded::Cuboid    (cuboid_a)),
        shape::Variant::Unbounded (shape::Unbounded::Orthant (orthant_b))
      ) => Self::query_cuboid_orthant (position_a, cuboid_a, position_b, orthant_b),
      ( shape::Variant::Unbounded (shape::Unbounded::Orthant (orthant_a)),
        shape::Variant::Bounded   (shape::Bounded::Cuboid    (cuboid_b))
      ) => Self::query_cuboid_orthant (position_b, cuboid_b, position_a, orthant_a)
        .flip(),
      //
      // hull vs. any
      //
      (shape::Variant::Bounded (shape::Bounded::Hull (_)), _) =>
        Self::query_gjk (object_a, object_b),
      (_, shape::Variant::Bounded (shape::Bounded::Hull (_))) =>
        Self::query_gjk (object_b, object_a).flip(),
      _ => unimplemented!("proximity query not implemented for (A, B): {:?}",
        (shape_a, shape_b))
    }
  }

  pub fn query_sphere_sphere (
    position_a : &Point3 <f64>, sphere_a : &shape::Sphere <f64>,
    position_b : &Point3 <f64>, sphere_b : &shape::Sphere <f64>
  ) -> Self {
    use num::Zero;
    let (radius_a, radius_b) = (*sphere_a.radius, *sphere_b.radius);
    let distance_normal_half_axis = |axis : &Vector3 <f64>| {
      let axis_distance = axis.magnitude();
      let distance      = axis_distance - radius_a - radius_b;
      let normal        = if axis.is_zero() {
        Unit3::axis_x()
      } else {
        Unit3::unchecked_approx (-(axis / axis_distance))
      };
      let half_axis     = -0.5 * distance * *normal;
      (distance, normal, half_axis)
    };
    let axis = position_b - position_a;
    let (distance, normal, half_axis) = distance_normal_half_axis (&axis);
    let midpoint = position_a - radius_a * *normal + half_axis;
    Proximity { distance, half_axis, midpoint, normal }
  }

  pub fn query_capsule_capsule (
    position_a : &Point3 <f64>, capsule_a : &shape::Capsule <f64>,
    position_b : &Point3 <f64>, capsule_b : &shape::Capsule <f64>
  ) -> Self {
    use num::Zero;
    let (half_height_a, half_height_b) =
      (*capsule_a.half_height, *capsule_b.half_height);
    let (radius_a, radius_b) = (*capsule_a.radius, *capsule_b.radius);
    let shaft_max_a = position_a + vector3 (0.0, 0.0, half_height_a);
    let shaft_min_a = position_a - vector3 (0.0, 0.0, half_height_a);
    let shaft_max_b = position_b + vector3 (0.0, 0.0, half_height_b);
    let shaft_min_b = position_b - vector3 (0.0, 0.0, half_height_b);
    let distance_normal_half_axis = |axis : &Vector3 <f64>| {
      let axis_distance = axis.magnitude();
      let distance      = axis_distance - radius_a - radius_b;
      let normal        = if axis.is_zero() {
        Unit3::axis_x()
      } else {
        Unit3::unchecked_approx (-(axis / axis_distance))
      };
      let half_axis     = -0.5 * distance * *normal;
      (distance, normal, half_axis)
    };
    if shaft_max_a.0.z <= shaft_min_b.0.z {
      // cylinder shaft of B is above cylinder shaft of A
      let axis     = shaft_min_b - shaft_max_a;
      let (distance, normal, half_axis) = distance_normal_half_axis (&axis);
      let midpoint = shaft_max_a - radius_a * *normal + half_axis;
      Proximity { distance, half_axis, midpoint, normal }
    } else if shaft_max_b.0.z <= shaft_min_a.0.z {
      // cylinder shaft of B is below cylinder shaft of A
      let axis     = shaft_max_b.0 - shaft_min_a.0;
      let (distance, normal, half_axis) = distance_normal_half_axis (&axis);
      let midpoint = shaft_min_a - radius_a * *normal + half_axis;
      Proximity { distance, half_axis, midpoint, normal }
    } else {
      // cylinder shafts are adjacent
      let axis =
        Point3::from ([position_b.0.x, position_b.0.y, 0.0]) -
        Point3::from ([position_a.0.x, position_a.0.y, 0.0]);
      let (distance, normal, half_axis) = distance_normal_half_axis (&axis);
      let midpoint = {
        let max_z  = f64::min (shaft_max_a.0.z, shaft_max_b.0.z);
        let min_z  = f64::max (shaft_min_a.0.z, shaft_min_b.0.z);
        let mid_z  = 0.5 * (max_z + min_z);
        let mid_a  = position_a - radius_a * *normal + half_axis;
        [mid_a.0.x, mid_a.0.y, mid_z].into()
      };
      Proximity { distance, half_axis, midpoint, normal }
    }
  }

  pub fn query_cuboid_cuboid (
    position_a : &Point3 <f64>, cuboid_a : &shape::Cuboid <f64>,
    position_b : &Point3 <f64>, cuboid_b : &shape::Cuboid <f64>
  ) -> Self {
    use VectorSpace;
    use num::Zero;
    fn nearest_ab_to_distance (
      nearest_a      : &Point3  <f64>,
      nearest_b      : &Point3  <f64>,
      default_normal : &Unit3 <f64>
    ) -> Proximity {
      let axis      = nearest_b - nearest_a;
      let distance  = axis.magnitude();
      let half_axis = 0.5 * (nearest_b - nearest_a);
      let midpoint  = nearest_a + half_axis;
      let normal    = if axis.is_zero() {
        *default_normal
      } else {
        Unit3::unchecked_approx (-axis / distance)
      };
      Proximity { distance, half_axis, midpoint, normal }
    }

    let max_a = position_a + cuboid_a.max().0;
    let min_a = position_a + cuboid_a.min().0;
    let max_b = position_b + cuboid_b.max().0;
    let min_b = position_b + cuboid_b.min().0;
    let overlap_x = min_b.0.x < max_a.0.x && min_a.0.x < max_b.0.x;
    let overlap_y = min_b.0.y < max_a.0.y && min_a.0.y < max_b.0.y;
    let overlap_z = min_b.0.z < max_a.0.z && min_a.0.z < max_b.0.z;
    match (overlap_x, overlap_y, overlap_z) {
      (false, false, false) => {
        // no overlaps: nearest features are corners
        let a_to_b = position_b - position_a;
        debug_assert!(!a_to_b.x.is_zero());
        debug_assert!(!a_to_b.y.is_zero());
        debug_assert!(!a_to_b.z.is_zero());
        let a_to_b_sigvec = a_to_b.map (f64::signum);
        let nearest_a = position_a + cuboid_a.max().0 * a_to_b_sigvec;
        let nearest_b = position_b + cuboid_b.max().0 * (-a_to_b_sigvec);
        nearest_ab_to_distance (&nearest_a, &nearest_b,
          &Unit3::normalize (-a_to_b_sigvec))
      }
      (true, false, false)  => {
        // overlap x only: nearest features are edges
        let max_x  = f64::min (max_a.0.x, max_b.0.x);
        let min_x  = f64::max (min_a.0.x, min_b.0.x);
        debug_assert_ne!(min_x, max_x);
        let mid_x  = 0.5 * (min_x + max_x);
        let a_to_b = position_b - position_a;
        let a_to_b_yz = vector3 (0.0, a_to_b.y, a_to_b.z);
        let a_to_b_yz_sigvec = Vector3::sigvec (a_to_b_yz);
        debug_assert!(a_to_b_yz_sigvec.x.is_zero());
        debug_assert!(!a_to_b_yz_sigvec.y.is_zero());
        debug_assert!(!a_to_b_yz_sigvec.z.is_zero());
        let nearest_a = {
          let y = *cuboid_a.extents.y * a_to_b_yz_sigvec.y;
          let z = *cuboid_a.extents.z * a_to_b_yz_sigvec.z;
          [mid_x, position_a.0.y + y, position_a.0.z + z].into()
        };
        let nearest_b = {
          let y = -*cuboid_b.extents.y * a_to_b_yz_sigvec.y;
          let z = -*cuboid_b.extents.z * a_to_b_yz_sigvec.z;
          [mid_x, position_b.0.y + y, position_b.0.z + z].into()
        };
        nearest_ab_to_distance (&nearest_a, &nearest_b,
          &Unit3::normalize (-a_to_b_yz_sigvec))
      }
      (false, true, false)  => {
        // overlap y only: nearest features are edges
        let max_y  = f64::min (max_a.0.y, max_b.0.y);
        let min_y  = f64::max (min_a.0.y, min_b.0.y);
        debug_assert_ne!(min_y, max_y);
        let mid_y  = 0.5 * (min_y + max_y);
        let a_to_b = position_b - position_a;
        let a_to_b_xz = vector3 (a_to_b.x, 0.0, a_to_b.z);
        let a_to_b_xz_sigvec = Vector3::sigvec (a_to_b_xz);
        debug_assert!(!a_to_b_xz_sigvec.x.is_zero());
        debug_assert!(a_to_b_xz_sigvec.y.is_zero());
        debug_assert!(!a_to_b_xz_sigvec.z.is_zero());
        let nearest_a = {
          let x = *cuboid_a.extents.x * a_to_b_xz_sigvec.x;
          let z = *cuboid_a.extents.z * a_to_b_xz_sigvec.z;
          [position_a.0.x + x, mid_y, position_a.0.z + z].into()
        };
        let nearest_b = {
          let x = -*cuboid_b.extents.x * a_to_b_xz_sigvec.x;
          let z = -*cuboid_b.extents.z * a_to_b_xz_sigvec.z;
          [position_b.0.x + x, mid_y, position_b.0.z + z].into()
        };
        nearest_ab_to_distance (&nearest_a, &nearest_b,
          &Unit3::normalize (-a_to_b_xz_sigvec))
      }
      (false, false, true)  => {
        // overlap z only: nearest features are edges
        let max_z  = f64::min (max_a.0.z, max_b.0.z);
        let min_z  = f64::max (min_a.0.z, min_b.0.z);
        debug_assert_ne!(min_z, max_z);
        let mid_z  = 0.5 * (min_z + max_z);
        let a_to_b = position_b - position_a;
        let a_to_b_xy = vector3 (a_to_b.x, a_to_b.y, 0.0);
        let a_to_b_xy_sigvec = Vector3::sigvec (a_to_b_xy);
        debug_assert!(!a_to_b_xy_sigvec.x.is_zero());
        debug_assert!(!a_to_b_xy_sigvec.y.is_zero());
        debug_assert!(a_to_b_xy_sigvec.z.is_zero());
        let nearest_a = {
          let x = *cuboid_a.extents.x * a_to_b_xy_sigvec.x;
          let y = *cuboid_a.extents.y * a_to_b_xy_sigvec.y;
          [position_a.0.x + x, position_a.0.y + y, mid_z].into()
        };
        let nearest_b = {
          let x = -*cuboid_b.extents.x * a_to_b_xy_sigvec.x;
          let y = -*cuboid_b.extents.y * a_to_b_xy_sigvec.y;
          [position_b.0.x + x, position_b.0.y + y, mid_z].into()
        };
        nearest_ab_to_distance (&nearest_a, &nearest_b,
          &Unit3::normalize (-a_to_b_xy_sigvec))
      }
      (true, true, false)   => {
        // overlap x and y: nearest features are faces
        let max_x    = f64::min (max_a.0.x, max_b.0.x);
        let min_x    = f64::max (min_a.0.x, min_b.0.x);
        let mid_x    = 0.5 * (min_x + max_x);
        let max_y    = f64::min (max_a.0.y, max_b.0.y);
        let min_y    = f64::max (min_a.0.y, min_b.0.y);
        let mid_y    = 0.5 * (min_y + max_y);
        let a_to_b_z = position_b.0.z - position_a.0.z;
        let a_to_b_signum_z = a_to_b_z.signum();
        debug_assert_ne!(a_to_b_signum_z, 0.0);
        let nearest_a = point3 (
          mid_x, mid_y,
          //position_a.0.z + a_to_b_signum_z * *cuboid_a.extents.z
          a_to_b_signum_z.mul_add (*cuboid_a.extents.z, position_a.0.z)
        );
        let nearest_b = point3 (
          mid_x, mid_y,
          //position_b.0.z - a_to_b_signum_z * *cuboid_b.extents.z
          a_to_b_signum_z.mul_add (-*cuboid_b.extents.z, position_b.0.z)
        );
        let nearest_a_to_b_z = nearest_b.0.z - nearest_a.0.z;
        let distance  = nearest_a_to_b_z.abs();
        let half_axis = [0.0, 0.0, 0.5 * nearest_a_to_b_z].into();
        let normal    = Unit3::unchecked ([0.0, 0.0, -a_to_b_signum_z].into());
        let midpoint  = nearest_a + half_axis;
        Proximity { distance, half_axis, midpoint, normal }
      }
      (true, false, true)   => {
        // overlap x and z: nearest features are faces
        let max_x    = f64::min (max_a.0.x, max_b.0.x);
        let min_x    = f64::max (min_a.0.x, min_b.0.x);
        let mid_x    = 0.5 * (min_x + max_x);
        let max_z    = f64::min (max_a.0.z, max_b.0.z);
        let min_z    = f64::max (min_a.0.z, min_b.0.z);
        let mid_z    = 0.5 * (min_z + max_z);
        let a_to_b_y = position_b.0.y - position_a.0.y;
        let a_to_b_signum_y = a_to_b_y.signum();
        debug_assert_ne!(a_to_b_signum_y, 0.0);
        let nearest_a = point3 (
          mid_x,
          //position_a.0.y + a_to_b_signum_y * *cuboid_a.extents.y,
          a_to_b_signum_y.mul_add (*cuboid_a.extents.y, position_a.0.y),
          mid_z
        );
        let nearest_b = point3 (
          mid_x,
          //position_b.0.y - a_to_b_signum_y * *cuboid_b.extents.y,
          a_to_b_signum_y.mul_add (-*cuboid_b.extents.y, position_b.0.y),
          mid_z
        );
        let nearest_a_to_b_y = nearest_b.0.y - nearest_a.0.y;
        let distance  = nearest_a_to_b_y.abs();
        let half_axis = [0.0, 0.5 * nearest_a_to_b_y, 0.0].into();
        let normal    = Unit3::unchecked ([0.0, -a_to_b_signum_y, 0.0].into());
        let midpoint  = nearest_a + half_axis;
        Proximity { distance, half_axis, midpoint, normal }
      }
      (false, true, true)   => {
        // overlap y and z: nearest features are faces
        let max_y    = f64::min (max_a.0.y, max_b.0.y);
        let min_y    = f64::max (min_a.0.y, min_b.0.y);
        let mid_y    = 0.5 * (min_y + max_y);
        let max_z    = f64::min (max_a.0.z, max_b.0.z);
        let min_z    = f64::max (min_a.0.z, min_b.0.z);
        let mid_z    = 0.5 * (min_z + max_z);
        let a_to_b_x = position_b.0.x - position_a.0.x;
        let a_to_b_signum_x = a_to_b_x.signum();
        debug_assert_ne!(a_to_b_signum_x, 0.0);
        let nearest_a = point3 (
          //position_a.0.x + a_to_b_signum_x * *cuboid_a.extents.x,
          a_to_b_signum_x.mul_add (*cuboid_a.extents.x, position_a.0.x),
          mid_y,
          mid_z
        );
        let nearest_b = point3 (
          //position_b.0.x - a_to_b_signum_x * *cuboid_b.extents.x,
          a_to_b_signum_x.mul_add (-*cuboid_b.extents.x, position_b.0.x),
          mid_y,
          mid_z
        );
        let nearest_a_to_b_x = nearest_b.0.x - nearest_a.0.x;
        let distance  = nearest_a_to_b_x.abs();
        let half_axis = [0.5 * nearest_a_to_b_x, 0.0, 0.0].into();
        let normal    = Unit3::unchecked ([-a_to_b_signum_x, 0.0, 0.0].into());
        let midpoint  = nearest_a + half_axis;
        Proximity { distance, half_axis, midpoint, normal }
      }
      (true, true, true)    => {
        // overlap on all axes: intersection
        let max_x     = f64::min (max_a.0.x, max_b.0.x);
        let min_x     = f64::max (min_a.0.x, min_b.0.x);
        let mid_x     = 0.5 * (min_x + max_x);
        let max_y     = f64::min (max_a.0.y, max_b.0.y);
        let min_y     = f64::max (min_a.0.y, min_b.0.y);
        let mid_y     = 0.5 * (min_y + max_y);
        let max_z     = f64::min (max_a.0.z, max_b.0.z);
        let min_z     = f64::max (min_a.0.z, min_b.0.z);
        let mid_z     = 0.5 * (min_z + max_z);
        let midpoint  = [mid_x, mid_y, mid_z].into();
        let resolve_x = if position_a.0.x < position_b.0.x {
          min_b.0.x - max_a.0.x
        } else {
          debug_assert!(position_a.0.x >= position_b.0.x);
          max_b.0.x - min_a.0.x
        };
        let resolve_y = if position_a.0.y < position_b.0.y {
          min_b.0.y - max_a.0.y
        } else {
          debug_assert!(position_a.0.y >= position_b.0.y);
          max_b.0.y - min_a.0.y
        };
        let resolve_z = if position_a.0.z < position_b.0.z {
          min_b.0.z - max_a.0.z
        } else {
          debug_assert!(position_a.0.z >= position_b.0.z);
          max_b.0.z - min_a.0.z
        };
        let resolve_x_abs = resolve_x.abs();
        let resolve_y_abs = resolve_y.abs();
        let resolve_z_abs = resolve_z.abs();
        let (distance, half_axis, normal) : (
          f64, Vector3 <f64>, Unit3 <f64>
        ) = if resolve_x_abs <= resolve_y_abs && resolve_x_abs <= resolve_z_abs {
          ( -resolve_x_abs,
            0.5 * vector3 (resolve_x, 0.0, 0.0),
            Unit3::unchecked ([resolve_x.signum(), 0.0, 0.0].into())
          )
        } else if resolve_y_abs <= resolve_x_abs && resolve_y_abs <= resolve_z_abs {
          ( -resolve_y_abs,
            0.5 * vector3 (0.0, resolve_y, 0.0),
            Unit3::unchecked ([0.0, resolve_y.signum(), 0.0].into())
          )
        } else {
          debug_assert!(resolve_z_abs <= resolve_x_abs);
          debug_assert!(resolve_z_abs <= resolve_y_abs);
          ( -resolve_z_abs,
            0.5 * vector3 (0.0, 0.0, resolve_z),
            Unit3::unchecked ([0.0, 0.0, resolve_z.signum()].into())
          )
        };
        debug_assert!(!half_axis.is_zero());
        Proximity { distance, half_axis, midpoint, normal }
      }
    }
  }

  /*
  /// Uses pre-computed sigvec normals with array lookup.
  ///
  /// This benchmarked 2-3ns slower than the version using `normalize`.
  pub fn query_cuboid_cuboid_array (
    position_a : &Point3 <f64>, cuboid_a : &shape::Cuboid <f64>,
    position_b : &Point3 <f64>, cuboid_b : &shape::Cuboid <f64>
  ) -> Self {
    use {EuclideanSpace, InnerSpace, Zero};

    fn nearest_ab_to_distance (
      nearest_a      : &Point3  <f64>,
      nearest_b      : &Point3  <f64>,
      default_normal : &Vector3 <f64>
    ) -> Proximity {
      let axis      = nearest_b - nearest_a;
      let distance  = axis.magnitude();
      let half_axis = 0.5 * (nearest_b - nearest_a);
      let midpoint  = nearest_a + half_axis;
      let normal    = if axis.is_zero() {
        *default_normal
      } else {
        -axis / distance
      };
      Proximity { distance, half_axis, midpoint, normal }
    }

    let max_a = position_a + cuboid_a.max().0;
    let min_a = position_a + cuboid_a.min().0;
    let max_b = position_b + cuboid_b.max().0;
    let min_b = position_b + cuboid_b.min().0;
    let overlap_x = min_b.0.x < max_a.0.x && min_a.0.x < max_b.0.x;
    let overlap_y = min_b.0.y < max_a.0.y && min_a.0.y < max_b.0.y;
    let overlap_z = min_b.0.z < max_a.0.z && min_a.0.z < max_b.0.z;
    match (overlap_x, overlap_y, overlap_z) {
      (false, false, false) => {
        use ElementWise;
        // no overlaps: nearest features are corners
        let a_to_b = position_b - position_a;
        debug_assert!(!a_to_b.x.is_zero());
        debug_assert!(!a_to_b.y.is_zero());
        debug_assert!(!a_to_b.z.is_zero());
        let a_to_b_sigvec = a_to_b.map (f64::signum);
        let nearest_a = position_a +
          cuboid_a.half_extents().mul_element_wise (a_to_b_sigvec);
        let nearest_b = position_b +
          cuboid_b.half_extents().mul_element_wise (-a_to_b_sigvec);
        nearest_ab_to_distance (&nearest_a, &nearest_b,
          &-sigvec_unit_f64_unchecked (a_to_b_sigvec))
      }
      (true, false, false)  => {
        // overlap x only: nearest features are edges
        let max_x  = f64::min (max_a.0.x, max_b.0.x);
        let min_x  = f64::max (min_a.0.x, min_b.0.x);
        debug_assert_ne!(min_x, max_x);
        let mid_x  = 0.5 * (min_x + max_x);
        let a_to_b = position_b - position_a;
        let a_to_b_yz = vector3 (0.0, a_to_b.y, a_to_b.z);
        let a_to_b_yz_sigvec = sigvec (a_to_b_yz);
        debug_assert!(a_to_b_yz_sigvec.x.is_zero());
        debug_assert!(!a_to_b_yz_sigvec.y.is_zero());
        debug_assert!(!a_to_b_yz_sigvec.z.is_zero());
        let nearest_a = {
          let y = cuboid_a.extents.y * a_to_b_yz_sigvec.y;
          let z = cuboid_a.extents.z * a_to_b_yz_sigvec.z;
          [mid_x, position_a.0.y + y, position_a.0.z + z].into()
        };
        let nearest_b = {
          let y = -cuboid_b.extents.y * a_to_b_yz_sigvec.y;
          let z = -cuboid_b.extents.z * a_to_b_yz_sigvec.z;
          [mid_x, position_b.0.y + y, position_b.0.z + z].into()
        };
        nearest_ab_to_distance (&nearest_a, &nearest_b,
          &-sigvec_unit_f64_unchecked (a_to_b_yz_sigvec))
      }
      (false, true, false)  => {
        // overlap y only: nearest features are edges
        let max_y  = f64::min (max_a.0.y, max_b.0.y);
        let min_y  = f64::max (min_a.0.y, min_b.0.y);
        debug_assert_ne!(min_y, max_y);
        let mid_y  = 0.5 * (min_y + max_y);
        let a_to_b = position_b - position_a;
        let a_to_b_xz = vector3 (a_to_b.x, 0.0, a_to_b.z);
        let a_to_b_xz_sigvec = sigvec (a_to_b_xz);
        debug_assert!(!a_to_b_xz_sigvec.x.is_zero());
        debug_assert!(a_to_b_xz_sigvec.y.is_zero());
        debug_assert!(!a_to_b_xz_sigvec.z.is_zero());
        let nearest_a = {
          let x = cuboid_a.extents.x * a_to_b_xz_sigvec.x;
          let z = cuboid_a.extents.z * a_to_b_xz_sigvec.z;
          [position_a.0.x + x, mid_y, position_a.0.z + z].into()
        };
        let nearest_b = {
          let x = -cuboid_b.extents.x * a_to_b_xz_sigvec.x;
          let z = -cuboid_b.extents.z * a_to_b_xz_sigvec.z;
          [position_b.0.x + x, mid_y, position_b.0.z + z].into()
        };
        nearest_ab_to_distance (&nearest_a, &nearest_b,
          &-sigvec_unit_f64_unchecked (a_to_b_xz_sigvec))
      }
      (false, false, true)  => {
        // overlap z only: nearest features are edges
        let max_z  = f64::min (max_a.0.z, max_b.0.z);
        let min_z  = f64::max (min_a.0.z, min_b.0.z);
        debug_assert_ne!(min_z, max_z);
        let mid_z  = 0.5 * (min_z + max_z);
        let a_to_b = position_b - position_a;
        let a_to_b_xy = vector3 (a_to_b.x, a_to_b.y, 0.0);
        let a_to_b_xy_sigvec = sigvec (a_to_b_xy);
        debug_assert!(!a_to_b_xy_sigvec.x.is_zero());
        debug_assert!(!a_to_b_xy_sigvec.y.is_zero());
        debug_assert!(a_to_b_xy_sigvec.z.is_zero());
        let nearest_a = {
          let x = cuboid_a.extents.x * a_to_b_xy_sigvec.x;
          let y = cuboid_a.extents.y * a_to_b_xy_sigvec.y;
          [position_a.0.x + x, position_a.0.y + y, mid_z].into()
        };
        let nearest_b = {
          let x = -cuboid_b.extents.x * a_to_b_xy_sigvec.x;
          let y = -cuboid_b.extents.y * a_to_b_xy_sigvec.y;
          [position_b.0.x + x, position_b.0.y + y, mid_z].into()
        };
        nearest_ab_to_distance (&nearest_a, &nearest_b,
          &-sigvec_unit_f64_unchecked (a_to_b_xy_sigvec))
      }
      (true, true, false)   => {
        // overlap x and y: nearest features are faces
        let max_x    = f64::min (max_a.0.x, max_b.0.x);
        let min_x    = f64::max (min_a.0.x, min_b.0.x);
        let mid_x    = 0.5 * (min_x + max_x);
        let max_y    = f64::min (max_a.0.y, max_b.0.y);
        let min_y    = f64::max (min_a.0.y, min_b.0.y);
        let mid_y    = 0.5 * (min_y + max_y);
        let a_to_b_z = position_b.0.z - position_a.0.z;
        let a_to_b_signum_z = a_to_b_z.signum();
        debug_assert_ne!(a_to_b_signum_z, 0.0);
        let nearest_a = point3 (
          mid_x, mid_y,
          position_a.0.z + a_to_b_signum_z * cuboid_a.extents.z
        );
        let nearest_b = point3 (
          mid_x, mid_y,
          position_b.0.z - a_to_b_signum_z * cuboid_b.extents.z
        );
        let nearest_a_to_b_z = nearest_b.0.z - nearest_a.0.z;
        let distance  = nearest_a_to_b_z.abs();
        let half_axis = [0.0, 0.0, 0.5 * nearest_a_to_b_z].into();
        let normal    = [0.0, 0.0, -a_to_b_signum_z].into();
        let midpoint  = nearest_a + half_axis;
        Proximity { distance, half_axis, midpoint, normal }
      }
      (true, false, true)   => {
        // overlap x and z: nearest features are faces
        let max_x    = f64::min (max_a.0.x, max_b.0.x);
        let min_x    = f64::max (min_a.0.x, min_b.0.x);
        let mid_x    = 0.5 * (min_x + max_x);
        let max_z    = f64::min (max_a.0.z, max_b.0.z);
        let min_z    = f64::max (min_a.0.z, min_b.0.z);
        let mid_z    = 0.5 * (min_z + max_z);
        let a_to_b_y = position_b.0.y - position_a.0.y;
        let a_to_b_signum_y = a_to_b_y.signum();
        debug_assert_ne!(a_to_b_signum_y, 0.0);
        let nearest_a = point3 (
          mid_x,
          position_a.0.y + a_to_b_signum_y * cuboid_a.extents.y,
          mid_z
        );
        let nearest_b = point3 (
          mid_x,
          position_b.0.y - a_to_b_signum_y * cuboid_b.extents.y,
          mid_z
        );
        let nearest_a_to_b_y = nearest_b.0.y - nearest_a.0.y;
        let distance  = nearest_a_to_b_y.abs();
        let half_axis = [0.0, 0.5 * nearest_a_to_b_y, 0.0].into();
        let normal    = [0.0, -a_to_b_signum_y, 0.0].into();
        let midpoint  = nearest_a + half_axis;
        Proximity { distance, half_axis, midpoint, normal }
      }
      (false, true, true)   => {
        // overlap y and z: nearest features are faces
        let max_y    = f64::min (max_a.0.y, max_b.0.y);
        let min_y    = f64::max (min_a.0.y, min_b.0.y);
        let mid_y    = 0.5 * (min_y + max_y);
        let max_z    = f64::min (max_a.0.z, max_b.0.z);
        let min_z    = f64::max (min_a.0.z, min_b.0.z);
        let mid_z    = 0.5 * (min_z + max_z);
        let a_to_b_x = position_b.0.x - position_a.0.x;
        let a_to_b_signum_x = a_to_b_x.signum();
        debug_assert_ne!(a_to_b_signum_x, 0.0);
        let nearest_a = point3 (
          position_a.0.x + a_to_b_signum_x * cuboid_a.extents.x,
          mid_y,
          mid_z
        );
        let nearest_b = point3 (
          position_b.0.x - a_to_b_signum_x * cuboid_b.extents.x,
          mid_y,
          mid_z
        );
        let nearest_a_to_b_x = nearest_b.0.x - nearest_a.0.x;
        let distance  = nearest_a_to_b_x.abs();
        let half_axis = [0.5 * nearest_a_to_b_x, 0.0, 0.0].into();
        let normal    = [-a_to_b_signum_x, 0.0, 0.0].into();
        let midpoint  = nearest_a + half_axis;
        Proximity { distance, half_axis, midpoint, normal }
      }
      (true, true, true)    => {
        // overlap on all axes: intersection
        let max_x     = f64::min (max_a.0.x, max_b.0.x);
        let min_x     = f64::max (min_a.0.x, min_b.0.x);
        let mid_x     = 0.5 * (min_x + max_x);
        let max_y     = f64::min (max_a.0.y, max_b.0.y);
        let min_y     = f64::max (min_a.0.y, min_b.0.y);
        let mid_y     = 0.5 * (min_y + max_y);
        let max_z     = f64::min (max_a.0.z, max_b.0.z);
        let min_z     = f64::max (min_a.0.z, min_b.0.z);
        let mid_z     = 0.5 * (min_z + max_z);
        let midpoint  = [mid_x, mid_y, mid_z].into();
        let resolve_x = if position_a.0.x < position_b.0.x {
          min_b.0.x - max_a.0.x
        } else {
          debug_assert!(position_a.0.x >= position_b.0.x);
          max_b.0.x - min_a.0.x
        };
        let resolve_y = if position_a.0.y < position_b.0.y {
          min_b.0.y - max_a.0.y
        } else {
          debug_assert!(position_a.0.y >= position_b.0.y);
          max_b.0.y - min_a.0.y
        };
        let resolve_z = if position_a.0.z < position_b.0.z {
          min_b.0.z - max_a.0.z
        } else {
          debug_assert!(position_a.0.z >= position_b.0.z);
          max_b.0.z - min_a.0.z
        };
        let resolve_x_abs = resolve_x.abs();
        let resolve_y_abs = resolve_y.abs();
        let resolve_z_abs = resolve_z.abs();
        let (distance, half_axis, normal) : (
          f64, Vector3 <f64>, Vector3 <f64>
        ) = if resolve_x_abs <= resolve_y_abs && resolve_x_abs <= resolve_z_abs {
          ( -resolve_x_abs,
            0.5 * vector3 (resolve_x, 0.0, 0.0),
            [resolve_x.signum(), 0.0, 0.0].into()
          )
        } else if
          resolve_y_abs <= resolve_x_abs && resolve_y_abs <= resolve_z_abs
        {
          ( -resolve_y_abs,
            0.5 * vector3 (0.0, resolve_y, 0.0),
            [0.0, resolve_y.signum(), 0.0].into()
          )
        } else {
          debug_assert!(resolve_z_abs <= resolve_x_abs);
          debug_assert!(resolve_z_abs <= resolve_y_abs);
          ( -resolve_z_abs,
            0.5 * vector3 (0.0, 0.0, resolve_z),
            [0.0, 0.0, resolve_z.signum()].into()
          )
        };
        debug_assert!(!half_axis.is_zero());
        Proximity { distance, half_axis, midpoint, normal }
      }
    }
  } // end query_cuboid_cuboid_array

  /// This is a version of the cuboid distance query with the three edge
  /// cases factored out into a function taking vector indices. This turns out
  /// to be slower in benchmarks (190ns vs. the 180ns inline versions above).
  pub fn query_cuboid_cuboid_refactor (
    position_a : &Point3 <f64>, cuboid_a : &shape::Cuboid <f64>,
    position_b : &Point3 <f64>, cuboid_b : &shape::Cuboid <f64>
  ) -> Self {
    use {EuclideanSpace, InnerSpace, Zero};

    #[inline]
    fn distance_edge (
      edge : usize, other1 : usize, other2 : usize,
      position_a : &Point3 <f64>, cuboid_a : &shape::Cuboid <f64>,
      position_b : &Point3 <f64>, cuboid_b : &shape::Cuboid <f64>,
      max_a      : &Point3 <f64>, min_a    : &Point3 <f64>,
      max_b      : &Point3 <f64>, min_b    : &Point3 <f64>
    ) -> Proximity {
      let max_edge  = f64::min (max_a[edge], max_b[edge]);
      let min_edge  = f64::max (min_a[edge], min_b[edge]);
      debug_assert_ne!(min_edge, max_edge);
      let mid_edge  = 0.5 * (min_edge + max_edge);
      let a_to_b = position_b - position_a;
      let a_to_b_others = {
        let mut a_to_b_others = a_to_b;
        a_to_b_others[edge] = 0.0;
        a_to_b_others
      };
      let a_to_b_others_signum = a_to_b_others.map (signum_or_zero);
      debug_assert!(a_to_b_others_signum[edge].is_zero());
      debug_assert!(!a_to_b_others_signum[other1].is_zero());
      debug_assert!(!a_to_b_others_signum[other2].is_zero());
      let nearest_a = {
        let other_component1 =
          cuboid_a.half_extents()[other1] * a_to_b_others_signum[other1];
        let other_component2 =
          cuboid_a.half_extents()[other2] * a_to_b_others_signum[other2];
        let mut nearest_a = [0.0; 3];
        nearest_a[edge]   = mid_edge;
        nearest_a[other1] = position_a[other1] + other_component1;
        nearest_a[other2] = position_a[other2] + other_component2;
        nearest_a.into()
      };
      let nearest_b = {
        let other_component1 =
          -cuboid_b.half_extents()[other1] * a_to_b_others_signum[other1];
        let other_component2 =
          -cuboid_b.half_extents()[other2] * a_to_b_others_signum[other2];
        let mut nearest_b = [0.0; 3];
        nearest_b[edge]   = mid_edge;
        nearest_b[other1] = position_b[other1] + other_component1;
        nearest_b[other2] = position_b[other2] + other_component2;
        nearest_b.into()
      };
      nearest_ab_to_distance (&nearest_a, &nearest_b, &a_to_b_others_signum)
    }

    fn nearest_ab_to_distance (
      nearest_a     : &Point3 <f64>,
      nearest_b     : &Point3 <f64>,
      a_to_b_signum : &Vector3 <f64>
    ) -> Proximity {
      let axis      = nearest_b - nearest_a;
      let distance  = axis.magnitude();
      let half_axis = 0.5 * (nearest_b - nearest_a);
      let midpoint  = nearest_a + half_axis;
      let normal    = if axis.is_zero() {
        -a_to_b_signum.normalize()
      } else {
        -axis / distance
      };
      Proximity { distance, half_axis, midpoint, normal }
    }

    let max_a = position_a + cuboid_a.max().0;
    let min_a = position_a + cuboid_a.min().0;
    let max_b = position_b + cuboid_b.max().0;
    let min_b = position_b + cuboid_b.min().0;

    let overlap_x = min_b.0.x < max_a.0.x && min_a.0.x < max_b.0.x;
    let overlap_y = min_b.0.y < max_a.0.y && min_a.0.y < max_b.0.y;
    let overlap_z = min_b.0.z < max_a.0.z && min_a.0.z < max_b.0.z;
    match (overlap_x, overlap_y, overlap_z) {
      (false, false, false) => {
        use ElementWise;
        // no overlaps: nearest features are corners
        let a_to_b = position_b - position_a;
        debug_assert!(!a_to_b.x.is_zero());
        debug_assert!(!a_to_b.y.is_zero());
        debug_assert!(!a_to_b.z.is_zero());
        let a_to_b_signum = a_to_b.map (f64::signum);
        let nearest_a = position_a +
          cuboid_a.half_extents().mul_element_wise (a_to_b_signum);
        let nearest_b = position_b +
          cuboid_b.half_extents().mul_element_wise (-a_to_b_signum);
        nearest_ab_to_distance (&nearest_a, &nearest_b, &a_to_b_signum)
      }
      (true, false, false)  => {
        // overlap x only: nearest features are edges
        distance_edge (0, 1, 2,
          position_a, cuboid_a, position_b, cuboid_b,
          &max_a, &min_a, &max_b, &min_b
        )
      }
      (false, true, false)  => {
        // overlap y only: nearest features are edges
        distance_edge (1, 0, 2,
          position_a, cuboid_a, position_b, cuboid_b,
          &max_a, &min_a, &max_b, &min_b
        )
      }
      (false, false, true)  => {
        // overlap z only: nearest features are edges
        distance_edge (2, 1, 2,
          position_a, cuboid_a, position_b, cuboid_b,
          &max_a, &min_a, &max_b, &min_b
        )
      }
      (true, true, false)   => {
        // overlap x and y: nearest features are faces
        unimplemented!()
      }
      (true, false, true)   => {
        // overlap x and z: nearest features are faces
        unimplemented!()
      }
      (false, true, true)   => {
        // overlap y and z: nearest features are faces
        unimplemented!()
      }
      (true, true, true)    => {
        // overlap on all axes: intersection
        unimplemented!()
      }
    }
  } // end query_cuboid_cuboid_refactor
  */

  pub fn query_sphere_capsule (
    position_a : &Point3 <f64>, sphere_a  : &shape::Sphere <f64>,
    position_b : &Point3 <f64>, capsule_b : &shape::Capsule <f64>
  ) -> Self {
    use num::Zero;
    let half_height_b = *capsule_b.half_height;
    let (radius_a, radius_b) = (*sphere_a.radius, *capsule_b.radius);
    let shaft_max_b = position_b + vector3 (0.0, 0.0, half_height_b);
    let shaft_min_b = position_b - vector3 (0.0, 0.0, half_height_b);
    let distance_normal_half_axis = |axis : &Vector3 <f64>| {
      let axis_distance = axis.magnitude();
      let distance      = axis_distance - radius_a - radius_b;
      let normal        = if axis.is_zero() {
        Unit3::axis_x()
      } else {
        Unit3::unchecked_approx (-(axis / axis_distance))
      };
      let half_axis     = -0.5 * distance * *normal;
      (distance, normal, half_axis)
    };
    if position_a.0.z <= shaft_min_b.0.z {
      // cylinder shaft of B is above center of A
      let axis     = shaft_min_b - position_a;
      let (distance, normal, half_axis) = distance_normal_half_axis (&axis);
      let midpoint = position_a - radius_a * *normal + half_axis;
      Proximity { distance, half_axis, midpoint, normal }
    } else if shaft_max_b.0.z <= position_a.0.z {
      // cylinder shaft of B is below center of A
      let axis     = shaft_max_b - position_a;
      let (distance, normal, half_axis) = distance_normal_half_axis (&axis);
      let midpoint = position_a - radius_a * *normal + half_axis;
      Proximity { distance, half_axis, midpoint, normal }
    } else {
      // cylinder shaft of B is adjacent to center of A
      let axis     =
        Point3::from ([position_b.0.x, position_b.0.y, 0.0]) -
        Point3::from ([position_a.0.x, position_a.0.y, 0.0]);
      let (distance, normal, half_axis) = distance_normal_half_axis (&axis);
      let midpoint = {
        let max_z  = f64::min (position_a.0.z, shaft_max_b.0.z);
        let min_z  = f64::max (position_a.0.z, shaft_min_b.0.z);
        let mid_z  = 0.5 * (max_z + min_z);
        let mid_a  = position_a - radius_a * *normal + half_axis;
        [mid_a.0.x, mid_a.0.y, mid_z].into()
      };
      Proximity { distance, half_axis, midpoint, normal }
    }
  }

  pub fn query_sphere_cuboid (
    position_a : &Point3 <f64>, sphere_a : &shape::Sphere <f64>,
    position_b : &Point3 <f64>, cuboid_b : &shape::Cuboid <f64>
  ) -> Self {
    use num::Zero;
    let cuboid_max_b = position_b + cuboid_b.max().0;
    let cuboid_min_b = position_b - cuboid_b.max().0;

    let overlap_x = cuboid_min_b.0.x < position_a.0.x &&
      position_a.0.x < cuboid_max_b.0.x;
    let overlap_y = cuboid_min_b.0.y < position_a.0.y &&
      position_a.0.y < cuboid_max_b.0.y;

    let sphere_a_radius   = *sphere_a.radius;
    let cuboid_b_extent_x = *cuboid_b.extents.x;
    let cuboid_b_extent_y = *cuboid_b.extents.y;
    let cuboid_b_extent_z = *cuboid_b.extents.z;
    if position_a.0.z >= cuboid_max_b.0.z {
      // sphere above cuboid
      match (overlap_x, overlap_y) {
        (false, false) => {
          // nearest/deepest point is a corner of the cuboid
          let b_to_a        = position_a - position_b;
          debug_assert!(!b_to_a.x.is_zero());
          debug_assert!(!b_to_a.y.is_zero());
          debug_assert!(!b_to_a.z.is_zero());
          let b_to_a_sigvec = b_to_a.map (f64::signum);
          let nearest_b     = position_b + cuboid_b.max().0 * b_to_a_sigvec;
          let a_to_b        = nearest_b - position_a;
          let nearest_a     = position_a + if !a_to_b.is_zero() {
            sphere_a_radius * a_to_b.normalized()
          } else {
            [-b_to_a_sigvec.x * sphere_a_radius, 0.0, 0.0].into()
          };
          let axis          = nearest_b - nearest_a;
          let distance      = axis.magnitude() * if !a_to_b.is_zero() {
            a_to_b.dot (axis).signum()
          } else {
            -1.0
          };
          let half_axis     = 0.5 * axis;
          let normal        = if axis.is_zero() {
            Unit3::normalize (b_to_a_sigvec)
          } else {
            Unit3::unchecked_approx (-axis / distance)
          };
          let midpoint      = nearest_a + half_axis;
          Proximity { distance, half_axis, midpoint, normal }
        }
        (true, false) => {
          // nearest/deepest point is on an edge parallel to the X axis
          let signum_y  = (position_a.0.y - position_b.0.y).signum();
          let nearest_b = point3 (
            position_a.0.x,
            //position_b.0.y + signum_y * cuboid_b_extent_y,
            signum_y.mul_add (cuboid_b_extent_y, position_b.0.y),
            cuboid_max_b.0.z
          );
          let a_to_b    = nearest_b - position_a;
          let nearest_a = position_a + if a_to_b.is_zero() {
            -signum_y * vector3 (0.0, sphere_a_radius, 0.0)
          } else {
            sphere_a_radius * a_to_b.normalized()
          };
          let axis      = nearest_b - nearest_a;
          let distance  = axis.magnitude() * if !a_to_b.is_zero() {
            a_to_b.dot (axis).signum()
          } else {
            -1.0
          };
          let half_axis = 0.5 * axis;
          let normal    = if axis.is_zero() {
            Unit3::normalize (vector3 (0.0, signum_y, 1.0))
          } else {
            Unit3::unchecked_approx (-axis / distance)
          };
          let midpoint  = nearest_a + half_axis;
          Proximity { distance, half_axis, midpoint, normal }
        }
        (false, true) => {
          // nearest/deepest point is on an edge parallel to the Y axis
          let signum_x  = (position_a.0.x - position_b.0.x).signum();
          let nearest_b = point3 (
            //position_b.0.x + signum_x * cuboid_b_extent_x,
            signum_x.mul_add (cuboid_b_extent_x, position_b.0.x),
            position_a.0.y,
            cuboid_max_b.0.z
          );
          let a_to_b    = nearest_b - position_a;
          let nearest_a = position_a + if a_to_b.is_zero() {
            -signum_x * vector3 (sphere_a_radius, 0.0, 0.0)
          } else {
            sphere_a_radius * a_to_b.normalized()
          };
          let axis      = nearest_b - nearest_a;
          let distance  = axis.magnitude() * if !a_to_b.is_zero() {
            a_to_b.dot (axis).signum()
          } else {
            -1.0
          };
          let half_axis = 0.5 * axis;
          let normal    = if axis.is_zero() {
            Unit3::normalize (vector3 (signum_x, 0.0, 1.0))
          } else {
            Unit3::unchecked_approx (-axis / distance)
          };
          let midpoint  = nearest_a + half_axis;
          Proximity { distance, half_axis, midpoint, normal }
        }
        (true, true) => {
          // nearest/deepest point is in the -Z direction
          let nearest_a : Point3 <f64> =
            position_a - vector3 (0.0, 0.0, sphere_a_radius);
          let nearest_b : Point3 <f64> =
            [position_a.0.x, position_a.0.y, cuboid_max_b.0.z].into();
          let axis      = nearest_b - nearest_a;
          let normal    = Unit3::axis_z();
          let half_axis = 0.5 * axis;
          let distance  = nearest_a.0.z - nearest_b.0.z;
          let midpoint  = nearest_a + half_axis;
          Proximity { distance, half_axis, midpoint, normal }
        }
      }
    } else if position_a.0.z <= cuboid_min_b.0.z {
      // sphere below cuboid
      match (overlap_x, overlap_y) {
        (false, false) => {
          // nearest/deepest point is a corner of the cuboid
          let b_to_a        = position_a - position_b;
          debug_assert!(!b_to_a.x.is_zero());
          debug_assert!(!b_to_a.y.is_zero());
          debug_assert!(!b_to_a.z.is_zero());
          let b_to_a_sigvec = b_to_a.map (f64::signum);
          let nearest_b     = position_b + cuboid_b.max().0 * b_to_a_sigvec;
          let a_to_b        = nearest_b - position_a;
          let nearest_a     = position_a + if !a_to_b.is_zero() {
            sphere_a_radius * a_to_b.normalized()
          } else {
            [-b_to_a_sigvec.x * sphere_a_radius, 0.0, 0.0].into()
          };
          let axis          = nearest_b - nearest_a;
          let distance      = axis.magnitude() * if !a_to_b.is_zero() {
            a_to_b.dot (axis).signum()
          } else {
            -1.0
          };
          let half_axis     = 0.5 * axis;
          let normal        = if axis.is_zero() {
            Unit3::normalize (b_to_a_sigvec)
          } else {
            Unit3::unchecked_approx (-axis / distance)
          };
          let midpoint      = nearest_a + half_axis;
          Proximity { distance, half_axis, midpoint, normal }
        }
        (true, false) => {
          // nearest/deepest point is on an edge parallel to the X axis
          let signum_y  = (position_a.0.y - position_b.0.y).signum();
          let nearest_b = point3 (
            position_a.0.x,
            //position_b.0.y + signum_y * cuboid_b_extent_y,
            signum_y.mul_add (cuboid_b_extent_y, position_b.0.y),
            cuboid_min_b.0.z
          );
          let a_to_b    = nearest_b - position_a;
          let nearest_a = position_a + if a_to_b.is_zero() {
            -signum_y * vector3 (0.0, sphere_a_radius, 0.0)
          } else {
            sphere_a_radius * a_to_b.normalized()
          };
          let axis      = nearest_b - nearest_a;
          let distance  = axis.magnitude() * if !a_to_b.is_zero() {
            a_to_b.dot (axis).signum()
          } else {
            -1.0
          };
          let half_axis = 0.5 * axis;
          let normal    = if axis.is_zero() {
            Unit3::normalize (vector3 (0.0, signum_y, -1.0))
          } else {
            Unit3::unchecked_approx (-axis / distance)
          };
          let midpoint  = nearest_a + half_axis;
          Proximity { distance, half_axis, midpoint, normal }
        }
        (false, true) => {
          // nearest/deepest point is on an edge parallel to the Y axis
          let signum_x  = (position_a.0.x - position_b.0.x).signum();
          let nearest_b = point3 (
            //position_b.0.x + signum_x * cuboid_b_extent_x,
            signum_x.mul_add (cuboid_b_extent_x, position_b.0.x),
            position_a.0.y,
            cuboid_min_b.0.z
          );
          let a_to_b    = nearest_b - position_a;
          let nearest_a = position_a + if a_to_b.is_zero() {
            -signum_x * vector3 (sphere_a_radius, 0.0, 0.0)
          } else {
            sphere_a_radius * a_to_b.normalized()
          };
          let axis      = nearest_b - nearest_a;
          let distance  = axis.magnitude() * if !a_to_b.is_zero() {
            a_to_b.dot (axis).signum()
          } else {
            -1.0
          };
          let half_axis = 0.5 * axis;
          let normal    = if axis.is_zero() {
            Unit3::normalize (vector3 (signum_x, 0.0, -1.0))
          } else {
            Unit3::unchecked_approx (-axis / distance)
          };
          let midpoint  = nearest_a + half_axis;
          Proximity { distance, half_axis, midpoint, normal }
        }
        (true, true) => {
          // nearest/deepest point is in the -Z direction
          let nearest_a : Point3 <f64> =
            position_a + vector3 (0.0, 0.0, sphere_a_radius);
          let nearest_b : Point3 <f64> =
            [position_a.0.x, position_a.0.y, cuboid_min_b.0.z].into();
          let axis      = nearest_b - nearest_a;
          let normal    = Unit3::axis_z().invert();
          let half_axis = 0.5 * axis;
          let distance  = nearest_b.0.z - nearest_a.0.z;
          let midpoint  = nearest_a + half_axis;
          Proximity { distance, half_axis, midpoint, normal }
        }
      }
    } else {
      // sphere center overlaps on z axis
      match (overlap_x, overlap_y) {
        (false, false) => {
          // nearest point is on a Z parallel edge of the cuboid
          let b_to_a_signum_x = (position_a.0.x - position_b.0.x).signum();
          let b_to_a_signum_y = (position_a.0.y - position_b.0.y).signum();
          let b_to_a_unit_sigvec = vector3 (b_to_a_signum_x, b_to_a_signum_y, 0.0)
            .normalized();
          let nearest_b = Point3::from ([
            //position_b.0.x + b_to_a_signum_x * cuboid_b_extent_x,
            b_to_a_signum_x.mul_add (cuboid_b_extent_x, position_b.0.x),
            //position_b.0.y + b_to_a_signum_y * cuboid_b_extent_y,
            b_to_a_signum_y.mul_add (cuboid_b_extent_y, position_b.0.y),
            position_a.0.z
          ]);
          let a_to_b    = vector3 (
            nearest_b.0.x - position_a.0.x,
            nearest_b.0.y - position_a.0.y,
            0.0
          );
          let nearest_a = point3 (position_a.0.x, position_a.0.y, position_a.0.z)
            + sphere_a_radius * if !a_to_b.is_zero() {
              a_to_b.normalized()
            } else {
              -b_to_a_unit_sigvec
            };
          let axis      = nearest_b - nearest_a;
          let distance  = axis.magnitude() * if !a_to_b.is_zero() {
            a_to_b.dot (axis).signum()
          } else {
            -1.0
          };
          let half_axis = 0.5 * axis;
          let normal    = Unit3::unchecked_approx (if axis.is_zero() {
            b_to_a_unit_sigvec
          } else {
            -axis / distance
          });
          let midpoint  = nearest_a + half_axis;
          Proximity { distance, half_axis, midpoint, normal }
        }
        (true, false)  => {
          // nearest point is on a Z/X parallel face of the cuboid
          let b_to_a_signum_y    = (position_a.0.y - position_b.0.y).signum();
          let b_to_a_unit_sigvec = vector3 (0.0, b_to_a_signum_y, 0.0);
          let nearest_b = Point3::from ([
            position_a.0.x,
            //position_b.0.y + b_to_a_signum_y * cuboid_b_extent_y,
            b_to_a_signum_y.mul_add (cuboid_b_extent_y, position_b.0.y),
            position_a.0.z
          ]);
          let a_to_b    = vector3 (0.0, nearest_b.0.y - position_a.0.y, 0.0);
          let nearest_a =
            point3 (position_a.0.x, position_a.0.y, position_a.0.z) +
            vector3 (0.0, -b_to_a_signum_y * sphere_a_radius, 0.0);
          let axis      = nearest_b - nearest_a;
          let distance  = axis.magnitude() * if !a_to_b.is_zero() {
            a_to_b.dot (axis).signum()
          } else {
            -1.0
          };
          let half_axis = 0.5 * axis;
          let normal    = Unit3::unchecked_approx (if axis.is_zero() {
            b_to_a_unit_sigvec
          } else {
            -axis / distance
          });
          let midpoint  = nearest_a + half_axis;
          Proximity { distance, half_axis, midpoint, normal }
        }
        (false, true)  => {
          // nearest point is on a Z/Y parallel face of the cuboid
          let b_to_a_signum_x    = (position_a.0.x - position_b.0.x).signum();
          let b_to_a_unit_sigvec = vector3 (b_to_a_signum_x, 0.0, 0.0);
          let nearest_b = Point3::from ([
            //position_b.0.x + b_to_a_signum_x * cuboid_b_extent_x,
            b_to_a_signum_x.mul_add (cuboid_b_extent_x, position_b.0.x),
            position_a.0.y,
            position_a.0.z
          ]);
          let a_to_b    = vector3 (nearest_b.0.x - position_a.0.x, 0.0, 0.0);
          let nearest_a =
            point3 (position_a.0.x, position_a.0.y, position_a.0.z) +
            vector3 (-b_to_a_signum_x * sphere_a_radius, 0.0, 0.0);
          let axis      = nearest_b - nearest_a;
          let distance  = axis.magnitude() * if !a_to_b.is_zero() {
            a_to_b.dot (axis).signum()
          } else {
            -1.0
          };
          let half_axis = 0.5 * axis;
          let normal    = Unit3::unchecked_approx (if axis.is_zero() {
            b_to_a_unit_sigvec
          } else {
            -axis / distance
          });
          let midpoint  = nearest_a + half_axis;
          Proximity { distance, half_axis, midpoint, normal }
        }
        (true, true) => {
          // center inside X/Y bounds
          if position_a == position_b {
            let depth_x = cuboid_b_extent_x + sphere_a_radius;
            let depth_y = cuboid_b_extent_y + sphere_a_radius;
            let depth_z = cuboid_b_extent_z + sphere_a_radius;
            if depth_x <= depth_y && depth_x <= depth_z {
              let distance  = -depth_x;
              let half_axis = [0.5 * depth_x, 0.0, 0.0].into();
              let normal    = Unit3::axis_x();
              let midpoint  = position_a + half_axis -
                vector3 (sphere_a_radius, 0.0, 0.0);
              Proximity { distance, half_axis, midpoint, normal }
            } else if depth_y <= depth_x && depth_y <= depth_z {
              let distance = -depth_y;
              let half_axis = [0.0, 0.5 * depth_y, 0.0].into();
              let normal    = Unit3::axis_y();
              let midpoint  = position_a + half_axis -
                vector3 (0.0, sphere_a_radius, 0.0);
              Proximity { distance, half_axis, midpoint, normal }
            } else {
              debug_assert!(depth_z < depth_x);
              debug_assert!(depth_z < depth_y);
              let distance = -depth_z;
              let half_axis = [0.0, 0.0, 0.5 * depth_z].into();
              let normal    = Unit3::axis_z();
              let midpoint  = position_a + half_axis -
                vector3 (0.0, 0.0, sphere_a_radius);
              Proximity { distance, half_axis, midpoint, normal }
            }
          } else {
            // NOTE: we want the following to map zero values to positive 1.0
            // so we use f64::signum instead of sigvec
            let b_to_a_sigvec = (position_a - position_b).map (f64::signum);
            let corner_b      = position_b + cuboid_b.max().0 * b_to_a_sigvec;
            let depth_x       = (corner_b.0.x - position_a.0.x).abs() + sphere_a_radius;
            let depth_y       = (corner_b.0.y - position_a.0.y).abs() + sphere_a_radius;
            let depth_z       = (corner_b.0.z - position_a.0.z).abs() + sphere_a_radius
              - if position_a.0.z.abs() > corner_b.0.z.abs() {
                2.0 * (position_a.0.z.abs() - corner_b.0.z.abs())
              } else {
                0.0
              };

            if depth_x <= depth_y && depth_x <= depth_z {
              let distance  = -depth_x;
              let half_axis = [0.5 * b_to_a_sigvec.x * depth_x, 0.0, 0.0].into();
              let normal = Unit3::unchecked ([b_to_a_sigvec.x, 0.0, 0.0].into());
              let midpoint  = position_a + half_axis +
                vector3 (-b_to_a_sigvec.x * sphere_a_radius, 0.0, 0.0);
              Proximity { distance, half_axis, midpoint, normal }
            } else if depth_y <= depth_x && depth_y <= depth_z {
              let distance  = -depth_y;
              let half_axis = [0.0, 0.5 * b_to_a_sigvec.y * depth_y, 0.0].into();
              let normal = Unit3::unchecked ([0.0, b_to_a_sigvec.y, 0.0].into());
              let midpoint  = position_a + half_axis +
                vector3 (0.0, -b_to_a_sigvec.y * sphere_a_radius, 0.0);
              Proximity { distance, half_axis, midpoint, normal }
            } else {
              debug_assert!(depth_z < depth_x);
              debug_assert!(depth_z < depth_y);
              let nearest_a = point3 (
                position_a.0.x,
                position_a.0.y,
                //position_a.0.z - b_to_a_sigvec.z * sphere_a_radius
                b_to_a_sigvec.z.mul_add (-sphere_a_radius, position_a.0.z)
              );
              let nearest_b = point3 (
                position_a.0.x,
                position_a.0.y,
                //position_b.0.z + b_to_a_sigvec.z * cuboid_b_extent_z
                b_to_a_sigvec.z.mul_add (cuboid_b_extent_z, position_b.0.z)
              );
              let axis      = nearest_b - nearest_a;
              debug_assert_eq!(axis.x, 0.0);
              debug_assert_eq!(axis.y, 0.0);
              let distance  = -axis.z.abs();
              let half_axis = 0.5 * axis;
              let normal = Unit3::unchecked ([0.0, 0.0, b_to_a_sigvec.z].into());
              let midpoint  = nearest_a + half_axis;
              Proximity { distance, half_axis, midpoint, normal }
            }
          }
        }
      }
    }
  } // end query_sphere_cuboid

  pub fn query_sphere_orthant (
    position_a : &Point3 <f64>, sphere_a : &shape::Sphere <f64>,
    position_b : &Point3 <f64>, orthant_b : &shape::Orthant
  ) -> Self {
    let radius_a    = *sphere_a.radius;
    let normal      = orthant_b.normal_axis.to_vec::<f64>();
    let normal_abs  = normal.map (f64::abs);
    let surface_sigvec = vector3 (1.0, 1.0, 1.0) - normal_abs;
    let nearest_a   = position_a - vector3 (radius_a, radius_a, radius_a) * normal;
    let nearest_b   =
      Point3 (position_a.0 * surface_sigvec + position_b.0 * normal_abs);
    let axis        = nearest_b - nearest_a;
    debug_assert_eq!(axis.sum().abs(), axis.magnitude());
    let distance    = -normal.dot (axis).signum() * axis.sum().abs();
    let half_axis   = 0.5 * axis;
    let midpoint    = nearest_a + half_axis;
    let normal      = Unit3::unchecked (normal);
    Proximity { distance, half_axis, midpoint, normal }
  } // end query_sphere_orthant

  pub fn query_capsule_cuboid (
    position_a : &Point3 <f64>, capsule_a : &shape::Capsule <f64>,
    position_b : &Point3 <f64>, cuboid_b  : &shape::Cuboid  <f64>
  ) -> Self {
    use num::Zero;
    let shaft_max_a  = position_a + vector3 (0.0, 0.0, *capsule_a.half_height);
    let shaft_min_a  = position_a - vector3 (0.0, 0.0, *capsule_a.half_height);
    let cuboid_max_b = position_b + cuboid_b.max().0;
    let cuboid_min_b = position_b - cuboid_b.max().0;

    let overlap_x = cuboid_min_b.0.x < position_a.0.x &&
      position_a.0.x < cuboid_max_b.0.x;
    let overlap_y = cuboid_min_b.0.y < position_a.0.y &&
      position_a.0.y < cuboid_max_b.0.y;

    let capsule_a_radius      = *capsule_a.radius;
    let capsule_a_half_height = *capsule_a.half_height;
    let cuboid_b_extent_x     = *cuboid_b.extents.x;
    let cuboid_b_extent_y     = *cuboid_b.extents.y;
    let cuboid_b_extent_z     = *cuboid_b.extents.z;
    if shaft_min_a.0.z >= cuboid_max_b.0.z {
      // capsule above cuboid
      match (overlap_x, overlap_y) {
        (false, false) => {
          // nearest/deepest point is a corner of the cuboid
          let b_to_a        = shaft_min_a - position_b;
          debug_assert!(!b_to_a.x.is_zero());
          debug_assert!(!b_to_a.y.is_zero());
          debug_assert!(!b_to_a.z.is_zero());
          let b_to_a_sigvec = b_to_a.map (f64::signum);
          let nearest_b     = position_b + cuboid_b.max().0 * b_to_a_sigvec;
          let a_to_b        = nearest_b - shaft_min_a;
          let nearest_a     = shaft_min_a + if !a_to_b.is_zero() {
            capsule_a_radius * a_to_b.normalized()
          } else {
            [-b_to_a_sigvec.x * capsule_a_radius, 0.0, 0.0].into()
          };
          let axis          = nearest_b - nearest_a;
          let distance      = axis.magnitude() * if !a_to_b.is_zero() {
            a_to_b.dot (axis).signum()
          } else {
            -1.0
          };
          let half_axis     = 0.5 * axis;
          let normal        = if axis.is_zero() {
            Unit3::normalize (b_to_a_sigvec)
          } else {
            Unit3::unchecked_approx (-axis / distance)
          };
          let midpoint      = nearest_a + half_axis;
          Proximity { distance, half_axis, midpoint, normal }
        }
        (true, false) => {
          // nearest/deepest point is on an edge parallel to the X axis
          let signum_y  = (shaft_min_a.0.y - position_b.0.y).signum();
          let nearest_b = point3 (
            shaft_min_a.0.x,
            position_b.0.y + signum_y * cuboid_b_extent_y,
            cuboid_max_b.0.z
          );
          let a_to_b    = nearest_b - shaft_min_a;
          let nearest_a = shaft_min_a + if a_to_b.is_zero() {
            -signum_y * vector3 (0.0, capsule_a_radius, 0.0)
          } else {
            capsule_a_radius * a_to_b.normalized()
          };
          let axis      = nearest_b - nearest_a;
          let distance  = axis.magnitude() * if !a_to_b.is_zero() {
            a_to_b.dot (axis).signum()
          } else {
            -1.0
          };
          let half_axis = 0.5 * axis;
          let normal    = if axis.is_zero() {
            Unit3::normalize (vector3 (0.0, signum_y, 1.0))
          } else {
            Unit3::unchecked_approx (-axis / distance)
          };
          let midpoint  = nearest_a + half_axis;
          Proximity { distance, half_axis, midpoint, normal }
        }
        (false, true) => {
          // nearest/deepest point is on an edge parallel to the Y axis
          let signum_x  = (shaft_min_a.0.x - position_b.0.x).signum();
          let nearest_b = point3 (
            position_b.0.x + signum_x * cuboid_b_extent_x,
            shaft_min_a.0.y,
            cuboid_max_b.0.z
          );
          let a_to_b    = nearest_b - shaft_min_a;
          let nearest_a = shaft_min_a + if a_to_b.is_zero() {
            -signum_x * vector3 (capsule_a_radius, 0.0, 0.0)
          } else {
            capsule_a_radius * a_to_b.normalized()
          };
          let axis      = nearest_b - nearest_a;
          let distance  = axis.magnitude() * if !a_to_b.is_zero() {
            a_to_b.dot (axis).signum()
          } else {
            -1.0
          };
          let half_axis = 0.5 * axis;
          let normal    = if axis.is_zero() {
            Unit3::normalize (vector3 (signum_x, 0.0, 1.0))
          } else {
            Unit3::unchecked_approx (-axis / distance)
          };
          let midpoint  = nearest_a + half_axis;
          Proximity { distance, half_axis, midpoint, normal }
        }
        (true, true) => {
          // nearest/deepest point is in the -Z direction
          let nearest_a : Point3 <f64> =
            shaft_min_a - vector3 (0.0, 0.0, capsule_a_radius);
          let nearest_b : Point3 <f64> =
            [position_a.0.x, position_a.0.y, cuboid_max_b.0.z].into();
          let axis      = nearest_b - nearest_a;
          let normal    = Unit3::axis_z();
          let half_axis = 0.5 * axis;
          let distance  = nearest_a.0.z - nearest_b.0.z;
          let midpoint  = nearest_a + half_axis;
          Proximity { distance, half_axis, midpoint, normal }
        }
      }
    } else if shaft_max_a.0.z <= cuboid_min_b.0.z {
      // capsule below cuboid
      match (overlap_x, overlap_y) {
        (false, false) => {
          // nearest/deepest point is a corner of the cuboid
          let b_to_a        = shaft_max_a - position_b;
          debug_assert!(!b_to_a.x.is_zero());
          debug_assert!(!b_to_a.y.is_zero());
          debug_assert!(!b_to_a.z.is_zero());
          let b_to_a_sigvec = b_to_a.map (f64::signum);
          let nearest_b     = position_b + cuboid_b.max().0 * b_to_a_sigvec;
          let a_to_b        = nearest_b - shaft_max_a;
          let nearest_a     = shaft_max_a + if !a_to_b.is_zero() {
            capsule_a_radius * a_to_b.normalized()
          } else {
            [-b_to_a_sigvec.x * capsule_a_radius, 0.0, 0.0].into()
          };
          let axis          = nearest_b - nearest_a;
          let distance      = axis.magnitude() * if !a_to_b.is_zero() {
            a_to_b.dot (axis).signum()
          } else {
            -1.0
          };
          let half_axis     = 0.5 * axis;
          let normal        = if axis.is_zero() {
            Unit3::normalize (b_to_a_sigvec)
          } else {
            Unit3::unchecked_approx (-axis / distance)
          };
          let midpoint      = nearest_a + half_axis;
          Proximity { distance, half_axis, midpoint, normal }
        }
        (true, false) => {
          // nearest/deepest point is on an edge parallel to the X axis
          let signum_y  = (shaft_max_a.0.y - position_b.0.y).signum();
          let nearest_b = point3 (
            shaft_max_a.0.x,
            position_b.0.y + signum_y * cuboid_b_extent_y,
            cuboid_min_b.0.z
          );
          let a_to_b    = nearest_b - shaft_max_a;
          let nearest_a = shaft_max_a + if a_to_b.is_zero() {
            -signum_y * vector3 (0.0, capsule_a_radius, 0.0)
          } else {
            capsule_a_radius * a_to_b.normalized()
          };
          let axis      = nearest_b - nearest_a;
          let distance  = axis.magnitude() * if !a_to_b.is_zero() {
            a_to_b.dot (axis).signum()
          } else {
            -1.0
          };
          let half_axis = 0.5 * axis;
          let normal    = if axis.is_zero() {
            Unit3::normalize (vector3 (0.0, signum_y, -1.0))
          } else {
            Unit3::unchecked_approx (-axis / distance)
          };
          let midpoint  = nearest_a + half_axis;
          Proximity { distance, half_axis, midpoint, normal }
        }
        (false, true) => {
          // nearest/deepest point is on an edge parallel to the Y axis
          let signum_x  = (shaft_max_a.0.x - position_b.0.x).signum();
          let nearest_b = point3 (
            position_b.0.x + signum_x * cuboid_b_extent_x,
            shaft_max_a.0.y,
            cuboid_min_b.0.z
          );
          let a_to_b    = nearest_b - shaft_max_a;
          let nearest_a = shaft_max_a + if a_to_b.is_zero() {
            -signum_x * vector3 (capsule_a_radius, 0.0, 0.0)
          } else {
            capsule_a_radius * a_to_b.normalized()
          };
          let axis      = nearest_b - nearest_a;
          let distance  = axis.magnitude() * if !a_to_b.is_zero() {
            a_to_b.dot (axis).signum()
          } else {
            -1.0
          };
          let half_axis = 0.5 * axis;
          let normal    = if axis.is_zero() {
            Unit3::normalize (vector3 (signum_x, 0.0, -1.0))
          } else {
            Unit3::unchecked_approx (-axis / distance)
          };
          let midpoint  = nearest_a + half_axis;
          Proximity { distance, half_axis, midpoint, normal }
        }
        (true, true) => {
          // nearest/deepest point is in the -Z direction
          let nearest_a : Point3 <f64> =
            shaft_max_a + vector3 (0.0, 0.0, capsule_a_radius);
          let nearest_b : Point3 <f64> =
            [position_a.0.x, position_a.0.y, cuboid_min_b.0.z].into();
          let axis      = nearest_b - nearest_a;
          let normal    = Unit3::axis_z().invert();
          let half_axis = 0.5 * axis;
          let distance  = nearest_b.0.z - nearest_a.0.z;
          let midpoint  = nearest_a + half_axis;
          Proximity { distance, half_axis, midpoint, normal }
        }
      }
    } else {
      // capsule cylinder shaft overlaps on z axis
      let max_z = f64::min (shaft_max_a.0.z, cuboid_max_b.0.z);
      let min_z = f64::max (shaft_min_a.0.z, cuboid_min_b.0.z);
      debug_assert_ne!(min_z, max_z);
      let mid_z = 0.5 * (max_z + min_z);
      match (overlap_x, overlap_y) {
        (false, false) => {
          // nearest point is on a Z parallel edge of the cuboid
          let b_to_a_signum_x = (position_a.0.x - position_b.0.x).signum();
          let b_to_a_signum_y = (position_a.0.y - position_b.0.y).signum();
          let b_to_a_unit_sigvec = vector3 (b_to_a_signum_x, b_to_a_signum_y, 0.0)
            .normalized();
          let nearest_b = Point3::from ([
            //position_b.0.x + b_to_a_signum_x * cuboid_b_extent_x,
            b_to_a_signum_x.mul_add (cuboid_b_extent_x, position_b.0.x),
            //position_b.0.y + b_to_a_signum_y * cuboid_b_extent_y,
            b_to_a_signum_y.mul_add (cuboid_b_extent_y, position_b.0.y),
            mid_z
          ]);
          let a_to_b    = vector3 (
            nearest_b.0.x - position_a.0.x,
            nearest_b.0.y - position_a.0.y,
            0.0
          );
          let nearest_a = point3 (position_a.0.x, position_a.0.y, mid_z)
            + capsule_a_radius * if !a_to_b.is_zero() {
              a_to_b.normalized()
            } else {
              -b_to_a_unit_sigvec
            };
          let axis      = nearest_b - nearest_a;
          let distance  = axis.magnitude() * if !a_to_b.is_zero() {
            a_to_b.dot (axis).signum()
          } else {
            -1.0
          };
          let half_axis = 0.5 * axis;
          let normal    = Unit3::unchecked_approx (if axis.is_zero() {
            b_to_a_unit_sigvec
          } else {
            -axis / distance
          });
          let midpoint  = nearest_a + half_axis;
          Proximity { distance, half_axis, midpoint, normal }
        }
        (true, false)  => {
          // nearest point is on a Z/X parallel face of the cuboid
          let b_to_a_signum_y    = (position_a.0.y - position_b.0.y).signum();
          let b_to_a_unit_sigvec = vector3 (0.0, b_to_a_signum_y, 0.0);
          let nearest_b = Point3::from ([
            position_a.0.x,
            //position_b.0.y + b_to_a_signum_y * cuboid_b_extent_y,
            b_to_a_signum_y.mul_add (cuboid_b_extent_y, position_b.0.y),
            mid_z
          ]);
          let a_to_b    = vector3 (0.0, nearest_b.0.y - position_a.0.y, 0.0);
          let nearest_a = point3 (position_a.0.x, position_a.0.y, mid_z) +
            vector3 (0.0, -b_to_a_signum_y * capsule_a_radius, 0.0);
          let axis      = nearest_b - nearest_a;
          let distance  = axis.magnitude() * if !a_to_b.is_zero() {
            a_to_b.dot (axis).signum()
          } else {
            -1.0
          };
          let half_axis = 0.5 * axis;
          let normal    = Unit3::unchecked_approx (if axis.is_zero() {
            b_to_a_unit_sigvec
          } else {
            -axis / distance
          });
          let midpoint  = nearest_a + half_axis;
          Proximity { distance, half_axis, midpoint, normal }
        }
        (false, true)  => {
          // nearest point is on a Z/Y parallel face of the cuboid
          let b_to_a_signum_x    = (position_a.0.x - position_b.0.x).signum();
          let b_to_a_unit_sigvec = vector3 (b_to_a_signum_x, 0.0, 0.0);
          let nearest_b = Point3::from ([
            //position_b.0.x + b_to_a_signum_x * cuboid_b_extent_x,
            b_to_a_signum_x.mul_add (cuboid_b_extent_x, position_b.0.x),
            position_a.0.y,
            mid_z
          ]);
          let a_to_b    = vector3 (nearest_b.0.x - position_a.0.x, 0.0, 0.0);
          let nearest_a = point3 (position_a.0.x, position_a.0.y, mid_z) +
            vector3 (-b_to_a_signum_x * capsule_a_radius, 0.0, 0.0);
          let axis      = nearest_b - nearest_a;
          let distance  = axis.magnitude() * if !a_to_b.is_zero() {
            a_to_b.dot (axis).signum()
          } else {
            -1.0
          };
          let half_axis = 0.5 * axis;
          let normal    = Unit3::unchecked_approx (if axis.is_zero() {
            b_to_a_unit_sigvec
          } else {
            -axis / distance
          });
          let midpoint  = nearest_a + half_axis;
          Proximity { distance, half_axis, midpoint, normal }
        }
        (true, true) => {
          // cylinder shaft inside X/Y bounds
          if position_a == position_b {
            let depth_x = cuboid_b_extent_x + capsule_a_radius;
            let depth_y = cuboid_b_extent_y + capsule_a_radius;
            let depth_z = cuboid_b_extent_z + capsule_a_half_height +
              capsule_a_radius;
            if depth_x <= depth_y && depth_x <= depth_z {
              let distance  = -depth_x;
              let half_axis = [0.5 * depth_x, 0.0, 0.0].into();
              let normal    = Unit3::axis_x();
              let midpoint  = position_a + half_axis -
                vector3 (capsule_a_radius, 0.0, 0.0);
              Proximity { distance, half_axis, midpoint, normal }
            } else if depth_y <= depth_x && depth_y <= depth_z {
              let distance = -depth_y;
              let half_axis = [0.0, 0.5 * depth_y, 0.0].into();
              let normal    = Unit3::axis_y();
              let midpoint  = position_a + half_axis -
                vector3 (0.0, capsule_a_radius, 0.0);
              Proximity { distance, half_axis, midpoint, normal }
            } else {
              debug_assert!(depth_z < depth_x);
              debug_assert!(depth_z < depth_y);
              let distance = -depth_z;
              let half_axis = [0.0, 0.0, 0.5 * depth_z].into();
              let normal    = Unit3::axis_z();
              let midpoint  = position_a + half_axis -
                vector3 (0.0, 0.0, capsule_a_radius + capsule_a_half_height);
              Proximity { distance, half_axis, midpoint, normal }
            }
          } else {
            // NOTE: we want the following to map zero values to positive 1.0
            // so we use f64::signum instead of sigvec
            let b_to_a_sigvec = (position_a - position_b).map (f64::signum);
            let corner_b      = position_b + cuboid_b.max().0 * b_to_a_sigvec;
            let depth_x       = (corner_b.0.x - position_a.0.x).abs() + capsule_a_radius;
            let depth_y       = (corner_b.0.y - position_a.0.y).abs() + capsule_a_radius;
            let depth_z       = (corner_b.0.z - position_a.0.z).abs() +
              capsule_a_radius + capsule_a_half_height -
              if position_a.0.z.abs() > corner_b.0.z.abs() {
                2.0 * (position_a.0.z.abs() - corner_b.0.z.abs())
              } else {
                0.0
              };

            if depth_x <= depth_y && depth_x <= depth_z {
              let distance  = -depth_x;
              let half_axis = [0.5 * b_to_a_sigvec.x * depth_x, 0.0, 0.0].into();
              let normal = Unit3::unchecked ([b_to_a_sigvec.x, 0.0, 0.0].into());
              let midpoint  = position_a + half_axis +
                vector3 (-b_to_a_sigvec.x * capsule_a_radius, 0.0, 0.0);
              Proximity { distance, half_axis, midpoint, normal }
            } else if depth_y <= depth_x && depth_y <= depth_z {
              let distance  = -depth_y;
              let half_axis = [0.0, 0.5 * b_to_a_sigvec.y * depth_y, 0.0].into();
              let normal = Unit3::unchecked ([0.0, b_to_a_sigvec.y, 0.0].into());
              let midpoint  = position_a + half_axis +
                vector3 (0.0, -b_to_a_sigvec.y * capsule_a_radius, 0.0);
              Proximity { distance, half_axis, midpoint, normal }
            } else {
              debug_assert!(depth_z < depth_x);
              debug_assert!(depth_z < depth_y);
              let nearest_a = point3 (
                position_a.0.x,
                position_a.0.y,
                //position_a.0.z - b_to_a_sigvec.z * (capsule_a_radius + capsule_a_half_height)
                b_to_a_sigvec.z
                  .mul_add (-(capsule_a_radius + capsule_a_half_height), position_a.0.z)
              );
              let nearest_b = point3 (
                position_a.0.x,
                position_a.0.y,
                //position_b.0.z + b_to_a_sigvec.z * cuboid_b_extent_z
                b_to_a_sigvec.z.mul_add (cuboid_b_extent_z, position_b.0.z)
              );
              let axis      = nearest_b - nearest_a;
              debug_assert_eq!(axis.x, 0.0);
              debug_assert_eq!(axis.y, 0.0);
              let distance  = -axis.z.abs();
              let half_axis = 0.5 * axis;
              let normal = Unit3::unchecked ([0.0, 0.0, b_to_a_sigvec.z].into());
              let midpoint  = nearest_a + half_axis;
              Proximity { distance, half_axis, midpoint, normal }
            }
          }
        }
      }
    }
  } // end query_capsule_cuboid

  pub fn query_capsule_orthant (
    position_a : &Point3 <f64>, capsule_a : &shape::Capsule <f64>,
    position_b : &Point3 <f64>, orthant_b : &shape::Orthant
  ) -> Self {
    let radius_a    = *capsule_a.radius;
    let half_height = *capsule_a.half_height;
    let normal      = orthant_b.normal_axis.to_vec::<f64>();
    let normal_abs  = normal.map (f64::abs);
    let surface_sigvec = vector3 (1.0, 1.0, 1.0) - normal_abs;
    let nearest_a   = position_a -
      vector3 (radius_a, radius_a, radius_a + half_height) * normal;
    let nearest_b   =
      Point3 (position_a.0 * surface_sigvec + position_b.0 * normal_abs);
    let axis        = nearest_b - nearest_a;
    debug_assert_eq!(axis.sum().abs(), axis.magnitude());
    let distance    = -normal.dot (axis).signum() * axis.sum().abs();
    let half_axis   = 0.5 * axis;
    let midpoint    = nearest_a + half_axis;
    let normal      = Unit3::unchecked (normal);
    Proximity { distance, half_axis, midpoint, normal }
  } // end query_capsule_orthant

  pub fn query_cuboid_orthant (
    position_a : &Point3 <f64>, cuboid_a  : &shape::Cuboid <f64>,
    position_b : &Point3 <f64>, orthant_b : &shape::Orthant
  ) -> Self {
    let normal      = orthant_b.normal_axis.to_vec::<f64>();
    let normal_abs  = normal.map (f64::abs);
    let surface_sigvec = vector3 (1.0, 1.0, 1.0) - normal_abs;
    let nearest_a   = position_a - cuboid_a.extents.map (|e| *e) * normal;
    let nearest_b   =
      Point3 (position_a.0 * surface_sigvec + position_b.0 * normal_abs);
    let axis        = nearest_b - nearest_a;
    debug_assert_eq!(axis.sum().abs(), axis.magnitude());
    let distance    = -normal.dot (axis).signum() * axis.sum().abs();
    let half_axis   = 0.5 * axis;
    let midpoint    = nearest_a + half_axis;
    let normal      = Unit3::unchecked (normal);
    Proximity { distance, half_axis, midpoint, normal }
  } // end query_cuboid_orthant

  pub fn query_triangle_triangle (
    triangle_a : &geometry::Triangle3 <f64>,
    triangle_b : &geometry::Triangle3 <f64>
  ) -> Self {
    fn query_triangle_triangle_recursive (
      triangle_a : &geometry::Triangle3 <f64>,
      triangle_b : &geometry::Triangle3 <f64>,
      is_recurse : bool
    ) -> Proximity {
      let mut out = None;
      let mut min_distance = f64::MAX;
      for edge in triangle_a.edges() {
        let proximity = Proximity::query_triangle_segment (triangle_b, &edge);
        if proximity.distance < -f64::EPSILON * 2.0.powi (12) {
          out = None;
          break
        }
        if proximity.distance < min_distance {
          min_distance = proximity.distance;
          // we swap the roles of object A and B
          out = Some (proximity.flip());
        }
      }
      if out.is_some() {
        // only check the other triangle if the first triangle didn't overlap
        for edge in triangle_b.edges() {
          let proximity = Proximity::query_triangle_segment (triangle_a, &edge);
          if proximity.distance < -f64::EPSILON * 2.0.powi (12) {
            out = None;
            break
          }
          if proximity.distance < min_distance {
            min_distance = proximity.distance;
            out = Some (proximity);
          }
        }
      }
      if out.is_none() {
        // triangles overlap
        if is_recurse {
          log::error!("triangle/triangle proximity query can't recurse more than once");
          panic!("triangle/triangle proximity query can't recurse more than once");
        }
        // edge intersected: construct minkowski hull A - B
        let points : [Point3 <f64>; 9] = std::array::from_fn (|i|{
          let p = triangle_a.points()[i / 3];
          let q = triangle_b.points()[i % 3];
          Point3 (p.0 - q.0)
        });
        // brute force search: there are 9 choose 3 = 84 unique triangles, for each
        // triangle, calculate the distance and whether or not it is exterior; discard
        // any triangles that are further away than the nearest exterior triangle
        let get_triangle = |a, b, c|
          geometry::Triangle3::new (points[a], points[b], points[c]);
        let mut nearest_distance : Option <geometry::Distance3 <f64>> = None;
        let tolerance = f64::EPSILON * 2.0.powi (10);
        for i in 0..9 {
          for j in (i + 1)..9 {
            for k in (j + 1)..9 {
              if let Some (triangle) = get_triangle (i, j, k) {
                let distance =
                  geometry::Distance3::triangle_point (triangle, Point3::origin());
                if distance.distance_squared() < nearest_distance.as_ref().map_or (
                  NonNegative::unchecked (f64::MAX),
                  geometry::Distance3::distance_squared)
                {
                  // check if it is exterior
                  let normal = {
                    let mut normal = triangle.normal();
                    if normal.dot (triangle.point_a().0) < 0.0 {
                      normal = -normal
                    }
                    normal
                  };
                  let mut exterior = true;
                  for point in points.iter() {
                    if normal.dot (point - triangle.point_a()) > tolerance {
                      exterior = false;
                      break
                    }
                  }
                  if exterior {
                    nearest_distance = Some (distance);
                  }
                }
              }
            }
          }
        }
        let mut nearest = nearest_distance.unwrap();
        let distance    = -*nearest.distance();
        let half_axis   = -nearest.nearest_a().0 * 0.5;
        let resolved_proximity = {
          let translate = |mut triangle : geometry::Triangle3 <f64>, displacement| {
            triangle.translate (displacement);
            triangle
          };
          let resolved_a = translate (*triangle_a, half_axis);
          let resolved_b = translate (*triangle_b, -half_axis);
          // TODO: instead of recursing can we query nearest features ?
          query_triangle_triangle_recursive (&resolved_a, &resolved_b, true)
        };
        out = Some(Proximity { distance, half_axis, .. resolved_proximity });
      }
      out.unwrap()
    }
    query_triangle_triangle_recursive (triangle_a, triangle_b, false)
  }

  pub fn query_triangle_line (
    triangle : &geometry::Triangle3 <f64>,
    line     : &geometry::Segment3 <f64>
  ) -> Self {
    use approx::{AbsDiffEq, RelativeEq};
    let (([s, t], nearest_a), (_, nearest_b)) =
      geometry::distance::nearest_triangle3_line3 (*triangle, *line);
    if approx::relative_eq!(nearest_a, nearest_b,
      max_relative = f64::default_max_relative() * 256.0,
      epsilon      = f64::default_epsilon() * 256.0)
    {
      // line and triangle intersect
      let distance;
      let half_axis;
      let midpoint;
      let normal;
      let triangle_normal = triangle.normal();
      if approx::abs_diff_eq!(triangle_normal.dot (*line.vector()), 0.0) {
        // line and triangle are parallel:
        // normal is perpendicular to triangle face
        let a_to_b = nearest_b - nearest_a;
        distance   = a_to_b.magnitude();
        half_axis  = a_to_b * 0.5;
        midpoint   = nearest_a + half_axis;
        normal     = triangle_normal.normalize();
      } else {
        let y = nearest_a;
        let (d0, d1, d2) = triangle_distance_to_edges (*triangle, *s, *t);
        // project to nearest edge
        let edge0 = triangle.point_c() - triangle.point_b();
        let edge1 = triangle.point_b() - triangle.point_a();
        let edge2 = triangle.point_c() - triangle.point_a();
        let len2_e0 = edge0.magnitude_squared();
        let len2_e1 = edge1.magnitude_squared();
        let len2_e2 = edge2.magnitude_squared();
        let mut e0_is_nearest = false;
        let mut e1_is_nearest = false;
        let mut e2_is_nearest = false;
        let nearest_a = if d0 <= d1 && d0 <= d2 {
          e0_is_nearest = true;
          let t = (y - triangle.point_b()).dot (edge0) / len2_e0;
          triangle.point_b() + edge0 * t
        } else if d1 <= d2 {
          e1_is_nearest = true;
          let t = (y - triangle.point_a()).dot (edge1) / len2_e1;
          triangle.point_a() + edge1 * t
        } else {
          debug_assert!(d2 < d0);
          debug_assert!(d2 < d1);
          e2_is_nearest = true;
          let t = (y - triangle.point_a()).dot (edge2) / len2_e2;
          triangle.point_a() + edge2 * t
        };
        let nearest_b = y;
        let a_to_b    = nearest_b - nearest_a;
        distance      = -a_to_b.magnitude();
        half_axis     = a_to_b * 0.5;
        midpoint      = nearest_a + half_axis;
        // TODO: would it be better to check distance approx eq 0.0 here? in theory
        // there are less comparisons, but if nearest_a and nearest_b are far from the
        // origin the distance may be large
        normal        = if approx::relative_eq!(nearest_a, nearest_b,
          max_relative = f64::default_max_relative() * 256.0,
          epsilon      = f64::default_epsilon() * 256.0)
        {
          let perp = if e0_is_nearest {
            // e0 is edge BC
            -triangle.perpendicular_bc()
          } else if e1_is_nearest {
            // e1 is edge AB
            -triangle.perpendicular_ab()
          } else {
            debug_assert!(e2_is_nearest);
            // e2 is edge AC
            -triangle.perpendicular_ca()
          };
          *perp
        } else {
          -a_to_b
        }.normalize();
      }
      return Proximity { distance, half_axis, midpoint, normal }
    }
    let a_to_b    = nearest_b - nearest_a;
    let distance  = a_to_b.magnitude();
    let half_axis = a_to_b * 0.5;
    let midpoint  = nearest_a + half_axis;
    let normal    = Unit3::normalize (-a_to_b);
    Proximity { distance, half_axis, midpoint, normal }
  }

  pub fn query_triangle_segment (
    triangle : &geometry::Triangle3 <f64>,
    segment  : &geometry::Segment3 <f64>
  ) -> Self {
    use approx::{AbsDiffEq, RelativeEq};
    let (([s, t], nearest_a), (_r, nearest_b)) =
      geometry::distance::nearest_triangle3_segment3 (*triangle, *segment);
    if approx::relative_eq!(nearest_a, nearest_b,
      max_relative = f64::default_max_relative() * 2.0.powi (17),
      epsilon      = f64::default_epsilon() * 2.0.powi (17))
    {
      // segment and triangle intersect
      // we can check edge intersection by looking at s and t
      if let Some((mut proximity, edge_perp)) = if approx::abs_diff_eq!(*s, 0.0,
          epsilon = f64::default_epsilon() * 2.0.powi (23))
        {
          // intersect edge AC
          Some (( Proximity::query_segment_segment (&triangle.edge_ca(), segment),
            triangle.perpendicular_ca() ))
        } else if approx::abs_diff_eq!(*t, 0.0,
          epsilon = f64::default_epsilon() * 2.0.powi (23))
        {
          // intersect edge AB
          Some (( Proximity::query_segment_segment (&triangle.edge_ab(), segment),
            triangle.perpendicular_ab() ))
        } else if approx::relative_eq!(*s + *t, 1.0,
          max_relative = f64::default_max_relative() * 2.0.powi (23),
          epsilon = f64::default_epsilon() * 2.0.powi (23))
        {
          // intersect edge BC
          Some (( Proximity::query_segment_segment (&triangle.edge_bc(), segment),
            triangle.perpendicular_bc() ))
        } else {
          None
        }
      {
        if edge_perp.dot (*proximity.normal) > 0.0 {
          // ensure normal points from segment towards triangle
          proximity.normal = -proximity.normal;
        }
        return proximity
      }
      // the penetration depth is the shortest distance to move the segment to one of
      // the triangle edges, or else to move the segment perpendicularly to the triangle
      // surface
      let distance;
      let half_axis;
      let midpoint;
      let normal;
      let triangle_normal          = triangle.normal();
      let mut triangle_normal_unit = None;
      let segment_vec              = segment.vector();
      if approx::abs_diff_eq!(triangle_normal.dot (*segment_vec), 0.0) {
        // segment and triangle are parallel:
        // penetration depth is zero and direction is perpendicular to triangle face
        distance  = 0.0;
        half_axis = Vector3::zero();
        midpoint  = nearest_a;
        normal    = triangle_normal.normalize();
      } else {
        // segment intersects triangle at a single point
        let (seg_depth_sq, seg_depth_vec) = if triangle_normal.cross (*segment_vec)
          == Vector3::zero()
        {
          // segment is perpendicular to triangle: just use endpoint distances
          let depth_vec_a = segment.point_a() - nearest_b;
          let depth_vec_b = segment.point_b() - nearest_b;
          let dsq_a = depth_vec_a.magnitude_squared();
          let dsq_b = depth_vec_b.magnitude_squared();
          if dsq_a < dsq_b {
            (dsq_a, depth_vec_a)
          } else {
            (dsq_b, depth_vec_b)
          }
        } else {
          // segment intersects at an angle
          // project onto triangle plane
          let normal_unit      = triangle_normal.normalize();
          triangle_normal_unit = Some (normal_unit);
          let tri0_to_seg0     = segment.point_a() - triangle.point_a();
          let tri0_to_seg1     = segment.point_b() - triangle.point_a();
          let normal_dot_seg0  = normal_unit.dot (tri0_to_seg0);
          let normal_dot_seg1  = normal_unit.dot (tri0_to_seg1);
          let proj_seg0        = segment.point_a() - *normal_unit * normal_dot_seg0;
          let proj_seg1        = segment.point_b() - *normal_unit * normal_dot_seg1;
          // convert projected points to barycentric coordinates
          let ab               = triangle.point_b() - triangle.point_a();
          let ac               = triangle.point_c() - triangle.point_a();
          let ab_mag_sq        = ab.magnitude_squared();
          let ab_dot_ac        = ab.dot (ac);
          let ac_mag_sq        = ac.magnitude_squared();
          //let determinant     = ab_mag_sq * ac_mag_sq - ab_dot_ac.squared();
          let determinant      = ab_mag_sq.mul_add (ac_mag_sq,  -ab_dot_ac.squared());
          // barycentric coords: b0 = 1 - s - t, b1 = s, b2 = t
          let barycentric_to_cartesian = |s, t| triangle.point_a() + ab * s + ac * t;
          let barycentric_coords = |p_proj : Point3 <f64>| {
            let to_p = p_proj - triangle.point_a();
            let to_p_dot_ab = to_p.dot (ab);
            let to_p_dot_ac = to_p.dot (ac);
            //let s = (ac_mag_sq * to_p_dot_ab - ab_dot_ac * to_p_dot_ac) / determinant;
            let s = ac_mag_sq.mul_add (to_p_dot_ab,  -ab_dot_ac * to_p_dot_ac)
              / determinant;
            //let t = (ab_mag_sq * to_p_dot_ac - ab_dot_ac * to_p_dot_ab) / determinant;
            let t = ab_mag_sq.mul_add (to_p_dot_ac, -ab_dot_ac * to_p_dot_ab)
              / determinant;
            if cfg!(debug_assertions) {
              approx::assert_relative_eq!(p_proj, barycentric_to_cartesian (s, t),
                epsilon = f64::default_epsilon() * 2.0.powi (40),
                max_relative = f64::default_max_relative() * 2.0.powi (40));
            }
            [s, t]
          };
          let [s0, t0] = barycentric_coords (proj_seg0);
          let [s1, t1] = barycentric_coords (proj_seg1);
          let (base_s, base_t) = (s0, t0);
          let (vec_s, vec_t)   = (s1 - s0, t1 - t0);
          // returns clipped barycentric coordinates and parameter along segment
          let clip_seg_endpoint = |s, t| {
            // clip against each edge in barycentric space
            // describe the projected line segment in barycentric space as parameterized
            // line `base + u * vec`
            // edge AB: t
            // edge BC: 1 - s - t
            // edge CA: s
            // clip s == 0
            let (mut s, mut t) = (s, t);
            let mut u = 0.0;
            #[expect(clippy::useless_let_if_seq)]
            let mut clipped = false;
            if s < 0.0 {
              // solve: 0 = base_s + u * vec_s
              u = -base_s / vec_s;
              s = 0.0;
              t = base_t + u * vec_t;
              clipped = true;
            }
            if t < 0.0 {
              // solve: 0 = base_t + u * vec_t
              u = -base_t / vec_t;
              s = base_s + u * vec_s;
              t = 0.0;
              clipped = true;
            }
            if s + t > 1.0 {
              // solve: for u such that s + t == 1
              u = (1.0 - base_s - base_t) / (vec_s + vec_t);
              s = base_s + u * vec_s;
              t = base_t + u * vec_t;
              if cfg!(debug_assertions) {
                approx::assert_relative_eq!(s + t, 1.0,
                  max_relative = f64::default_max_relative() * 2.0.powi (16),
                  epsilon = f64::default_epsilon() * 2.0.powi (16));
              }
              debug_assert!(s >= 0.0);
              debug_assert!(t >= 0.0);
              clipped = true;
            }
            if clipped {
              Some (([s, t], u))
            } else {
              None
            }
          };
          let (dsq_a, depth_vec_a) =
            if let Some (([s, t], u)) = clip_seg_endpoint (s0, t0) {
              let clipped   = segment.point_a() + *segment_vec * u;
              let proj      = barycentric_to_cartesian (s, t);
              let depth_vec = clipped - proj;
              (depth_vec.magnitude_squared(), depth_vec)
            } else {
              let depth_vec = segment.point_a() - proj_seg0;
              (depth_vec.magnitude_squared(), depth_vec)
            };
          let (dsq_b, depth_vec_b) =
            if let Some (([s, t], u)) = clip_seg_endpoint (s1, t1) {
              let clipped   = segment.point_a() + *segment_vec * u;
              let proj      = barycentric_to_cartesian (s, t);
              let depth_vec = clipped - proj;
              (depth_vec.magnitude_squared(), depth_vec)
            } else {
              let depth_vec = segment.point_b() - proj_seg1;
              (depth_vec.magnitude_squared(), depth_vec)
            };
          if dsq_a < dsq_b {
            (dsq_a, depth_vec_a)
          } else {
            (dsq_b, depth_vec_b)
          }
        };
        // now compute the distance to nearest edge
        let (d0, d1, d2) = triangle_distance_to_edges (*triangle, *s, *t);
        let edge_dist_sq = f64::min (d0 * d0, f64::min (d1 * d1, d2 * d2));
        // choose minimum between perpendicular depth and edge distance
        if seg_depth_sq <= edge_dist_sq {
          // perpendicular movement is shorter
          let unit_normal = triangle_normal_unit
            .unwrap_or_else (|| triangle_normal.normalize());
          if unit_normal.dot (seg_depth_vec) < 0.0 {
            normal = -unit_normal;
          } else {
            normal = unit_normal;
          }
          distance  = -seg_depth_sq.sqrt();
          half_axis = *normal * -distance * 0.5;
          midpoint  = nearest_a + half_axis;
        } else {
          // moving to edge is shorter
          let edge_dist = edge_dist_sq.sqrt();
          distance      = -edge_dist;
          // find which edge is nearest
          let edge0     = triangle.point_b() - triangle.point_a();
          let edge1     = triangle.point_c() - triangle.point_a();
          let edge2     = triangle.point_c() - triangle.point_b();
          let len2_e0   = edge0.magnitude_squared();
          let len2_e1   = edge1.magnitude_squared();
          let len2_e2   = edge2.magnitude_squared();
          let (nearest_edge_point, nearest_edge_perp) = if d0 <= d1 && d0 <= d2 {
            let t = (nearest_b - triangle.point_b()).dot (edge0) / len2_e0;
            (triangle.point_b() + edge0 * t, -triangle.perpendicular_ab())
          } else if d1 <= d2 {
            let t = (nearest_b - triangle.point_a()).dot (edge1) / len2_e1;
            (triangle.point_a() + edge1 * t, -triangle.perpendicular_ca())
          } else {
            let t = (nearest_b - triangle.point_a()).dot (edge2) / len2_e2;
            (triangle.point_a() + edge2 * t, -triangle.perpendicular_bc())
          };
          half_axis = (nearest_b - nearest_edge_point) * 0.5;
          midpoint  = nearest_a + half_axis;
          normal    = nearest_edge_perp.normalize();
        }
      }
      return Proximity { distance, half_axis, midpoint, normal }
    }
    let a_to_b    = nearest_b - nearest_a;
    let distance  = a_to_b.magnitude();
    let half_axis = a_to_b * 0.5;
    let midpoint  = nearest_a + half_axis;
    let normal    = -a_to_b.normalize();
    Proximity { distance, half_axis, midpoint, normal }
  }

  pub fn query_triangle_point (
    triangle : &geometry::Triangle3 <f64>,
    point    : &Point3 <f64>
  ) -> Self {
    let ([s, t], nearest_a) = triangle.nearest_point (*point);
    let a_to_b    = point - nearest_a;
    let distance  = a_to_b.magnitude();
    let half_axis = a_to_b * 0.5;
    let midpoint  = nearest_a + half_axis;
    let normal    = if *point == nearest_a {
      // TODO: for efficiency should we just return the triangle normal in all cases ?
      let perp = if *t == 0.0 {
        // on edge ab
        -triangle.perpendicular_ab()
      } else if *s == 0.0 {
        // on edge ac
        -triangle.perpendicular_ca()
      } else if *s + *t == 1.0 {
        // on edge bc
        -triangle.perpendicular_bc()
      } else {
        // use the triangle normal
        triangle.normal()
      };
      *perp
    } else {
      -a_to_b
    }.normalize();
    Proximity { distance, half_axis, midpoint, normal }
  }

  pub fn query_line_line (
    line_a : &geometry::Segment3 <f64>,
    line_b : &geometry::Segment3 <f64>
  ) -> Self {
    use approx::{AbsDiffEq, RelativeEq};
    use num::Zero;
    let ((_, nearest_a), (_, nearest_b)) =
      geometry::distance::nearest_line3_line3 (*line_a, *line_b);
    let a_to_b    = nearest_b - nearest_a;
    let distance  = a_to_b.magnitude();
    let half_axis = a_to_b * 0.5;
    let midpoint  = nearest_a + half_axis;
    let normal    = if approx::relative_eq!(nearest_a, nearest_b,
      max_relative = f64::default_max_relative() * 256.0,
      epsilon      = f64::default_epsilon() * 256.0)
    {
      // intersection
      let line_a_dir = line_a.vector();
      let line_b_dir = line_b.vector();
      let mut perp   = line_a_dir.cross (*line_b_dir);
      if perp.is_zero() {
        perp = line_a_dir.orthogonal();
      }
      perp
    } else {
      -a_to_b
    }.normalize();
    Proximity { distance, half_axis, midpoint, normal }
  }

  pub fn query_line_segment (
    line    : &geometry::Segment3 <f64>,
    segment : &geometry::Segment3 <f64>
  ) -> Self {
    use approx::{AbsDiffEq, RelativeEq};
    let seg_dir   = segment.vector();
    let line_dir  = line.vector();
    let ((_, nearest_a), (s1, nearest_b)) =
      geometry::distance::nearest_line3_segment3 (*line, *segment);
    let a_to_b    = nearest_b - nearest_a;
    let distance  = a_to_b.magnitude();
    let half_axis = a_to_b * 0.5;
    let midpoint  = nearest_a + half_axis;
    let normal    = if approx::relative_eq!(nearest_a, nearest_b,
      max_relative = f64::default_max_relative() * 256.0,
      epsilon      = f64::default_epsilon() * 256.0)
    {
      // intersection
      if 0.0 < *s1 && *s1 < 1.0 {
        // intersection internal to segment
        let mut perp = seg_dir.cross (*line_dir);
        if perp.is_approx_zero() {
          // line and segment are parallel
          perp = seg_dir.orthogonal();
        }
        debug_assert!(!perp.is_approx_zero());
        perp
      } else {
        // when an endpoint is touching, choose the normal to point away from the
        // segment
        debug_assert!(*s1 == 0.0 || *s1 == 1.0);
        let mut perp = seg_dir.cross (*line_dir).cross (*line_dir);
        if perp.is_approx_zero() {
          // parallel
          perp = seg_dir.orthogonal()
        } else if *s1 == 1.0 {
          perp = -perp
        }
        debug_assert!((nearest_a - segment.point_a()).dot (perp) >= 0.0);
        debug_assert!((nearest_a - segment.point_b()).dot (perp) >= 0.0);
        perp
      }
    } else {
      -a_to_b
    }.normalize();
    Proximity { distance, half_axis, midpoint, normal }
  }

  pub fn query_line_point (
    line  : &geometry::Segment3 <f64>,
    point : &Point3 <f64>
  ) -> Self {
    let (_, nearest_a) = line.nearest_point (*point);
    let nearest_b      = *point;
    let a_to_b         = nearest_b - nearest_a;
    let distance       = a_to_b.magnitude();
    let half_axis      = a_to_b * 0.5;
    let midpoint       = nearest_a + half_axis;
    let normal         = if approx::relative_eq!(nearest_a, nearest_b) {
      *line.perpendicular()
    } else {
      -a_to_b
    }.normalize();
    Proximity { distance, half_axis, midpoint, normal }
  }

  pub fn query_segment_segment (
    segment_a : &geometry::Segment3 <f64>,
    segment_b : &geometry::Segment3 <f64>
  ) -> Self {
    use approx::{AbsDiffEq, RelativeEq};
    // method from Eberly GeometricTools:
    // <https://github.com/davideberly/GeometricTools/blob/a9a2b149485857273ea1f1cfa5416b392f495882/GTE/Mathematics/DistSegmentSegment.h>
    let vec_a = segment_a.vector();
    let vec_b = segment_b.vector();
    let vec_b0a0 = segment_a.point_a() - segment_b.point_a();
    let a = Positive::unchecked (*vec_a.self_dot());
    let b = vec_a.dot (*vec_b);
    let c = Positive::unchecked (*vec_b.self_dot());
    let d = vec_a.dot (vec_b0a0);
    let e = vec_b.dot (vec_b0a0);
    let f00 = d;
    let f10 = f00 + *a;
    let f01 = f00 - b;
    let f11 = f10 - b;
    let g00 = -e;
    let g10 = g00 - b;
    let g01 = g00 + *c;
    let g11 = g10 + *c;
    let svalue = [
      clamped_root (*a, f00, f10),
      clamped_root (*a, f01, f11)
    ];
    let mut classify = [i8::MAX; 2];
    for i in 0..2 {
      if svalue[i] <= 0.0 {
        classify[i] = -1
      } else if svalue[i] >= 1.0 {
        classify[i] = 1
      } else {
        classify[i] = 0
      }
    }
    let param_a;
    let param_b;
    if classify[0] == -1 && classify[1] == -1 {
      param_a = 0.0;
      param_b = clamped_root (*c, g00, g01);
    } else if classify[0] == 1 && classify[1] == 1 {
      param_a = 1.0;
      param_b = clamped_root (*c, g10, g11);
    } else {
      let mut edge = [u8::MAX; 2];
      let mut end  = [[f64::NAN; 2]; 2];
      if classify[0] == -1 {
        edge[0] = 0;
        end[0][0] = 0.0;
        end[0][1] = f00 / b;
        if end[0][1] < 0.0 || end[0][1] > 1.0 {
          end[0][1] = 0.5
        }
        if classify[1] == 0 {
          edge[1] = 3;
          end[1][0] = svalue[1];
          end[1][1] = 1.0;
        } else {
          debug_assert_eq!(classify[1], 1);
          edge[1] = 1;
          end[1][0] = 1.0;
          end[1][1] = f10 / b;
          if end[1][1] < 0.0 || end[1][1] > 1.0 {
            end[1][1] = 0.5
          }
        }
      } else if classify[0] == 0 {
        edge[0] = 2;
        end[0][0] = svalue[0];
        end[0][1] = 0.0;
        if classify[1] == -1 {
          edge[1] = 0;
          end[1][0] = 0.0;
          end[1][1] = f00 / b;
          if end[1][1] < 0.0 || end[1][1] > 1.0 {
            end[1][1] = 0.5
          }
        } else if classify[1] == 0 {
          edge[1] = 3;
          end[1][0] = svalue[1];
          end[1][1] = 1.0;
        } else {
          debug_assert_eq!(classify[1], 1);
          edge[1] = 1;
          end[1][0] = 1.0;
          end[1][1] = f10 / b;
          if end[1][1] < 0.0 || end[1][1] > 1.0 {
            end[1][1] = 0.5
          }
        }
      } else {
        debug_assert_eq!(classify[0], 1);
        edge[0] = 1;
        end[0][0] = 1.0;
        end[0][1] = f10 / b;
        if end[0][1] < 0.0 || end[0][1] > 1.0 {
          end[0][1] = 0.5
        }
        if classify[1] == 0 {
          edge[1] = 3;
          end[1][0] = svalue[1];
          end[1][1] = 1.0;
        } else {
          edge[1] = 0;
          end[1][0] = 0.0;
          end[1][1] = f00 / b;
          if end[1][1] < 0.0 || end[1][1] > 1.0 {
            end[1][1] = 0.5
          }
        }
      }
      // compute minimum parameters
      let delta = end[1][1] - end[0][1];
      //let h0 = delta * (-b * end[0][0] + c * end[0][1] - e);
      let h0 = delta * ((-b).mul_add (end[0][0], c.mul_add (end[0][1], -e)));
      if h0 >= 0.0 {
        param_a = (end[0][0] + end[1][0]) / 2.0;
        param_b = (end[0][1] + end[1][1]) / 2.0;
      } else {
        //let h1 = delta * (-b * end[1][0] + c * end[1][1] - e);
        let h1 = delta * ((-b).mul_add (end[1][0], c.mul_add (end[1][1], -e)));
        if h1 <= 0.0 {
          if edge[1] == 0 {
            param_a = 0.0;
            param_b = clamped_root (*c, g00, g01);
          } else if edge[1] == 1 {
            param_a = 1.0;
            param_b = clamped_root (*c, g10, g11);
          } else {
            param_a = end[1][0];
            param_b = end[1][1];
          }
        } else {
          debug_assert!(h0 < 0.0);
          debug_assert!(h1 > 0.0);
          let z = (h0 / (h0 - h1)).clamp (0.0, 1.0);
          let omz = 1.0 - z;
          param_a = omz * end[0][0] + z * end[1][0];
          param_b = omz * end[0][1] + z * end[1][1];
        }
      }
      debug_assert_ne!(edge[0], u8::MAX);
      debug_assert_ne!(edge[1], u8::MAX);
      debug_assert_ne!(end[0][0], f64::NAN);
      debug_assert_ne!(end[0][1], f64::NAN);
      debug_assert_ne!(end[1][0], f64::NAN);
      debug_assert_ne!(end[1][1], f64::NAN);
    }
    let nearest_a = Point3::from (segment_a.point_a().0 * (1.0 - param_a)
      + segment_a.point_b().0 * param_a);
    let nearest_b = Point3::from (segment_b.point_a().0 * (1.0 - param_b)
      + segment_b.point_b().0 * param_b);
    let a_to_b    = nearest_b - nearest_a;
    let distance  = a_to_b.magnitude();
    let half_axis = a_to_b * 0.5;
    let midpoint  = nearest_a + half_axis;
    let normal    = if approx::relative_eq!(nearest_a, nearest_b,
      max_relative = f64::default_max_relative() * 2.0.powi (20),
      epsilon      = f64::default_epsilon() * 2.0.powi (20))
    {
      // generate a normal vector perpendicular to the edges
      let mut perp = vec_a.cross (*vec_b);
      if perp == Vector3::zero() {
        // edges are parallel
        perp = vec_a.cross (Vector3::unit_z());
        if perp == Vector3::zero() {
          // edges are parallel to the Z axis, use the X axis instead
          perp = vec_a.cross (Vector3::unit_x());
        }
      }
      Unit3::normalize (perp)
    } else {
      Unit3::normalize (-a_to_b)
    };
    Proximity { distance, half_axis, midpoint, normal }
  }

  pub fn query_segment_point (
    segment : &geometry::Segment3 <f64>,
    point   : &Point3 <f64>
  ) -> Self {
    let (s, nearest_a) = segment.nearest_point (*point);
    let nearest_b      = *point;
    let a_to_b         = nearest_b - nearest_a;
    let distance       = a_to_b.magnitude();
    let half_axis      = a_to_b * 0.5;
    let midpoint       = nearest_a + half_axis;
    let normal         = if approx::relative_eq!(nearest_a, nearest_b) {
      let perp = if *s == 0.0 {
        debug_assert_eq!(nearest_a, segment.point_a());
        -segment.vector()
      } else if *s == 1.0 {
        segment.vector()
      } else {
        segment.perpendicular()
      };
      *perp
    } else {
      -a_to_b
    }.normalize();
    Proximity { distance, half_axis, midpoint, normal }
  }
}

/// Given a point on a triangle defined as s(b - a) + t(c - a), compute the distance to
/// each edge of the triangle
fn triangle_distance_to_edges (triangle : geometry::Triangle3 <f64>, s : f64, t : f64)
  -> (f64, f64, f64)
{
  use approx::{AbsDiffEq, RelativeEq};
  debug_assert!(s + t <= 1.0);
  // barycentric coordinates of each vertex
  let b0        = 1.0 - s - t;
  let b1        = s;
  let b2        = t;
  // edges opposite to each vertex
  let edge0     = triangle.point_c() - triangle.point_b();
  let edge1     = triangle.point_c() - triangle.point_a();
  let edge2     = triangle.point_b() - triangle.point_a();
  let len_e0_sq = edge0.magnitude_squared();
  let len_e1_sq = edge1.magnitude_squared();
  let len_e2_sq = edge2.magnitude_squared();
  // length of triangle normal is 2 * area
  let area_2_sq = triangle.normal().magnitude_squared();
  // convert barycentric to trilinear coordinates
  let d0        = b0 * (area_2_sq / len_e0_sq).sqrt();
  let d1        = b1 * (area_2_sq / len_e1_sq).sqrt();
  let d2        = b2 * (area_2_sq / len_e2_sq).sqrt();
  approx::assert_relative_eq!(
    triangle.barycentric ([b0, b1, b2].into()),
    triangle.trilinear ([d0, d1, d2].into()),
    max_relative = f64::default_max_relative() * 32.0,
    epsilon = f64::default_epsilon() * 32.0);
  (d0, d1, d2)
}

fn clamped_root (slope : f64, h0 : f64, h1 : f64) -> f64 {
  // helper function for Eberly GeometricTools methods
  if h0 < 0.0 {
    if h1 > 0.0 {
      let mut root = -h0 / slope;
      if root > 1.0 {
        root = 0.5
      }
      root
    } else {
      1.0
    }
  } else {
    0.0
  }
}

#[cfg(test)]
mod tests {
  use approx::AbsDiffEq;
  use rand;
  use rand_distr;
  use rand_xorshift::XorShiftRng;
  use std::f64::consts::*;
  use super::*;

  #[test]
  fn segment_point() {
    use rand::SeedableRng;
    use rand_distr::Distribution;
    // point above segment
    let segment = geometry::Segment3::noisy (
      [-1.0, 0.0, 0.0].into(),
      [ 1.0, 0.0, 0.0].into());
    let point = point3 (0.0, 1.0, 0.0);
    let proximity = Proximity::query_segment_point (&segment, &point);
    assert_eq!(proximity.nearest_a(), [0.0, 0.0, 0.0].into());
    assert_eq!(proximity.nearest_b(), [0.0, 1.0, 0.0].into());
    assert_eq!(proximity.distance, 1.0);
    assert_eq!(proximity.half_axis, [0.0, 0.5, 0.0].into());
    assert_eq!(proximity.midpoint, [0.0, 0.5, 0.0].into());
    assert_eq!(proximity.normal, -Unit3::axis_y());
    // point on segment: middle
    let segment = geometry::Segment3::noisy (
      [-1.0, 0.0, 0.0].into(),
      [ 1.0, 0.0, 0.0].into());
    let point = point3 (0.0, 0.0, 0.0);
    let proximity = Proximity::query_segment_point (&segment, &point);
    assert_eq!(proximity.nearest_a(), [0.0, 0.0, 0.0].into());
    assert_eq!(proximity.nearest_b(), [0.0, 0.0, 0.0].into());
    assert_eq!(proximity.distance, 0.0);
    assert_eq!(proximity.half_axis, [0.0, 0.0, 0.0].into());
    assert_eq!(proximity.midpoint, [0.0, 0.0, 0.0].into());
    assert_eq!(proximity.normal, Unit3::axis_z());
    // point on segment: end
    let segment = geometry::Segment3::noisy (
      [-1.0, 0.0, 0.0].into(),
      [ 1.0, 0.0, 0.0].into());
    let point = point3 (1.0, 0.0, 0.0);
    let proximity = Proximity::query_segment_point (&segment, &point);
    assert_eq!(proximity.nearest_a(), [1.0, 0.0, 0.0].into());
    assert_eq!(proximity.nearest_b(), [1.0, 0.0, 0.0].into());
    assert_eq!(proximity.distance, 0.0);
    assert_eq!(proximity.half_axis, [0.0, 0.0, 0.0].into());
    assert_eq!(proximity.midpoint, [1.0, 0.0, 0.0].into());
    assert_eq!(proximity.normal, segment.vector().normalize());
    // test that normal is invariant under translation
    let segment = geometry::Segment3::noisy (
      [-1.0, 0.0, 0.0].into(),
      [ 1.0, 0.0, 0.0].into());
    let point = point3 (0.0, 0.0, 0.0);
    let mut rng = XorShiftRng::seed_from_u64 (0);
    let std_normal = rand_distr::StandardNormal;
    // cauchy sample
    let randf = |rng : &mut XorShiftRng| {
      let x : f64 = std_normal.sample (rng);
      let y : f64 = std_normal.sample (rng);
      x / y
    };
    for _ in 0..10000 {
      let translation = vector3 (randf (&mut rng), randf (&mut rng), randf (&mut rng));
      let mut segment = segment;
      segment.translate (translation);
      let point = point + translation;
      let proximity = Proximity::query_segment_point (&segment, &point);
      assert_eq!(proximity.normal, Unit3::axis_z());
    }
  }

  #[test]
  fn segment_segment() {
    use approx::RelativeEq;
    use rand::SeedableRng;
    use rand_distr::Distribution;
    // parallel
    let segment1 = geometry::Segment3::noisy (
      [-2.0, 2.0, 0.0].into(),
      [ 2.0, 2.0, 0.0].into());
    let segment2 = geometry::Segment3::noisy (
      [-1.0, 0.0, 0.0].into(),
      [ 0.0, 0.0, 0.0].into());
    let proximity = Proximity::query_segment_segment (&segment1, &segment2);
    assert_eq!(proximity.nearest_a(), [-0.5, 2.0, 0.0].into());
    assert_eq!(proximity.nearest_b(), [-0.5, 0.0, 0.0].into());
    assert_eq!(proximity.distance, 2.0);
    assert_eq!(proximity.half_axis, [0.0, -1.0, 0.0].into());
    assert_eq!(proximity.midpoint, [-0.5, 1.0, 0.0].into());
    assert_eq!(proximity.normal, Unit3::noisy ([0.0, 1.0, 0.0].into()));
    let segment1 = geometry::Segment3::noisy (
      [-2.0, 2.0, 0.0].into(),
      [ 2.0, 2.0, 0.0].into());
    let segment2 = geometry::Segment3::noisy (
      [ 0.0, 0.0, 0.0].into(),
      [-3.0, 0.0, 0.0].into());
    let proximity = Proximity::query_segment_segment (&segment1, &segment2);
    assert_eq!(proximity.nearest_a(), [-1.0, 2.0, 0.0].into());
    assert_eq!(proximity.nearest_b(), [-1.0, 0.0, 0.0].into());
    assert_eq!(proximity.distance, 2.0);
    assert_eq!(proximity.half_axis, [0.0, -1.0, 0.0].into());
    assert_eq!(proximity.midpoint, [-1.0, 1.0, 0.0].into());
    assert_eq!(proximity.normal, Unit3::noisy ([0.0, 1.0, 0.0].into()));
    let segment1 = geometry::Segment3::noisy (
      [-2.0, 2.0, 0.0].into(),
      [ 2.0, 2.0, 0.0].into());
    let segment2 = geometry::Segment3::noisy (
      [ 0.0, 0.0, 0.0].into(),
      [ 3.0, 0.0, 0.0].into());
    let proximity = Proximity::query_segment_segment (&segment1, &segment2);
    assert_eq!(proximity.nearest_a(), [ 1.0, 2.0, 0.0].into());
    assert_eq!(proximity.nearest_b(), [ 1.0, 0.0, 0.0].into());
    assert_eq!(proximity.distance, 2.0);
    assert_eq!(proximity.half_axis, [0.0, -1.0, 0.0].into());
    assert_eq!(proximity.midpoint, [ 1.0, 1.0, 0.0].into());
    assert_eq!(proximity.normal, Unit3::noisy ([0.0, 1.0, 0.0].into()));
    let segment1 = geometry::Segment3::noisy (
      [-2.0, 2.0, 0.0].into(),
      [ 2.0, 2.0, 0.0].into());
    let segment2 = geometry::Segment3::noisy (
      [ 3.0, 0.0, 0.0].into(),
      [ 0.0, 0.0, 0.0].into());
    let proximity = Proximity::query_segment_segment (&segment1, &segment2);
    assert_eq!(proximity.nearest_a(), [ 1.0, 2.0, 0.0].into());
    assert_eq!(proximity.nearest_b(), [ 1.0, 0.0, 0.0].into());
    assert_eq!(proximity.distance, 2.0);
    assert_eq!(proximity.half_axis, [0.0, -1.0, 0.0].into());
    assert_eq!(proximity.midpoint, [ 1.0, 1.0, 0.0].into());
    assert_eq!(proximity.normal, Unit3::noisy ([0.0, 1.0, 0.0].into()));
    let segment1 = geometry::Segment3::noisy (
      [-2.0, 2.0, 0.0].into(),
      [ 2.0, 2.0, 0.0].into());
    let segment2 = geometry::Segment3::noisy (
      [-3.0, 0.0, 0.0].into(),
      [ 0.0, 0.0, 0.0].into());
    let proximity = Proximity::query_segment_segment (&segment1, &segment2);
    assert_eq!(proximity.nearest_a(), [0.0, 2.0, 0.0].into());
    assert_eq!(proximity.nearest_b(), [0.0, 0.0, 0.0].into());
    assert_eq!(proximity.distance, 2.0);
    assert_eq!(proximity.half_axis, [0.0, -1.0, 0.0].into());
    assert_eq!(proximity.midpoint, [0.0, 1.0, 0.0].into());
    assert_eq!(proximity.normal, Unit3::noisy ([0.0, 1.0, 0.0].into()));
    let segment1 = geometry::Segment3::noisy (
      [-2.0, 2.0, 0.0].into(),
      [ 2.0, 2.0, 0.0].into());
    let segment2 = geometry::Segment3::noisy (
      [ 1.0, -1.0, 0.0].into(),
      [ 1.0,  1.0, 0.0].into());
    let proximity = Proximity::query_segment_segment (&segment1, &segment2);
    assert_eq!(proximity.nearest_a(), [ 1.0, 2.0, 0.0].into());
    assert_eq!(proximity.nearest_b(), [ 1.0, 1.0, 0.0].into());
    assert_eq!(proximity.distance, 1.0);
    assert_eq!(proximity.half_axis, [0.0, -0.5, 0.0].into());
    assert_eq!(proximity.midpoint, [ 1.0, 1.5, 0.0].into());
    assert_eq!(proximity.normal, Unit3::noisy ([0.0, 1.0, 0.0].into()));
    let segment1 = geometry::Segment3::noisy (
      [-2.0, 2.0, 0.0].into(),
      [ 2.0, 2.0, 0.0].into());
    let segment2 = geometry::Segment3::noisy (
      [-3.0, 0.0, 0.0].into(),
      [-3.0, 3.0, 0.0].into());
    let proximity = Proximity::query_segment_segment (&segment1, &segment2);
    assert_eq!(proximity.nearest_a(), [-2.0, 2.0, 0.0].into());
    assert_eq!(proximity.nearest_b(), [-3.0, 2.0, 0.0].into());
    assert_eq!(proximity.distance, 1.0);
    assert_eq!(proximity.half_axis, [-0.5, 0.0, 0.0].into());
    assert_eq!(proximity.midpoint, [-2.5, 2.0, 0.0].into());
    assert_eq!(proximity.normal, Unit3::noisy ([1.0, 0.0, 0.0].into()));
    let segment1 = geometry::Segment3::noisy (
      [-2.0, 2.0, 0.0].into(),
      [ 2.0, 2.0, 0.0].into());
    let segment2 = geometry::Segment3::noisy (
      [ 3.0, 0.0, 0.0].into(),
      [ 3.0, 3.0, 0.0].into());
    let proximity = Proximity::query_segment_segment (&segment1, &segment2);
    assert_eq!(proximity.nearest_a(), [2.0, 2.0, 0.0].into());
    assert_eq!(proximity.nearest_b(), [3.0, 2.0, 0.0].into());
    assert_eq!(proximity.distance, 1.0);
    assert_eq!(proximity.half_axis, [0.5, 0.0, 0.0].into());
    assert_eq!(proximity.midpoint, [2.5, 2.0, 0.0].into());
    assert_eq!(proximity.normal, Unit3::noisy ([-1.0, 0.0, 0.0].into()));
    // normal should be invariant under translation
    let segment1 = geometry::Segment3::noisy (
      [ 0.0,  1.0, 0.0].into(),
      [ 0.0, -1.0, 0.0].into());
    let segment2 = geometry::Segment3::noisy (
      [ -0.5, -0.5, -0.5].into(),
      [  0.5,  0.5,  0.5].into());
    let proximity1 = Proximity::query_segment_segment (&segment1, &segment2);
    let mut rng = XorShiftRng::seed_from_u64 (0);
    let std_normal = rand_distr::StandardNormal;
    // cauchy sample
    let randf = |rng : &mut XorShiftRng| {
      let x : f64 = std_normal.sample (rng);
      let y : f64 = std_normal.sample (rng);
      x / y
    };
    for _ in 0..10000 {
      let translation = vector3 (randf (&mut rng), randf (&mut rng), randf (&mut rng));
      let mut segment1 = segment1;
      let mut segment2 = segment2;
      segment1.translate (translation);
      segment2.translate (translation);
      let proximity = Proximity::query_segment_segment (&segment1, &segment2);
      approx::assert_relative_eq!(*proximity.normal, *proximity1.normal,
        max_relative = f64::default_max_relative() * 2.0.powi (10))
    }
  }

  #[test]
  fn line_point() {
    use rand::SeedableRng;
    use rand_distr::Distribution;
    // point above line
    let line = geometry::Segment3::noisy (
      [-1.0, 0.0, 0.0].into(),
      [ 1.0, 0.0, 0.0].into());
    let point = point3 (0.0, 1.0, 0.0);
    let proximity = Proximity::query_line_point (&line, &point);
    assert_eq!(proximity.nearest_a(), [0.0, 0.0, 0.0].into());
    assert_eq!(proximity.nearest_b(), [0.0, 1.0, 0.0].into());
    assert_eq!(proximity.distance, 1.0);
    assert_eq!(proximity.half_axis, [0.0, 0.5, 0.0].into());
    assert_eq!(proximity.midpoint, [0.0, 0.5, 0.0].into());
    assert_eq!(proximity.normal, -Unit3::axis_y());
    // point on line: middle
    let line = geometry::Segment3::noisy (
      [-1.0, 0.0, 0.0].into(),
      [ 1.0, 0.0, 0.0].into());
    let point = point3 (0.0, 0.0, 0.0);
    let proximity = Proximity::query_line_point (&line, &point);
    assert_eq!(proximity.nearest_a(), [0.0, 0.0, 0.0].into());
    assert_eq!(proximity.nearest_b(), [0.0, 0.0, 0.0].into());
    assert_eq!(proximity.distance, 0.0);
    assert_eq!(proximity.half_axis, [0.0, 0.0, 0.0].into());
    assert_eq!(proximity.midpoint, [0.0, 0.0, 0.0].into());
    assert_eq!(proximity.normal, Unit3::axis_z());
    // point on line: end
    let line = geometry::Segment3::noisy (
      [-1.0, 0.0, 0.0].into(),
      [ 1.0, 0.0, 0.0].into());
    let point = point3 (1.0, 0.0, 0.0);
    let proximity = Proximity::query_line_point (&line, &point);
    assert_eq!(proximity.nearest_a(), [1.0, 0.0, 0.0].into());
    assert_eq!(proximity.nearest_b(), [1.0, 0.0, 0.0].into());
    assert_eq!(proximity.distance, 0.0);
    assert_eq!(proximity.half_axis, [0.0, 0.0, 0.0].into());
    assert_eq!(proximity.midpoint, [1.0, 0.0, 0.0].into());
    assert_eq!(proximity.normal, Unit3::axis_z());
    // test that normal is invariant under translation
    let line = geometry::Segment3::noisy (
      [-1.0, 0.0, 0.0].into(),
      [ 1.0, 0.0, 0.0].into());
    let point = point3 (0.0, 0.0, 0.0);
    let mut rng = XorShiftRng::seed_from_u64 (0);
    let std_normal = rand_distr::StandardNormal;
    // cauchy sample
    let randf = |rng : &mut XorShiftRng| {
      let x : f64 = std_normal.sample (rng);
      let y : f64 = std_normal.sample (rng);
      x / y
    };
    for _ in 0..10000 {
      let translation = vector3 (randf (&mut rng), randf (&mut rng), randf (&mut rng));
      let mut line = line;
      line.translate (translation);
      let point = point + translation;
      let proximity = Proximity::query_line_point (&line, &point);
      assert_eq!(proximity.normal, Unit3::axis_z());
    }
  }

  #[test]
  fn line_segment() {
    use approx::RelativeEq;
    use rand::SeedableRng;
    use rand_distr::Distribution;
    // parallel
    let line    = geometry::Segment3::noisy (
      [-2.0, 0.0, 0.0].into(),
      [ 2.0, 0.0, 0.0].into());
    let segment = geometry::Segment3::noisy (
      [-2.0, 2.0, 0.0].into(),
      [ 2.0, 2.0, 0.0].into());
    let proximity = Proximity::query_line_segment (&line, &segment);
    assert_eq!(proximity.nearest_a(), [-2.0, 0.0, 0.0].into());
    assert_eq!(proximity.nearest_b(), [-2.0, 2.0, 0.0].into());
    assert_eq!(proximity.distance, 2.0);
    assert_eq!(proximity.half_axis, [0.0, 1.0, 0.0].into());
    assert_eq!(proximity.midpoint, [-2.0, 1.0, 0.0].into());
    assert_eq!(proximity.normal, -Unit3::axis_y());
    let line    = geometry::Segment3::noisy (
      [-5.0, 0.0, 0.0].into(),
      [-2.0, 0.0, 0.0].into());
    let segment = geometry::Segment3::noisy (
      [-2.0, 2.0, 0.0].into(),
      [ 2.0, 2.0, 0.0].into());
    let proximity = Proximity::query_line_segment (&line, &segment);
    assert_eq!(proximity.nearest_a(), [-2.0, 0.0, 0.0].into());
    assert_eq!(proximity.nearest_b(), [-2.0, 2.0, 0.0].into());
    assert_eq!(proximity.distance, 2.0);
    assert_eq!(proximity.half_axis, [0.0, 1.0, 0.0].into());
    assert_eq!(proximity.midpoint, [-2.0, 1.0, 0.0].into());
    assert_eq!(proximity.normal, -Unit3::axis_y());
    // perpendicular
    let line    = geometry::Segment3::noisy (
      [0.0, -1.0, 0.0].into(),
      [0.0,  1.0, 0.0].into());
    let segment = geometry::Segment3::noisy (
      [1.0,  1.0, 0.0].into(),
      [3.0,  1.0, 0.0].into());
    let proximity = Proximity::query_line_segment (&line, &segment);
    assert_eq!(proximity.nearest_a(), [0.0, 1.0, 0.0].into());
    assert_eq!(proximity.nearest_b(), [1.0, 1.0, 0.0].into());
    assert_eq!(proximity.distance, 1.0);
    assert_eq!(proximity.half_axis, [0.5, 0.0, 0.0].into());
    assert_eq!(proximity.midpoint, [0.5, 1.0, 0.0].into());
    assert_eq!(proximity.normal, -Unit3::axis_x());
    let line    = geometry::Segment3::noisy (
      [ 0.0, -1.0, 0.0].into(),
      [ 0.0,  1.0, 0.0].into());
    let segment = geometry::Segment3::noisy (
      [-3.0,  1.0, 0.0].into(),
      [-1.0,  1.0, 0.0].into());
    let proximity = Proximity::query_line_segment (&line, &segment);
    assert_eq!(proximity.nearest_a(), [ 0.0, 1.0, 0.0].into());
    assert_eq!(proximity.nearest_b(), [-1.0, 1.0, 0.0].into());
    assert_eq!(proximity.distance, 1.0);
    assert_eq!(proximity.half_axis, [-0.5, 0.0, 0.0].into());
    assert_eq!(proximity.midpoint, [-0.5, 1.0, 0.0].into());
    assert_eq!(proximity.normal, Unit3::axis_x());
    // in front
    let line    = geometry::Segment3::noisy (
      [ 0.0, -1.0, 0.0].into(),
      [ 0.0,  1.0, 0.0].into());
    let segment = geometry::Segment3::noisy (
      [-2.0,  2.0, 1.0].into(),
      [ 2.0,  2.0, 1.0].into());
    let proximity = Proximity::query_line_segment (&line, &segment);
    assert_eq!(proximity.nearest_a(), [ 0.0, 2.0, 0.0].into());
    assert_eq!(proximity.nearest_b(), [ 0.0, 2.0, 1.0].into());
    assert_eq!(proximity.distance, 1.0);
    assert_eq!(proximity.half_axis, [0.0, 0.0, 0.5].into());
    assert_eq!(proximity.midpoint, [0.0, 2.0, 0.5].into());
    assert_eq!(proximity.normal, -Unit3::axis_z());
    // intersection: middle
    let line    = geometry::Segment3::noisy (
      [ 0.0, -1.0, 0.0].into(),
      [ 0.0,  1.0, 0.0].into());
    let segment = geometry::Segment3::noisy (
      [-2.0,  2.0, 0.0].into(),
      [ 2.0,  2.0, 0.0].into());
    let proximity = Proximity::query_line_segment (&line, &segment);
    assert_eq!(proximity.nearest_a(), [ 0.0, 2.0, 0.0].into());
    assert_eq!(proximity.nearest_b(), [ 0.0, 2.0, 0.0].into());
    assert_eq!(proximity.distance, 0.0);
    assert_eq!(proximity.half_axis, [0.0, 0.0, 0.0].into());
    assert_eq!(proximity.midpoint, [0.0, 2.0, 0.0].into());
    assert_eq!(proximity.normal, Unit3::axis_z());
    // intersection: endpoint A
    let line    = geometry::Segment3::noisy (
      [-2.0, -1.0, 0.0].into(),
      [-2.0,  1.0, 0.0].into());
    let segment = geometry::Segment3::noisy (
      [-2.0,  2.0, 0.0].into(),
      [ 2.0,  2.0, 0.0].into());
    let proximity = Proximity::query_line_segment (&line, &segment);
    assert_eq!(proximity.nearest_a(), [-2.0, 2.0, 0.0].into());
    assert_eq!(proximity.nearest_b(), [-2.0, 2.0, 0.0].into());
    assert_eq!(proximity.distance, 0.0);
    assert_eq!(proximity.half_axis, [0.0, 0.0, 0.0].into());
    assert_eq!(proximity.midpoint, [-2.0, 2.0, 0.0].into());
    assert_eq!(proximity.normal, -Unit3::axis_x());
    let line    = geometry::Segment3::noisy (
      [-2.0, -1.0, 0.0].into(),
      [-2.0,  1.0, 0.0].into());
    let segment = geometry::Segment3::noisy (
      [-2.0,   2.0, 0.0].into(),
      [ 2.0,  -2.0, 0.0].into());
    let proximity = Proximity::query_line_segment (&line, &segment);
    assert_eq!(proximity.nearest_a(), [-2.0, 2.0, 0.0].into());
    assert_eq!(proximity.nearest_b(), [-2.0, 2.0, 0.0].into());
    assert_eq!(proximity.distance, 0.0);
    assert_eq!(proximity.half_axis, [0.0, 0.0, 0.0].into());
    assert_eq!(proximity.midpoint, [-2.0, 2.0, 0.0].into());
    assert_eq!(proximity.normal, -Unit3::axis_x());
    // intersection: endpoint B
    let line    = geometry::Segment3::noisy (
      [ 2.0, -1.0, 0.0].into(),
      [ 2.0,  1.0, 0.0].into());
    let segment = geometry::Segment3::noisy (
      [-2.0,  2.0, 0.0].into(),
      [ 2.0,  2.0, 0.0].into());
    let proximity = Proximity::query_line_segment (&line, &segment);
    assert_eq!(proximity.nearest_a(), [2.0, 2.0, 0.0].into());
    assert_eq!(proximity.nearest_b(), [2.0, 2.0, 0.0].into());
    assert_eq!(proximity.distance, 0.0);
    assert_eq!(proximity.half_axis, [0.0, 0.0, 0.0].into());
    assert_eq!(proximity.midpoint, [2.0, 2.0, 0.0].into());
    assert_eq!(proximity.normal, Unit3::axis_x());
    let line    = geometry::Segment3::noisy (
      [ 2.0, -1.0, 0.0].into(),
      [ 2.0,  1.0, 0.0].into());
    let segment = geometry::Segment3::noisy (
      [-2.0,  2.0, 0.0].into(),
      [ 2.0,  0.0, 0.0].into());
    let proximity = Proximity::query_line_segment (&line, &segment);
    assert_eq!(proximity.nearest_a(), [2.0, 0.0, 0.0].into());
    assert_eq!(proximity.nearest_b(), [2.0, 0.0, 0.0].into());
    assert_eq!(proximity.distance, 0.0);
    assert_eq!(proximity.half_axis, [0.0, 0.0, 0.0].into());
    assert_eq!(proximity.midpoint, [2.0, 0.0, 0.0].into());
    assert_eq!(proximity.normal, Unit3::axis_x());
    // intersection: parallel
    let line    = geometry::Segment3::noisy (
      [-1.0, 2.0, 0.0].into(),
      [ 0.0, 2.0, 0.0].into());
    let segment = geometry::Segment3::noisy (
      [-2.0, 2.0, 0.0].into(),
      [ 2.0, 2.0, 0.0].into());
    let proximity = Proximity::query_line_segment (&line, &segment);
    assert_eq!(proximity.nearest_a(), [-2.0, 2.0, 0.0].into());
    assert_eq!(proximity.nearest_b(), [-2.0, 2.0, 0.0].into());
    assert_eq!(proximity.distance, 0.0);
    assert_eq!(proximity.half_axis, [0.0, 0.0, 0.0].into());
    assert_eq!(proximity.midpoint, [-2.0, 2.0, 0.0].into());
    assert_eq!(proximity.normal, Unit3::axis_z());
    // normal should be invariant under translation
    let line = geometry::Segment3::noisy (
      [ -0.5, -0.5, -0.5].into(),
      [  0.5,  0.5,  0.5].into());
    let segment = geometry::Segment3::noisy (
      [ 0.0,  1.0, 0.0].into(),
      [ 0.0, -1.0, 0.0].into());
    let proximity1 = Proximity::query_line_segment (&line, &segment);
    let mut rng = XorShiftRng::seed_from_u64 (0);
    let std_normal = rand_distr::StandardNormal;
    // cauchy sample
    let randf = |rng : &mut XorShiftRng| {
      let x : f64 = std_normal.sample (rng);
      let y : f64 = std_normal.sample (rng);
      x / y
    };
    for _ in 0..10000 {
      let translation = vector3 (randf (&mut rng), randf (&mut rng), randf (&mut rng));
      let mut line    = line;
      let mut segment = segment;
      line.translate (translation);
      segment.translate (translation);
      let proximity = Proximity::query_line_segment (&line, &segment);
      approx::assert_relative_eq!(*proximity.normal, *proximity1.normal,
        max_relative = f64::default_max_relative() * 2.0.powi (10))
    }
  }

  #[test]
  fn line_line() {
    use approx::RelativeEq;
    use rand::SeedableRng;
    use rand_distr::Distribution;
    // parallel
    let line_a    = geometry::Segment3::noisy (
      [-2.0, 0.0, 0.0].into(),
      [ 2.0, 0.0, 0.0].into());
    let line_b = geometry::Segment3::noisy (
      [-2.0, 2.0, 0.0].into(),
      [ 2.0, 2.0, 0.0].into());
    let proximity = Proximity::query_line_line (&line_a, &line_b);
    assert_eq!(proximity.nearest_a(), [-2.0, 0.0, 0.0].into());
    assert_eq!(proximity.nearest_b(), [-2.0, 2.0, 0.0].into());
    assert_eq!(proximity.distance, 2.0);
    assert_eq!(proximity.half_axis, [0.0, 1.0, 0.0].into());
    assert_eq!(proximity.midpoint, [-2.0, 1.0, 0.0].into());
    assert_eq!(proximity.normal, -Unit3::axis_y());
    let line_a    = geometry::Segment3::noisy (
      [-5.0, 0.0, 0.0].into(),
      [-2.0, 0.0, 0.0].into());
    let line_b = geometry::Segment3::noisy (
      [-2.0, 2.0, 0.0].into(),
      [ 2.0, 2.0, 0.0].into());
    let proximity = Proximity::query_line_line (&line_a, &line_b);
    assert_eq!(proximity.nearest_a(), [-2.0, 0.0, 0.0].into());
    assert_eq!(proximity.nearest_b(), [-2.0, 2.0, 0.0].into());
    assert_eq!(proximity.distance, 2.0);
    assert_eq!(proximity.half_axis, [0.0, 1.0, 0.0].into());
    assert_eq!(proximity.midpoint, [-2.0, 1.0, 0.0].into());
    assert_eq!(proximity.normal, -Unit3::axis_y());
    // perpendicular
    let line_a    = geometry::Segment3::noisy (
      [0.0, -1.0, 0.0].into(),
      [0.0,  1.0, 0.0].into());
    let line_b = geometry::Segment3::noisy (
      [1.0,  1.0, 0.0].into(),
      [3.0,  1.0, 0.0].into());
    let proximity = Proximity::query_line_line (&line_a, &line_b);
    assert_eq!(proximity.nearest_a(), [0.0, 1.0, 0.0].into());
    assert_eq!(proximity.nearest_b(), [0.0, 1.0, 0.0].into());
    assert_eq!(proximity.distance, 0.0);
    assert_eq!(proximity.half_axis, [0.0, 0.0, 0.0].into());
    assert_eq!(proximity.midpoint, [0.0, 1.0, 0.0].into());
    assert_eq!(proximity.normal, -Unit3::axis_z());
    let line_a    = geometry::Segment3::noisy (
      [ 0.0, -1.0, 0.0].into(),
      [ 0.0,  1.0, 0.0].into());
    let line_b = geometry::Segment3::noisy (
      [-3.0,  1.0, 0.0].into(),
      [-1.0,  1.0, 0.0].into());
    let proximity = Proximity::query_line_line (&line_a, &line_b);
    assert_eq!(proximity.nearest_a(), [0.0, 1.0, 0.0].into());
    assert_eq!(proximity.nearest_b(), [0.0, 1.0, 0.0].into());
    assert_eq!(proximity.distance, 0.0);
    assert_eq!(proximity.half_axis, [0.0, 0.0, 0.0].into());
    assert_eq!(proximity.midpoint, [0.0, 1.0, 0.0].into());
    assert_eq!(proximity.normal, -Unit3::axis_z());
    // in front
    let line_a    = geometry::Segment3::noisy (
      [ 0.0, -1.0, 0.0].into(),
      [ 0.0,  1.0, 0.0].into());
    let line_b = geometry::Segment3::noisy (
      [-2.0,  2.0, 1.0].into(),
      [ 2.0,  2.0, 1.0].into());
    let proximity = Proximity::query_line_line (&line_a, &line_b);
    assert_eq!(proximity.nearest_a(), [ 0.0, 2.0, 0.0].into());
    assert_eq!(proximity.nearest_b(), [ 0.0, 2.0, 1.0].into());
    assert_eq!(proximity.distance, 1.0);
    assert_eq!(proximity.half_axis, [0.0, 0.0, 0.5].into());
    assert_eq!(proximity.midpoint, [0.0, 2.0, 0.5].into());
    assert_eq!(proximity.normal, -Unit3::axis_z());
    // intersection: middle
    let line_a    = geometry::Segment3::noisy (
      [ 0.0, -1.0, 0.0].into(),
      [ 0.0,  1.0, 0.0].into());
    let line_b = geometry::Segment3::noisy (
      [-2.0,  2.0, 0.0].into(),
      [ 2.0,  2.0, 0.0].into());
    let proximity = Proximity::query_line_line (&line_a, &line_b);
    assert_eq!(proximity.nearest_a(), [ 0.0, 2.0, 0.0].into());
    assert_eq!(proximity.nearest_b(), [ 0.0, 2.0, 0.0].into());
    assert_eq!(proximity.distance, 0.0);
    assert_eq!(proximity.half_axis, [0.0, 0.0, 0.0].into());
    assert_eq!(proximity.midpoint, [0.0, 2.0, 0.0].into());
    assert_eq!(proximity.normal, -Unit3::axis_z());
    // intersection: endpoint A
    let line_a    = geometry::Segment3::noisy (
      [-2.0, -1.0, 0.0].into(),
      [-2.0,  1.0, 0.0].into());
    let line_b = geometry::Segment3::noisy (
      [-2.0,  2.0, 0.0].into(),
      [ 2.0,  2.0, 0.0].into());
    let proximity = Proximity::query_line_line (&line_a, &line_b);
    assert_eq!(proximity.nearest_a(), [-2.0, 2.0, 0.0].into());
    assert_eq!(proximity.nearest_b(), [-2.0, 2.0, 0.0].into());
    assert_eq!(proximity.distance, 0.0);
    assert_eq!(proximity.half_axis, [0.0, 0.0, 0.0].into());
    assert_eq!(proximity.midpoint, [-2.0, 2.0, 0.0].into());
    assert_eq!(proximity.normal, -Unit3::axis_z());
    let line_a    = geometry::Segment3::noisy (
      [-2.0, -1.0, 0.0].into(),
      [-2.0,  1.0, 0.0].into());
    let line_b = geometry::Segment3::noisy (
      [-2.0,   2.0, 0.0].into(),
      [ 2.0,  -2.0, 0.0].into());
    let proximity = Proximity::query_line_line (&line_a, &line_b);
    assert_eq!(proximity.nearest_a(), [-2.0, 2.0, 0.0].into());
    assert_eq!(proximity.nearest_b(), [-2.0, 2.0, 0.0].into());
    assert_eq!(proximity.distance, 0.0);
    assert_eq!(proximity.half_axis, [0.0, 0.0, 0.0].into());
    assert_eq!(proximity.midpoint, [-2.0, 2.0, 0.0].into());
    assert_eq!(proximity.normal, -Unit3::axis_z());
    // intersection: endpoint B
    let line_a    = geometry::Segment3::noisy (
      [ 2.0, -1.0, 0.0].into(),
      [ 2.0,  1.0, 0.0].into());
    let line_b = geometry::Segment3::noisy (
      [-2.0,  2.0, 0.0].into(),
      [ 2.0,  2.0, 0.0].into());
    let proximity = Proximity::query_line_line (&line_a, &line_b);
    assert_eq!(proximity.nearest_a(), [2.0, 2.0, 0.0].into());
    assert_eq!(proximity.nearest_b(), [2.0, 2.0, 0.0].into());
    assert_eq!(proximity.distance, 0.0);
    assert_eq!(proximity.half_axis, [0.0, 0.0, 0.0].into());
    assert_eq!(proximity.midpoint, [2.0, 2.0, 0.0].into());
    assert_eq!(proximity.normal, -Unit3::axis_z());
    let line_a    = geometry::Segment3::noisy (
      [ 2.0, -1.0, 0.0].into(),
      [ 2.0,  1.0, 0.0].into());
    let line_b = geometry::Segment3::noisy (
      [-2.0,  2.0, 0.0].into(),
      [ 2.0,  0.0, 0.0].into());
    let proximity = Proximity::query_line_line (&line_a, &line_b);
    assert_eq!(proximity.nearest_a(), [2.0, 0.0, 0.0].into());
    assert_eq!(proximity.nearest_b(), [2.0, 0.0, 0.0].into());
    assert_eq!(proximity.distance, 0.0);
    assert_eq!(proximity.half_axis, [0.0, 0.0, 0.0].into());
    assert_eq!(proximity.midpoint, [2.0, 0.0, 0.0].into());
    assert_eq!(proximity.normal, -Unit3::axis_z());
    // intersection: parallel
    let line_a    = geometry::Segment3::noisy (
      [-1.0, 2.0, 0.0].into(),
      [ 0.0, 2.0, 0.0].into());
    let line_b = geometry::Segment3::noisy (
      [-2.0, 2.0, 0.0].into(),
      [ 2.0, 2.0, 0.0].into());
    let proximity = Proximity::query_line_line (&line_a, &line_b);
    assert_eq!(proximity.nearest_a(), [-2.0, 2.0, 0.0].into());
    assert_eq!(proximity.nearest_b(), [-2.0, 2.0, 0.0].into());
    assert_eq!(proximity.distance, 0.0);
    assert_eq!(proximity.half_axis, [0.0, 0.0, 0.0].into());
    assert_eq!(proximity.midpoint, [-2.0, 2.0, 0.0].into());
    assert_eq!(proximity.normal, Unit3::axis_z());
    // normal should be invariant under translation
    let line_a = geometry::Segment3::noisy (
      [ -0.5, -0.5, -0.5].into(),
      [  0.5,  0.5,  0.5].into());
    let line_b = geometry::Segment3::noisy (
      [ 0.0,  1.0, 0.0].into(),
      [ 0.0, -1.0, 0.0].into());
    let proximity1 = Proximity::query_line_line (&line_a, &line_b);
    let mut rng = XorShiftRng::seed_from_u64 (0);
    let std_normal = rand_distr::StandardNormal;
    // cauchy sample
    let randf = |rng : &mut XorShiftRng| {
      let x : f64 = std_normal.sample (rng);
      let y : f64 = std_normal.sample (rng);
      x / y
    };
    for _ in 0..10000 {
      let translation = vector3 (randf (&mut rng), randf (&mut rng), randf (&mut rng));
      let mut line_a = line_a;
      let mut line_b = line_b;
      line_a.translate (translation);
      line_b.translate (translation);
      let proximity = Proximity::query_line_line (&line_a, &line_b);
      approx::assert_relative_eq!(*proximity.normal, *proximity1.normal,
        max_relative = f64::default_max_relative() * 2.0.powi (10))
    }
  }

  #[test]
  fn triangle_point() {
    use approx::RelativeEq;
    use rand::SeedableRng;
    use rand_distr::Distribution;
    let triangle = geometry::Triangle3::noisy (
      [-1.0, -1.0, 0.0].into(),
      [ 1.0, -1.0, 0.0].into(),
      [ 0.0,  1.0, 0.0].into());
    let point = point3 (0.0, 0.0, 0.0);
    let proximity = Proximity::query_triangle_point (&triangle, &point);
    assert_eq!(proximity.nearest_a(), [0.0, 0.0, 0.0].into());
    assert_eq!(proximity.nearest_b(), [0.0, 0.0, 0.0].into());
    assert_eq!(proximity.distance, 0.0);
    assert_eq!(proximity.half_axis, [0.0, 0.0, 0.0].into());
    assert_eq!(proximity.midpoint, [0.0, 0.0, 0.0].into());
    assert_eq!(proximity.normal, Unit3::axis_z());
    let triangle = geometry::Triangle3::noisy (
      [-1.0, -1.0, 0.0].into(),
      [ 1.0, -1.0, 0.0].into(),
      [ 0.0,  1.0, 0.0].into());
    let point = point3 (0.5, 0.0, 0.0);
    let proximity = Proximity::query_triangle_point (&triangle, &point);
    assert_eq!(proximity.nearest_a(), [0.5, 0.0, 0.0].into());
    assert_eq!(proximity.nearest_b(), [0.5, 0.0, 0.0].into());
    assert_eq!(proximity.distance, 0.0);
    assert_eq!(proximity.half_axis, [0.0, 0.0, 0.0].into());
    assert_eq!(proximity.midpoint, [0.5, 0.0, 0.0].into());
    assert_eq!(proximity.normal, Unit3::normalize ([-2.0, -1.0, 0.0].into()));
    let triangle = geometry::Triangle3::noisy (
      [ 0.0, -2.0, 3.5].into(),
      [ 2.0,  2.0, 3.5].into(),
      [-2.0,  2.0, 3.5].into());
    let point = point3 (0.0, 0.0, 2.0);
    let proximity = Proximity::query_triangle_point (&triangle, &point);
    assert_eq!(proximity.distance, 1.5);
    assert_eq!(proximity.half_axis, [0.0, 0.0, -0.75].into());
    assert_eq!(proximity.midpoint, [0.0, 0.0, 2.75].into());
    assert_eq!(proximity.normal, Unit3::axis_z());
    assert_eq!(proximity.nearest_a(), [0.0, 0.0, 3.5].into());
    assert_eq!(proximity.nearest_b(), [0.0, 0.0, 2.0].into());
    // test that normal is invariant under translation
    let triangle = geometry::Triangle3::noisy (
      [-1.0, -1.0, 0.0].into(),
      [ 1.0, -1.0, 0.0].into(),
      [ 0.0,  1.0, 0.0].into());
    let point = point3 (1.0, 0.0, 0.0);
    let proximity1 = Proximity::query_triangle_point (&triangle, &point);
    let mut rng = XorShiftRng::seed_from_u64 (0);
    let std_normal = rand_distr::StandardNormal;
    // cauchy sample
    let randf = |rng : &mut XorShiftRng| {
      let x : f64 = std_normal.sample (rng);
      let y : f64 = std_normal.sample (rng);
      x / y
    };
    for _ in 0..10000 {
      let translation = vector3 (randf (&mut rng), randf (&mut rng), randf (&mut rng));
      let mut triangle = triangle;
      triangle.translate (translation);
      let point = point + translation;
      let proximity = Proximity::query_triangle_point (&triangle, &point);
      approx::assert_relative_eq!(*proximity.normal, *proximity1.normal,
        max_relative = f64::default_max_relative() * 2.0.powi (12));
    }
  }

  #[test]
  fn triangle_segment() {
    use approx::RelativeEq;
    use rand::SeedableRng;
    use rand_distr::Distribution;
    let mut triangle = geometry::Triangle3::noisy (
      [-1.0, -1.0, 0.0].into(),
      [ 1.0, -1.0, 0.0].into(),
      [ 0.0,  1.0, 0.0].into());
    let mut segment = geometry::Segment3::noisy (
      [-2.0, -2.0, -0.01].into(),
      [ 2.0,  2.0,  0.01].into());
    let proximity = Proximity::query_triangle_segment (&triangle, &segment);
    triangle.translate (proximity.half_axis);
    segment.translate (-proximity.half_axis);
    let proximity = Proximity::query_triangle_segment (&triangle, &segment);
    approx::assert_abs_diff_eq!(proximity.distance, 0.0,
      epsilon = f64::default_epsilon() * 2.0);
    let mut triangle = geometry::Triangle3::noisy (
      [-1.0, -1.0, 0.0].into(),
      [ 1.0, -1.0, 0.0].into(),
      [ 0.0,  1.0, 0.0].into());
    let mut segment = geometry::Segment3::noisy (
      [-2.0,  2.0, -10.0].into(),
      [ 2.0, -2.0,  10.0].into());
    let proximity = Proximity::query_triangle_segment (&triangle, &segment);
    triangle.translate (proximity.half_axis);
    segment.translate (-proximity.half_axis);
    let proximity = Proximity::query_triangle_segment (&triangle, &segment);
    approx::assert_abs_diff_eq!(proximity.distance, 0.0,
      epsilon = f64::default_epsilon() * 2.0);
    // normal should be invariant under translation
    let triangle = geometry::Triangle3::noisy (
      [-1.0, -1.0, 0.0].into(),
      [ 1.0, -1.0, 0.0].into(),
      [ 0.0,  1.0, 0.0].into());
    let segment = geometry::Segment3::noisy (
      [-1.5, 0.0, -2.0].into(),
      [ 2.0, 0.0,  2.0].into());
    let proximity1 = Proximity::query_triangle_segment (&triangle, &segment);
    println!("normal: {}", *proximity1.normal);
    let mut rng = XorShiftRng::seed_from_u64 (0);
    let std_normal = rand_distr::StandardNormal;
    // cauchy sample
    let randf = |rng : &mut XorShiftRng| {
      let x : f64 = std_normal.sample (rng);
      let y : f64 = std_normal.sample (rng);
      x / y
    };
    for _ in 0..10000 {
      let translation  = vector3 (randf (&mut rng), randf (&mut rng), randf (&mut rng));
      let mut triangle = triangle;
      let mut segment  = segment;
      triangle.translate (translation);
      segment.translate (translation);
      let proximity = Proximity::query_triangle_segment (&triangle, &segment);
      approx::assert_relative_eq!(*proximity.normal, *proximity1.normal,
        max_relative = f64::default_max_relative() * 2.0.powi (10))
    }
  }

  #[test]
  fn triangle_line() {
    use approx::RelativeEq;
    use rand::SeedableRng;
    use rand_distr::Distribution;
    let triangle = geometry::Triangle3::noisy (
      [-1.0, -1.0, 0.0].into(),
      [ 1.0, -1.0, 0.0].into(),
      [ 0.0,  1.0, 0.0].into()
    );
    let line = geometry::Segment3::noisy (
      [-2.0, 2.0, 0.0].into(),
      [ 2.0, 2.0, 0.0].into());
    let proximity = Proximity::query_triangle_line (&triangle, &line);
    assert_eq!(proximity.nearest_a(), [0.0, 1.0, 0.0].into());
    assert_eq!(proximity.nearest_b(), [0.0, 2.0, 0.0].into());
    assert_eq!(proximity.distance, 1.0);
    assert_eq!(proximity.half_axis, [0.0, 0.5, 0.0].into());
    assert_eq!(proximity.midpoint, [0.0, 1.5, 0.0].into());
    assert_eq!(proximity.normal, -Unit3::axis_y());
    // normal should be invariant under translation
    let triangle = geometry::Triangle3::noisy (
      [-1.0, -1.0, 0.0].into(),
      [ 1.0, -1.0, 0.0].into(),
      [ 0.0,  1.0, 0.0].into());
    let line = geometry::Segment3::noisy (
      [-2.0, 0.0, 0.0].into(),
      [ 2.0, 0.0, 0.0].into());
    let proximity1 = Proximity::query_triangle_line (&triangle, &line);
    let mut rng = XorShiftRng::seed_from_u64 (0);
    let std_normal = rand_distr::StandardNormal;
    // cauchy sample
    let randf = |rng : &mut XorShiftRng| {
      let x : f64 = std_normal.sample (rng);
      let y : f64 = std_normal.sample (rng);
      x / y
    };
    for _ in 0..10000 {
      let translation  = vector3 (randf (&mut rng), randf (&mut rng), randf (&mut rng));
      let mut triangle = triangle;
      let mut line     = line;
      triangle.translate (translation);
      line.translate (translation);
      let proximity = Proximity::query_triangle_line (&triangle, &line);
      approx::assert_relative_eq!(*proximity.normal, *proximity1.normal,
        max_relative = f64::default_max_relative() * 2.0.powi (10))
    }
  }

  #[test]
  #[expect(clippy::unreadable_literal)]
  fn triangle_triangle() {
    use approx::RelativeEq;
    use rand::SeedableRng;
    use rand_distr::Distribution;
    let triangle_a = geometry::Triangle3::noisy (
      [-3.0, 0.0, 0.0].into(),
      [-1.0, 0.0, 0.0].into(),
      [-2.0, 2.0, 0.0].into());
    let triangle_b = geometry::Triangle3::noisy (
      [ 3.0,  0.0, 0.0].into(),
      [ 1.0,  0.0, 0.0].into(),
      [ 2.0,  2.0, 0.0].into());
    let proximity = Proximity::query_triangle_triangle (&triangle_a, &triangle_b);
    assert_eq!(proximity, Proximity {
      distance:  2.0,
      half_axis: [1.0, 0.0, 0.0].into(),
      midpoint:  [0.0, 0.0, 0.0].into(),
      normal:    -Unit3::axis_x()
    });
    let triangle_a = geometry::Triangle3::noisy (
      [-2.0, -2.0, 0.0].into(),
      [ 2.0, -2.0, 0.0].into(),
      [ 0.0,  2.0, 0.0].into());
    let triangle_b = geometry::Triangle3::noisy (
      [-0.5, 0.0, -0.1].into(),
      [ 0.5, 0.0, -0.1].into(),
      [ 0.0, 0.0,  2.0].into());
    let proximity = Proximity::query_triangle_triangle (&triangle_a, &triangle_b);
    assert_eq!(proximity.distance, -0.1);
    approx::assert_relative_eq!(proximity.half_axis, [0.0, 0.0, -0.05].into());
    let triangle_a = geometry::Triangle3::noisy (
      [-2.0, -2.0, 0.0].into(),
      [ 2.0, -2.0, 0.0].into(),
      [ 0.0,  2.0, 0.0].into());
    let triangle_b = geometry::Triangle3::noisy (
      [-1.0, 0.0, -0.1].into(),
      [ 1.0, 0.0, -0.1].into(),
      [ 0.0, 0.0,  2.0].into());
    let proximity = Proximity::query_triangle_triangle (&triangle_a, &triangle_b);
    assert_eq!(proximity.distance, -0.1);
    assert_eq!(proximity.half_axis, [0.0, 0.0, -0.05].into());
    let mut triangle_a = geometry::Triangle3::noisy (
      [-2.0, -2.0, 0.0].into(),
      [ 2.0, -2.0, 0.0].into(),
      [ 0.0,  2.0, 0.0].into());
    let mut triangle_b = geometry::Triangle3::noisy (
      [ 0.0, 0.0, -2.0].into(),
      [-1.0, 0.0,  2.0].into(),
      [ 1.0, 0.0,  2.0].into());
    let proximity = Proximity::query_triangle_triangle (&triangle_a, &triangle_b);
    assert_eq!(proximity, Proximity {
      distance:  -1.3093073414159542,
      half_axis: [0.5714285714285714, -0.2857142857142857, -0.1428571428571429].into(),
      midpoint:  [-0.14285714285714285, 0.28571428571428564, -0.1428571428571429].into(),
      normal:    Unit3::noisy (
        [0.8728715609439696, -0.4364357804719848, -0.2182178902359924].into())
    });
    // make the arrangement asymmetrical so that the result is consistent
    triangle_a.translate (proximity.half_axis * 0.5);
    triangle_b.translate (-proximity.half_axis * 0.5);
    let proximity1 = Proximity::query_triangle_triangle (&triangle_a, &triangle_b);
    approx::assert_relative_eq!(proximity1.distance * 2.0, proximity.distance,
      max_relative = f64::default_max_relative() * 2.0);
    approx::assert_relative_eq!(proximity1.half_axis * 2.0, proximity.half_axis,
      max_relative = f64::default_max_relative() * 16.0);
    approx::assert_relative_eq!(proximity1.midpoint, proximity.midpoint,
      max_relative = f64::default_max_relative() * 8.0);
    approx::assert_relative_eq!(*proximity1.normal, *proximity.normal);
    // normal should be invariant under translation
    let mut rng = XorShiftRng::seed_from_u64 (0);
    let std_normal = rand_distr::StandardNormal;
    // cauchy sample
    let randf = |rng : &mut XorShiftRng| {
      let x : f64 = std_normal.sample (rng);
      let y : f64 = std_normal.sample (rng);
      x / y
    };
    for _i in 0..5000 {
      let translation  = vector3 (randf (&mut rng), randf (&mut rng), randf (&mut rng));
      let mut triangle_a = triangle_a;
      let mut triangle_b = triangle_b;
      triangle_a.translate (translation);
      triangle_b.translate (translation);
      let proximity = Proximity::query_triangle_triangle (&triangle_a, &triangle_b);
      approx::assert_relative_eq!(*proximity.normal, *proximity1.normal,
        max_relative = f64::default_max_relative() * 2.0.powi (10))
    }
  }

  #[test]
  #[expect(clippy::suboptimal_flops)]
  fn distance_query() {
    use rand::{RngExt, SeedableRng};
    use approx::{assert_relative_eq, assert_ulps_eq};

    //
    //  capsule v. capsule
    //

    // intersection: positions of capsules A and B are identical-- half_axis and
    // normal parallel to +X
    let a = object::Static {
      position: component::Position ([0.0, 0.0, 0.0].into()),
      bound:    component::Bound::from (shape::Capsule::noisy (2.0, 3.0)),
      material: component::MATERIAL_STONE,
      collidable: true
    };
    let b = object::Static {
      position: component::Position ([0.0, 0.0, 0.0].into()),
      bound:    component::Bound::from (shape::Capsule::noisy (1.0, 2.0)),
      material: component::MATERIAL_STONE,
      collidable: true
    };
    assert_eq!(
      Proximity::query (&a, &b),
      Proximity {
        distance:  -3.0,
        half_axis: [ 1.5, 0.0, 0.0].into(),
        normal:    Unit3::axis_x(),
        midpoint:  [-0.5, 0.0, 0.0].into()
      }
    );

    // intersection: cylinder shafts of A and B overlap-- half_axis and normal
    // parallel to +X
    let a = object::Static {
      position: component::Position ([0.0, 0.0, 1.0].into()),  .. a
    };
    let b = object::Static {
      position: component::Position ([0.0, 0.0, -1.0].into()), .. b
    };
    assert_eq!(
      Proximity::query (&a, &b),
      Proximity {
        distance:  -3.0,
        half_axis: [ 1.5, 0.0,  0.0].into(),
        normal:    Unit3::axis_x(),
        midpoint:  [-0.5, 0.0, -0.5].into()
      }
    );

    // intersection: cylinder shaft max of A coincides with cylinder shaft min
    // of B
    let a = object::Static {
      position: component::Position ([0.0, 0.0, -3.0].into()), .. a
    };
    let b = object::Static {
      position: component::Position ([0.0, 0.0, 2.0].into()),  .. b
    };
    assert_eq!(
      Proximity::query (&a, &b),
      Proximity {
        distance:  -3.0,
        half_axis: [ 1.5, 0.0, 0.0].into(),
        normal:    Unit3::axis_x(),
        midpoint:  [-0.5, 0.0, 0.0].into()
      }
    );

    // intersection: cylinder shaft min of A coincides with cylinder shaft max of B
    let a = object::Static {
      position: component::Position ([0.0, 0.0, 3.0].into()),  .. a
    };
    let b = object::Static {
      position: component::Position ([0.0, 0.0, -2.0].into()), .. b
    };
    assert_eq!(
      Proximity::query (&a, &b),
      Proximity {
        distance:  -3.0,
        half_axis: [ 1.5, 0.0, 0.0].into(),
        normal:    Unit3::axis_x(),
        midpoint:  [-0.5, 0.0, 0.0].into()
      }
    );

    // intersection: cylinder shaft of A adjacent to cylinder shaft of B
    let a = object::Static {
      position: component::Position ([-1.0, 0.0, 1.0].into()), .. a
    };
    let b = object::Static {
      position: component::Position ([0.0, 0.0, -1.0].into()), .. b
    };
    assert_eq!(
      Proximity::query (&a, &b),
      Proximity {
        distance:  -2.0,
        half_axis: [-1.0, 0.0,  0.0].into(),
        normal:    Unit3::axis_x().invert(),
        midpoint:  [ 0.0, 0.0, -0.5].into()
      }
    );

    // intersection: cylinder shaft of B above cylinder shaft of A
    let a = object::Static {
      position: component::Position ([0.0, 0.0, -3.0].into()), .. a
    };
    let b = object::Static {
      position: component::Position ([0.0, 0.0,  3.0].into()), .. b
    };
    assert_eq!(
      Proximity::query (&a, &b),
      Proximity {
        distance:  -2.0,
        half_axis: [ 0.0, 0.0, -1.0].into(),
        normal:    Unit3::axis_z().invert(),
        midpoint:  [ 0.0, 0.0,  1.0].into()
      }
    );

    // intersection: cylinder shaft of A above cylinder shaft of B
    let a = object::Static {
      position: component::Position ([0.0, 0.0,  3.0].into()), .. a
    };
    let b = object::Static {
      position: component::Position ([0.0, 0.0, -3.0].into()), .. b
    };
    assert_eq!(
      Proximity::query (&a, &b),
      Proximity {
        distance:  -2.0,
        half_axis: [ 0.0, 0.0,  1.0].into(),
        normal:    Unit3::axis_z(),
        midpoint:  [ 0.0, 0.0, -1.0].into()
      }
    );

    // intersection: cylinder shafts of A and B are adjacent
    let a = object::Static {
      position: component::Position ([ 1.0, 0.0, 0.0].into()), .. a
    };
    let b = object::Static {
      position: component::Position ([-1.0, 0.0, 0.0].into()), .. b
    };
    assert_eq!(
      Proximity::query (&a, &b),
      Proximity {
        distance:  -1.0,
        half_axis: [ 0.5, 0.0, 0.0].into(),
        normal:    Unit3::axis_x(),
        midpoint:  [-0.5, 0.0, 0.0].into()
      }
    );

    // disjoint: capsule A above and to the right (+X) of capsule B
    let a = object::Static {
      position: component::Position ([ 2.0, 0.0,  5.0].into()), .. a
    };
    let b = object::Static {
      position: component::Position ([-1.0, 0.0, -3.0].into()), .. b
    };

    let proximity = Proximity::query (&a, &b);
    let distance_manual = {
      let distance  = vector3 (3.0, 0.0, 3.0).magnitude() - 3.0;
      let half_axis : Vector3 <_> =
        0.5 * distance * vector3 (-1.0, 0.0, -1.0).normalized();
      let normal    = Unit3::normalize ([1.0, 0.0, 1.0].into());
      let midpoint  = Point3::from ([-1.0, 0.0, -1.0]) + *normal - half_axis;
      Proximity { distance, half_axis, midpoint, normal }
    };
    assert_eq!(proximity.distance, distance_manual.distance);
    assert_relative_eq!(proximity.half_axis, distance_manual.half_axis);
    assert_relative_eq!(*proximity.normal,   *distance_manual.normal);
    assert_relative_eq!(proximity.midpoint,  distance_manual.midpoint,
      epsilon = 2.0 * f64::EPSILON);

    // disjoint: cylinder shafts of A and B are adjacent
    let a = object::Static {
      position: component::Position ([ 3.0, 0.0, 3.0].into()), .. a
    };
    let b = object::Static {
      position: component::Position ([-1.0, 0.0, 0.0].into()), .. b
    };
    assert_eq!(
      Proximity::query (&a, &b),
      Proximity {
        distance:  1.0,
        half_axis: [-0.5, 0.0, 0.0].into(),
        normal:    Unit3::axis_x(),
        midpoint:  [ 0.5, 0.0, 1.0].into()
      }
    );

    //
    //  cuboid v. cuboid
    //

    // disjoint: nearest points are corners
    let a = object::Static {
      position: component::Position ([-5.0, -5.0, -5.0].into()),
      bound:    component::Bound::from (shape::Cuboid::noisy ([1.0, 2.0, 3.0])),
      .. a
    };
    let b = object::Static {
      position: component::Position ([5.0, 5.0, 5.0].into()),
      bound:    component::Bound::from (shape::Cuboid::noisy ([3.0, 1.0, 2.0])),
      .. b
    };
    let proximity = Proximity::query (&a, &b);
    assert_eq!(proximity.distance, f64::sqrt (110.0));
    assert_eq!(proximity.half_axis, [3.0, 3.5, 2.5].into());
    assert_ulps_eq!(*proximity.normal, -vector3 (3.0, 3.5, 2.5).normalized());
    assert_eq!(proximity.midpoint,  [-1.0, 0.5, 0.5].into());

    // disjoint: overlap on X axis only
    let a = object::Static {
      position: component::Position ([-5.0, 0.0, -5.0].into()), .. a
    };
    let b = object::Static {
      position: component::Position ([5.0, 0.0, 5.0].into()), .. b
    };
    let proximity = Proximity::query (&a, &b);
    assert_eq!(proximity.distance, f64::sqrt (61.0));
    assert_eq!(proximity.half_axis, [3.0, 0.0, 2.5].into());
    assert_ulps_eq!(*proximity.normal, -vector3 (3.0, 0.0, 2.5).normalized());
    assert_eq!(proximity.midpoint,  [-1.0, 0.0, 0.5].into());

    //
    //  capsule v. cuboid
    //

    // disjoint: capsule nearest cuboid max corner
    let a = object::Static {
      position: component::Position ([5.0, 5.0, 3.0].into()),
      bound:    component::Bound::from (shape::Capsule::noisy (SQRT_2, 2.0)),
      .. a
    };
    let b = object::Static {
      position: component::Position ([0.0, 0.0, 0.0].into()),
      bound:    component::Bound::from (shape::Cuboid::noisy ([1.0, 1.0, 1.0])),
      .. b
    };
    let proximity = Proximity::query (&a, &b);
    assert_ulps_eq!(proximity.distance, 3.0 * SQRT_2);
    assert_eq!(proximity.half_axis, [-1.5, -1.5, 0.0].into());
    assert_ulps_eq!(*proximity.normal, vector3 (1.0, 1.0, 0.0).normalized());
    assert_eq!(proximity.midpoint,  [ 2.5, 2.5, 1.0].into());

    // disjoint: capsule nearest cuboid min corner
    let a = object::Static {
      position: component::Position ([-5.0, -5.0, -3.0].into()),
      bound:    component::Bound::from (shape::Capsule::noisy (SQRT_2, 2.0)),
      .. a
    };
    let b = object::Static {
      position: component::Position ([0.0, 0.0, 0.0].into()),
      bound:    component::Bound::from (shape::Cuboid::noisy ([1.0, 1.0, 1.0])),
      .. b
    };
    let proximity = Proximity::query (&a, &b);
    assert_ulps_eq!(proximity.distance, 3.0 * SQRT_2);
    assert_eq!(proximity.half_axis, [ 1.5,  1.5, 0.0].into());
    assert_ulps_eq!(*proximity.normal, vector3 (-1.0, -1.0, 0.0).normalized());
    assert_eq!(proximity.midpoint,  [-2.5,-2.5,-1.0].into());

    // disjoint: capsule nearest cuboid max Y edge
    let a = object::Static {
      position: component::Position ([4.0, 0.5, 4.0].into()),
      bound:    component::Bound::from (shape::Capsule::noisy (1.0, 2.0)),
      .. a
    };
    let b = object::Static {
      position: component::Position ([0.0, 0.0, 0.0].into()),
      bound:    component::Bound::from (shape::Cuboid::noisy ([2.0, 2.0, 2.0])),
      .. b
    };
    assert_eq!(
      Proximity::query (&a, &b),
      Proximity {
        distance:  1.0,
        half_axis: [-0.5, 0.0, 0.0].into(),
        normal:    Unit3::axis_x(),
        midpoint:  [ 2.5, 0.5, 2.0].into()
      }
    );

    // disjoint: capsule nearest cuboid min Y edge
    let a = object::Static {
      position: component::Position ([4.0, 0.5, -4.0].into()), .. a
    };
    let b = object::Static {
      position: component::Position ([0.0, 0.0, 0.0].into()), .. b
    };
    assert_eq!(
      Proximity::query (&a, &b),
      Proximity {
        distance:  1.0,
        half_axis: [-0.5, 0.0,  0.0].into(),
        normal:    Unit3::axis_x(),
        midpoint:  [ 2.5, 0.5, -2.0].into()
      }
    );

    // disjoint: capsule nearest cuboid max X edge
    let a = object::Static {
      position: component::Position ([0.5, 4.0, 4.0].into()), .. a
    };
    let b = object::Static {
      position: component::Position ([0.0, 0.0, 0.0].into()), .. b
    };
    assert_eq!(
      Proximity::query (&a, &b),
      Proximity {
        distance:  1.0,
        half_axis: [ 0.0, -0.5, 0.0].into(),
        normal:    Unit3::axis_y(),
        midpoint:  [ 0.5,  2.5, 2.0].into()
      }
    );

    // disjoint: capsule nearest cuboid min X edge
    let a = object::Static {
      position: component::Position ([0.5, 4.0, -4.0].into()), .. a
    };
    let b = object::Static {
      position: component::Position ([0.0, 0.0, 0.0].into()), .. b
    };
    assert_eq!(
      Proximity::query (&a, &b),
      Proximity {
        distance:  1.0,
        half_axis: [ 0.0, -0.5,  0.0].into(),
        normal:    Unit3::axis_y(),
        midpoint:  [ 0.5,  2.5, -2.0].into()
      }
    );

    // disjoint: capsule A above cuboid B
    let a = object::Static {
      position: component::Position ([0.5, 0.5, 6.0].into()), .. a
    };
    let b = object::Static {
      position: component::Position ([0.0, 0.0, 0.0].into()), .. b
    };
    assert_eq!(
      Proximity::query (&a, &b),
      Proximity {
        distance:  1.0,
        half_axis: [ 0.0, 0.0, -0.5].into(),
        normal:    Unit3::axis_z(),
        midpoint:  [ 0.5, 0.5,  2.5].into()
      }
    );

    // disjoint: capsule A below cuboid B
    let a = object::Static {
      position: component::Position ([0.5, 0.5, -6.0].into()), .. a
    };
    let b = object::Static {
      position: component::Position ([0.0, 0.0, 0.0].into()), .. b
    };
    assert_eq!(
      Proximity::query (&a, &b),
      Proximity {
        distance:  1.0,
        half_axis: [ 0.0, 0.0,  0.5].into(),
        normal:    Unit3::axis_z().invert(),
        midpoint:  [ 0.5, 0.5, -2.5].into()
      }
    );

    // disjoint: capsule A cylinder nearest max Z edge of cuboid B
    let a = object::Static {
      position: component::Position ([3.0, 3.0, 2.0].into()), .. a
    };
    let b = object::Static {
      position: component::Position ([0.0, 0.0, 0.0].into()), .. b
    };
    let proximity = Proximity::query (&a, &b);
    assert_eq!(proximity.distance, SQRT_2 - 1.0);
    assert_eq!(proximity.half_axis, [
      -0.5 * f64::sqrt (0.5 * (SQRT_2 - 1.0).powi (2)),
      -0.5 * f64::sqrt (0.5 * (SQRT_2 - 1.0).powi (2)),
       0.0
    ].into());
    assert_ulps_eq!(*proximity.normal,
      [FRAC_1_SQRT_2, FRAC_1_SQRT_2, 0.0].into());
    assert_eq!(proximity.midpoint, [
      2.0 + 0.5 * f64::sqrt (0.5 * (SQRT_2 - 1.0).powi (2)),
      2.0 + 0.5 * f64::sqrt (0.5 * (SQRT_2 - 1.0).powi (2)),
      1.0
    ].into());

    // disjoint: capsule A cylinder nearest min Z edge of cuboid B
    let a = object::Static {
      position: component::Position ([-3.0, -3.0, -2.0].into()), .. a
    };
    let b = object::Static {
      position: component::Position ([0.0, 0.0, 0.0].into()), .. b
    };
    let proximity = Proximity::query (&a, &b);
    assert_eq!(proximity.distance, SQRT_2 - 1.0);
    assert_eq!(proximity.half_axis, [
       0.5 * f64::sqrt (0.5 * (SQRT_2 - 1.0).powi (2)),
       0.5 * f64::sqrt (0.5 * (SQRT_2 - 1.0).powi (2)),
       0.0
    ].into());
    assert_ulps_eq!(*proximity.normal, [-FRAC_1_SQRT_2, -FRAC_1_SQRT_2, 0.0].into());
    assert_eq!(proximity.midpoint, [
      -2.0 - 0.5 * f64::sqrt (0.5 * (SQRT_2 - 1.0).powi (2)),
      -2.0 - 0.5 * f64::sqrt (0.5 * (SQRT_2 - 1.0).powi (2)),
      -1.0
    ].into());

    // disjoint: capsule A cylinder nearest max X/Z face of cuboid B
    let a = object::Static {
      position: component::Position ([1.0, 4.0, 2.0].into()), .. a
    };
    let b = object::Static {
      position: component::Position ([0.0, 0.0, 0.0].into()), .. b
    };
    assert_eq!(
      Proximity::query (&a, &b),
      Proximity {
        distance:  1.0,
        half_axis: [ 0.0, -0.5,  0.0].into(),
        normal:    Unit3::axis_y(),
        midpoint:  [ 1.0,  2.5,  1.0].into()
      }
    );

    // disjoint: capsule A cylinder nearest min X/Z face of cuboid B
    let a = object::Static {
      position: component::Position ([-1.0, -4.0, -2.0].into()), .. a
    };
    let b = object::Static {
      position: component::Position ([0.0, 0.0, 0.0].into()),    .. b
    };
    assert_eq!(
      Proximity::query (&a, &b),
      Proximity {
        distance:  1.0,
        half_axis: [ 0.0,  0.5,  0.0].into(),
        normal:    Unit3::axis_y().invert(),
        midpoint:  [-1.0, -2.5, -1.0].into()
      }
    );

    // disjoint: capsule A cylinder nearest max Y/Z face of cuboid B
    let a = object::Static {
      position: component::Position ([4.0, 1.0, 2.0].into()), .. a
    };
    let b = object::Static {
      position: component::Position ([0.0, 0.0, 0.0].into()), .. b
    };
    assert_eq!(
      Proximity::query (&a, &b),
      Proximity {
        distance:  1.0,
        half_axis: [-0.5, 0.0, 0.0].into(),
        normal:    Unit3::axis_x(),
        midpoint:  [ 2.5, 1.0, 1.0].into()
      }
    );

    // disjoint: capsule A cylinder nearest min Y/Z face of cuboid B
    let a = object::Static {
      position: component::Position ([-4.0, -1.0, -2.0].into()), .. a
    };
    let b = object::Static {
      position: component::Position ([0.0, 0.0, 0.0].into()),    .. b
    };
    assert_eq!(
      Proximity::query (&a, &b),
      Proximity {
        distance:  1.0,
        half_axis: [ 0.5,  0.0,  0.0].into(),
        normal:    Unit3::axis_x().invert(),
        midpoint:  [-2.5, -1.0, -1.0].into()
      }
    );

    // intersect: capsule cylinder shaft min equal to cuboid max corner
    let a = object::Static {
      position: component::Position ([2.0, 2.0, 4.0].into()),
      bound:    component::Bound::from (shape::Capsule::noisy (1.0, 2.0)),
      .. a
    };
    let b = object::Static {
      position: component::Position ([0.0, 0.0, 0.0].into()),
      bound:    component::Bound::from (shape::Cuboid::noisy ([2.0, 2.0, 2.0])),
      .. b
    };
    assert_eq!(
      Proximity::query (&a, &b),
      Proximity {
        distance:  -1.0,
        half_axis: [ 0.5, 0.0, 0.0].into(),
        normal:    Unit3::axis_x(),
        midpoint:  [ 1.5, 2.0, 2.0].into()
      }
    );

    // intersect: capsule cylinder shaft max equal to cuboid min corner
    let a = object::Static {
      position: component::Position ([-2.0, -2.0, -4.0].into()), .. a
    };
    let b = object::Static {
      position: component::Position ([0.0, 0.0, 0.0].into()),  .. b
    };
    assert_eq!(
      Proximity::query (&a, &b),
      Proximity {
        distance:  -1.0,
        half_axis: [-0.5,  0.0, 0.0].into(),
        normal:    Unit3::axis_x().invert(),
        midpoint:  [-1.5, -2.0,-2.0].into()
      }
    );

    // intersect: capsule cylinder shaft min equal to cuboid max Y edge
    let a = object::Static {
      position: component::Position ([2.0, 1.0, 4.0].into()), .. a
    };
    let b = object::Static {
      position: component::Position ([0.0, 0.0, 0.0].into()), .. b
    };
    assert_eq!(
      Proximity::query (&a, &b),
      Proximity {
        distance:  -1.0,
        half_axis: [ 0.5, 0.0, 0.0].into(),
        normal:    Unit3::axis_x(),
        midpoint:  [ 1.5, 1.0, 2.0].into()
      }
    );

    // intersect: capsule cylinder shaft max equal to cuboid min Y edge
    let a = object::Static {
      position: component::Position ([-2.0, 1.0, -4.0].into()), .. a
    };
    let b = object::Static {
      position: component::Position ([0.0, 0.0, 0.0].into()),  .. b
    };
    assert_eq!(
      Proximity::query (&a, &b),
      Proximity {
        distance:  -1.0,
        half_axis: [-0.5, 0.0,  0.0].into(),
        normal:    Unit3::axis_x().invert(),
        midpoint:  [-1.5, 1.0, -2.0].into()
      }
    );

    // intersect: capsule cylinder shaft min equal to cuboid max X edge
    let a = object::Static {
      position: component::Position ([1.0, 2.0, 4.0].into()), .. a
    };
    let b = object::Static {
      position: component::Position ([0.0, 0.0, 0.0].into()), .. b
    };
    assert_eq!(
      Proximity::query (&a, &b),
      Proximity {
        distance:  -1.0,
        half_axis: [ 0.0, 0.5, 0.0].into(),
        normal:    Unit3::axis_y(),
        midpoint:  [ 1.0, 1.5, 2.0].into()
      }
    );

    // intersect: capsule cylinder shaft max equal to cuboid min X edge
    let a = object::Static {
      position: component::Position ([ 1.0, -2.0, -4.0].into()), .. a
    };
    let b = object::Static {
      position: component::Position ([0.0, 0.0, 0.0].into()),  .. b
    };
    assert_eq!(
      Proximity::query (&a, &b),
      Proximity {
        distance:  -1.0,
        half_axis: [0.0, -0.5,  0.0].into(),
        normal:    Unit3::axis_y().invert(),
        midpoint:  [1.0, -1.5, -2.0].into()
      }
    );

    // intersect: capsule cylinder shaft min equal to cuboid max Z face
    let a = object::Static {
      position: component::Position ([1.0, 1.0, 4.0].into()), .. a
    };
    let b = object::Static {
      position: component::Position ([0.0, 0.0, 0.0].into()), .. b
    };
    assert_eq!(
      Proximity::query (&a, &b),
      Proximity {
        distance:  -1.0,
        half_axis: [ 0.0, 0.0, 0.5].into(),
        normal:    Unit3::axis_z(),
        midpoint:  [ 1.0, 1.0, 1.5].into()
      }
    );

    // intersect: capsule cylinder shaft max equal to cuboid min Z face
    let a = object::Static {
      position: component::Position ([-1.0, -1.0, -4.0].into()), .. a
    };
    let b = object::Static {
      position: component::Position ([0.0, 0.0, 0.0].into()), .. b
    };
    assert_eq!(
      Proximity::query (&a, &b),
      Proximity {
        distance:  -1.0,
        half_axis: [ 0.0,  0.0, -0.5].into(),
        normal:    Unit3::axis_z().invert(),
        midpoint:  [-1.0, -1.0, -1.5].into()
      }
    );

    // intersect: capsule and cuboid position coincide
    let a = object::Static {
      position: component::Position ([0.0, 0.0, 0.0].into()), .. a
    };
    let b = object::Static {
      position: component::Position ([0.0, 0.0, 0.0].into()), .. b
    };
    assert_eq!(
      Proximity::query (&a, &b),
      Proximity {
        distance:  -3.0,
        half_axis: [ 1.5,  0.0,  0.0].into(),
        normal:    Unit3::axis_x(),
        midpoint:  [ 0.5,  0.0,  0.0].into()
      }
    );

    // intersect: capsule nearest positive X face of cuboid
    let a = object::Static {
      position: component::Position ([1.0, 0.0, 0.0].into()), .. a
    };
    let b = object::Static {
      position: component::Position ([0.0, 0.0, 0.0].into()), .. b
    };
    assert_eq!(
      Proximity::query (&a, &b),
      Proximity {
        distance:  -2.0,
        half_axis: [ 1.0,  0.0,  0.0].into(),
        normal:    Unit3::axis_x(),
        midpoint:  [ 1.0,  0.0,  0.0].into()
      }
    );

    // intersect: capsule nearest negative X face of cuboid
    let a = object::Static {
      position: component::Position ([-1.0, 0.0, 0.0].into()), .. a
    };
    let b = object::Static {
      position: component::Position ([0.0, 0.0, 0.0].into()), .. b
    };
    assert_eq!(
      Proximity::query (&a, &b),
      Proximity {
        distance:  -2.0,
        half_axis: [-1.0,  0.0,  0.0].into(),
        normal:    Unit3::axis_x().invert(),
        midpoint:  [-1.0,  0.0,  0.0].into()
      }
    );

    // intersect: capsule nearest positive Y face of cuboid
    let a = object::Static {
      position: component::Position ([0.0, 1.0, 0.0].into()), .. a
    };
    let b = object::Static {
      position: component::Position ([0.0, 0.0, 0.0].into()), .. b
    };
    assert_eq!(
      Proximity::query (&a, &b),
      Proximity {
        distance:  -2.0,
        half_axis: [0.0, 1.0, 0.0].into(),
        normal:    Unit3::axis_y(),
        midpoint:  [0.0, 1.0, 0.0].into()
      }
    );

    // intersect: capsule nearest negative Y face of cuboid
    let a = object::Static {
      position: component::Position ([0.0, -1.0, 0.0].into()), .. a
    };
    let b = object::Static {
      position: component::Position ([0.0, 0.0, 0.0].into()), .. b
    };
    assert_eq!(
      Proximity::query (&a, &b),
      Proximity {
        distance:  -2.0,
        half_axis: [0.0, -1.0, 0.0].into(),
        normal:    Unit3::axis_y().invert(),
        midpoint:  [0.0, -1.0, 0.0].into()
      }
    );

    // intersect: capsule nearest positive Z face of cuboid
    let a = object::Static {
      position: component::Position ([0.0, 0.0, 3.0].into()),
      bound:    component::Bound::from (shape::Capsule::noisy (1.0, 2.0)),
      .. a
    };
    let b = object::Static {
      position: component::Position ([0.0, 0.0, 0.0].into()),
      bound:    component::Bound::from (shape::Cuboid::noisy ([4.0, 4.0, 2.0])),
      .. b
    };
    assert_eq!(
      Proximity::query (&a, &b),
      Proximity {
        distance:  -2.0,
        half_axis: [0.0, 0.0, 1.0].into(),
        normal:    Unit3::axis_z(),
        midpoint:  [0.0, 0.0, 1.0].into()
      }
    );

    // intersect: capsule nearest positive Z face of cuboid
    let a = object::Static {
      position: component::Position ([0.0, 0.0, 1.0].into()),
      bound:    component::Bound::from (shape::Capsule::noisy (1.0, 2.0)),
      .. a
    };
    let b = object::Static {
      position: component::Position ([0.0, 0.0, 0.0].into()),
      bound:    component::Bound::from (shape::Cuboid::noisy ([4.0, 4.0, 2.0])),
      .. b
    };
    assert_eq!(
      Proximity::query (&a, &b),
      Proximity {
        distance:  -4.0,
        half_axis: [0.0, 0.0, 2.0].into(),
        normal:    Unit3::axis_z(),
        midpoint:  [0.0, 0.0, 0.0].into()
      }
    );

    // intersect: capsule nearest negative Z face of cuboid
    let a = object::Static {
      position: component::Position ([0.0, 0.0, -3.0].into()),
      bound:    component::Bound::from (shape::Capsule::noisy (1.0, 2.0)),
      .. a
    };
    let b = object::Static {
      position: component::Position ([0.0, 0.0, 0.0].into()),
      bound:    component::Bound::from (shape::Cuboid::noisy ([4.0, 4.0, 2.0])),
      .. b
    };
    assert_eq!(
      Proximity::query (&a, &b),
      Proximity {
        distance:  -2.0,
        half_axis: [0.0, 0.0, -1.0].into(),
        normal:    Unit3::axis_z().invert(),
        midpoint:  [0.0, 0.0, -1.0].into()
      }
    );

    // intersect: capsule nearest negative Z face of cuboid
    let a = object::Static {
      position: component::Position ([0.0, 0.0, -1.0].into()),
      bound:    component::Bound::from (shape::Capsule::noisy (1.0, 2.0)),
      .. a
    };
    let b = object::Static {
      position: component::Position ([0.0, 0.0, 0.0].into()),
      bound:    component::Bound::from (shape::Cuboid::noisy ([4.0, 4.0, 2.0])),
      .. b
    };
    assert_eq!(
      Proximity::query (&a, &b),
      Proximity {
        distance:  -4.0,
        half_axis: [0.0, 0.0, -2.0].into(),
        normal:    Unit3::axis_z().invert(),
        midpoint:  [0.0, 0.0,  0.0].into()
      }
    );

    // random queries
    let mut rng = XorShiftRng::seed_from_u64 (0);
    let a = object::Static {
      position: component::Position ([0.0, 0.0, 0.0].into()),
      bound:    component::Bound::from (shape::Capsule::noisy (1.5, 1.5)),
      .. a
    };
    for _ in 0..100 {
      let position = component::Position ([
        rng.random_range (-40.0..40.0),
        rng.random_range (-40.0..40.0),
        rng.random_range (-40.0..40.0)
      ].into());
      let bound    = shape::Cuboid::noisy ([
        rng.random_range (0.5..5.0),
        rng.random_range (0.5..5.0),
        rng.random_range (0.5..5.0)
      ]).into();
      let b = object::Static {
        position, bound, material: component::MATERIAL_STONE, collidable: true
      };
      let _ = Proximity::query (&a, &b);
    }

    //
    //  capsule v. orthant
    //
    use SignedAxis3;

    // disjoint: capsule outside of orthant +Z
    let a = object::Static {
      position: component::Position ([0.0, 0.0, 4.0].into()),
      bound:    component::Bound::from (shape::Capsule::noisy (1.0, 2.0)),
      .. a
    };
    let b = object::Static {
      position: component::Position ([0.0, 0.0, 0.0].into()),
      bound:    shape::Orthant { normal_axis: SignedAxis3::PosZ }.into(),
      .. b
    };
    assert_eq!(
      Proximity::query (&a, &b),
      Proximity {
        distance:  1.0,
        half_axis: [0.0, 0.0, -0.5].into(),
        normal:    Unit3::axis_z(),
        midpoint:  [0.0, 0.0,  0.5].into()
      }
    );

    // disjoint: capsule outside of orthant -Z
    let a = object::Static {
      position: component::Position ([0.0, 0.0, -4.0].into()),
      bound:    component::Bound::from (shape::Capsule::noisy (1.0, 2.0)),
      .. a
    };
    let b = object::Static {
      position: component::Position ([0.0, 0.0, 0.0].into()),
      bound:    shape::Orthant { normal_axis: SignedAxis3::NegZ }.into(),
      .. b
    };
    assert_eq!(
      Proximity::query (&a, &b),
      Proximity {
        distance:  1.0,
        half_axis: [0.0, 0.0,  0.5].into(),
        normal:    Unit3::axis_z().invert(),
        midpoint:  [0.0, 0.0, -0.5].into()
      }
    );

    // disjoint: capsule outside orthant +Y
    let a = object::Static {
      position: component::Position ([0.0, 2.0, 0.0].into()),
      bound:    component::Bound::from (shape::Capsule::noisy (1.0, 2.0)),
      .. a
    };
    let b = object::Static {
      position: component::Position ([0.0, 0.0, 0.0].into()),
      bound:    shape::Orthant { normal_axis: SignedAxis3::PosY }.into(),
      .. b
    };
    assert_eq!(
      Proximity::query (&a, &b),
      Proximity {
        distance:  1.0,
        half_axis: [0.0, -0.5, 0.0].into(),
        normal:    Unit3::axis_y(),
        midpoint:  [0.0,  0.5, 0.0].into()
      }
    );

    // disjoint: capsule outside orthant -Y
    let a = object::Static {
      position: component::Position ([0.0, -2.0, 0.0].into()),
      bound:    component::Bound::from (shape::Capsule::noisy (1.0, 2.0)),
      .. a
    };
    let b = object::Static {
      position: component::Position ([0.0, 0.0, 0.0].into()),
      bound:    shape::Orthant { normal_axis: SignedAxis3::NegY }.into(),
      .. b
    };
    assert_eq!(
      Proximity::query (&a, &b),
      Proximity {
        distance:  1.0,
        half_axis: [0.0,  0.5, 0.0].into(),
        normal:    Unit3::axis_y().invert(),
        midpoint:  [0.0, -0.5, 0.0].into()
      }
    );

    // disjoint: capsule outside orthant +X
    let a = object::Static {
      position: component::Position ([2.0, 0.0, 0.0].into()),
      bound:    component::Bound::from (shape::Capsule::noisy (1.0, 2.0)),
      .. a
    };
    let b = object::Static {
      position: component::Position ([0.0, 0.0, 0.0].into()),
      bound:    shape::Orthant { normal_axis: SignedAxis3::PosX }.into(),
      .. b
    };
    assert_eq!(
      Proximity::query (&a, &b),
      Proximity {
        distance:  1.0,
        half_axis: [-0.5, 0.0, 0.0].into(),
        normal:    Unit3::axis_x(),
        midpoint:  [ 0.5, 0.0, 0.0].into()
      }
    );

    // disjoint: outside orthant -X
    let a = object::Static {
      position: component::Position ([-2.0, 0.0, 0.0].into()),
      bound:    component::Bound::from (shape::Capsule::noisy (1.0, 2.0)),
      .. a
    };
    let b = object::Static {
      position: component::Position ([0.0, 0.0, 0.0].into()),
      bound:    shape::Orthant { normal_axis: SignedAxis3::NegX }.into(),
      .. b
    };
    assert_eq!(
      Proximity::query (&a, &b),
      Proximity {
        distance:  1.0,
        half_axis: [ 0.5, 0.0, 0.0].into(),
        normal:    Unit3::axis_x().invert(),
        midpoint:  [-0.5, 0.0, 0.0].into()
      }
    );

    // intersect: capsule and orthant position coincide; normal +Z
    let a = object::Static {
      position: component::Position ([0.0, 0.0, 0.0].into()),
      bound:    component::Bound::from (shape::Capsule::noisy (1.0, 2.0)),
      .. a
    };
    let b = object::Static {
      position: component::Position ([0.0, 0.0, 0.0].into()),
      bound:    shape::Orthant { normal_axis: SignedAxis3::PosZ }.into(),
      .. b
    };
    assert_eq!(
      Proximity::query (&a, &b),
      Proximity {
        distance:  -3.0,
        half_axis: [0.0, 0.0,  1.5].into(),
        normal:    Unit3::axis_z(),
        midpoint:  [0.0, 0.0, -1.5].into()
      }
    );

    // intersect: capsule and orthant position coincide; normal -Z
    let a = object::Static {
      position: component::Position ([0.0, 0.0, 0.0].into()),
      bound:    component::Bound::from (shape::Capsule::noisy (1.0, 2.0)),
      .. a
    };
    let b = object::Static {
      position: component::Position ([0.0, 0.0, 0.0].into()),
      bound:    shape::Orthant { normal_axis: SignedAxis3::NegZ }.into(),
      .. b
    };
    assert_eq!(
      Proximity::query (&a, &b),
      Proximity {
        distance:  -3.0,
        half_axis: [0.0, 0.0, -1.5].into(),
        normal:    Unit3::axis_z().invert(),
        midpoint:  [0.0, 0.0,  1.5].into()
      }
    );

    // intersect: capsule and orthant position coincide; normal +Y
    let a = object::Static {
      position: component::Position ([0.0, 0.0, 0.0].into()),
      bound:    component::Bound::from (shape::Capsule::noisy (1.0, 2.0)),
      .. a
    };
    let b = object::Static {
      position: component::Position ([0.0, 0.0, 0.0].into()),
      bound:    shape::Orthant { normal_axis: SignedAxis3::PosY }.into(),
      .. b
    };
    assert_eq!(
      Proximity::query (&a, &b),
      Proximity {
        distance:  -1.0,
        half_axis: [0.0,  0.5, 0.0].into(),
        normal:    Unit3::axis_y(),
        midpoint:  [0.0, -0.5, 0.0].into()
      }
    );

    // intersect: capsule and orthant position coincide; normal -Y
    let a = object::Static {
      position: component::Position ([0.0, 0.0, 0.0].into()),
      bound:    component::Bound::from (shape::Capsule::noisy (1.0, 2.0)),
      .. a
    };
    let b = object::Static {
      position: component::Position ([0.0, 0.0, 0.0].into()),
      bound:    shape::Orthant { normal_axis: SignedAxis3::NegY }.into(),
      .. b
    };
    assert_eq!(
      Proximity::query (&a, &b),
      Proximity {
        distance:  -1.0,
        half_axis: [0.0, -0.5, 0.0].into(),
        normal:    Unit3::axis_y().invert(),
        midpoint:  [0.0,  0.5, 0.0].into()
      }
    );

    // intersect: capsule and orthant position coincide; normal +X
    let a = object::Static {
      position: component::Position ([0.0, 0.0, 0.0].into()),
      bound:    component::Bound::from (shape::Capsule::noisy (1.0, 2.0)),
      .. a
    };
    let b = object::Static {
      position: component::Position ([0.0, 0.0, 0.0].into()),
      bound:    shape::Orthant { normal_axis: SignedAxis3::PosX }.into(),
      .. b
    };
    assert_eq!(
      Proximity::query (&a, &b),
      Proximity {
        distance:  -1.0,
        half_axis: [ 0.5, 0.0, 0.0].into(),
        normal:    Unit3::axis_x(),
        midpoint:  [-0.5, 0.0, 0.0].into()
      }
    );

    // intersect: capsule and orthant position coincide; normal -X
    let a = object::Static {
      position: component::Position ([0.0, 0.0, 0.0].into()),
      bound:    component::Bound::from (shape::Capsule::noisy (1.0, 2.0)),
      .. a
    };
    let b = object::Static {
      position: component::Position ([0.0, 0.0, 0.0].into()),
      bound:    shape::Orthant { normal_axis: SignedAxis3::NegX }.into(),
      .. b
    };
    assert_eq!(
      Proximity::query (&a, &b),
      Proximity {
        distance:  -1.0,
        half_axis: [-0.5, 0.0, 0.0].into(),
        normal:    Unit3::axis_x().invert(),
        midpoint:  [ 0.5, 0.0, 0.0].into()
      }
    );

    //
    //  sphere v. orthant
    //

    // disjoint: sphere outside of orthant +Z
    let a = object::Static {
      position: component::Position ([0.0, 0.0, 4.0].into()),
      bound:    component::Bound::from (shape::Sphere::noisy (1.0)),
      .. a
    };
    let b = object::Static {
      position: component::Position ([0.0, 0.0, 0.0].into()),
      bound:    shape::Orthant { normal_axis: SignedAxis3::PosZ }.into(),
      .. b
    };
    assert_eq!(
      Proximity::query (&a, &b),
      Proximity {
        distance:  3.0,
        half_axis: [0.0, 0.0, -1.5].into(),
        normal:    Unit3::axis_z(),
        midpoint:  [0.0, 0.0,  1.5].into()
      }
    );

    // disjoint: sphere outside of orthant -Z
    let a = object::Static {
      position: component::Position ([0.0, 0.0, -4.0].into()),
      bound:    component::Bound::from (shape::Sphere::noisy (1.0)),
      .. a
    };
    let b = object::Static {
      position: component::Position ([0.0, 0.0, 0.0].into()),
      bound:    shape::Orthant { normal_axis: SignedAxis3::NegZ }.into(),
      .. b
    };
    assert_eq!(
      Proximity::query (&a, &b),
      Proximity {
        distance:  3.0,
        half_axis: [0.0, 0.0,  1.5].into(),
        normal:    Unit3::axis_z().invert(),
        midpoint:  [0.0, 0.0, -1.5].into()
      }
    );

    // disjoint: sphere outside orthant +Y
    let a = object::Static {
      position: component::Position ([0.0, 2.0, 0.0].into()),
      bound:    component::Bound::from (shape::Sphere::noisy (1.0)),
      .. a
    };
    let b = object::Static {
      position: component::Position ([0.0, 0.0, 0.0].into()),
      bound:    shape::Orthant { normal_axis: SignedAxis3::PosY }.into(),
      .. b
    };
    assert_eq!(
      Proximity::query (&a, &b),
      Proximity {
        distance:  1.0,
        half_axis: [0.0, -0.5, 0.0].into(),
        normal:    Unit3::axis_y(),
        midpoint:  [0.0,  0.5, 0.0].into()
      }
    );

    // disjoint: sphere outside orthant -Y
    let a = object::Static {
      position: component::Position ([0.0, -2.0, 0.0].into()),
      bound:    component::Bound::from (shape::Sphere::noisy (1.0)),
      .. a
    };
    let b = object::Static {
      position: component::Position ([0.0, 0.0, 0.0].into()),
      bound:    shape::Orthant { normal_axis: SignedAxis3::NegY }.into(),
      .. b
    };
    assert_eq!(
      Proximity::query (&a, &b),
      Proximity {
        distance:  1.0,
        half_axis: [0.0,  0.5, 0.0].into(),
        normal:    Unit3::axis_y().invert(),
        midpoint:  [0.0, -0.5, 0.0].into()
      }
    );

    // disjoint: sphere outside orthant +X
    let a = object::Static {
      position: component::Position ([2.0, 0.0, 0.0].into()),
      bound:    component::Bound::from (shape::Sphere::noisy (1.0)),
      .. a
    };
    let b = object::Static {
      position: component::Position ([0.0, 0.0, 0.0].into()),
      bound:    shape::Orthant { normal_axis: SignedAxis3::PosX }.into(),
      .. b
    };
    assert_eq!(
      Proximity::query (&a, &b),
      Proximity {
        distance:  1.0,
        half_axis: [-0.5, 0.0, 0.0].into(),
        normal:    Unit3::axis_x(),
        midpoint:  [ 0.5, 0.0, 0.0].into()
      }
    );

    // disjoint: outside orthant -X
    let a = object::Static {
      position: component::Position ([-2.0, 0.0, 0.0].into()),
      bound:    component::Bound::from (shape::Sphere::noisy (1.0)),
      .. a
    };
    let b = object::Static {
      position: component::Position ([0.0, 0.0, 0.0].into()),
      bound:    shape::Orthant { normal_axis: SignedAxis3::NegX }.into(),
      .. b
    };
    assert_eq!(
      Proximity::query (&a, &b),
      Proximity {
        distance:  1.0,
        half_axis: [ 0.5, 0.0, 0.0].into(),
        normal:    Unit3::axis_x().invert(),
        midpoint:  [-0.5, 0.0, 0.0].into()
      }
    );

    // intersect: sphere and orthant position coincide; normal +Z
    let a = object::Static {
      position: component::Position ([0.0, 0.0, 0.0].into()),
      bound:    component::Bound::from (shape::Sphere::noisy (1.0)),
      .. a
    };
    let b = object::Static {
      position: component::Position ([0.0, 0.0, 0.0].into()),
      bound:    shape::Orthant { normal_axis: SignedAxis3::PosZ }.into(),
      .. b
    };
    assert_eq!(
      Proximity::query (&a, &b),
      Proximity {
        distance:  -1.0,
        half_axis: [0.0, 0.0,  0.5].into(),
        normal:    Unit3::axis_z(),
        midpoint:  [0.0, 0.0, -0.5].into()
      }
    );

    // intersect: sphere and orthant position coincide; normal -Z
    let a = object::Static {
      position: component::Position ([0.0, 0.0, 0.0].into()),
      bound:    component::Bound::from (shape::Sphere::noisy (1.0)),
      .. a
    };
    let b = object::Static {
      position: component::Position ([0.0, 0.0, 0.0].into()),
      bound:    shape::Orthant { normal_axis: SignedAxis3::NegZ }.into(),
      .. b
    };
    assert_eq!(
      Proximity::query (&a, &b),
      Proximity {
        distance:  -1.0,
        half_axis: [0.0, 0.0, -0.5].into(),
        normal:    Unit3::axis_z().invert(),
        midpoint:  [0.0, 0.0,  0.5].into()
      }
    );

    // intersect: sphere and orthant position coincide; normal +Y
    let a = object::Static {
      position: component::Position ([0.0, 0.0, 0.0].into()),
      bound:    component::Bound::from (shape::Sphere::noisy (1.0)),
      .. a
    };
    let b = object::Static {
      position: component::Position ([0.0, 0.0, 0.0].into()),
      bound:    shape::Orthant { normal_axis: SignedAxis3::PosY }.into(),
      .. b
    };
    assert_eq!(
      Proximity::query (&a, &b),
      Proximity {
        distance:  -1.0,
        half_axis: [0.0,  0.5, 0.0].into(),
        normal:    Unit3::axis_y(),
        midpoint:  [0.0, -0.5, 0.0].into()
      }
    );

    // intersect: sphere and orthant position coincide; normal -Y
    let a = object::Static {
      position: component::Position ([0.0, 0.0, 0.0].into()),
      bound:    component::Bound::from (shape::Sphere::noisy (1.0)),
      .. a
    };
    let b = object::Static {
      position: component::Position ([0.0, 0.0, 0.0].into()),
      bound:    shape::Orthant { normal_axis: SignedAxis3::NegY }.into(),
      .. b
    };
    assert_eq!(
      Proximity::query (&a, &b),
      Proximity {
        distance:  -1.0,
        half_axis: [0.0, -0.5, 0.0].into(),
        normal:    Unit3::axis_y().invert(),
        midpoint:  [0.0,  0.5, 0.0].into()
      }
    );

    // intersect: sphere and orthant position coincide; normal +X
    let a = object::Static {
      position: component::Position ([0.0, 0.0, 0.0].into()),
      bound:    component::Bound::from (shape::Sphere::noisy (1.0)),
      .. a
    };
    let b = object::Static {
      position: component::Position ([0.0, 0.0, 0.0].into()),
      bound:    shape::Orthant { normal_axis: SignedAxis3::PosX }.into(),
      .. b
    };
    assert_eq!(
      Proximity::query (&a, &b),
      Proximity {
        distance:  -1.0,
        half_axis: [ 0.5, 0.0, 0.0].into(),
        normal:    Unit3::axis_x(),
        midpoint:  [-0.5, 0.0, 0.0].into()
      }
    );

    // intersect: sphere and orthant position coincide; normal -X
    let a = object::Static {
      position: component::Position ([0.0, 0.0, 0.0].into()),
      bound:    component::Bound::from (shape::Sphere::noisy (1.0)),
      .. a
    };
    let b = object::Static {
      position: component::Position ([0.0, 0.0, 0.0].into()),
      bound:    shape::Orthant { normal_axis: SignedAxis3::NegX }.into(),
      .. b
    };
    assert_eq!(
      Proximity::query (&a, &b),
      Proximity {
        distance:  -1.0,
        half_axis: [-0.5, 0.0, 0.0].into(),
        normal:    Unit3::axis_x().invert(),
        midpoint:  [ 0.5, 0.0, 0.0].into()
      }
    );

    //
    //  cuboid v. orthant
    //

    // disjoint: cuboid outside of orthant +Z
    let a = object::Static {
      position: component::Position ([0.0, 0.0, 4.0].into()),
      bound:    component::Bound::from (shape::Cuboid::noisy ([1.0, 1.0, 1.0])),
      .. a
    };
    let b = object::Static {
      position: component::Position ([0.0, 0.0, 0.0].into()),
      bound:    shape::Orthant { normal_axis: SignedAxis3::PosZ }.into(),
      .. b
    };
    assert_eq!(
      Proximity::query (&a, &b),
      Proximity {
        distance:  3.0,
        half_axis: [0.0, 0.0, -1.5].into(),
        normal:    Unit3::axis_z(),
        midpoint:  [0.0, 0.0,  1.5].into()
      }
    );

    // disjoint: cuboid outside of orthant -Z
    let a = object::Static {
      position: component::Position ([0.0, 0.0, -4.0].into()),
      bound:    component::Bound::from (shape::Cuboid::noisy ([1.0, 1.0, 1.0])),
      .. a
    };
    let b = object::Static {
      position: component::Position ([0.0, 0.0, 0.0].into()),
      bound:    shape::Orthant { normal_axis: SignedAxis3::NegZ }.into(),
      .. b
    };
    assert_eq!(
      Proximity::query (&a, &b),
      Proximity {
        distance:  3.0,
        half_axis: [0.0, 0.0,  1.5].into(),
        normal:    Unit3::axis_z().invert(),
        midpoint:  [0.0, 0.0, -1.5].into()
      }
    );

    // disjoint: cuboid outside orthant +Y
    let a = object::Static {
      position: component::Position ([0.0, 2.0, 0.0].into()),
      bound:    component::Bound::from (shape::Cuboid::noisy ([1.0, 1.0, 1.0])),
      .. a
    };
    let b = object::Static {
      position: component::Position ([0.0, 0.0, 0.0].into()),
      bound:    shape::Orthant { normal_axis: SignedAxis3::PosY }.into(),
      .. b
    };
    assert_eq!(
      Proximity::query (&a, &b),
      Proximity {
        distance:  1.0,
        half_axis: [0.0, -0.5, 0.0].into(),
        normal:    Unit3::axis_y(),
        midpoint:  [0.0,  0.5, 0.0].into()
      }
    );

    // disjoint: cuboid outside orthant -Y
    let a = object::Static {
      position: component::Position ([0.0, -2.0, 0.0].into()),
      bound:    component::Bound::from (shape::Cuboid::noisy ([1.0, 1.0, 1.0])),
      .. a
    };
    let b = object::Static {
      position: component::Position ([0.0, 0.0, 0.0].into()),
      bound:    shape::Orthant { normal_axis: SignedAxis3::NegY }.into(),
      .. b
    };
    assert_eq!(
      Proximity::query (&a, &b),
      Proximity {
        distance:  1.0,
        half_axis: [0.0,  0.5, 0.0].into(),
        normal:    Unit3::axis_y().invert(),
        midpoint:  [0.0, -0.5, 0.0].into()
      }
    );

    // disjoint: cuboid outside orthant +X
    let a = object::Static {
      position: component::Position ([2.0, 0.0, 0.0].into()),
      bound:    component::Bound::from (shape::Cuboid::noisy ([1.0, 1.0, 1.0])),
      .. a
    };
    let b = object::Static {
      position: component::Position ([0.0, 0.0, 0.0].into()),
      bound:    shape::Orthant { normal_axis: SignedAxis3::PosX }.into(),
      .. b
    };
    assert_eq!(
      Proximity::query (&a, &b),
      Proximity {
        distance:  1.0,
        half_axis: [-0.5, 0.0, 0.0].into(),
        normal:    Unit3::axis_x(),
        midpoint:  [ 0.5, 0.0, 0.0].into()
      }
    );

    // disjoint: outside orthant -X
    let a = object::Static {
      position: component::Position ([-2.0, 0.0, 0.0].into()),
      bound:    component::Bound::from (shape::Cuboid::noisy ([1.0, 1.0, 1.0])),
      .. a
    };
    let b = object::Static {
      position: component::Position ([0.0, 0.0, 0.0].into()),
      bound:    shape::Orthant { normal_axis: SignedAxis3::NegX }.into(),
      .. b
    };
    assert_eq!(
      Proximity::query (&a, &b),
      Proximity {
        distance:  1.0,
        half_axis: [ 0.5, 0.0, 0.0].into(),
        normal:    Unit3::axis_x().invert(),
        midpoint:  [-0.5, 0.0, 0.0].into()
      }
    );

    // intersect: cuboid and orthant position coincide; normal +Z
    let a = object::Static {
      position: component::Position ([0.0, 0.0, 0.0].into()),
      bound:    component::Bound::from (shape::Cuboid::noisy ([1.0, 1.0, 1.0])),
      .. a
    };
    let b = object::Static {
      position: component::Position ([0.0, 0.0, 0.0].into()),
      bound:    shape::Orthant { normal_axis: SignedAxis3::PosZ }.into(),
      .. b
    };
    assert_eq!(
      Proximity::query (&a, &b),
      Proximity {
        distance:  -1.0,
        half_axis: [0.0, 0.0,  0.5].into(),
        normal:    Unit3::axis_z(),
        midpoint:  [0.0, 0.0, -0.5].into()
      }
    );

    // intersect: cuboid and orthant position coincide; normal -Z
    let a = object::Static {
      position: component::Position ([0.0, 0.0, 0.0].into()),
      bound:    component::Bound::from (shape::Cuboid::noisy ([1.0, 1.0, 1.0])),
      .. a
    };
    let b = object::Static {
      position: component::Position ([0.0, 0.0, 0.0].into()),
      bound:    shape::Orthant { normal_axis: SignedAxis3::NegZ }.into(),
      .. b
    };
    assert_eq!(
      Proximity::query (&a, &b),
      Proximity {
        distance:  -1.0,
        half_axis: [0.0, 0.0, -0.5].into(),
        normal:    Unit3::axis_z().invert(),
        midpoint:  [0.0, 0.0,  0.5].into()
      }
    );

    // intersect: cuboid and orthant position coincide; normal +Y
    let a = object::Static {
      position: component::Position ([0.0, 0.0, 0.0].into()),
      bound:    component::Bound::from (shape::Cuboid::noisy ([1.0, 1.0, 1.0])),
      .. a
    };
    let b = object::Static {
      position: component::Position ([0.0, 0.0, 0.0].into()),
      bound:    shape::Orthant { normal_axis: SignedAxis3::PosY }.into(),
      .. b
    };
    assert_eq!(
      Proximity::query (&a, &b),
      Proximity {
        distance:  -1.0,
        half_axis: [0.0,  0.5, 0.0].into(),
        normal:    Unit3::axis_y(),
        midpoint:  [0.0, -0.5, 0.0].into()
      }
    );

    // intersect: cuboid and orthant position coincide; normal -Y
    let a = object::Static {
      position: component::Position ([0.0, 0.0, 0.0].into()),
      bound:    component::Bound::from (shape::Cuboid::noisy ([1.0, 1.0, 1.0])),
      .. a
    };
    let b = object::Static {
      position: component::Position ([0.0, 0.0, 0.0].into()),
      bound:    shape::Orthant { normal_axis: SignedAxis3::NegY }.into(),
      .. b
    };
    assert_eq!(
      Proximity::query (&a, &b),
      Proximity {
        distance:  -1.0,
        half_axis: [0.0, -0.5, 0.0].into(),
        normal:    Unit3::axis_y().invert(),
        midpoint:  [0.0,  0.5, 0.0].into()
      }
    );

    // intersect: cuboid and orthant position coincide; normal +X
    let a = object::Static {
      position: component::Position ([0.0, 0.0, 0.0].into()),
      bound:    component::Bound::from (shape::Cuboid::noisy ([1.0, 1.0, 1.0])),
      .. a
    };
    let b = object::Static {
      position: component::Position ([0.0, 0.0, 0.0].into()),
      bound:    shape::Orthant { normal_axis: SignedAxis3::PosX }.into(),
      .. b
    };
    assert_eq!(
      Proximity::query (&a, &b),
      Proximity {
        distance:  -1.0,
        half_axis: [ 0.5, 0.0, 0.0].into(),
        normal:    Unit3::axis_x(),
        midpoint:  [-0.5, 0.0, 0.0].into()
      }
    );

    // intersect: cuboid and orthant position coincide; normal -X
    let a = object::Static {
      position: component::Position ([0.0, 0.0, 0.0].into()),
      bound:    component::Bound::from (shape::Cuboid::noisy ([1.0, 1.0, 1.0])),
      .. a
    };
    let b = object::Static {
      position: component::Position ([0.0, 0.0, 0.0].into()),
      bound:    shape::Orthant { normal_axis: SignedAxis3::NegX }.into(),
      .. b
    };
    assert_eq!(
      Proximity::query (&a, &b),
      Proximity {
        distance:  -1.0,
        half_axis: [-0.5, 0.0, 0.0].into(),
        normal:    Unit3::axis_x().invert(),
        midpoint:  [ 0.5, 0.0, 0.0].into()
      }
    );

  } // end test_distance_query

  #[test]
  #[expect(clippy::approx_constant)]
  fn query_hull() {
    // cuboid/hull
    let object_a = (
      component::Position ([0.0, 1.0, 1.6].into()),
      component::Bound::from (shape::Cuboid::noisy ([2.2, 2.2, 0.2]))
    );
    let object_b = {
      let hull = {
        let points = [
          [0.7333798838354574, 0.07387687099658402, -0.4254516842354885],
          [0.28934048979234384, -0.41943067939022566, 0.20302532977593885],
          [0.006562630606422983, 0.7345982182090985, 0.6784706255126608],
          [-0.16787537854609894, -0.7656664645854194, 0.004114727975213302],
          [-0.6400977561459589, -0.062079241531162725, -0.33763483879239625],
          [-0.3844054415702367, 0.6216012595122049, -0.04739397447380115]
        ].map (Point3::from);
        geometry::Hull3::from_points (&points).unwrap()
      };
      ( component::Position (
          [-0.021233076872148704, 1.4862487906056179, 2.2279516842354887].into()),
        component::Bound::from (hull)
      )
    };
    let proximity = Proximity::query (&object_a, &object_b);
    assert_eq!(proximity.distance, 0.0025000000000001688);

    let object_a = {
      let hull = {
        let points = [
          [-1.6970562748477143, -2.2,  1.4142135623730951],
          [-1.4142135623730954, -2.2,  1.697056274847714],
          [-1.6970562748477143,  2.2,  1.4142135623730951],
          [-1.4142135623730954,  2.2,  1.697056274847714],
          [ 1.4142135623730954, -2.2, -1.697056274847714],
          [ 1.6970562748477143, -2.2, -1.4142135623730951],
          [ 1.4142135623730954,  2.2, -1.697056274847714],
          [ 1.6970562748477143,  2.2, -1.4142135623730951]
        ].map (Point3::from);
        geometry::Hull3::from_points (&points).unwrap()
      };
      ( component::Position ([-3.9, 1.0, 11.2125].into()),
        component::Bound::from (hull)
      )
    };
    let object_b = (
      component::Position (
        [-4.004858348343004, -1.2161640324940572, 11.869089780656527].into()),
      component::Bound::from (shape::Sphere::noisy (0.5))
    );
    let proximity = Proximity::query (&object_a, &object_b);
    assert_eq!(proximity.distance, -0.3091811079468084);

    let object_a = {
      let hull = {
        let points = [
          [-2.2, -1.4142135623730954, -1.697056274847714],
          [-2.2, -1.6970562748477143, -1.4142135623730951],
          [-2.2,  1.6970562748477143,  1.4142135623730951],
          [-2.2,  1.4142135623730954,  1.697056274847714],
          [ 2.2, -1.4142135623730954, -1.697056274847714],
          [ 2.2, -1.6970562748477143, -1.4142135623730951],
          [ 2.2,  1.6970562748477143,  1.4142135623730951],
          [ 2.2,  1.4142135623730954,  1.69705627484771]
        ].map (Point3::from);
        geometry::Hull3::from_points (&points).unwrap()
      };
      ( component::Position ([0.0, 4.9, 11.2125].into()),
        component::Bound::from (hull)
      )
    };
    let object_b = (
      component::Position (
        [1.3113904930865399, 3.1609676773871804, 9.819471286970987].into()),
      component::Bound::from (shape::Sphere::noisy (0.5))
    );
    let proximity = Proximity::query (&object_a, &object_b);
    assert_eq!(proximity.distance, -0.4529809974587258);

    let object_a = {
      let hull = {
        let points = [
          [-2.2, -1.6970562748477138,   1.4142135623730956],
          [-2.2, -1.414213562373095,    1.6970562748477145],
          [-2.2,  1.414213562373095,   -1.6970562748477145],
          [-2.2,  1.6970562748477138,  -1.4142135623730956],
          [ 2.2, -1.6970562748477138,   1.4142135623730956],
          [ 2.2, -1.414213562373095,    1.6970562748477145],
          [ 2.2,  1.414213562373095,   -1.6970562748477145],
          [ 2.2,  1.6970562748477138,  -1.4142135623730956]
        ].map (Point3::from);
        geometry::Hull3::from_points (&points).unwrap()
      };
      ( component::Position ([0.0, -2.9, 11.2125].into()),
        component::Bound::from (hull)
      )
    };
    let object_b = (
      component::Position (
        [2.256458188002921, -0.9674933916701636, 10.065811222437437].into()),
      component::Bound::from (shape::Sphere::noisy (0.5))
    );
    let proximity = Proximity::query (&object_a, &object_b);
    assert_eq!(proximity.distance, -0.13988958350691275);

    let object_a = {
      let hull = {
        let points = [
          [-1.8000000000000003, -1.4142135623730954, -1.697056274847714 ],
          [-1.8000000000000003, -1.6970562748477143, -1.4142135623730951],
          [-1.8000000000000003,  1.6970562748477143,  1.4142135623730951],
          [-1.8000000000000003,  1.4142135623730954,  1.697056274847714 ],
          [ 1.8000000000000003, -1.4142135623730954, -1.697056274847714 ],
          [ 1.8000000000000003, -1.6970562748477143, -1.4142135623730951],
          [ 1.8000000000000003,  1.6970562748477143,  1.4142135623730951],
          [ 1.8000000000000003,  1.4142135623730954,  1.697056274847714 ]
        ].map (Point3::from);
        geometry::Hull3::from_points (&points).unwrap()
      };
      ( component::Position ([0.0, 4.5, 11.2125].into()),
        component::Bound::from (hull)
      )
    };
    let object_b = (
      component::Position (
        [-1.8483359321030037, 3.6499052772728637, 11.360649949330503].into()),
      component::Bound::from (shape::Sphere::noisy (0.5))
    );
    let proximity = Proximity::query (&object_a, &object_b);
    assert_eq!(proximity.distance, 0.008169602403322415);
  }

  #[test]
  fn distance_to_triangle_edges() {
    let triangle = geometry::Triangle3::noisy (
      [-2.0, -1.0, 0.0].into(),
      [ 1.0, -1.0, 0.0].into(),
      [ 1.0,  5.0, 0.0].into()
    );
    // origin
    let s = 0.5;
    let t = 1.0 / 6.0;
    let (d0, d1, d2) = triangle_distance_to_edges (triangle, s, t);
    approx::assert_relative_eq!(d0, 1.0);
    assert_eq!(d1, 3.0 / 5.0.sqrt());
    assert_eq!(d2, 1.0);
  }

} // end tests