math_utils/geometry/primitive/
mod.rs

1//! Affine, convex, and combinatorial spaces
2
3use approx;
4use num_traits as num;
5use rand;
6
7use crate::*;
8use super::{intersect, shape};
9
10#[cfg(feature = "derive_serdes")]
11use serde::{Deserialize, Serialize};
12
13/// Primitives over signed integers
14pub mod integer;
15
16mod simplex;
17pub use self::simplex::{simplex2, simplex3, Simplex2, Simplex3};
18pub use simplex2::Segment     as Segment2;
19pub use simplex3::Segment     as Segment3;
20pub use simplex2::Triangle    as Triangle2;
21pub use simplex3::Triangle    as Triangle3;
22pub use simplex3::Tetrahedron as Tetrahedron3;
23
24pub trait Primitive <S, V> where
25  S : Field,
26  V : VectorSpace <S>
27{
28  fn translate (self, displacement : V) -> Self;
29  fn scale     (self, scale : NonZero <S>) -> Self;
30}
31
32////////////////////////////////////////////////////////////////////////////////
33//  structs                                                                   //
34////////////////////////////////////////////////////////////////////////////////
35
36/// 1D interval
37#[cfg_attr(feature = "derive_serdes", derive(Serialize, Deserialize))]
38#[derive(Clone, Copy, Debug, PartialEq)]
39pub struct Interval <S> {
40  min : S,
41  max : S
42}
43
44/// 2D axis-aligned bounding box.
45///
46/// Min/max can be co-linear:
47/// ```
48/// # use math_utils::geometry::Aabb2;
49/// let x = Aabb2::with_minmax (
50///   [0.0, 0.0].into(),
51///   [0.0, 1.0].into()
52/// );
53/// ```
54/// But not identical:
55/// ```should_panic
56/// # use math_utils::geometry::Aabb2;
57/// let x = Aabb2::with_minmax (
58///   [0.0, 0.0].into(),
59///   [0.0, 0.0].into()
60/// );
61/// ```
62#[cfg_attr(feature = "derive_serdes", derive(Serialize, Deserialize))]
63#[derive(Clone, Copy, Debug, PartialEq)]
64pub struct Aabb2 <S> {
65  min : Point2 <S>,
66  max : Point2 <S>
67}
68
69/// 3D axis-aligned bounding box.
70///
71/// See also `shape::Aabb` trait for primitive and shape types that can compute
72/// a 3D AABB.
73///
74/// Min/max can be co-planar or co-linear:
75/// ```
76/// # use math_utils::geometry::Aabb3;
77/// let x = Aabb3::with_minmax (
78///   [0.0, 0.0, 0.0].into(),
79///   [0.0, 1.0, 1.0].into()
80/// );
81/// let y = Aabb3::with_minmax (
82///   [0.0, 0.0, 0.0].into(),
83///   [0.0, 1.0, 0.0].into()
84/// );
85/// ```
86/// But not identical:
87/// ```should_panic
88/// # use math_utils::geometry::Aabb3;
89/// let x = Aabb3::with_minmax (
90///   [0.0, 0.0, 0.0].into(),
91///   [0.0, 0.0, 0.0].into()
92/// );
93/// ```
94#[cfg_attr(feature = "derive_serdes", derive(Serialize, Deserialize))]
95#[derive(Clone, Copy, Debug, PartialEq)]
96pub struct Aabb3 <S> {
97  min : Point3 <S>,
98  max : Point3 <S>
99}
100
101/// 3D Z-axis-aligned cylinder
102#[cfg_attr(feature = "derive_serdes", derive(Serialize, Deserialize))]
103#[derive(Clone, Copy, Debug, PartialEq)]
104pub struct Cylinder3 <S> {
105  pub center      : Point3   <S>,
106  pub half_height : Positive <S>,
107  pub radius      : Positive <S>
108}
109
110/// 3D Z-axis-aligned capsule
111#[cfg_attr(feature = "derive_serdes", derive(Serialize, Deserialize))]
112#[derive(Clone, Copy, Debug, PartialEq)]
113pub struct Capsule3 <S> {
114  pub center      : Point3   <S>,
115  pub half_height : Positive <S>,
116  pub radius      : Positive <S>
117}
118
119/// An infinitely extended line in 2D space defined by a base point and
120/// normalized direction
121#[cfg_attr(feature = "derive_serdes", derive(Serialize, Deserialize))]
122#[derive(Clone, Copy, Debug, PartialEq)]
123pub struct Line2 <S> {
124  pub base      : Point2 <S>,
125  pub direction : Unit2  <S>
126}
127
128/// An infinitely extended line in 3D space defined by a base point and
129/// normalized direction
130#[cfg_attr(feature = "derive_serdes", derive(Serialize, Deserialize))]
131#[derive(Clone, Copy, Debug, PartialEq)]
132pub struct Line3 <S> {
133  pub base      : Point3 <S>,
134  pub direction : Unit3  <S>
135}
136
137/// A plane in 3D space defined by a base point and (unit) normal vector
138#[cfg_attr(feature = "derive_serdes", derive(Serialize, Deserialize))]
139#[derive(Clone, Copy, Debug, PartialEq)]
140pub struct Plane3 <S> {
141  pub base   : Point3 <S>,
142  pub normal : Unit3  <S>
143}
144
145/// Sphere in 2D space (a circle)
146#[cfg_attr(feature = "derive_serdes", derive(Serialize, Deserialize))]
147#[derive(Clone, Copy, Debug, PartialEq)]
148pub struct Sphere2 <S> {
149  pub center : Point2   <S>,
150  pub radius : Positive <S>
151}
152
153/// Sphere in 3D space
154#[cfg_attr(feature = "derive_serdes", derive(Serialize, Deserialize))]
155#[derive(Clone, Copy, Debug, PartialEq)]
156pub struct Sphere3 <S> {
157  pub center : Point3   <S>,
158  pub radius : Positive <S>
159}
160
161////////////////////////////////////////////////////////////////////////////////
162//  functions                                                                 //
163////////////////////////////////////////////////////////////////////////////////
164
165/// Returns true when three points lie on the same line in 2D space.
166///
167/// ```
168/// # use math_utils::geometry::colinear_2d;
169/// assert!(colinear_2d (
170///   &[-1.0, -1.0].into(),
171///   &[ 0.0,  0.0].into(),
172///   &[ 1.0,  1.0].into())
173/// );
174/// assert!(!colinear_2d (
175///   &[-1.0, -1.0].into(),
176///   &[ 0.0,  1.0].into(),
177///   &[ 1.0, -1.0].into())
178/// );
179/// ```
180pub fn colinear_2d <S> (
181  a : &Point2 <S>,
182  b : &Point2 <S>,
183  c : &Point2 <S>
184) -> bool where
185  S : Field + approx::RelativeEq
186{
187  // early exit: a pair of points are equal
188  if a == b || a == c || b == c {
189    true
190  } else {
191    let a = Vector3::from_point_2d (a.0).into();
192    let b = Vector3::from_point_2d (b.0).into();
193    let c = Vector3::from_point_2d (c.0).into();
194    let determinant = determinant_3d (&a, &b, &c);
195    // a zero determinant indicates that the points are colinear
196    approx::relative_eq!(S::zero(), determinant)
197  }
198}
199
200/// Returns true when three points lie on the same line in 3D space.
201///
202/// ```
203/// # use math_utils::geometry::colinear_3d;
204/// assert!(colinear_3d (
205///   &[-1.0, -1.0, -1.0].into(),
206///   &[ 0.0,  0.0,  0.0].into(),
207///   &[ 1.0,  1.0,  1.0].into())
208/// );
209/// ```
210pub fn colinear_3d <S> (a : &Point3 <S>, b : &Point3 <S>, c : &Point3 <S>)
211  -> bool
212where
213  S : Field + Sqrt + approx::RelativeEq
214{
215  // early exit: a pair of points are equal
216  if a == b || a == c || b == c {
217    true
218  } else {
219    let determinant = determinant_3d (a, b, c);
220    // a zero determinant indicates that the points are coplanar
221    if approx::relative_eq!(S::zero(), determinant) {
222      approx::relative_eq!(S::zero(), triangle_3d_area2 (a, b, c))
223    } else {
224      false
225    }
226  }
227}
228
229/// Returns true when four points lie on the same plane in 3D space.
230///
231/// ```
232/// # use math_utils::geometry::coplanar_3d;
233/// assert!(coplanar_3d (
234///   &[-1.0, -1.0, -1.0].into(),
235///   &[ 1.0,  1.0,  1.0].into(),
236///   &[-1.0,  1.0,  0.0].into(),
237///   &[ 1.0, -1.0,  0.0].into()
238/// ));
239/// ```
240pub fn coplanar_3d <S> (
241  a : &Point3 <S>,
242  b : &Point3 <S>,
243  c : &Point3 <S>,
244  d : &Point3 <S>
245) -> bool where
246  S : Ring + approx::RelativeEq
247{
248  // early exit: a pair of points are equal
249  if a == b || a == c || a == d || b == c || b == d || c == d {
250    true
251  } else {
252    let ab_cross_ac_dot_ad = (*b-a).cross (*c-a).dot (*d-a);
253    approx::relative_eq!(S::zero(), ab_cross_ac_dot_ad)
254  }
255}
256
257/// Square area of three points in 3D space.
258///
259/// Uses a numerically stable Heron's formula:
260/// <https://en.wikipedia.org/wiki/Heron%27s_formula#Numerical_stability>
261///
262/// ```
263/// # use math_utils::approx::assert_relative_eq;
264/// # use math_utils::geometry::triangle_3d_area2;
265/// assert_relative_eq!(
266///   3.0/4.0,
267///   triangle_3d_area2 (
268///     &[-1.0,  0.0,  0.0].into(),
269///     &[ 0.0,  0.0,  1.0].into(),
270///     &[ 0.0,  1.0,  0.0].into())
271/// );
272/// ```
273///
274/// If the area squared is zero then the points are colinear:
275///
276/// ```
277/// # use math_utils::geometry::triangle_3d_area2;
278/// assert_eq!(
279///   0.0,
280///   triangle_3d_area2 (
281///     &[-1.0, -1.0, -1.0].into(),
282///     &[ 0.0,  0.0,  0.0].into(),
283///     &[ 1.0,  1.0,  1.0].into())
284/// );
285/// ```
286pub fn triangle_3d_area2 <S> (a : &Point3 <S>, b : &Point3 <S>, c : &Point3 <S>)
287  -> S
288where
289  S : Field + Sqrt
290{
291  if a == b || a == c || b == c {
292    return S::zero()
293  }
294  // compute the length of each side
295  let ab_mag = (*b-a).norm();
296  let ac_mag = (*c-a).norm();
297  let bc_mag = (*c-b).norm();
298  // order as max >= mid >= min
299  let max = S::max (S::max (ab_mag, ac_mag), bc_mag);
300  let min = S::min (S::min (ab_mag, ac_mag), bc_mag);
301  let mid = if min <= ab_mag && ab_mag <= max {
302    ab_mag
303  } else if min <= ac_mag && ac_mag <= max {
304    ac_mag
305  } else {
306    bc_mag
307  };
308  let two = S::one() + S::one();
309  let sixteen = two * two * two * two;
310  // 1.0/16.0
311  let frac = S::one() / sixteen;
312  frac
313    * (max + (mid + min))
314    * (min - (max - mid))
315    * (min + (max - mid))
316    * (max + (mid - min))
317}
318
319/// Computes the determinant of a matrix formed by the three points as columns
320#[inline]
321pub fn determinant_3d <S : Ring> (
322  a : &Point3 <S>,
323  b : &Point3 <S>,
324  c : &Point3 <S>
325) -> S {
326  Matrix3::from_col_arrays ([
327    a.0.into_array(), b.0.into_array(), c.0.into_array()
328  ]).determinant()
329}
330
331/// Coordinate-wise min
332pub fn point2_min <S : Ring> (a : &Point2 <S>, b : &Point2 <S>) -> Point2 <S> {
333  Vector2::partial_min (a.0, b.0).into()
334}
335
336/// Coordinate-wise max
337pub fn point2_max <S : Ring> (a : &Point2 <S>, b : &Point2 <S>) -> Point2 <S> {
338  Vector2::partial_max (a.0, b.0).into()
339}
340
341/// Coordinate-wise min
342pub fn point3_min <S : Ring> (a : &Point3 <S>, b : &Point3 <S>) -> Point3 <S> {
343  Vector3::partial_min (a.0, b.0).into()
344}
345
346/// Coordinate-wise max
347pub fn point3_max <S : Ring> (a : &Point3 <S>, b : &Point3 <S>) -> Point3 <S> {
348  Vector3::partial_max (a.0, b.0).into()
349}
350
351/// Given a 2D point and a 2D line, returns the nearest point on the line to
352/// the given point
353pub fn project_2d_point_on_line <S : Real>
354  (point : &Point2 <S>, line : &Line2 <S>) -> Point2 <S>
355{
356  let dot_dir = line.direction.dot (*point - line.base);
357  line.point (dot_dir)
358}
359
360/// Given a 3D point and a 3D line, returns the nearest point on the line to
361/// the given point
362pub fn project_3d_point_on_line <S : Real>
363  (point : &Point3 <S>, line : &Line3 <S>) -> Point3 <S>
364{
365  let dot_dir = line.direction.dot (*point - line.base);
366  line.point (dot_dir)
367}
368
369////////////////////////////////////////////////////////////////////////////////
370//  impls                                                                     //
371////////////////////////////////////////////////////////////////////////////////
372
373impl <S : Ring> Interval <S> {
374  /// Construct a new AABB with given min and max points.
375  ///
376  /// Debug panic if points are not min/max:
377  ///
378  /// ```should_panic
379  /// # use math_utils::geometry::*;
380  /// let b = Interval::with_minmax (1.0, 0.0);  // panic!
381  /// ```
382  ///
383  /// Debug panic if points are identical:
384  ///
385  /// ```should_panic
386  /// # use math_utils::geometry::*;
387  /// let b = Interval::with_minmax (0.0, 0.0);  // panic!
388  /// ```
389  #[inline]
390  pub fn with_minmax (min : S, max : S) -> Self where S : std::fmt::Debug {
391    debug_assert_ne!(min, max);
392    debug_assert_eq!(min, S::min (min, max));
393    debug_assert_eq!(max, S::max (min, max));
394    Interval { min, max }
395  }
396  /// Construct a new AABB using the two given points to determine min/max.
397  ///
398  /// Panic if points are identical:
399  ///
400  /// ```should_panic
401  /// # use math_utils::geometry::*;
402  /// let b = Interval::from_points (0.0, 0.0);  // panic!
403  /// ```
404  #[inline]
405  pub fn from_points (a : S, b : S) -> Self where S : std::fmt::Debug {
406    debug_assert_ne!(a, b);
407    let min = S::min (a, b);
408    let max = S::max (a, b);
409    Interval { min, max }
410  }
411  /// Construct the minimum AABB containing the given set of points.
412  ///
413  /// Debug panic if fewer than 2 points are given:
414  ///
415  /// ```should_panic
416  /// # use math_utils::geometry::*;
417  /// let b = Interval::<f32>::containing (&[]);  // panic!
418  /// ```
419  ///
420  /// ```should_panic
421  /// # use math_utils::geometry::*;
422  /// let b = Interval::containing (&[0.0]);  // panic!
423  /// ```
424  #[inline]
425  pub fn containing (points : &[S]) -> Self where S : std::fmt::Debug {
426    debug_assert!(points.len() >= 2);
427    let mut min = S::min (points[0], points[1]);
428    let mut max = S::max (points[0], points[1]);
429    for point in points.iter().skip (2) {
430      min = S::min (min, *point);
431      max = S::max (max, *point);
432    }
433    Interval::with_minmax (min, max)
434  }
435  #[inline]
436  pub fn numcast <T> (self) -> Option <Interval <T>> where
437    S : num::ToPrimitive,
438    T : num::NumCast
439  {
440    Some (Interval {
441      min: T::from (self.min)?,
442      max: T::from (self.max)?
443    })
444  }
445  /// Create a new AABB that is the union of the two input AABBs
446  #[inline]
447  pub fn union (a : &Interval <S>, b : &Interval <S>) -> Self where
448    S : std::fmt::Debug
449  {
450    Interval::with_minmax (
451      S::min (a.min(), b.min()),
452      S::max (a.max(), b.max()))
453  }
454  #[inline]
455  pub fn min (&self) -> S {
456    self.min
457  }
458  #[inline]
459  pub fn max (&self) -> S {
460    self.max
461  }
462  #[inline]
463  pub fn width (&self) -> NonNegative <S> {
464    NonNegative::unchecked (self.max - self.min)
465  }
466  #[inline]
467  pub fn contains (&self, point : &S) -> bool {
468    self.min < *point && *point < self.max
469  }
470  /// Clamp a given point to the AABB interval.
471  ///
472  /// ```
473  /// # use math_utils::geometry::*;
474  /// let b = Interval::from_points (-1.0, 1.0);
475  /// assert_eq!(b.clamp (&-2.0), -1.0);
476  /// assert_eq!(b.clamp ( &2.0),  1.0);
477  /// assert_eq!(b.clamp ( &0.0),  0.0);
478  /// ```
479  pub fn clamp (&self, point : &S) -> S {
480    S::max (S::min (self.max, *point), self.min)
481  }
482  /// Generate a random point contained in the AABB
483  ///
484  /// ```
485  /// # use rand;
486  /// # use rand_xorshift;
487  /// # use math_utils::geometry::*;
488  /// # fn main () {
489  /// use rand_xorshift;
490  /// use rand::SeedableRng;
491  /// // random sequence will be the same each time this is run
492  /// let mut rng = rand_xorshift::XorShiftRng::seed_from_u64 (0);
493  /// let aabb = Interval::<f32>::with_minmax (-10.0, 10.0);
494  /// let point = aabb.rand_point (&mut rng);
495  /// assert!(aabb.contains (&point));
496  /// let point = aabb.rand_point (&mut rng);
497  /// assert!(aabb.contains (&point));
498  /// let point = aabb.rand_point (&mut rng);
499  /// assert!(aabb.contains (&point));
500  /// let point = aabb.rand_point (&mut rng);
501  /// assert!(aabb.contains (&point));
502  /// let point = aabb.rand_point (&mut rng);
503  /// assert!(aabb.contains (&point));
504  /// # }
505  /// ```
506  #[inline]
507  pub fn rand_point <R> (&self, rng : &mut R) -> S where
508    S : rand::distr::uniform::SampleUniform,
509    R : rand::Rng
510  {
511    rng.random_range (self.min..self.max)
512  }
513  #[inline]
514  pub fn intersects (&self, other : &Interval <S>) -> bool {
515    intersect::discrete_interval (self, other)
516  }
517  #[inline]
518  pub fn intersection (self, other : Interval <S>) -> Option <Interval <S>>
519    where S : std::fmt::Debug
520  {
521    intersect::continuous_interval (&self, &other)
522  }
523
524  /// Increase or decrease each endpoint by the given amount.
525  ///
526  /// Debug panic if the interval width is less than or equal to twice the given
527  /// amount:
528  /// ```should_panic
529  /// # use math_utils::geometry::*;
530  /// let x = Interval::<f32>::with_minmax (-1.0, 1.0);
531  /// let y = x.dilate (-1.0); // panic!
532  /// ```
533  pub fn dilate (mut self, amount : S) -> Self {
534    self.min -= amount;
535    self.max += amount;
536    debug_assert!(self.min < self.max);
537    self
538  }
539}
540
541impl <S> Primitive <S, S> for Interval <S> where
542  S : Field + VectorSpace <S>
543{
544  fn translate (mut self, displacement : S) -> Self {
545    self.min += displacement;
546    self.max += displacement;
547    self
548  }
549  fn scale (mut self, scale : NonZero <S>) -> Self {
550    let old_width = *self.width();
551    let new_width = old_width * *scale;
552    let half_difference = (new_width - old_width) / S::two();
553    self.min -= half_difference;
554    self.max += half_difference;
555    self
556  }
557}
558
559impl <S> std::ops::Add <S> for Interval <S> where
560  S    : Field + VectorSpace <S>,
561  Self : Primitive <S, S>
562{
563  type Output = Self;
564  fn add (self, displacement : S) -> Self {
565    self.translate (displacement)
566  }
567}
568
569impl <S> std::ops::Sub <S> for Interval <S> where
570  S    : Field + VectorSpace <S>,
571  Self : Primitive <S, S>
572{
573  type Output = Self;
574  fn sub (self, displacement : S) -> Self {
575    self.translate (-displacement)
576  }
577}
578
579impl_numcast_fields!(Aabb2, min, max);
580impl <S : Ring> Aabb2 <S> {
581  /// Construct a new AABB with given min and max points.
582  ///
583  /// Debug panic if points are not min/max:
584  ///
585  /// ```should_panic
586  /// # use math_utils::geometry::*;
587  /// let b = Aabb2::with_minmax ([1.0, 1.0].into(), [0.0, 0.0].into());  // panic!
588  /// ```
589  ///
590  /// Debug panic if points are identical:
591  ///
592  /// ```should_panic
593  /// # use math_utils::geometry::*;
594  /// let b = Aabb2::with_minmax ([0.0, 0.0].into(), [0.0, 0.0].into());  // panic!
595  /// ```
596  #[inline]
597  pub fn with_minmax (min : Point2 <S>, max : Point2 <S>) -> Self where
598    S : std::fmt::Debug
599  {
600    debug_assert_ne!(min, max);
601    debug_assert_eq!(min, point2_min (&min, &max));
602    debug_assert_eq!(max, point2_max (&min, &max));
603    Aabb2 { min, max }
604  }
605  /// Construct a new AABB using the two given points to determine min/max.
606  ///
607  /// Debug panic if points are identical:
608  ///
609  /// ```should_panic
610  /// # use math_utils::geometry::*;
611  /// let b = Aabb2::from_points ([0.0, 0.0].into(), [0.0, 0.0].into());  // panic!
612  /// ```
613  #[inline]
614  pub fn from_points (a : Point2 <S>, b : Point2 <S>) -> Self where
615    S : std::fmt::Debug
616  {
617    debug_assert_ne!(a, b);
618    let min = point2_min (&a, &b);
619    let max = point2_max (&a, &b);
620    Aabb2 { min, max }
621  }
622  /// Construct the minimum AABB containing the given set of points.
623  ///
624  /// Debug panic if fewer than 2 points are given:
625  ///
626  /// ```should_panic
627  /// # use math_utils::geometry::*;
628  /// let b = Aabb2::<f32>::containing (&[]);  // panic!
629  /// ```
630  ///
631  /// ```should_panic
632  /// # use math_utils::geometry::*;
633  /// let b = Aabb2::containing (&[[0.0, 0.0].into()]);  // panic!
634  /// ```
635  #[inline]
636  pub fn containing (points : &[Point2 <S>]) -> Self where S : std::fmt::Debug {
637    debug_assert!(points.len() >= 2);
638    let mut min = point2_min (&points[0], &points[1]);
639    let mut max = point2_max (&points[0], &points[1]);
640    for point in points.iter().skip (2) {
641      min = point2_min (&min, point);
642      max = point2_max (&max, point);
643    }
644    Aabb2::with_minmax (min, max)
645  }
646  /// Create a new AABB that is the union of the two input AABBs
647  #[inline]
648  pub fn union (a : &Aabb2 <S>, b : &Aabb2 <S>) -> Self where
649    S : std::fmt::Debug
650  {
651    Aabb2::with_minmax (
652      point2_min (a.min(), b.min()),
653      point2_max (a.max(), b.max())
654    )
655  }
656  #[inline]
657  pub fn min (&self) -> &Point2 <S> {
658    &self.min
659  }
660  #[inline]
661  pub fn max (&self) -> &Point2 <S> {
662    &self.max
663  }
664  #[inline]
665  pub fn center (&self) -> Point2 <S> where S : Field {
666    Point2 ((self.min.0 + self.max.0) / S::two())
667  }
668  #[inline]
669  pub fn dimensions (&self) -> Vector2 <S> {
670    self.max.0 - self.min.0
671  }
672  #[inline]
673  pub fn width (&self) -> NonNegative <S> {
674    NonNegative::unchecked (self.max.0.x - self.min.0.x)
675  }
676  #[inline]
677  pub fn height (&self) -> NonNegative <S> {
678    NonNegative::unchecked (self.max.0.y - self.min.0.y)
679  }
680  #[inline]
681  pub fn contains (&self, point : &Point2 <S>) -> bool {
682    self.min.0.x < point.0.x && point.0.x < self.max.0.x &&
683    self.min.0.y < point.0.y && point.0.y < self.max.0.y
684  }
685  /// Clamp a given point to the AABB.
686  ///
687  /// ```
688  /// # use math_utils::geometry::*;
689  /// let b = Aabb2::from_points ([-1.0, -1.0].into(), [1.0, 1.0].into());
690  /// assert_eq!(b.clamp (&[-2.0, 0.0].into()), [-1.0, 0.0].into());
691  /// assert_eq!(b.clamp (&[ 2.0, 2.0].into()), [ 1.0, 1.0].into());
692  /// assert_eq!(b.clamp (&[-0.5, 0.5].into()), [-0.5, 0.5].into());
693  /// ```
694  pub fn clamp (&self, point : &Point2 <S>) -> Point2 <S> {
695    [ S::max (S::min (self.max.0.x, point.0.x), self.min.0.x),
696      S::max (S::min (self.max.0.y, point.0.y), self.min.0.y),
697    ].into()
698  }
699  /// Generate a random point contained in the AABB
700  ///
701  /// ```
702  /// # use rand;
703  /// # use rand_xorshift;
704  /// # use math_utils::geometry::*;
705  /// # fn main () {
706  /// use rand_xorshift;
707  /// use rand::SeedableRng;
708  /// // random sequence will be the same each time this is run
709  /// let mut rng = rand_xorshift::XorShiftRng::seed_from_u64 (0);
710  /// let aabb = Aabb2::<f32>::with_minmax (
711  ///   [-10.0, -10.0].into(),
712  ///   [ 10.0,  10.0].into());
713  /// let point = aabb.rand_point (&mut rng);
714  /// assert!(aabb.contains (&point));
715  /// let point = aabb.rand_point (&mut rng);
716  /// assert!(aabb.contains (&point));
717  /// let point = aabb.rand_point (&mut rng);
718  /// assert!(aabb.contains (&point));
719  /// let point = aabb.rand_point (&mut rng);
720  /// assert!(aabb.contains (&point));
721  /// let point = aabb.rand_point (&mut rng);
722  /// assert!(aabb.contains (&point));
723  /// # }
724  /// ```
725  #[inline]
726  pub fn rand_point <R> (&self, rng : &mut R) -> Point2 <S> where
727    S : rand::distr::uniform::SampleUniform,
728    R : rand::Rng
729  {
730    let x_range = self.min.0.x..self.max.0.x;
731    let y_range = self.min.0.y..self.max.0.y;
732    let x = if x_range.is_empty() {
733      self.min.0.x
734    } else {
735      rng.random_range (x_range)
736    };
737    let y = if y_range.is_empty() {
738      self.min.0.y
739    } else {
740      rng.random_range (y_range)
741    };
742    [x, y].into()
743  }
744  #[inline]
745  pub fn intersects (&self, other : &Aabb2 <S>) -> bool {
746    intersect::discrete_aabb2_aabb2 (self, other)
747  }
748  #[inline]
749  pub fn intersection (self, other : Aabb2 <S>) -> Option <Aabb2 <S>> where
750    S : std::fmt::Debug
751  {
752    intersect::continuous_aabb2_aabb2 (&self, &other)
753  }
754
755  pub fn corner (&self, quadrant : Quadrant) -> Point2 <S> {
756    match quadrant {
757      // upper-right
758      Quadrant::PosPos => self.max,
759      // lower-right
760      Quadrant::PosNeg => [self.max.0.x, self.min.0.y].into(),
761      // upper-left
762      Quadrant::NegPos => [self.min.0.x, self.max.0.y].into(),
763      // lower-left
764      Quadrant::NegNeg => self.min
765    }
766  }
767
768  pub fn edge (&self, direction : SignedAxis2) -> Self where
769    S : std::fmt::Debug
770  {
771    let (min, max) = match direction {
772      SignedAxis2::PosX => (
773        self.corner (Quadrant::PosNeg),
774        self.corner (Quadrant::PosPos) ),
775      SignedAxis2::NegX => (
776        self.corner (Quadrant::NegNeg),
777        self.corner (Quadrant::NegPos) ),
778      SignedAxis2::PosY => (
779        self.corner (Quadrant::NegPos),
780        self.corner (Quadrant::PosPos) ),
781      SignedAxis2::NegY => (
782        self.corner (Quadrant::NegNeg),
783        self.corner (Quadrant::PosNeg) )
784    };
785    Aabb2::with_minmax (min, max)
786  }
787
788  pub fn extrude (&self, axis : SignedAxis2, amount : Positive <S>) -> Self {
789    let (min, max) = match axis {
790      SignedAxis2::PosX => (
791        self.min + Vector2::zero().with_x (*self.width()),
792        self.max + Vector2::zero().with_x (*amount)
793      ),
794      SignedAxis2::NegX => (
795        self.min - Vector2::zero().with_x (*amount),
796        self.max - Vector2::zero().with_x (*self.height())
797      ),
798      SignedAxis2::PosY => (
799        self.min + Vector2::zero().with_y (*self.height()),
800        self.max + Vector2::zero().with_y (*amount)
801      ),
802      SignedAxis2::NegY => (
803        self.min - Vector2::zero().with_y (*amount),
804        self.max - Vector2::zero().with_y (*self.height())
805      )
806    };
807    Aabb2 { min, max }
808  }
809
810  pub fn extend (&mut self, axis : SignedAxis2, amount : Positive <S>) {
811    match axis {
812      SignedAxis2::PosX => self.max.0.x += *amount,
813      SignedAxis2::NegX => self.min.0.x -= *amount,
814      SignedAxis2::PosY => self.max.0.y += *amount,
815      SignedAxis2::NegY => self.min.0.y -= *amount
816    }
817  }
818
819  pub fn with_z (self, z : Interval <S>) -> Aabb3 <S> where S : std::fmt::Debug {
820    Aabb3::with_minmax (
821      self.min.0.with_z (z.min()).into(),
822      self.max.0.with_z (z.max()).into()
823    )
824  }
825
826  /// Increase or decrease each dimension by the given amount.
827  ///
828  /// Debug panic if any dimension would become negative:
829  /// ```should_panic
830  /// # use math_utils::geometry::*;
831  /// let x = Aabb2::with_minmax ([0.0, 0.0].into(), [1.0, 2.0].into());
832  /// let y = x.dilate (-1.0); // panic!
833  /// ```
834  pub fn dilate (mut self, amount : S) -> Self where S : std::fmt::Debug {
835    self.min -= Vector2::broadcast (amount);
836    self.max += Vector2::broadcast (amount);
837    debug_assert_ne!(self.min, self.max);
838    debug_assert_eq!(self.min, point2_min (&self.min, &self.max));
839    debug_assert_eq!(self.max, point2_max (&self.min, &self.max));
840    self
841  }
842
843  /// Project *along* the given axis.
844  ///
845  /// For example, projecting *along* the X-axis projects *onto* the Y-axis.
846  pub fn project (self, axis : Axis2) -> Interval <S> where S : std::fmt::Debug {
847    let (min, max) = match axis {
848      Axis2::X => (self.min.0.y, self.max.0.y),
849      Axis2::Y => (self.min.0.x, self.max.0.x)
850    };
851    Interval::with_minmax (min, max)
852  }
853}
854
855impl <S : Field> Primitive <S, Vector2 <S>> for Aabb2 <S> {
856  fn translate (mut self, displacement : Vector2 <S>) -> Self {
857    self.min.0 += displacement;
858    self.max.0 += displacement;
859    self
860  }
861  fn scale (mut self, scale : NonZero <S>) -> Self {
862    let center = self.center();
863    self = self.translate (-center.0);
864    self.min.0 *= *scale;
865    self.max.0 *= *scale;
866    self.translate (center.0)
867  }
868}
869
870impl <S> std::ops::Add <Vector2 <S>> for Aabb2 <S> where
871  S    : Field,
872  Self : Primitive <S, Vector2 <S>>
873{
874  type Output = Self;
875  fn add (self, displacement : Vector2 <S>) -> Self {
876    self.translate (displacement)
877  }
878}
879
880impl <S> std::ops::Sub <Vector2 <S>> for Aabb2 <S> where
881  S    : Field,
882  Self : Primitive <S, Vector2 <S>>
883{
884  type Output = Self;
885  fn sub (self, displacement : Vector2 <S>) -> Self {
886    self.translate (-displacement)
887  }
888}
889
890impl_numcast_fields!(Aabb3, min, max);
891impl <S : Ring> Aabb3 <S> {
892  /// Construct a new AABB with given min and max points.
893  ///
894  /// Debug panic if points are not min/max:
895  ///
896  /// ```should_panic
897  /// # use math_utils::geometry::*;
898  /// let b = Aabb3::with_minmax ([1.0, 1.0, 1.0].into(), [0.0, 0.0, 0.0].into());  // panic!
899  /// ```
900  ///
901  /// Debug panic if points are identical:
902  ///
903  /// ```should_panic
904  /// # use math_utils::geometry::*;
905  /// let b = Aabb3::with_minmax ([0.0, 0.0, 0.0].into(), [0.0, 0.0, 0.0].into());  // panic!
906  /// ```
907  #[inline]
908  pub fn with_minmax (min : Point3 <S>, max : Point3 <S>) -> Self where
909    S : std::fmt::Debug
910  {
911    debug_assert_ne!(min, max);
912    debug_assert_eq!(min, point3_min (&min, &max));
913    debug_assert_eq!(max, point3_max (&min, &max));
914    Aabb3 { min, max }
915  }
916  /// Construct a new AABB using the two given points to determine min/max.
917  ///
918  /// Panic if points are identical:
919  ///
920  /// ```should_panic
921  /// # use math_utils::geometry::*;
922  /// let b = Aabb3::from_points (
923  ///   [0.0, 0.0, 0.0].into(),
924  ///   [0.0, 0.0, 0.0].into());  // panic!
925  /// ```
926  #[inline]
927  pub fn from_points (a : Point3 <S>, b : Point3 <S>) -> Self where
928    S : std::fmt::Debug
929  {
930    debug_assert_ne!(a, b);
931    let min = point3_min (&a, &b);
932    let max = point3_max (&a, &b);
933    Aabb3 { min, max }
934  }
935  /// Construct the minimum AABB containing the given set of points.
936  ///
937  /// Debug panic if fewer than 2 points are given:
938  ///
939  /// ```should_panic
940  /// # use math_utils::geometry::*;
941  /// let b = Aabb3::<f32>::containing (&[]);  // panic!
942  /// ```
943  ///
944  /// ```should_panic
945  /// # use math_utils::geometry::*;
946  /// let b = Aabb3::containing (&[[0.0, 0.0, 0.0].into()]);  // panic!
947  /// ```
948  #[inline]
949  pub fn containing (points : &[Point3 <S>]) -> Self where S : std::fmt::Debug {
950    debug_assert!(points.len() >= 2);
951    let mut min = point3_min (&points[0], &points[1]);
952    let mut max = point3_max (&points[0], &points[1]);
953    for point in points.iter().skip (2) {
954      min = point3_min (&min, point);
955      max = point3_max (&max, point);
956    }
957    Aabb3::with_minmax (min, max)
958  }
959  /// Create a new AABB that is the union of the two input AABBs
960  #[inline]
961  pub fn union (a : &Aabb3 <S>, b : &Aabb3 <S>) -> Self where
962    S : std::fmt::Debug
963  {
964    Aabb3::with_minmax (
965      point3_min (a.min(), b.min()),
966      point3_max (a.max(), b.max())
967    )
968  }
969  #[inline]
970  pub fn min (&self) -> &Point3 <S> {
971    &self.min
972  }
973  #[inline]
974  pub fn max (&self) -> &Point3 <S> {
975    &self.max
976  }
977  #[inline]
978  pub fn center (&self) -> Point3 <S> where S : Field {
979    Point3 ((self.min.0 + self.max.0) / S::two())
980  }
981  #[inline]
982  pub fn dimensions (&self) -> Vector3 <S> {
983    self.max.0 - self.min.0
984  }
985  /// X extent
986  #[inline]
987  pub fn width (&self) -> NonNegative <S> {
988    NonNegative::unchecked (self.max.0.x - self.min.0.x)
989  }
990  /// Y extent
991  #[inline]
992  pub fn depth (&self) -> NonNegative <S> {
993    NonNegative::unchecked (self.max.0.y - self.min.0.y)
994  }
995  /// Z extent
996  #[inline]
997  pub fn height (&self) -> NonNegative <S> {
998    NonNegative::unchecked (self.max.0.z - self.min.0.z)
999  }
1000  #[inline]
1001  pub fn contains (&self, point : &Point3 <S>) -> bool {
1002    self.min.0.x < point.0.x && point.0.x < self.max.0.x &&
1003    self.min.0.y < point.0.y && point.0.y < self.max.0.y &&
1004    self.min.0.z < point.0.z && point.0.z < self.max.0.z
1005  }
1006  /// Clamp a given point to the AABB.
1007  ///
1008  /// ```
1009  /// # use math_utils::geometry::*;
1010  /// let b = Aabb3::with_minmax ([-1.0, -1.0, -1.0].into(), [1.0, 1.0, 1.0].into());
1011  /// assert_eq!(b.clamp (&[-2.0, 0.0, 0.0].into()), [-1.0, 0.0, 0.0].into());
1012  /// assert_eq!(b.clamp (&[ 2.0, 2.0, 0.0].into()), [ 1.0, 1.0, 0.0].into());
1013  /// assert_eq!(b.clamp (&[-1.0, 2.0, 3.0].into()), [-1.0, 1.0, 1.0].into());
1014  /// assert_eq!(b.clamp (&[-0.5, 0.5, 0.0].into()), [-0.5, 0.5, 0.0].into());
1015  /// ```
1016  pub fn clamp (&self, point : &Point3 <S>) -> Point3 <S> {
1017    [ S::max (S::min (self.max.0.x, point.0.x), self.min.0.x),
1018      S::max (S::min (self.max.0.y, point.0.y), self.min.0.y),
1019      S::max (S::min (self.max.0.z, point.0.z), self.min.0.z)
1020    ].into()
1021  }
1022  /// Generate a random point contained in the AABB
1023  ///
1024  /// ```
1025  /// # use rand;
1026  /// # use rand_xorshift;
1027  /// # use math_utils::geometry::*;
1028  /// # fn main () {
1029  /// use rand_xorshift;
1030  /// use rand::SeedableRng;
1031  /// // random sequence will be the same each time this is run
1032  /// let mut rng = rand_xorshift::XorShiftRng::seed_from_u64 (0);
1033  /// let aabb = Aabb3::<f32>::with_minmax (
1034  ///   [-10.0, -10.0, -10.0].into(),
1035  ///   [ 10.0,  10.0,  10.0].into());
1036  /// let point = aabb.rand_point (&mut rng);
1037  /// assert!(aabb.contains (&point));
1038  /// let point = aabb.rand_point (&mut rng);
1039  /// assert!(aabb.contains (&point));
1040  /// let point = aabb.rand_point (&mut rng);
1041  /// assert!(aabb.contains (&point));
1042  /// let point = aabb.rand_point (&mut rng);
1043  /// assert!(aabb.contains (&point));
1044  /// let point = aabb.rand_point (&mut rng);
1045  /// assert!(aabb.contains (&point));
1046  /// # }
1047  /// ```
1048  #[inline]
1049  pub fn rand_point <R> (&self, rng : &mut R) -> Point3 <S> where
1050    S : rand::distr::uniform::SampleUniform,
1051    R : rand::Rng
1052  {
1053    let x_range = self.min.0.x..self.max.0.x;
1054    let y_range = self.min.0.y..self.max.0.y;
1055    let z_range = self.min.0.z..self.max.0.z;
1056    let x = if x_range.is_empty() {
1057      self.min.0.x
1058    } else {
1059      rng.random_range (x_range)
1060    };
1061    let y = if y_range.is_empty() {
1062      self.min.0.y
1063    } else {
1064      rng.random_range (y_range)
1065    };
1066    let z = if z_range.is_empty() {
1067      self.min.0.z
1068    } else {
1069      rng.random_range (z_range)
1070    };
1071    [x, y, z].into()
1072  }
1073  #[inline]
1074  pub fn intersects (&self, other : &Aabb3 <S>) -> bool {
1075    intersect::discrete_aabb3_aabb3 (self, other)
1076  }
1077  #[inline]
1078  pub fn intersection (self, other : Aabb3 <S>) -> Option <Aabb3 <S>> where
1079    S : std::fmt::Debug
1080  {
1081    intersect::continuous_aabb3_aabb3 (&self, &other)
1082  }
1083
1084  pub fn corner (&self, octant : Octant) -> Point3 <S> {
1085    match octant {
1086      Octant::PosPosPos => self.max,
1087      Octant::NegPosPos => [self.min.0.x, self.max.0.y, self.max.0.z].into(),
1088      Octant::PosNegPos => [self.max.0.x, self.min.0.y, self.max.0.z].into(),
1089      Octant::NegNegPos => [self.min.0.x, self.min.0.y, self.max.0.z].into(),
1090      Octant::PosPosNeg => [self.max.0.x, self.max.0.y, self.min.0.z].into(),
1091      Octant::NegPosNeg => [self.min.0.x, self.max.0.y, self.min.0.z].into(),
1092      Octant::PosNegNeg => [self.max.0.x, self.min.0.y, self.min.0.z].into(),
1093      Octant::NegNegNeg => self.min
1094    }
1095  }
1096
1097  pub fn face (&self, direction : SignedAxis3) -> Self where
1098    S : std::fmt::Debug
1099  {
1100    let (min, max) = match direction {
1101      SignedAxis3::PosX => (
1102        self.corner (Octant::PosNegNeg),
1103        self.corner (Octant::PosPosPos) ),
1104      SignedAxis3::NegX => (
1105        self.corner (Octant::NegNegNeg),
1106        self.corner (Octant::NegPosPos) ),
1107      SignedAxis3::PosY => (
1108        self.corner (Octant::NegPosNeg),
1109        self.corner (Octant::PosPosPos) ),
1110      SignedAxis3::NegY => (
1111        self.corner (Octant::NegNegNeg),
1112        self.corner (Octant::PosNegPos) ),
1113      SignedAxis3::PosZ => (
1114        self.corner (Octant::NegNegPos),
1115        self.corner (Octant::PosPosPos) ),
1116      SignedAxis3::NegZ => (
1117        self.corner (Octant::NegNegNeg),
1118        self.corner (Octant::PosPosNeg) )
1119    };
1120    Aabb3::with_minmax (min, max)
1121  }
1122
1123  pub fn extrude (&self, axis : SignedAxis3, amount : Positive <S>) -> Self {
1124    let (min, max) = match axis {
1125      SignedAxis3::PosX => (
1126        self.min + Vector3::zero().with_x (*self.width()),
1127        self.max + Vector3::zero().with_x (*amount)
1128      ),
1129      SignedAxis3::NegX => (
1130        self.min - Vector3::zero().with_x (*amount),
1131        self.max - Vector3::zero().with_x (*self.width())
1132      ),
1133      SignedAxis3::PosY => (
1134        self.min + Vector3::zero().with_y (*self.depth()),
1135        self.max + Vector3::zero().with_y (*amount)
1136      ),
1137      SignedAxis3::NegY => (
1138        self.min - Vector3::zero().with_y (*amount),
1139        self.max - Vector3::zero().with_y (*self.depth())
1140      ),
1141      SignedAxis3::PosZ => (
1142        self.min + Vector3::zero().with_z (*self.height()),
1143        self.max + Vector3::zero().with_z (*amount)
1144      ),
1145      SignedAxis3::NegZ => (
1146        self.min - Vector3::zero().with_z (*amount),
1147        self.max - Vector3::zero().with_z (*self.height())
1148      )
1149    };
1150    Aabb3 { min, max }
1151  }
1152
1153  pub fn extend (&mut self, axis : SignedAxis3, amount : Positive <S>) {
1154    match axis {
1155      SignedAxis3::PosX => self.max.0.x += *amount,
1156      SignedAxis3::NegX => self.min.0.x -= *amount,
1157      SignedAxis3::PosY => self.max.0.y += *amount,
1158      SignedAxis3::NegY => self.min.0.y -= *amount,
1159      SignedAxis3::PosZ => self.max.0.z += *amount,
1160      SignedAxis3::NegZ => self.min.0.z -= *amount
1161    }
1162  }
1163
1164  /// Increase or decrease each dimension by the given amount.
1165  ///
1166  /// Debug panic if any dimension would become negative:
1167  /// ```should_panic
1168  /// # use math_utils::geometry::*;
1169  /// let x = Aabb3::with_minmax ([0.0, 0.0, 0.0].into(), [1.0, 2.0, 3.0].into());
1170  /// let y = x.dilate (-1.0); // panic!
1171  /// ```
1172  pub fn dilate (mut self, amount : S) -> Self where S : std::fmt::Debug {
1173    self.min -= Vector3::broadcast (amount);
1174    self.max += Vector3::broadcast (amount);
1175    debug_assert_ne!(self.min, self.max);
1176    debug_assert_eq!(self.min, point3_min (&self.min, &self.max));
1177    debug_assert_eq!(self.max, point3_max (&self.min, &self.max));
1178    self
1179  }
1180
1181  /// Project *along* the given axis.
1182  ///
1183  /// For example, projecting *along* the X-axis projects *onto* the YZ-plane.
1184  pub fn project (self, axis : Axis3) -> Aabb2 <S> where S : std::fmt::Debug {
1185    let (min, max) = match axis {
1186      Axis3::X => ([self.min.0.y, self.min.0.z], [self.max.0.y, self.max.0.z]),
1187      Axis3::Y => ([self.min.0.x, self.min.0.z], [self.max.0.x, self.max.0.z]),
1188      Axis3::Z => ([self.min.0.x, self.min.0.y], [self.max.0.x, self.max.0.y]),
1189    };
1190    Aabb2::with_minmax (min.into(), max.into())
1191  }
1192}
1193
1194impl <S : Field> Primitive <S, Vector3 <S>> for Aabb3 <S> {
1195  fn translate (mut self, displacement : Vector3 <S>) -> Self {
1196    self.min.0 += displacement;
1197    self.max.0 += displacement;
1198    self
1199  }
1200  fn scale (mut self, scale : NonZero <S>) -> Self {
1201    let center = self.center();
1202    self = self.translate (-center.0);
1203    self.min.0 *= *scale;
1204    self.max.0 *= *scale;
1205    self.translate (center.0)
1206  }
1207}
1208
1209impl <S> std::ops::Add <Vector3 <S>> for Aabb3 <S> where
1210  S    : Field,
1211  Self : Primitive <S, Vector3 <S>>
1212{
1213  type Output = Self;
1214  fn add (self, displacement : Vector3 <S>) -> Self {
1215    self.translate (displacement)
1216  }
1217}
1218
1219impl <S> std::ops::Sub <Vector3 <S>> for Aabb3 <S> where
1220  S    : Field,
1221  Self : Primitive <S, Vector3 <S>>
1222{
1223  type Output = Self;
1224  fn sub (self, displacement : Vector3 <S>) -> Self {
1225    self.translate (-displacement)
1226  }
1227}
1228
1229impl_numcast_fields!(Capsule3, center, half_height, radius);
1230impl <S : Ring> Capsule3 <S> {
1231  pub fn aabb3 (&self) -> Aabb3 <S> where S : Real + std::fmt::Debug {
1232    use shape::Aabb;
1233    let shape_aabb = shape::Capsule {
1234      radius:      self.radius,
1235      half_height: self.half_height
1236    }.aabb();
1237    let center_vec = self.center.0;
1238    let min        = *shape_aabb.min() + center_vec;
1239    let max        = *shape_aabb.max() + center_vec;
1240    Aabb3::with_minmax (min, max)
1241  }
1242  /// Return the (upper sphere, cylinder, lower sphere) making up this capsule
1243  pub fn decompose (&self)
1244    -> (Sphere3 <S>, Option <Cylinder3 <S>>, Sphere3 <S>)
1245  {
1246    let cylinder = if *self.half_height > S::zero() {
1247      Some (Cylinder3 {
1248        center:      self.center,
1249        radius:      self.radius,
1250        half_height: Positive::unchecked (*self.half_height)
1251      })
1252    } else {
1253      None
1254    };
1255    let upper_sphere = Sphere3 {
1256      center: self.center + (Vector3::unit_z() * *self.half_height),
1257      radius: self.radius
1258    };
1259    let lower_sphere = Sphere3 {
1260      center: self.center - (Vector3::unit_z() * *self.half_height),
1261      radius: self.radius
1262    };
1263    (upper_sphere, cylinder, lower_sphere)
1264  }
1265}
1266
1267// TODO: impl Primitive for Capsule3
1268
1269impl_numcast_fields!(Cylinder3, center, half_height, radius);
1270impl <S : Ring> Cylinder3 <S> {
1271  pub fn aabb3 (&self) -> Aabb3 <S> where S : Real + std::fmt::Debug {
1272    use shape::Aabb;
1273    let shape_aabb = shape::Cylinder::unchecked (*self.radius, *self.half_height)
1274      .aabb();
1275    let center_vec = self.center.0;
1276    let min        = *shape_aabb.min() + center_vec;
1277    let max        = *shape_aabb.max() + center_vec;
1278    Aabb3::with_minmax (min, max)
1279  }
1280}
1281
1282// TODO: impl Primitive for Cylinder3
1283
1284impl_numcast_fields!(Line2, base, direction);
1285impl <S : Real> Line2 <S> {
1286  /// Construct a new 2D line
1287  #[inline]
1288  pub fn new (base : Point2 <S>, direction : Unit2 <S>) -> Self {
1289    Line2 { base, direction }
1290  }
1291  #[inline]
1292  pub fn point (&self, t : S) -> Point2 <S> {
1293    self.base + *self.direction * t
1294  }
1295  #[inline]
1296  pub fn intersect_aabb (&self, aabb : &Aabb2 <S>)
1297    -> Option <((S, Point2 <S>), (S, Point2 <S>))>
1298  where
1299    S : std::fmt::Debug
1300  {
1301    intersect::continuous_line2_aabb2 (self, aabb)
1302  }
1303  #[inline]
1304  pub fn intersect_sphere (&self, sphere : &Sphere2 <S>)
1305    -> Option <((S, Point2 <S>), (S, Point2 <S>))>
1306  {
1307    intersect::continuous_line2_sphere2 (self, sphere)
1308  }
1309}
1310impl <S : Real> Default for Line2 <S> {
1311  fn default() -> Self {
1312    Line2 {
1313      base:      Point2::origin(),
1314      direction: Unit2::axis_y()
1315    }
1316  }
1317}
1318
1319// TODO: impl Primitive for Line2
1320
1321impl_numcast_fields!(Line3, base, direction);
1322impl <S : Real> Line3 <S> {
1323  /// Consructs a new 3D line with given base point and direction vector
1324  #[inline]
1325  pub fn new (base : Point3 <S>, direction : Unit3 <S>) -> Self {
1326    Line3 { base, direction }
1327  }
1328  #[inline]
1329  pub fn point (&self, t : S) -> Point3 <S> {
1330    self.base + *self.direction * t
1331  }
1332  #[inline]
1333  pub fn intersect_plane (&self, plane : &Plane3 <S>)
1334    -> Option <(S, Point3 <S>)>
1335  where
1336    S : approx::RelativeEq
1337  {
1338    intersect::continuous_line3_plane3 (self, plane)
1339  }
1340  #[inline]
1341  pub fn intersect_aabb (&self, aabb : &Aabb3 <S>)
1342    -> Option <((S, Point3 <S>), (S, Point3 <S>))>
1343  where
1344    S : num_traits::Float + approx::RelativeEq <Epsilon=S> + std::fmt::Debug
1345  {
1346    intersect::continuous_line3_aabb3 (self, aabb)
1347  }
1348  #[inline]
1349  pub fn intersect_sphere (&self, sphere : &Sphere3 <S>)
1350    -> Option <((S, Point3 <S>), (S, Point3 <S>))>
1351  {
1352    intersect::continuous_line3_sphere3 (self, sphere)
1353  }
1354}
1355impl <S : Real> Default for Line3 <S> {
1356  fn default() -> Self {
1357    Line3 {
1358      base:      Point3::origin(),
1359      direction: Unit3::axis_z()
1360    }
1361  }
1362}
1363
1364// TODO: impl Primitive for Line3
1365
1366impl_numcast_fields!(Plane3, base, normal);
1367impl <S : Real> Plane3 <S> {
1368  /// Consructs a new 3D plane with given base point and normal vector
1369  #[inline]
1370  pub fn new (base : Point3 <S>, normal : Unit3 <S>) -> Self {
1371    Plane3 { base, normal }
1372  }
1373}
1374impl <S : Real> Default for Plane3 <S> {
1375  fn default() -> Self {
1376    Plane3 {
1377      base:   Point3::origin(),
1378      normal: Unit3::axis_z()
1379    }
1380  }
1381}
1382
1383// TODO: impl Primitive for Plane3
1384
1385impl_numcast_fields!(Sphere2, center, radius);
1386impl <S : Ring> Sphere2 <S> {
1387  /// Unit circle
1388  #[inline]
1389  pub fn unit() -> Self where S : Field {
1390    Sphere2 {
1391      center: Point2::origin(),
1392      radius: num::One::one()
1393    }
1394  }
1395  /// Discrete intersection test with another sphere
1396  #[inline]
1397  pub fn intersects (&self, other : &Self) -> bool {
1398    intersect::discrete_sphere2_sphere2 (self, other)
1399  }
1400}
1401
1402// TODO: impl Primitive for Sphere2
1403
1404impl_numcast_fields!(Sphere3, center, radius);
1405impl <S : Ring> Sphere3 <S> {
1406  /// Unit sphere
1407  #[inline]
1408  pub fn unit() -> Self where S : Field {
1409    Sphere3 {
1410      center: Point3::origin(),
1411      radius: num::One::one()
1412    }
1413  }
1414  /// Discrete intersection test with another sphere
1415  #[inline]
1416  pub fn intersects (&self, other : &Self) -> bool {
1417    intersect::discrete_sphere3_sphere3 (self, other)
1418  }
1419}
1420
1421// TODO: impl Primitive for Sphere3
1422
1423////////////////////////////////////////////////////////////////////////////////
1424//  tests                                                                     //
1425////////////////////////////////////////////////////////////////////////////////
1426
1427#[cfg(test)]
1428mod tests {
1429  use approx::{assert_relative_eq, assert_ulps_eq};
1430  use super::*;
1431
1432  #[test]
1433  fn test_project_2d_point_on_line() {
1434    use Unit2;
1435    let point : Point2 <f64> = [2.0, 2.0].into();
1436    let line = Line2::<f64>::new ([0.0, 0.0].into(), Unit2::axis_x());
1437    assert_eq!(project_2d_point_on_line (&point, &line), [2.0, 0.0].into());
1438    let line = Line2::<f64>::new ([0.0, 0.0].into(), Unit2::axis_y());
1439    assert_eq!(project_2d_point_on_line (&point, &line), [0.0, 2.0].into());
1440    let point : Point2 <f64> = [0.0, 1.0].into();
1441    let line = Line2::<f64>::new (
1442      [0.0, -1.0].into(), Unit2::normalize ([1.0, 1.0].into()));
1443    assert_relative_eq!(
1444      project_2d_point_on_line (&point, &line), [1.0, 0.0].into());
1445    // the answer should be the same for lines with equivalent definitions
1446    let point : Point2 <f64> = [1.0, 3.0].into();
1447    let line_a = Line2::<f64>::new (
1448      [0.0, -1.0].into(), Unit2::normalize ([2.0, 1.0].into()));
1449    let line_b = Line2::<f64>::new (
1450      [2.0, 0.0].into(),  Unit2::normalize ([2.0, 1.0].into()));
1451    assert_relative_eq!(
1452      project_2d_point_on_line (&point, &line_a),
1453      project_2d_point_on_line (&point, &line_b));
1454    let line_a = Line2::<f64>::new (
1455      [0.0, 0.0].into(),   Unit2::normalize ([1.0, 1.0].into()));
1456    let line_b = Line2::<f64>::new (
1457      [-2.0, -2.0].into(), Unit2::normalize ([1.0, 1.0].into()));
1458    assert_ulps_eq!(
1459      project_2d_point_on_line (&point, &line_a),
1460      project_2d_point_on_line (&point, &line_b)
1461    );
1462    assert_relative_eq!(
1463      project_2d_point_on_line (&point, &line_a),
1464      [2.0, 2.0].into());
1465  }
1466
1467  #[test]
1468  fn test_project_3d_point_on_line() {
1469    use Unit3;
1470    // all the tests from 2d projection with 0.0 for the Z component
1471    let point : Point3 <f64> = [2.0, 2.0, 0.0].into();
1472    let line = Line3::<f64>::new ([0.0, 0.0, 0.0].into(), Unit3::axis_x());
1473    assert_eq!(project_3d_point_on_line (&point, &line), [2.0, 0.0, 0.0].into());
1474    let line = Line3::<f64>::new ([0.0, 0.0, 0.0].into(), Unit3::axis_y());
1475    assert_eq!(project_3d_point_on_line (&point, &line), [0.0, 2.0, 0.0].into());
1476    let point : Point3 <f64> = [0.0, 1.0, 0.0].into();
1477    let line = Line3::<f64>::new (
1478      [0.0, -1.0, 0.0].into(), Unit3::normalize ([1.0, 1.0, 0.0].into()));
1479    assert_relative_eq!(
1480      project_3d_point_on_line (&point, &line), [1.0, 0.0, 0.0].into());
1481    // the answer should be the same for lines with equivalent definitions
1482    let point : Point3 <f64> = [1.0, 3.0, 0.0].into();
1483    let line_a = Line3::<f64>::new (
1484      [0.0, -1.0, 0.0].into(), Unit3::normalize ([2.0, 1.0, 0.0].into()));
1485    let line_b = Line3::<f64>::new (
1486      [2.0, 0.0, 0.0].into(),  Unit3::normalize ([2.0, 1.0, 0.0].into()));
1487    assert_relative_eq!(
1488      project_3d_point_on_line (&point, &line_a),
1489      project_3d_point_on_line (&point, &line_b));
1490    let line_a = Line3::<f64>::new (
1491      [0.0, 0.0, 0.0].into(), Unit3::normalize ([1.0, 1.0, 0.0].into()));
1492    let line_b = Line3::<f64>::new (
1493      [2.0, 2.0, 0.0].into(), Unit3::normalize ([1.0, 1.0, 0.0].into()));
1494    assert_relative_eq!(
1495      project_3d_point_on_line (&point, &line_a),
1496      project_3d_point_on_line (&point, &line_b));
1497    assert_relative_eq!(
1498      project_3d_point_on_line (&point, &line_a),
1499      [2.0, 2.0, 0.0].into());
1500    // more tests
1501    let point : Point3 <f64> = [0.0, 0.0, 2.0].into();
1502    let line_a = Line3::<f64>::new (
1503      [-4.0, -4.0, -4.0].into(), Unit3::normalize ([1.0, 1.0, 1.0].into()));
1504    let line_b = Line3::<f64>::new (
1505      [4.0, 4.0, 4.0].into(),    Unit3::normalize ([1.0, 1.0, 1.0].into()));
1506    assert_relative_eq!(
1507      project_3d_point_on_line (&point, &line_a),
1508      project_3d_point_on_line (&point, &line_b),
1509      epsilon = 0.00000000000001
1510    );
1511    assert_relative_eq!(
1512      project_3d_point_on_line (&point, &line_a),
1513      [2.0/3.0, 2.0/3.0, 2.0/3.0].into(),
1514      epsilon = 0.00000000000001
1515    );
1516  }
1517
1518} // end tests