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 hull;
17mod simplex;
18
19pub use self::hull::*;
20pub use self::simplex::{simplex2, simplex3, Simplex2, Simplex3};
21pub use simplex2::Segment     as Segment2;
22pub use simplex3::Segment     as Segment3;
23pub use simplex2::Triangle    as Triangle2;
24pub use simplex3::Triangle    as Triangle3;
25pub use simplex3::Tetrahedron as Tetrahedron3;
26
27pub trait Primitive <S, V> where
28  S : Field,
29  V : VectorSpace <S>
30{
31  fn translate (self, displacement : V) -> Self;
32  fn scale     (self, scale : NonZero <S>) -> Self;
33}
34
35////////////////////////////////////////////////////////////////////////////////
36//  structs                                                                   //
37////////////////////////////////////////////////////////////////////////////////
38
39/// 1D interval
40#[cfg_attr(feature = "derive_serdes", derive(Serialize, Deserialize))]
41#[derive(Clone, Copy, Debug, PartialEq)]
42pub struct Interval <S> {
43  min : S,
44  max : S
45}
46
47/// 2D axis-aligned bounding box.
48///
49/// Min/max can be co-linear:
50/// ```
51/// # use math_utils::geometry::Aabb2;
52/// let x = Aabb2::with_minmax (
53///   [0.0, 0.0].into(),
54///   [0.0, 1.0].into()
55/// );
56/// ```
57/// But not identical:
58/// ```should_panic
59/// # use math_utils::geometry::Aabb2;
60/// // panic!
61/// let x = Aabb2::with_minmax (
62///   [0.0, 0.0].into(),
63///   [0.0, 0.0].into()
64/// );
65/// ```
66#[cfg_attr(feature = "derive_serdes", derive(Serialize, Deserialize))]
67#[derive(Clone, Copy, Debug, PartialEq)]
68pub struct Aabb2 <S> {
69  min : Point2 <S>,
70  max : Point2 <S>
71}
72
73/// 3D axis-aligned bounding box.
74///
75/// See also `shape::Aabb` trait for primitive and shape types that can compute
76/// a 3D AABB.
77///
78/// Min/max can be co-planar or co-linear:
79/// ```
80/// # use math_utils::geometry::Aabb3;
81/// let x = Aabb3::with_minmax (
82///   [0.0, 0.0, 0.0].into(),
83///   [0.0, 1.0, 1.0].into()
84/// );
85/// let y = Aabb3::with_minmax (
86///   [0.0, 0.0, 0.0].into(),
87///   [0.0, 1.0, 0.0].into()
88/// );
89/// ```
90/// But not identical:
91/// ```should_panic
92/// # use math_utils::geometry::Aabb3;
93/// let x = Aabb3::with_minmax (
94///   [0.0, 0.0, 0.0].into(),
95///   [0.0, 0.0, 0.0].into()
96/// );
97/// ```
98#[cfg_attr(feature = "derive_serdes", derive(Serialize, Deserialize))]
99#[derive(Clone, Copy, Debug, PartialEq)]
100pub struct Aabb3 <S> {
101  min : Point3 <S>,
102  max : Point3 <S>
103}
104
105/// 2D oriented bounding box
106#[cfg_attr(feature = "derive_serdes", derive(Serialize, Deserialize))]
107#[derive(Clone, Copy, Debug, PartialEq)]
108pub struct Obb2 <S> {
109  center        : Point2 <S>,
110  half_extent_x : Positive <S>,
111  half_extent_y : Positive <S>,
112  orientation   : AngleWrapped <S>
113}
114
115/// 3D oriented bounding box
116#[cfg_attr(feature = "derive_serdes", derive(Serialize, Deserialize))]
117#[derive(Clone, Copy, Debug, PartialEq)]
118pub struct Obb3 <S> {
119  center        : Point3 <S>,
120  half_extent_x : Positive <S>,
121  half_extent_y : Positive <S>,
122  half_extent_z : Positive <S>,
123  orientation   : Angles3 <S>
124}
125
126/// 3D Z-axis-aligned cylinder
127#[cfg_attr(feature = "derive_serdes", derive(Serialize, Deserialize))]
128#[derive(Clone, Copy, Debug, PartialEq)]
129pub struct Cylinder3 <S> {
130  pub center      : Point3   <S>,
131  pub half_height : Positive <S>,
132  pub radius      : Positive <S>
133}
134
135/// 3D Z-axis-aligned capsule
136#[cfg_attr(feature = "derive_serdes", derive(Serialize, Deserialize))]
137#[derive(Clone, Copy, Debug, PartialEq)]
138pub struct Capsule3 <S> {
139  pub center      : Point3   <S>,
140  pub half_height : Positive <S>,
141  pub radius      : Positive <S>
142}
143
144/// An infinitely extended line in 2D space defined by a base point and
145/// normalized direction
146#[cfg_attr(feature = "derive_serdes", derive(Serialize, Deserialize))]
147#[derive(Clone, Copy, Debug, PartialEq)]
148pub struct Line2 <S> {
149  pub base      : Point2 <S>,
150  pub direction : Unit2  <S>
151}
152
153/// An infinitely extended line in 3D space defined by a base point and
154/// normalized direction
155#[cfg_attr(feature = "derive_serdes", derive(Serialize, Deserialize))]
156#[derive(Clone, Copy, Debug, PartialEq)]
157pub struct Line3 <S> {
158  pub base      : Point3 <S>,
159  pub direction : Unit3  <S>
160}
161
162/// A plane in 3D space defined by a base point and (unit) normal vector
163#[cfg_attr(feature = "derive_serdes", derive(Serialize, Deserialize))]
164#[derive(Clone, Copy, Debug, PartialEq)]
165pub struct Plane3 <S> {
166  pub base   : Point3 <S>,
167  pub normal : Unit3  <S>
168}
169
170/// Sphere in 2D space (a circle)
171#[cfg_attr(feature = "derive_serdes", derive(Serialize, Deserialize))]
172#[derive(Clone, Copy, Debug, PartialEq)]
173pub struct Sphere2 <S> {
174  pub center : Point2   <S>,
175  pub radius : Positive <S>
176}
177
178/// Sphere in 3D space
179#[cfg_attr(feature = "derive_serdes", derive(Serialize, Deserialize))]
180#[derive(Clone, Copy, Debug, PartialEq)]
181pub struct Sphere3 <S> {
182  pub center : Point3   <S>,
183  pub radius : Positive <S>
184}
185
186/// Describes points of a 2D convex hull that support a rectangle and the area
187/// of the rectangle
188#[derive(Clone, Debug)]
189struct Rect <S> {
190  pub edge    : Vector2 <S>,
191  pub perp    : Vector2 <S>,
192  /// bottom, right, top, left
193  pub indices : [u32; 4],
194  pub area    : S,
195  pub edge_len_squared : S
196}
197
198////////////////////////////////////////////////////////////////////////////////
199//  functions                                                                 //
200////////////////////////////////////////////////////////////////////////////////
201
202/// Returns true when three points lie on the same line in 2D space.
203///
204/// ```
205/// # use math_utils::geometry::colinear_2d;
206/// assert!(colinear_2d (
207///   &[-1.0, -1.0].into(),
208///   &[ 0.0,  0.0].into(),
209///   &[ 1.0,  1.0].into())
210/// );
211/// assert!(!colinear_2d (
212///   &[-1.0, -1.0].into(),
213///   &[ 0.0,  1.0].into(),
214///   &[ 1.0, -1.0].into())
215/// );
216/// ```
217pub fn colinear_2d <S> (
218  a : &Point2 <S>,
219  b : &Point2 <S>,
220  c : &Point2 <S>
221) -> bool where
222  S : Field + approx::RelativeEq
223{
224  // early exit: a pair of points are equal
225  if a == b || a == c || b == c {
226    true
227  } else {
228    let a = Vector3::from_point_2d (a.0).into();
229    let b = Vector3::from_point_2d (b.0).into();
230    let c = Vector3::from_point_2d (c.0).into();
231    let determinant = determinant_3d (&a, &b, &c);
232    // a zero determinant indicates that the points are colinear
233    approx::relative_eq!(S::zero(), determinant)
234  }
235}
236
237/// Returns true when three points lie on the same line in 3D space.
238///
239/// ```
240/// # use math_utils::geometry::colinear_3d;
241/// assert!(colinear_3d (
242///   &[-1.0, -1.0, -1.0].into(),
243///   &[ 0.0,  0.0,  0.0].into(),
244///   &[ 1.0,  1.0,  1.0].into())
245/// );
246/// ```
247pub fn colinear_3d <S> (a : &Point3 <S>, b : &Point3 <S>, c : &Point3 <S>)
248  -> bool
249where
250  S : Field + Sqrt + approx::RelativeEq
251{
252  // early exit: a pair of points are equal
253  if a == b || a == c || b == c {
254    true
255  } else {
256    let determinant = determinant_3d (a, b, c);
257    // a zero determinant indicates that the points are coplanar
258    if approx::relative_eq!(S::zero(), determinant) {
259      approx::relative_eq!(S::zero(), triangle_3d_area2 (a, b, c))
260    } else {
261      false
262    }
263  }
264}
265
266/// Returns true when four points lie on the same plane in 3D space.
267///
268/// ```
269/// # use math_utils::geometry::coplanar_3d;
270/// assert!(coplanar_3d (
271///   &[-1.0, -1.0, -1.0].into(),
272///   &[ 1.0,  1.0,  1.0].into(),
273///   &[-1.0,  1.0,  0.0].into(),
274///   &[ 1.0, -1.0,  0.0].into()
275/// ));
276/// ```
277pub fn coplanar_3d <S> (
278  a : &Point3 <S>,
279  b : &Point3 <S>,
280  c : &Point3 <S>,
281  d : &Point3 <S>
282) -> bool where
283  S : Ring + approx::RelativeEq
284{
285  // early exit: a pair of points are equal
286  if a == b || a == c || a == d || b == c || b == d || c == d {
287    true
288  } else {
289    let ab_cross_ac_dot_ad = (*b-a).cross (*c-a).dot (*d-a);
290    approx::relative_eq!(S::zero(), ab_cross_ac_dot_ad)
291  }
292}
293
294/// Square area of three points in 3D space.
295///
296/// Uses a numerically stable Heron's formula:
297/// <https://en.wikipedia.org/wiki/Heron%27s_formula#Numerical_stability>
298///
299/// ```
300/// # use math_utils::approx::assert_relative_eq;
301/// # use math_utils::geometry::triangle_3d_area2;
302/// assert_relative_eq!(
303///   3.0/4.0,
304///   triangle_3d_area2 (
305///     &[-1.0,  0.0,  0.0].into(),
306///     &[ 0.0,  0.0,  1.0].into(),
307///     &[ 0.0,  1.0,  0.0].into())
308/// );
309/// ```
310///
311/// If the area squared is zero then the points are colinear:
312///
313/// ```
314/// # use math_utils::geometry::triangle_3d_area2;
315/// assert_eq!(
316///   0.0,
317///   triangle_3d_area2 (
318///     &[-1.0, -1.0, -1.0].into(),
319///     &[ 0.0,  0.0,  0.0].into(),
320///     &[ 1.0,  1.0,  1.0].into())
321/// );
322/// ```
323pub fn triangle_3d_area2 <S> (a : &Point3 <S>, b : &Point3 <S>, c : &Point3 <S>)
324  -> S
325where
326  S : Field + Sqrt
327{
328  if a == b || a == c || b == c {
329    return S::zero()
330  }
331  // compute the length of each side
332  let ab_mag = (*b-a).norm();
333  let ac_mag = (*c-a).norm();
334  let bc_mag = (*c-b).norm();
335  // order as max >= mid >= min
336  let max = S::max (S::max (ab_mag, ac_mag), bc_mag);
337  let min = S::min (S::min (ab_mag, ac_mag), bc_mag);
338  let mid = if min <= ab_mag && ab_mag <= max {
339    ab_mag
340  } else if min <= ac_mag && ac_mag <= max {
341    ac_mag
342  } else {
343    bc_mag
344  };
345  let two = S::one() + S::one();
346  let sixteen = two * two * two * two;
347  // 1.0/16.0
348  let frac = S::one() / sixteen;
349  frac
350    * (max + (mid + min))
351    * (min - (max - mid))
352    * (min + (max - mid))
353    * (max + (mid - min))
354}
355
356/// Computes the determinant of a matrix formed by the three points as columns
357#[inline]
358pub fn determinant_3d <S : Ring> (
359  a : &Point3 <S>,
360  b : &Point3 <S>,
361  c : &Point3 <S>
362) -> S {
363  Matrix3::from_col_arrays ([
364    a.0.into_array(), b.0.into_array(), c.0.into_array()
365  ]).determinant()
366}
367
368/// Coordinate-wise min
369pub fn point2_min <S : Ring> (a : &Point2 <S>, b : &Point2 <S>) -> Point2 <S> {
370  Vector2::partial_min (a.0, b.0).into()
371}
372
373/// Coordinate-wise max
374pub fn point2_max <S : Ring> (a : &Point2 <S>, b : &Point2 <S>) -> Point2 <S> {
375  Vector2::partial_max (a.0, b.0).into()
376}
377
378/// Coordinate-wise min
379pub fn point3_min <S : Ring> (a : &Point3 <S>, b : &Point3 <S>) -> Point3 <S> {
380  Vector3::partial_min (a.0, b.0).into()
381}
382
383/// Coordinate-wise max
384pub fn point3_max <S : Ring> (a : &Point3 <S>, b : &Point3 <S>) -> Point3 <S> {
385  Vector3::partial_max (a.0, b.0).into()
386}
387
388/// Given a 2D point and a 2D line, returns the nearest point on the line to
389/// the given point
390pub fn project_2d_point_on_line <S : Real>
391  (point : &Point2 <S>, line : &Line2 <S>) -> Point2 <S>
392{
393  let dot_dir = line.direction.dot (*point - line.base);
394  line.point (dot_dir)
395}
396
397/// Given a 3D point and a 3D line, returns the nearest point on the line to
398/// the given point
399pub fn project_3d_point_on_line <S : Real>
400  (point : &Point3 <S>, line : &Line3 <S>) -> Point3 <S>
401{
402  let dot_dir = line.direction.dot (*point - line.base);
403  line.point (dot_dir)
404}
405
406// private
407
408/// Smallest rectangle for given edge indices
409fn smallest_rect <S : Real> (i0 : u32, i1 : u32, hull : &Hull2 <S>)
410  -> Rect <S>
411{
412  // NOTE: the following algorithm assumes that the points in the hull are
413  // counter-clockwise sorted and that the hull does not contain any colinear
414  // points. The construction of the hull should ensure these conditions.
415  let points = hull.points();
416  let edge             = points[i1 as usize] - points[i0 as usize];
417  let perp             = vector2 (-edge.y, edge.x);
418  let edge_len_squared = edge.norm_squared();
419  let origin           = points[i1 as usize];
420  let mut support      = [Point2::<S>::origin(); 4];
421  let mut rect_index   = [i1, i1, i1, i1];
422  for (i, point) in points.iter().enumerate() {
423    let diff = point - origin;
424    let v    = point2 (edge.dot (diff), perp.dot (diff));
425    if v.0.x > support[1].0.x ||
426      v.0.x == support[1].0.x && v.0.y > support[1].0.y
427    {
428      rect_index[1] = i as u32;
429      support[1]    = v;
430    }
431    if v.0.y > support[2].0.y ||
432      v.0.y == support[2].0.y && v.0.x < support[2].0.x
433    {
434      rect_index[2] = i as u32;
435      support[2]    = v;
436    }
437    if v.0.x < support[3].0.x ||
438      v.0.x == support[3].0.x && v.0.y < support[3].0.y
439    {
440      rect_index[3] = i as u32;
441      support[3]    = v;
442    }
443  }
444  let scaled_width  = support[1].0.x - support[3].0.x;
445  let scaled_height = support[2].0.y;
446  let area = scaled_width * scaled_height / edge_len_squared;
447
448  Rect {
449    edge, perp, area, edge_len_squared,
450    indices: [rect_index[0], rect_index[1], rect_index[2], rect_index[3]]
451  }
452}
453
454////////////////////////////////////////////////////////////////////////////////
455//  impls                                                                     //
456////////////////////////////////////////////////////////////////////////////////
457
458impl <S : Ring> Interval <S> {
459  /// Construct a new AABB with given min and max points.
460  ///
461  /// Debug panic if points are not min/max:
462  ///
463  /// ```should_panic
464  /// # use math_utils::geometry::*;
465  /// let b = Interval::with_minmax (1.0, 0.0);  // panic!
466  /// ```
467  ///
468  /// Debug panic if points are identical:
469  ///
470  /// ```should_panic
471  /// # use math_utils::geometry::*;
472  /// let b = Interval::with_minmax (0.0, 0.0);  // panic!
473  /// ```
474  #[inline]
475  pub fn with_minmax (min : S, max : S) -> Self where S : std::fmt::Debug {
476    debug_assert_ne!(min, max);
477    debug_assert_eq!(min, S::min (min, max));
478    debug_assert_eq!(max, S::max (min, max));
479    Interval { min, max }
480  }
481  /// Construct a new AABB using the two given points to determine min/max.
482  ///
483  /// Panic if points are identical:
484  ///
485  /// ```should_panic
486  /// # use math_utils::geometry::*;
487  /// let b = Interval::from_points (0.0, 0.0);  // panic!
488  /// ```
489  #[inline]
490  pub fn from_points (a : S, b : S) -> Self where S : std::fmt::Debug {
491    debug_assert_ne!(a, b);
492    let min = S::min (a, b);
493    let max = S::max (a, b);
494    Interval { min, max }
495  }
496  /// Construct the minimum AABB containing the given set of points.
497  ///
498  /// Debug panic if fewer than 2 points are given:
499  ///
500  /// ```should_panic
501  /// # use math_utils::geometry::*;
502  /// let b = Interval::<f32>::containing (&[]);  // panic!
503  /// ```
504  ///
505  /// ```should_panic
506  /// # use math_utils::geometry::*;
507  /// let b = Interval::containing (&[0.0]);  // panic!
508  /// ```
509  #[inline]
510  pub fn containing (points : &[S]) -> Self where S : std::fmt::Debug {
511    debug_assert!(points.len() >= 2);
512    let mut min = S::min (points[0], points[1]);
513    let mut max = S::max (points[0], points[1]);
514    for point in points.iter().skip (2) {
515      min = S::min (min, *point);
516      max = S::max (max, *point);
517    }
518    Interval::with_minmax (min, max)
519  }
520  #[inline]
521  pub fn numcast <T> (self) -> Option <Interval <T>> where
522    S : num::ToPrimitive,
523    T : num::NumCast
524  {
525    Some (Interval {
526      min: T::from (self.min)?,
527      max: T::from (self.max)?
528    })
529  }
530  /// Create a new AABB that is the union of the two input AABBs
531  #[inline]
532  pub fn union (a : &Interval <S>, b : &Interval <S>) -> Self where
533    S : std::fmt::Debug
534  {
535    Interval::with_minmax (
536      S::min (a.min(), b.min()),
537      S::max (a.max(), b.max()))
538  }
539  #[inline]
540  pub fn min (&self) -> S {
541    self.min
542  }
543  #[inline]
544  pub fn max (&self) -> S {
545    self.max
546  }
547  #[inline]
548  pub fn length (&self) -> NonNegative <S> {
549    self.width()
550  }
551  #[inline]
552  pub fn width (&self) -> NonNegative <S> {
553    NonNegative::unchecked (self.max - self.min)
554  }
555  #[inline]
556  pub fn contains (&self, point : &S) -> bool {
557    self.min < *point && *point < self.max
558  }
559  /// Clamp a given point to the AABB interval.
560  ///
561  /// ```
562  /// # use math_utils::geometry::*;
563  /// let b = Interval::from_points (-1.0, 1.0);
564  /// assert_eq!(b.clamp (&-2.0), -1.0);
565  /// assert_eq!(b.clamp ( &2.0),  1.0);
566  /// assert_eq!(b.clamp ( &0.0),  0.0);
567  /// ```
568  pub fn clamp (&self, point : &S) -> S {
569    S::max (S::min (self.max, *point), self.min)
570  }
571  /// Generate a random point contained in the AABB
572  ///
573  /// ```
574  /// # use rand;
575  /// # use rand_xorshift;
576  /// # use math_utils::geometry::*;
577  /// # fn main () {
578  /// use rand_xorshift;
579  /// use rand::SeedableRng;
580  /// // random sequence will be the same each time this is run
581  /// let mut rng = rand_xorshift::XorShiftRng::seed_from_u64 (0);
582  /// let aabb = Interval::<f32>::with_minmax (-10.0, 10.0);
583  /// let point = aabb.rand_point (&mut rng);
584  /// assert!(aabb.contains (&point));
585  /// let point = aabb.rand_point (&mut rng);
586  /// assert!(aabb.contains (&point));
587  /// let point = aabb.rand_point (&mut rng);
588  /// assert!(aabb.contains (&point));
589  /// let point = aabb.rand_point (&mut rng);
590  /// assert!(aabb.contains (&point));
591  /// let point = aabb.rand_point (&mut rng);
592  /// assert!(aabb.contains (&point));
593  /// # }
594  /// ```
595  #[inline]
596  pub fn rand_point <R> (&self, rng : &mut R) -> S where
597    S : rand::distr::uniform::SampleUniform,
598    R : rand::Rng
599  {
600    rng.random_range (self.min..self.max)
601  }
602  #[inline]
603  pub fn intersects (&self, other : &Interval <S>) -> bool {
604    intersect::discrete_interval (self, other)
605  }
606  #[inline]
607  pub fn intersection (self, other : Interval <S>) -> Option <Interval <S>>
608    where S : std::fmt::Debug
609  {
610    intersect::continuous_interval (&self, &other)
611  }
612
613  /// Increase or decrease each endpoint by the given amount.
614  ///
615  /// Debug panic if the interval width is less than or equal to twice the given
616  /// amount:
617  /// ```should_panic
618  /// # use math_utils::geometry::*;
619  /// let x = Interval::<f32>::with_minmax (-1.0, 1.0);
620  /// let y = x.dilate (-1.0); // panic!
621  /// ```
622  pub fn dilate (mut self, amount : S) -> Self {
623    self.min -= amount;
624    self.max += amount;
625    debug_assert!(self.min < self.max);
626    self
627  }
628}
629
630impl <S> Primitive <S, S> for Interval <S> where
631  S : Field + VectorSpace <S>
632{
633  fn translate (mut self, displacement : S) -> Self {
634    self.min += displacement;
635    self.max += displacement;
636    self
637  }
638  fn scale (mut self, scale : NonZero <S>) -> Self {
639    let old_width = *self.width();
640    let new_width = old_width * *scale;
641    let half_difference = (new_width - old_width) / S::two();
642    self.min -= half_difference;
643    self.max += half_difference;
644    self
645  }
646}
647
648impl <S> std::ops::Add <S> for Interval <S> where
649  S    : Field + VectorSpace <S>,
650  Self : Primitive <S, S>
651{
652  type Output = Self;
653  fn add (self, displacement : S) -> Self {
654    self.translate (displacement)
655  }
656}
657
658impl <S> std::ops::Sub <S> for Interval <S> where
659  S    : Field + VectorSpace <S>,
660  Self : Primitive <S, S>
661{
662  type Output = Self;
663  fn sub (self, displacement : S) -> Self {
664    self.translate (-displacement)
665  }
666}
667
668impl_numcast_fields!(Aabb2, min, max);
669impl <S : Ring> Aabb2 <S> {
670  /// Construct a new AABB with given min and max points.
671  ///
672  /// Debug panic if points are not min/max:
673  ///
674  /// ```should_panic
675  /// # use math_utils::geometry::*;
676  /// let b = Aabb2::with_minmax ([1.0, 1.0].into(), [0.0, 0.0].into());  // panic!
677  /// ```
678  ///
679  /// Debug panic if points are identical:
680  ///
681  /// ```should_panic
682  /// # use math_utils::geometry::*;
683  /// let b = Aabb2::with_minmax ([0.0, 0.0].into(), [0.0, 0.0].into());  // panic!
684  /// ```
685  #[inline]
686  pub fn with_minmax (min : Point2 <S>, max : Point2 <S>) -> Self where
687    S : std::fmt::Debug
688  {
689    debug_assert_ne!(min, max);
690    debug_assert_eq!(min, point2_min (&min, &max));
691    debug_assert_eq!(max, point2_max (&min, &max));
692    Aabb2 { min, max }
693  }
694  /// Construct a new AABB using the two given points to determine min/max.
695  ///
696  /// Debug panic if points are identical:
697  ///
698  /// ```should_panic
699  /// # use math_utils::geometry::*;
700  /// let b = Aabb2::from_points ([0.0, 0.0].into(), [0.0, 0.0].into());  // panic!
701  /// ```
702  #[inline]
703  pub fn from_points (a : Point2 <S>, b : Point2 <S>) -> Self where
704    S : std::fmt::Debug
705  {
706    debug_assert_ne!(a, b);
707    let min = point2_min (&a, &b);
708    let max = point2_max (&a, &b);
709    Aabb2 { min, max }
710  }
711  /// Construct the minimum AABB containing the given set of points.
712  ///
713  /// Debug panic if fewer than 2 points are given:
714  ///
715  /// ```should_panic
716  /// # use math_utils::geometry::*;
717  /// let b = Aabb2::<f32>::containing (&[]);  // panic!
718  /// ```
719  ///
720  /// ```should_panic
721  /// # use math_utils::geometry::*;
722  /// let b = Aabb2::containing (&[[0.0, 0.0].into()]);  // panic!
723  /// ```
724  #[inline]
725  pub fn containing (points : &[Point2 <S>]) -> Self where S : std::fmt::Debug {
726    debug_assert!(points.len() >= 2);
727    let mut min = point2_min (&points[0], &points[1]);
728    let mut max = point2_max (&points[0], &points[1]);
729    for point in points.iter().skip (2) {
730      min = point2_min (&min, point);
731      max = point2_max (&max, point);
732    }
733    Aabb2::with_minmax (min, max)
734  }
735  /// Create a new AABB that is the union of the two input AABBs
736  #[inline]
737  pub fn union (a : &Aabb2 <S>, b : &Aabb2 <S>) -> Self where
738    S : std::fmt::Debug
739  {
740    Aabb2::with_minmax (
741      point2_min (a.min(), b.min()),
742      point2_max (a.max(), b.max())
743    )
744  }
745  #[inline]
746  pub fn min (&self) -> &Point2 <S> {
747    &self.min
748  }
749  #[inline]
750  pub fn max (&self) -> &Point2 <S> {
751    &self.max
752  }
753  #[inline]
754  pub fn center (&self) -> Point2 <S> where S : Field {
755    Point2 ((self.min.0 + self.max.0) / S::two())
756  }
757  #[inline]
758  pub fn dimensions (&self) -> Vector2 <S> {
759    self.max.0 - self.min.0
760  }
761  #[inline]
762  pub fn extents (&self) -> Vector2 <S> where S : Field {
763    self.dimensions() * S::half()
764  }
765  /// X dimension
766  #[inline]
767  pub fn width (&self) -> NonNegative <S> {
768    NonNegative::unchecked (self.max.0.x - self.min.0.x)
769  }
770  /// Y dimension
771  #[inline]
772  pub fn height (&self) -> NonNegative <S> {
773    NonNegative::unchecked (self.max.0.y - self.min.0.y)
774  }
775  #[inline]
776  pub fn area (&self) -> NonNegative <S> {
777    self.width() * self.height()
778  }
779  #[inline]
780  pub fn contains (&self, point : &Point2 <S>) -> bool {
781    self.min.0.x < point.0.x && point.0.x < self.max.0.x &&
782    self.min.0.y < point.0.y && point.0.y < self.max.0.y
783  }
784  /// Clamp a given point to the AABB.
785  ///
786  /// ```
787  /// # use math_utils::geometry::*;
788  /// let b = Aabb2::from_points ([-1.0, -1.0].into(), [1.0, 1.0].into());
789  /// assert_eq!(b.clamp (&[-2.0, 0.0].into()), [-1.0, 0.0].into());
790  /// assert_eq!(b.clamp (&[ 2.0, 2.0].into()), [ 1.0, 1.0].into());
791  /// assert_eq!(b.clamp (&[-0.5, 0.5].into()), [-0.5, 0.5].into());
792  /// ```
793  pub fn clamp (&self, point : &Point2 <S>) -> Point2 <S> {
794    [ S::max (S::min (self.max.0.x, point.0.x), self.min.0.x),
795      S::max (S::min (self.max.0.y, point.0.y), self.min.0.y),
796    ].into()
797  }
798  /// Generate a random point contained in the AABB
799  ///
800  /// ```
801  /// # use rand;
802  /// # use rand_xorshift;
803  /// # use math_utils::geometry::*;
804  /// # fn main () {
805  /// use rand_xorshift;
806  /// use rand::SeedableRng;
807  /// // random sequence will be the same each time this is run
808  /// let mut rng = rand_xorshift::XorShiftRng::seed_from_u64 (0);
809  /// let aabb = Aabb2::<f32>::with_minmax (
810  ///   [-10.0, -10.0].into(),
811  ///   [ 10.0,  10.0].into());
812  /// let point = aabb.rand_point (&mut rng);
813  /// assert!(aabb.contains (&point));
814  /// let point = aabb.rand_point (&mut rng);
815  /// assert!(aabb.contains (&point));
816  /// let point = aabb.rand_point (&mut rng);
817  /// assert!(aabb.contains (&point));
818  /// let point = aabb.rand_point (&mut rng);
819  /// assert!(aabb.contains (&point));
820  /// let point = aabb.rand_point (&mut rng);
821  /// assert!(aabb.contains (&point));
822  /// # }
823  /// ```
824  #[inline]
825  pub fn rand_point <R> (&self, rng : &mut R) -> Point2 <S> where
826    S : rand::distr::uniform::SampleUniform,
827    R : rand::Rng
828  {
829    let x_range = self.min.0.x..self.max.0.x;
830    let y_range = self.min.0.y..self.max.0.y;
831    let x = if x_range.is_empty() {
832      self.min.0.x
833    } else {
834      rng.random_range (x_range)
835    };
836    let y = if y_range.is_empty() {
837      self.min.0.y
838    } else {
839      rng.random_range (y_range)
840    };
841    [x, y].into()
842  }
843  #[inline]
844  pub fn intersects (&self, other : &Aabb2 <S>) -> bool {
845    intersect::discrete_aabb2_aabb2 (self, other)
846  }
847  #[inline]
848  pub fn intersection (self, other : Aabb2 <S>) -> Option <Aabb2 <S>> where
849    S : std::fmt::Debug
850  {
851    intersect::continuous_aabb2_aabb2 (&self, &other)
852  }
853
854  pub fn corner (&self, quadrant : Quadrant) -> Point2 <S> {
855    match quadrant {
856      // upper-right
857      Quadrant::PosPos => self.max,
858      // lower-right
859      Quadrant::PosNeg => [self.max.0.x, self.min.0.y].into(),
860      // upper-left
861      Quadrant::NegPos => [self.min.0.x, self.max.0.y].into(),
862      // lower-left
863      Quadrant::NegNeg => self.min
864    }
865  }
866
867  pub fn edge (&self, direction : SignedAxis2) -> Self where
868    S : std::fmt::Debug
869  {
870    let (min, max) = match direction {
871      SignedAxis2::PosX => (
872        self.corner (Quadrant::PosNeg),
873        self.corner (Quadrant::PosPos) ),
874      SignedAxis2::NegX => (
875        self.corner (Quadrant::NegNeg),
876        self.corner (Quadrant::NegPos) ),
877      SignedAxis2::PosY => (
878        self.corner (Quadrant::NegPos),
879        self.corner (Quadrant::PosPos) ),
880      SignedAxis2::NegY => (
881        self.corner (Quadrant::NegNeg),
882        self.corner (Quadrant::PosNeg) )
883    };
884    Aabb2::with_minmax (min, max)
885  }
886
887  pub fn extrude (&self, axis : SignedAxis2, amount : Positive <S>) -> Self {
888    let (min, max) = match axis {
889      SignedAxis2::PosX => (
890        self.min + Vector2::zero().with_x (*self.width()),
891        self.max + Vector2::zero().with_x (*amount)
892      ),
893      SignedAxis2::NegX => (
894        self.min - Vector2::zero().with_x (*amount),
895        self.max - Vector2::zero().with_x (*self.height())
896      ),
897      SignedAxis2::PosY => (
898        self.min + Vector2::zero().with_y (*self.height()),
899        self.max + Vector2::zero().with_y (*amount)
900      ),
901      SignedAxis2::NegY => (
902        self.min - Vector2::zero().with_y (*amount),
903        self.max - Vector2::zero().with_y (*self.height())
904      )
905    };
906    Aabb2 { min, max }
907  }
908
909  pub fn extend (&mut self, axis : SignedAxis2, amount : Positive <S>) {
910    match axis {
911      SignedAxis2::PosX => self.max.0.x += *amount,
912      SignedAxis2::NegX => self.min.0.x -= *amount,
913      SignedAxis2::PosY => self.max.0.y += *amount,
914      SignedAxis2::NegY => self.min.0.y -= *amount
915    }
916  }
917
918  pub fn with_z (self, z : Interval <S>) -> Aabb3 <S> where S : std::fmt::Debug {
919    Aabb3::with_minmax (
920      self.min.0.with_z (z.min()).into(),
921      self.max.0.with_z (z.max()).into()
922    )
923  }
924
925  /// Increase or decrease each dimension by the given amount.
926  ///
927  /// Debug panic if any dimension would become negative:
928  /// ```should_panic
929  /// # use math_utils::geometry::*;
930  /// let x = Aabb2::with_minmax ([0.0, 0.0].into(), [1.0, 2.0].into());
931  /// let y = x.dilate (-1.0); // panic!
932  /// ```
933  pub fn dilate (mut self, amount : S) -> Self where S : std::fmt::Debug {
934    self.min -= Vector2::broadcast (amount);
935    self.max += Vector2::broadcast (amount);
936    debug_assert_ne!(self.min, self.max);
937    debug_assert_eq!(self.min, point2_min (&self.min, &self.max));
938    debug_assert_eq!(self.max, point2_max (&self.min, &self.max));
939    self
940  }
941
942  /// Project *along* the given axis.
943  ///
944  /// For example, projecting *along* the X-axis projects *onto* the Y-axis.
945  pub fn project (self, axis : Axis2) -> Interval <S> where S : std::fmt::Debug {
946    let (min, max) = match axis {
947      Axis2::X => (self.min.0.y, self.max.0.y),
948      Axis2::Y => (self.min.0.x, self.max.0.x)
949    };
950    Interval::with_minmax (min, max)
951  }
952}
953
954impl <S : Field> Primitive <S, Vector2 <S>> for Aabb2 <S> {
955  fn translate (mut self, displacement : Vector2 <S>) -> Self {
956    self.min.0 += displacement;
957    self.max.0 += displacement;
958    self
959  }
960  fn scale (mut self, scale : NonZero <S>) -> Self {
961    let center = self.center();
962    self = self.translate (-center.0);
963    self.min.0 *= *scale;
964    self.max.0 *= *scale;
965    self.translate (center.0)
966  }
967}
968
969impl <S> std::ops::Add <Vector2 <S>> for Aabb2 <S> where
970  S    : Field,
971  Self : Primitive <S, Vector2 <S>>
972{
973  type Output = Self;
974  fn add (self, displacement : Vector2 <S>) -> Self {
975    self.translate (displacement)
976  }
977}
978
979impl <S> std::ops::Sub <Vector2 <S>> for Aabb2 <S> where
980  S    : Field,
981  Self : Primitive <S, Vector2 <S>>
982{
983  type Output = Self;
984  fn sub (self, displacement : Vector2 <S>) -> Self {
985    self.translate (-displacement)
986  }
987}
988
989impl_numcast_fields!(Aabb3, min, max);
990impl <S : Ring> Aabb3 <S> {
991  /// Construct a new AABB with given min and max points.
992  ///
993  /// Debug panic if points are not min/max:
994  ///
995  /// ```should_panic
996  /// # use math_utils::geometry::*;
997  /// let b = Aabb3::with_minmax ([1.0, 1.0, 1.0].into(), [0.0, 0.0, 0.0].into());  // panic!
998  /// ```
999  ///
1000  /// Debug panic if points are identical:
1001  ///
1002  /// ```should_panic
1003  /// # use math_utils::geometry::*;
1004  /// let b = Aabb3::with_minmax ([0.0, 0.0, 0.0].into(), [0.0, 0.0, 0.0].into());  // panic!
1005  /// ```
1006  #[inline]
1007  pub fn with_minmax (min : Point3 <S>, max : Point3 <S>) -> Self where
1008    S : std::fmt::Debug
1009  {
1010    debug_assert_ne!(min, max);
1011    debug_assert_eq!(min, point3_min (&min, &max));
1012    debug_assert_eq!(max, point3_max (&min, &max));
1013    Aabb3 { min, max }
1014  }
1015  /// Construct a new AABB using the two given points to determine min/max.
1016  ///
1017  /// Panic if points are identical:
1018  ///
1019  /// ```should_panic
1020  /// # use math_utils::geometry::*;
1021  /// let b = Aabb3::from_points (
1022  ///   [0.0, 0.0, 0.0].into(),
1023  ///   [0.0, 0.0, 0.0].into());  // panic!
1024  /// ```
1025  #[inline]
1026  pub fn from_points (a : Point3 <S>, b : Point3 <S>) -> Self where
1027    S : std::fmt::Debug
1028  {
1029    debug_assert_ne!(a, b);
1030    let min = point3_min (&a, &b);
1031    let max = point3_max (&a, &b);
1032    Aabb3 { min, max }
1033  }
1034  /// Construct the minimum AABB containing the given set of points.
1035  ///
1036  /// Debug panic if fewer than 2 points are given:
1037  ///
1038  /// ```should_panic
1039  /// # use math_utils::geometry::*;
1040  /// let b = Aabb3::<f32>::containing (&[]);  // panic!
1041  /// ```
1042  ///
1043  /// ```should_panic
1044  /// # use math_utils::geometry::*;
1045  /// let b = Aabb3::containing (&[[0.0, 0.0, 0.0].into()]);  // panic!
1046  /// ```
1047  #[inline]
1048  pub fn containing (points : &[Point3 <S>]) -> Self where S : std::fmt::Debug {
1049    debug_assert!(points.len() >= 2);
1050    let mut min = point3_min (&points[0], &points[1]);
1051    let mut max = point3_max (&points[0], &points[1]);
1052    for point in points.iter().skip (2) {
1053      min = point3_min (&min, point);
1054      max = point3_max (&max, point);
1055    }
1056    Aabb3::with_minmax (min, max)
1057  }
1058  /// Create a new AABB that is the union of the two input AABBs
1059  #[inline]
1060  pub fn union (a : &Aabb3 <S>, b : &Aabb3 <S>) -> Self where
1061    S : std::fmt::Debug
1062  {
1063    Aabb3::with_minmax (
1064      point3_min (a.min(), b.min()),
1065      point3_max (a.max(), b.max())
1066    )
1067  }
1068  #[inline]
1069  pub fn min (&self) -> &Point3 <S> {
1070    &self.min
1071  }
1072  #[inline]
1073  pub fn max (&self) -> &Point3 <S> {
1074    &self.max
1075  }
1076  #[inline]
1077  pub fn center (&self) -> Point3 <S> where S : Field {
1078    Point3 ((self.min.0 + self.max.0) / S::two())
1079  }
1080  #[inline]
1081  pub fn dimensions (&self) -> Vector3 <S> {
1082    self.max.0 - self.min.0
1083  }
1084  #[inline]
1085  pub fn extents (&self) -> Vector3 <S> where S : Field {
1086    self.dimensions() * S::half()
1087  }
1088  /// X dimension
1089  #[inline]
1090  pub fn width (&self) -> NonNegative <S> {
1091    NonNegative::unchecked (self.max.0.x - self.min.0.x)
1092  }
1093  /// Y dimension
1094  #[inline]
1095  pub fn depth (&self) -> NonNegative <S> {
1096    NonNegative::unchecked (self.max.0.y - self.min.0.y)
1097  }
1098  /// Z dimension
1099  #[inline]
1100  pub fn height (&self) -> NonNegative <S> {
1101    NonNegative::unchecked (self.max.0.z - self.min.0.z)
1102  }
1103  #[inline]
1104  pub fn volume (&self) -> NonNegative <S> {
1105    self.width() * self.depth() * self.height()
1106  }
1107  #[inline]
1108  pub fn contains (&self, point : &Point3 <S>) -> bool {
1109    self.min.0.x < point.0.x && point.0.x < self.max.0.x &&
1110    self.min.0.y < point.0.y && point.0.y < self.max.0.y &&
1111    self.min.0.z < point.0.z && point.0.z < self.max.0.z
1112  }
1113  /// Clamp a given point to the AABB.
1114  ///
1115  /// ```
1116  /// # use math_utils::geometry::*;
1117  /// let b = Aabb3::with_minmax ([-1.0, -1.0, -1.0].into(), [1.0, 1.0, 1.0].into());
1118  /// assert_eq!(b.clamp (&[-2.0, 0.0, 0.0].into()), [-1.0, 0.0, 0.0].into());
1119  /// assert_eq!(b.clamp (&[ 2.0, 2.0, 0.0].into()), [ 1.0, 1.0, 0.0].into());
1120  /// assert_eq!(b.clamp (&[-1.0, 2.0, 3.0].into()), [-1.0, 1.0, 1.0].into());
1121  /// assert_eq!(b.clamp (&[-0.5, 0.5, 0.0].into()), [-0.5, 0.5, 0.0].into());
1122  /// ```
1123  pub fn clamp (&self, point : &Point3 <S>) -> Point3 <S> {
1124    [ S::max (S::min (self.max.0.x, point.0.x), self.min.0.x),
1125      S::max (S::min (self.max.0.y, point.0.y), self.min.0.y),
1126      S::max (S::min (self.max.0.z, point.0.z), self.min.0.z)
1127    ].into()
1128  }
1129  /// Generate a random point contained in the AABB
1130  ///
1131  /// ```
1132  /// # use rand;
1133  /// # use rand_xorshift;
1134  /// # use math_utils::geometry::*;
1135  /// # fn main () {
1136  /// use rand_xorshift;
1137  /// use rand::SeedableRng;
1138  /// // random sequence will be the same each time this is run
1139  /// let mut rng = rand_xorshift::XorShiftRng::seed_from_u64 (0);
1140  /// let aabb = Aabb3::<f32>::with_minmax (
1141  ///   [-10.0, -10.0, -10.0].into(),
1142  ///   [ 10.0,  10.0,  10.0].into());
1143  /// let point = aabb.rand_point (&mut rng);
1144  /// assert!(aabb.contains (&point));
1145  /// let point = aabb.rand_point (&mut rng);
1146  /// assert!(aabb.contains (&point));
1147  /// let point = aabb.rand_point (&mut rng);
1148  /// assert!(aabb.contains (&point));
1149  /// let point = aabb.rand_point (&mut rng);
1150  /// assert!(aabb.contains (&point));
1151  /// let point = aabb.rand_point (&mut rng);
1152  /// assert!(aabb.contains (&point));
1153  /// # }
1154  /// ```
1155  #[inline]
1156  pub fn rand_point <R> (&self, rng : &mut R) -> Point3 <S> where
1157    S : rand::distr::uniform::SampleUniform,
1158    R : rand::Rng
1159  {
1160    let x_range = self.min.0.x..self.max.0.x;
1161    let y_range = self.min.0.y..self.max.0.y;
1162    let z_range = self.min.0.z..self.max.0.z;
1163    let x = if x_range.is_empty() {
1164      self.min.0.x
1165    } else {
1166      rng.random_range (x_range)
1167    };
1168    let y = if y_range.is_empty() {
1169      self.min.0.y
1170    } else {
1171      rng.random_range (y_range)
1172    };
1173    let z = if z_range.is_empty() {
1174      self.min.0.z
1175    } else {
1176      rng.random_range (z_range)
1177    };
1178    [x, y, z].into()
1179  }
1180  #[inline]
1181  pub fn intersects (&self, other : &Aabb3 <S>) -> bool {
1182    intersect::discrete_aabb3_aabb3 (self, other)
1183  }
1184  #[inline]
1185  pub fn intersection (self, other : Aabb3 <S>) -> Option <Aabb3 <S>> where
1186    S : std::fmt::Debug
1187  {
1188    intersect::continuous_aabb3_aabb3 (&self, &other)
1189  }
1190
1191  pub fn corner (&self, octant : Octant) -> Point3 <S> {
1192    match octant {
1193      Octant::PosPosPos => self.max,
1194      Octant::NegPosPos => [self.min.0.x, self.max.0.y, self.max.0.z].into(),
1195      Octant::PosNegPos => [self.max.0.x, self.min.0.y, self.max.0.z].into(),
1196      Octant::NegNegPos => [self.min.0.x, self.min.0.y, self.max.0.z].into(),
1197      Octant::PosPosNeg => [self.max.0.x, self.max.0.y, self.min.0.z].into(),
1198      Octant::NegPosNeg => [self.min.0.x, self.max.0.y, self.min.0.z].into(),
1199      Octant::PosNegNeg => [self.max.0.x, self.min.0.y, self.min.0.z].into(),
1200      Octant::NegNegNeg => self.min
1201    }
1202  }
1203
1204  pub fn face (&self, direction : SignedAxis3) -> Self where
1205    S : std::fmt::Debug
1206  {
1207    let (min, max) = match direction {
1208      SignedAxis3::PosX => (
1209        self.corner (Octant::PosNegNeg),
1210        self.corner (Octant::PosPosPos) ),
1211      SignedAxis3::NegX => (
1212        self.corner (Octant::NegNegNeg),
1213        self.corner (Octant::NegPosPos) ),
1214      SignedAxis3::PosY => (
1215        self.corner (Octant::NegPosNeg),
1216        self.corner (Octant::PosPosPos) ),
1217      SignedAxis3::NegY => (
1218        self.corner (Octant::NegNegNeg),
1219        self.corner (Octant::PosNegPos) ),
1220      SignedAxis3::PosZ => (
1221        self.corner (Octant::NegNegPos),
1222        self.corner (Octant::PosPosPos) ),
1223      SignedAxis3::NegZ => (
1224        self.corner (Octant::NegNegNeg),
1225        self.corner (Octant::PosPosNeg) )
1226    };
1227    Aabb3::with_minmax (min, max)
1228  }
1229
1230  pub fn extrude (&self, axis : SignedAxis3, amount : Positive <S>) -> Self {
1231    let (min, max) = match axis {
1232      SignedAxis3::PosX => (
1233        self.min + Vector3::zero().with_x (*self.width()),
1234        self.max + Vector3::zero().with_x (*amount)
1235      ),
1236      SignedAxis3::NegX => (
1237        self.min - Vector3::zero().with_x (*amount),
1238        self.max - Vector3::zero().with_x (*self.width())
1239      ),
1240      SignedAxis3::PosY => (
1241        self.min + Vector3::zero().with_y (*self.depth()),
1242        self.max + Vector3::zero().with_y (*amount)
1243      ),
1244      SignedAxis3::NegY => (
1245        self.min - Vector3::zero().with_y (*amount),
1246        self.max - Vector3::zero().with_y (*self.depth())
1247      ),
1248      SignedAxis3::PosZ => (
1249        self.min + Vector3::zero().with_z (*self.height()),
1250        self.max + Vector3::zero().with_z (*amount)
1251      ),
1252      SignedAxis3::NegZ => (
1253        self.min - Vector3::zero().with_z (*amount),
1254        self.max - Vector3::zero().with_z (*self.height())
1255      )
1256    };
1257    Aabb3 { min, max }
1258  }
1259
1260  pub fn extend (&mut self, axis : SignedAxis3, amount : Positive <S>) {
1261    match axis {
1262      SignedAxis3::PosX => self.max.0.x += *amount,
1263      SignedAxis3::NegX => self.min.0.x -= *amount,
1264      SignedAxis3::PosY => self.max.0.y += *amount,
1265      SignedAxis3::NegY => self.min.0.y -= *amount,
1266      SignedAxis3::PosZ => self.max.0.z += *amount,
1267      SignedAxis3::NegZ => self.min.0.z -= *amount
1268    }
1269  }
1270
1271  /// Increase or decrease each dimension by the given amount.
1272  ///
1273  /// Debug panic if any dimension would become negative:
1274  /// ```should_panic
1275  /// # use math_utils::geometry::*;
1276  /// let x = Aabb3::with_minmax ([0.0, 0.0, 0.0].into(), [1.0, 2.0, 3.0].into());
1277  /// let y = x.dilate (-1.0); // panic!
1278  /// ```
1279  pub fn dilate (mut self, amount : S) -> Self where S : std::fmt::Debug {
1280    self.min -= Vector3::broadcast (amount);
1281    self.max += Vector3::broadcast (amount);
1282    debug_assert_ne!(self.min, self.max);
1283    debug_assert_eq!(self.min, point3_min (&self.min, &self.max));
1284    debug_assert_eq!(self.max, point3_max (&self.min, &self.max));
1285    self
1286  }
1287
1288  /// Project *along* the given axis.
1289  ///
1290  /// For example, projecting *along* the X-axis projects *onto* the YZ-plane.
1291  pub fn project (self, axis : Axis3) -> Aabb2 <S> where S : std::fmt::Debug {
1292    let (min, max) = match axis {
1293      Axis3::X => ([self.min.0.y, self.min.0.z], [self.max.0.y, self.max.0.z]),
1294      Axis3::Y => ([self.min.0.x, self.min.0.z], [self.max.0.x, self.max.0.z]),
1295      Axis3::Z => ([self.min.0.x, self.min.0.y], [self.max.0.x, self.max.0.y]),
1296    };
1297    Aabb2::with_minmax (min.into(), max.into())
1298  }
1299}
1300
1301impl <S : Field> Primitive <S, Vector3 <S>> for Aabb3 <S> {
1302  fn translate (mut self, displacement : Vector3 <S>) -> Self {
1303    self.min.0 += displacement;
1304    self.max.0 += displacement;
1305    self
1306  }
1307  fn scale (mut self, scale : NonZero <S>) -> Self {
1308    let center = self.center();
1309    self = self.translate (-center.0);
1310    self.min.0 *= *scale;
1311    self.max.0 *= *scale;
1312    self.translate (center.0)
1313  }
1314}
1315
1316impl <S> std::ops::Add <Vector3 <S>> for Aabb3 <S> where
1317  S    : Field,
1318  Self : Primitive <S, Vector3 <S>>
1319{
1320  type Output = Self;
1321  fn add (self, displacement : Vector3 <S>) -> Self {
1322    self.translate (displacement)
1323  }
1324}
1325
1326impl <S> std::ops::Sub <Vector3 <S>> for Aabb3 <S> where
1327  S    : Field,
1328  Self : Primitive <S, Vector3 <S>>
1329{
1330  type Output = Self;
1331  fn sub (self, displacement : Vector3 <S>) -> Self {
1332    self.translate (-displacement)
1333  }
1334}
1335
1336impl_numcast_fields!(Obb2, center, half_extent_x, half_extent_y, orientation);
1337impl <S : Ring> Obb2 <S> {
1338  /// Construct a new OBB
1339  #[inline]
1340  pub fn new (
1341    center        : Point2 <S>,
1342    half_extent_x : Positive <S>,
1343    half_extent_y : Positive <S>,
1344    orientation   : AngleWrapped <S>
1345  ) -> Self {
1346    Obb2 { center, half_extent_x, half_extent_y, orientation }
1347  }
1348  /// Construct the minimum OBB containing the convex hull of a given set of
1349  /// points using rotating calipers algorithm. To construct the OBB containing
1350  /// a convex hull directoy use the `Obb2::containing` method.
1351  ///
1352  /// Returns None if fewer than 3 points are given or if points span an empty
1353  /// set:
1354  ///
1355  /// ```
1356  /// # use math_utils::*;
1357  /// # use math_utils::geometry::*;
1358  /// assert!(Obb2::containing_points (
1359  ///   &[[0.0, 1.0], [1.0, 0.0]].map (Point2::from)
1360  /// ).is_none());
1361  ///
1362  /// assert!(Obb2::containing_points (
1363  ///   &[[1.0, 1.0], [0.0, 0.0], [-1.0, -1.0]].map (Point2::from)
1364  /// ).is_none());
1365  /// ```
1366  pub fn containing_points (points : &[Point2 <S>]) -> Option <Self> where
1367    S : Real
1368  {
1369    if points.len() < 3 {
1370      return None
1371    }
1372    let hull = Hull2::from_points (points)?;
1373    Self::containing (&hull)
1374  }
1375  /// Construct the minimum OBB containing the given convex hull of points using
1376  /// rotating calipers algorithm.
1377  ///
1378  /// Returns None if if points span an empty set:
1379  ///
1380  /// ```
1381  /// # use math_utils::*;
1382  /// # use math_utils::geometry::*;
1383  /// let hull = Hull2::from_points (
1384  ///   &[[1.0, 1.0], [0.0, 0.0], [-1.0, -1.0]].map (Point2::from)).unwrap();
1385  /// assert!(Obb2::containing (&hull).is_none());
1386  /// ```
1387  // code follows GeometricTools:
1388  // <https://github.com/davideberly/GeometricTools/blob/f78dd0b65bc3a0a4723586de6991dd2c339b08ad/GTE/Mathematics/MinimumAreaBox2.h>
1389  pub fn containing (hull : &Hull2 <S>) -> Option <Self> where S : Real {
1390    let points = hull.points();
1391    if points.len() < 3 {
1392      return None
1393    }
1394    let mut visited  = vec![false; points.len()];
1395    let mut min_rect = smallest_rect (points.len() as u32 - 1, 0, hull);
1396    visited[min_rect.indices[0] as usize] = true;
1397    // rotating calipers
1398    let mut rect = min_rect.clone();
1399    for _ in 0..points.len() {
1400      let mut angles  = [(S::zero(), u32::MAX); 4];
1401      let mut nangles = u32::MAX;
1402      if !compute_angles (hull, &rect, &mut angles, &mut nangles) {
1403        break
1404      }
1405      let sorted = sort_angles (&mut angles, nangles);
1406      if !update_support (
1407        &angles, nangles, &sorted, hull, &mut visited, &mut rect
1408      ) {
1409        break
1410      }
1411      if rect.area < min_rect.area {
1412        min_rect = rect.clone();
1413      }
1414    }
1415    // construct obb from min rect
1416    let sum = [
1417      points[min_rect.indices[1] as usize].0 +
1418        points[min_rect.indices[3] as usize].0,
1419      points[min_rect.indices[2] as usize].0 +
1420        points[min_rect.indices[0] as usize].0
1421    ];
1422    let diff = [
1423      points[min_rect.indices[1] as usize].0 -
1424        points[min_rect.indices[3] as usize].0,
1425      points[min_rect.indices[2] as usize].0 -
1426        points[min_rect.indices[0] as usize].0
1427    ];
1428    let center = Point2::from ((min_rect.edge * min_rect.edge.dot (sum[0]) +
1429      min_rect.perp * min_rect.perp.dot (sum[1]))
1430      * S::half() / min_rect.edge_len_squared);
1431    let half_extent_x = Positive::unchecked (S::sqrt (
1432      (S::half() * min_rect.edge.dot (diff[0])).squared()
1433      / min_rect.edge_len_squared));
1434    let half_extent_y = Positive::unchecked (S::sqrt (
1435      (S::half() * min_rect.perp.dot (diff[1])).squared()
1436      / min_rect.edge_len_squared));
1437    let orientation =
1438      AngleWrapped::wrap (Rad (min_rect.edge.y.atan2 (min_rect.edge.x)));
1439    let obb = Obb2 { center, half_extent_x, half_extent_y, orientation };
1440    fn compute_angles <S : Real> (
1441      hull    : &Hull2 <S>,
1442      rect    : &Rect <S>,
1443      angles  : &mut [(S, u32); 4],
1444      nangles : &mut u32
1445    ) -> bool {
1446      *nangles = 0;
1447      let points = hull.points();
1448      let mut k0 = 3;
1449      let mut k1 = 0;
1450      while k1 < 4 {
1451        if rect.indices[k0] != rect.indices[k1] {
1452          let d = {
1453            let mut d = if k0 & 1 != 0 {
1454              rect.perp
1455            } else {
1456              rect.edge
1457            };
1458            if k0 & 2 != 0 {
1459              d = -d;
1460            }
1461            d
1462          };
1463          let j0 = rect.indices[k0];
1464          let mut j1 = j0 + 1;
1465          if j1 == points.len() as u32 {
1466            j1 = 0;
1467          }
1468          let e = points[j1 as usize] - points[j0 as usize];
1469          let dp = d.dot ([e.y, -e.x].into());
1470          let e_lensq = e.norm_squared();
1471          let sin_theta_squared = (dp * dp) / e_lensq;
1472          angles[*nangles as usize] = (sin_theta_squared, k0 as u32);
1473          *nangles += 1;
1474        }
1475        k0 = k1;
1476        k1 += 1;
1477      }
1478      *nangles > 0
1479    }
1480    fn sort_angles <S : PartialOrd> (angles : &mut [(S, u32); 4], nangles : u32)
1481      -> [u32; 4]
1482    {
1483      let mut sorted = [0u32, 1, 2, 3];
1484      match nangles {
1485        0 | 1 => {}
1486        2 => if angles[sorted[0] as usize].0 > angles[sorted[1] as usize].0 {
1487          sorted.swap (0, 1)
1488        }
1489        3 => {
1490          if angles[sorted[0] as usize].0 > angles[sorted[1] as usize].0 {
1491            sorted.swap (0, 1)
1492          }
1493          if angles[sorted[0] as usize].0 > angles[sorted[2] as usize].0 {
1494            sorted.swap (0, 2)
1495          }
1496          if angles[sorted[1] as usize].0 > angles[sorted[2] as usize].0 {
1497            sorted.swap (1, 2)
1498          }
1499        }
1500        4 => {
1501          if angles[sorted[0] as usize].0 > angles[sorted[1] as usize].0 {
1502            sorted.swap (0, 1)
1503          }
1504          if angles[sorted[2] as usize].0 > angles[sorted[3] as usize].0 {
1505            sorted.swap (2, 3)
1506          }
1507          if angles[sorted[0] as usize].0 > angles[sorted[2] as usize].0 {
1508            sorted.swap (0, 2)
1509          }
1510          if angles[sorted[1] as usize].0 > angles[sorted[3] as usize].0 {
1511            sorted.swap (1, 3)
1512          }
1513          if angles[sorted[1] as usize].0 > angles[sorted[2] as usize].0 {
1514            sorted.swap (1, 2)
1515          }
1516        }
1517        _ => unreachable!()
1518      }
1519      sorted
1520    }
1521    fn update_support <S : Real> (
1522      angles  : &[(S, u32); 4],
1523      nangles : u32,
1524      sorted  : &[u32; 4],
1525      hull    : &Hull2 <S>,
1526      visited : &mut [bool],
1527      rect    : &mut Rect <S>
1528    ) -> bool {
1529      let points = hull.points();
1530      let min_angle = angles[sorted[0] as usize];
1531      for k in 0..nangles as usize {
1532        let (angle, index) = angles[sorted[k] as usize];
1533        if angle == min_angle.0 {
1534          rect.indices[index as usize] += 1;
1535          if rect.indices[index as usize] == points.len() as u32 {
1536            rect.indices[index as usize] = 0;
1537          }
1538        }
1539      }
1540      let bottom = rect.indices[min_angle.1 as usize];
1541      if visited[bottom as usize] {
1542        return false
1543      }
1544      visited[bottom as usize] = true;
1545      let mut next_index = [u32::MAX; 4];
1546      for k in 0u32..4 {
1547        next_index[k as usize] = rect.indices[((min_angle.1 + k) % 4) as usize];
1548      }
1549      rect.indices = next_index;
1550      let j1 = rect.indices[0] as usize;
1551      let j0 = if j1 == 0 {
1552        points.len() - 1
1553      } else {
1554        j1 - 1
1555      };
1556      rect.edge = points[j1] - points[j0];
1557      rect.perp = vector2 (-rect.edge.y, rect.edge.x);
1558      let edge_len_squared = rect.edge.norm_squared();
1559      let diff = [
1560        points[rect.indices[1] as usize] - points[rect.indices[3] as usize],
1561        points[rect.indices[2] as usize] - points[rect.indices[0] as usize]
1562      ];
1563      rect.area = rect.edge.dot (diff[0]) * rect.perp.dot (diff[1]) /
1564        edge_len_squared;
1565      true
1566    }
1567    Some (obb)
1568  }
1569  #[inline]
1570  pub fn center (&self) -> Point2 <S> {
1571    self.center
1572  }
1573  #[inline]
1574  pub fn half_extent_x (&self) -> Positive <S> {
1575    self.half_extent_x
1576  }
1577  #[inline]
1578  pub fn half_extent_y (&self) -> Positive <S> {
1579    self.half_extent_y
1580  }
1581  #[inline]
1582  pub fn orientation (&self) -> AngleWrapped <S> {
1583    self.orientation
1584  }
1585  pub fn dimensions (&self) -> Vector2 <S> where S : Field {
1586    self.extents() * S::two()
1587  }
1588  pub fn extents (&self) -> Vector2 <S> {
1589    [*self.half_extent_x, *self.half_extent_y].into()
1590  }
1591  /// X dimension
1592  #[inline]
1593  pub fn width (&self) -> Positive <S> where S : Field {
1594    self.half_extent_x * Positive::unchecked (S::two())
1595  }
1596  /// Y dimension
1597  #[inline]
1598  pub fn height (&self) -> Positive <S> where S : Field {
1599    self.half_extent_y * Positive::unchecked (S::two())
1600  }
1601  #[inline]
1602  pub fn area (&self) -> Positive <S> where S : Field {
1603    self.width() * self.height()
1604  }
1605  #[inline]
1606  pub fn corners (&self) -> [Point2 <S>; 4] where
1607    S : Real + num::real::Real + MaybeSerDes
1608  {
1609    let extents = self.extents();
1610    let mut points = [Point2::origin(); 4];
1611    for i in 0..2 {
1612      for j in 0..2 {
1613        let sign_x = if i % 2 == 0 { -S::one() } else { S::one() };
1614        let sign_y = if j % 2 == 0 { -S::one() } else { S::one() };
1615        points[i * 2 + j] = [extents.x * sign_x, extents.y * sign_y].into();
1616      }
1617    }
1618    let rotation = Rotation2::from_angle (self.orientation.angle());
1619    points.map (|point| rotation.rotate (point) + self.center.0)
1620  }
1621
1622  /// Check if point is contained by the bounding box.
1623  ///
1624  /// ```
1625  /// # use math_utils::*;
1626  /// # use math_utils::geometry::*;
1627  /// let x = Obb2::new (
1628  ///   [3.0, 6.0].into(),
1629  ///   Positive::unchecked (2.0),
1630  ///   Positive::unchecked (0.75),
1631  ///   AngleWrapped::wrap (Rad (std::f32::consts::FRAC_PI_4)));
1632  /// assert!(x.contains (&[2.0f32, 5.0].into()));
1633  /// assert!(!x.contains (&[3.0f32, 2.0].into()));
1634  /// ```
1635  pub fn contains (&self, point : &Point2 <S>) -> bool where
1636    S : Real + num::real::Real + MaybeSerDes
1637  {
1638    // re-center
1639    let p = point - self.center().0;
1640    // reverse rotation
1641    let p = Rotation2::from_angle (-self.orientation().angle()).rotate (p);
1642    p.0.x.abs() < *self.half_extent_x && p.0.y.abs() < *self.half_extent_y
1643  }
1644
1645  /// Increase or decrease each dimension by the given amount.
1646  ///
1647  /// Panics if any dimension would become zero or negative:
1648  ///
1649  /// ```should_panic
1650  /// # use math_utils::*;
1651  /// # use math_utils::geometry::*;
1652  /// # use math_utils::num_traits::Zero;
1653  /// let x = Obb2::new (
1654  ///   Point2::origin(),
1655  ///   Positive::unchecked (0.5),
1656  ///   Positive::unchecked (0.5),
1657  ///   AngleWrapped::zero()
1658  /// );
1659  /// let y = x.dilate (-1.0); // panic!
1660  /// ```
1661  pub fn dilate (mut self, amount : S) -> Self where S : Field {
1662    self.half_extent_x = Positive::noisy (*self.half_extent_x + amount);
1663    self.half_extent_y = Positive::noisy (*self.half_extent_y + amount);
1664    self
1665  }
1666}
1667impl <S : Real> From <Aabb2 <S>> for Obb2 <S> {
1668  /// Panics if AABB has a zero dimension:
1669  ///
1670  /// ```should_panic
1671  /// # use math_utils::*;
1672  /// # use math_utils::geometry::*;
1673  /// # use math_utils::num_traits::Zero;
1674  /// let x = Aabb2::with_minmax (
1675  ///   [0.0, 0.0].into(),
1676  ///   [0.0, 1.0].into()
1677  /// );
1678  /// let y = Obb2::from (x); // panic!
1679  /// ```
1680  fn from (aabb : Aabb2 <S>) -> Self {
1681    use num::Zero;
1682    let center        = aabb.center();
1683    let half_extent_x = Positive::noisy (*aabb.width()  / S::two());
1684    let half_extent_y = Positive::noisy (*aabb.height() / S::two());
1685    let orientation   = AngleWrapped::zero();
1686    Obb2::new (center, half_extent_x, half_extent_y, orientation)
1687  }
1688}
1689
1690impl_numcast_fields!(Obb3, center, half_extent_x, half_extent_y, half_extent_z,
1691  orientation);
1692impl <S : Ring> Obb3 <S> {
1693  /// Construct a new OBB
1694  #[inline]
1695  pub fn new (
1696    center        : Point3 <S>,
1697    half_extent_x : Positive <S>,
1698    half_extent_y : Positive <S>,
1699    half_extent_z : Positive <S>,
1700    orientation   : Angles3 <S>
1701  ) -> Self {
1702    Obb3 { center, half_extent_x, half_extent_y, half_extent_z, orientation }
1703  }
1704  /// Construct the minimum OBB containing the convex hull of a given set of
1705  /// points
1706  pub fn containing_points (_points : &[Point3 <S>]) -> Option <Self> where
1707    S : Real
1708  {
1709    unimplemented!("TODO")
1710  }
1711  /// Construct the minimum OBB containing the given convex hull of points
1712  pub fn containing (_hull : &Hull3 <S>) -> Option <Self> where S : Real {
1713    unimplemented!("TODO")
1714  }
1715  #[inline]
1716  pub fn center (&self) -> Point3 <S> {
1717    self.center
1718  }
1719  #[inline]
1720  pub fn half_extent_x (&self) -> Positive <S> {
1721    self.half_extent_x
1722  }
1723  #[inline]
1724  pub fn half_extent_y (&self) -> Positive <S> {
1725    self.half_extent_y
1726  }
1727  #[inline]
1728  pub fn half_extent_z (&self) -> Positive <S> {
1729    self.half_extent_z
1730  }
1731  #[inline]
1732  pub fn orientation (&self) -> Angles3 <S> {
1733    self.orientation
1734  }
1735  pub fn dimensions (&self) -> Vector3 <S> where S : Field {
1736    self.extents() * S::two()
1737  }
1738  pub fn extents (&self) -> Vector3 <S> where S : Field {
1739    [*self.half_extent_x, *self.half_extent_y, *self.half_extent_z].into()
1740  }
1741  /// X dimension
1742  #[inline]
1743  pub fn width (&self) -> Positive <S> where S : Field {
1744    self.half_extent_x * Positive::unchecked (S::two())
1745  }
1746  /// Y dimension
1747  #[inline]
1748  pub fn depth (&self) -> Positive <S> where S : Field {
1749    self.half_extent_y * Positive::unchecked (S::two())
1750  }
1751  /// Z dimension
1752  #[inline]
1753  pub fn height (&self) -> Positive <S> where S : Field {
1754    self.half_extent_z * Positive::unchecked (S::two())
1755  }
1756  #[inline]
1757  pub fn volume (&self) -> Positive <S> where S : Field {
1758    self.width() * self.depth() * self.height()
1759  }
1760  #[inline]
1761  pub fn corners (&self) -> [Point3 <S>; 8] where
1762    S : Real + num::real::Real + MaybeSerDes
1763  {
1764    let extents = self.extents();
1765    let mut points = [Point3::origin(); 8];
1766    for i in 0..2 {
1767      for j in 0..2 {
1768        for k in 0..2 {
1769          let sign_x = if i % 2 == 0 { -S::one() } else { S::one() };
1770          let sign_y = if j % 2 == 0 { -S::one() } else { S::one() };
1771          let sign_z = if k % 2 == 0 { -S::one() } else { S::one() };
1772          points[i * 4 + j * 2 + k] =
1773            [extents.x * sign_x, extents.y * sign_y, extents.z * sign_z].into();
1774        }
1775      }
1776    }
1777    let rotation = Rotation3::from (self.orientation);
1778    points.map (|point| rotation.rotate (point) + self.center.0)
1779  }
1780
1781  /// Check if point is contained by the bounding box
1782  pub fn contains (&self, _point : &Point2 <S>) -> bool where
1783    S : Real + num::real::Real + MaybeSerDes
1784  {
1785    unimplemented!("TODO")
1786  }
1787
1788  /// Increase or decrease each dimension by the given amount.
1789  ///
1790  /// Panics if any dimension would become zero or negative:
1791  ///
1792  /// ```should_panic
1793  /// # use math_utils::*;
1794  /// # use math_utils::geometry::*;
1795  /// let x = Obb3::new (
1796  ///   Point3::origin(),
1797  ///   Positive::unchecked (0.5),
1798  ///   Positive::unchecked (0.5),
1799  ///   Positive::unchecked (0.5),
1800  ///   Angles3::zero()
1801  /// );
1802  /// let y = x.dilate (-1.0); // panic!
1803  /// ```
1804  pub fn dilate (mut self, amount : S) -> Self where S : Field {
1805    self.half_extent_x = Positive::noisy (*self.half_extent_x + amount);
1806    self.half_extent_y = Positive::noisy (*self.half_extent_y + amount);
1807    self.half_extent_z = Positive::noisy (*self.half_extent_z + amount);
1808    self
1809  }
1810}
1811impl <S : Real> From <Aabb3 <S>> for Obb3 <S> {
1812  /// Panics if AABB has a zero dimension:
1813  ///
1814  /// ```should_panic
1815  /// # use math_utils::*;
1816  /// # use math_utils::geometry::*;
1817  /// # use math_utils::num_traits::Zero;
1818  /// let x = Aabb3::with_minmax (
1819  ///   [0.0, 0.0, 0.0].into(),
1820  ///   [0.0, 1.0, 1.0].into()
1821  /// );
1822  /// let y = Obb3::from (x); // panic!
1823  /// ```
1824  fn from (aabb : Aabb3 <S>) -> Self {
1825    let center = aabb.center();
1826    let [half_extent_x, half_extent_y, half_extent_z] =
1827      aabb.extents().into_array().map (Positive::noisy);
1828    let orientation = Angles3::zero();
1829    Obb3::new (center, half_extent_x, half_extent_y, half_extent_z, orientation)
1830  }
1831}
1832
1833impl_numcast_fields!(Capsule3, center, half_height, radius);
1834impl <S : Ring> Capsule3 <S> {
1835  pub fn aabb3 (&self) -> Aabb3 <S> where S : Real + std::fmt::Debug {
1836    use shape::Aabb;
1837    let shape_aabb = shape::Capsule {
1838      radius:      self.radius,
1839      half_height: self.half_height
1840    }.aabb();
1841    let center_vec = self.center.0;
1842    let min        = *shape_aabb.min() + center_vec;
1843    let max        = *shape_aabb.max() + center_vec;
1844    Aabb3::with_minmax (min, max)
1845  }
1846  /// Return the (upper sphere, cylinder, lower sphere) making up this capsule
1847  pub fn decompose (&self)
1848    -> (Sphere3 <S>, Option <Cylinder3 <S>>, Sphere3 <S>)
1849  {
1850    let cylinder = if *self.half_height > S::zero() {
1851      Some (Cylinder3 {
1852        center:      self.center,
1853        radius:      self.radius,
1854        half_height: Positive::unchecked (*self.half_height)
1855      })
1856    } else {
1857      None
1858    };
1859    let upper_sphere = Sphere3 {
1860      center: self.center + (Vector3::unit_z() * *self.half_height),
1861      radius: self.radius
1862    };
1863    let lower_sphere = Sphere3 {
1864      center: self.center - (Vector3::unit_z() * *self.half_height),
1865      radius: self.radius
1866    };
1867    (upper_sphere, cylinder, lower_sphere)
1868  }
1869}
1870
1871// TODO: impl Primitive for Capsule3
1872
1873impl_numcast_fields!(Cylinder3, center, half_height, radius);
1874impl <S : Ring> Cylinder3 <S> {
1875  pub fn aabb3 (&self) -> Aabb3 <S> where S : Real + std::fmt::Debug {
1876    use shape::Aabb;
1877    let shape_aabb = shape::Cylinder::unchecked (*self.radius, *self.half_height)
1878      .aabb();
1879    let center_vec = self.center.0;
1880    let min        = *shape_aabb.min() + center_vec;
1881    let max        = *shape_aabb.max() + center_vec;
1882    Aabb3::with_minmax (min, max)
1883  }
1884}
1885
1886// TODO: impl Primitive for Cylinder3
1887
1888impl_numcast_fields!(Line2, base, direction);
1889impl <S : Real> Line2 <S> {
1890  /// Construct a new 2D line
1891  #[inline]
1892  pub fn new (base : Point2 <S>, direction : Unit2 <S>) -> Self {
1893    Line2 { base, direction }
1894  }
1895  #[inline]
1896  pub fn point (&self, t : S) -> Point2 <S> {
1897    self.base + *self.direction * t
1898  }
1899  #[inline]
1900  #[allow(clippy::type_complexity)]
1901  pub fn intersect_aabb (&self, aabb : &Aabb2 <S>)
1902    -> Option <((S, Point2 <S>), (S, Point2 <S>))>
1903  where S : std::fmt::Debug {
1904    intersect::continuous_line2_aabb2 (self, aabb)
1905  }
1906  #[inline]
1907  #[allow(clippy::type_complexity)]
1908  pub fn intersect_sphere (&self, sphere : &Sphere2 <S>)
1909    -> Option <((S, Point2 <S>), (S, Point2 <S>))>
1910  {
1911    intersect::continuous_line2_sphere2 (self, sphere)
1912  }
1913}
1914impl <S : Real> Default for Line2 <S> {
1915  fn default() -> Self {
1916    Line2 {
1917      base:      Point2::origin(),
1918      direction: Unit2::axis_y()
1919    }
1920  }
1921}
1922
1923// TODO: impl Primitive for Line2
1924
1925impl_numcast_fields!(Line3, base, direction);
1926impl <S : Real> Line3 <S> {
1927  /// Consructs a new 3D line with given base point and direction vector
1928  #[inline]
1929  pub fn new (base : Point3 <S>, direction : Unit3 <S>) -> Self {
1930    Line3 { base, direction }
1931  }
1932  #[inline]
1933  pub fn point (&self, t : S) -> Point3 <S> {
1934    self.base + *self.direction * t
1935  }
1936  #[inline]
1937  pub fn intersect_plane (&self, plane : &Plane3 <S>)
1938    -> Option <(S, Point3 <S>)>
1939  where
1940    S : approx::RelativeEq
1941  {
1942    intersect::continuous_line3_plane3 (self, plane)
1943  }
1944  #[inline]
1945  #[allow(clippy::type_complexity)]
1946  pub fn intersect_aabb (&self, aabb : &Aabb3 <S>)
1947    -> Option <((S, Point3 <S>), (S, Point3 <S>))>
1948  where S : num::Float + approx::RelativeEq <Epsilon=S> + std::fmt::Debug {
1949    intersect::continuous_line3_aabb3 (self, aabb)
1950  }
1951  #[inline]
1952  #[allow(clippy::type_complexity)]
1953  pub fn intersect_sphere (&self, sphere : &Sphere3 <S>)
1954    -> Option <((S, Point3 <S>), (S, Point3 <S>))>
1955  {
1956    intersect::continuous_line3_sphere3 (self, sphere)
1957  }
1958}
1959impl <S : Real> Default for Line3 <S> {
1960  fn default() -> Self {
1961    Line3 {
1962      base:      Point3::origin(),
1963      direction: Unit3::axis_z()
1964    }
1965  }
1966}
1967
1968// TODO: impl Primitive for Line3
1969
1970impl_numcast_fields!(Plane3, base, normal);
1971impl <S : Real> Plane3 <S> {
1972  /// Consructs a new 3D plane with given base point and normal vector
1973  #[inline]
1974  pub fn new (base : Point3 <S>, normal : Unit3 <S>) -> Self {
1975    Plane3 { base, normal }
1976  }
1977}
1978impl <S : Real> Default for Plane3 <S> {
1979  fn default() -> Self {
1980    Plane3 {
1981      base:   Point3::origin(),
1982      normal: Unit3::axis_z()
1983    }
1984  }
1985}
1986
1987// TODO: impl Primitive for Plane3
1988
1989impl_numcast_fields!(Sphere2, center, radius);
1990impl <S : Ring> Sphere2 <S> {
1991  /// Unit circle
1992  #[inline]
1993  pub fn unit() -> Self where S : Field {
1994    Sphere2 {
1995      center: Point2::origin(),
1996      radius: num::One::one()
1997    }
1998  }
1999  /// Discrete intersection test with another sphere
2000  #[inline]
2001  pub fn intersects (&self, other : &Self) -> bool {
2002    intersect::discrete_sphere2_sphere2 (self, other)
2003  }
2004}
2005
2006// TODO: impl Primitive for Sphere2
2007
2008impl_numcast_fields!(Sphere3, center, radius);
2009impl <S : Ring> Sphere3 <S> {
2010  /// Unit sphere
2011  #[inline]
2012  pub fn unit() -> Self where S : Field {
2013    Sphere3 {
2014      center: Point3::origin(),
2015      radius: num::One::one()
2016    }
2017  }
2018  /// Discrete intersection test with another sphere
2019  #[inline]
2020  pub fn intersects (&self, other : &Self) -> bool {
2021    intersect::discrete_sphere3_sphere3 (self, other)
2022  }
2023}
2024
2025// TODO: impl Primitive for Sphere3
2026
2027////////////////////////////////////////////////////////////////////////////////
2028//  tests                                                                     //
2029////////////////////////////////////////////////////////////////////////////////
2030
2031#[cfg(test)]
2032mod tests {
2033  use approx::{assert_relative_eq, assert_ulps_eq};
2034  use super::*;
2035
2036  #[test]
2037  fn test_project_2d_point_on_line() {
2038    use Unit2;
2039    let point : Point2 <f64> = [2.0, 2.0].into();
2040    let line = Line2::<f64>::new ([0.0, 0.0].into(), Unit2::axis_x());
2041    assert_eq!(project_2d_point_on_line (&point, &line), [2.0, 0.0].into());
2042    let line = Line2::<f64>::new ([0.0, 0.0].into(), Unit2::axis_y());
2043    assert_eq!(project_2d_point_on_line (&point, &line), [0.0, 2.0].into());
2044    let point : Point2 <f64> = [0.0, 1.0].into();
2045    let line = Line2::<f64>::new (
2046      [0.0, -1.0].into(), Unit2::normalize ([1.0, 1.0].into()));
2047    assert_relative_eq!(
2048      project_2d_point_on_line (&point, &line), [1.0, 0.0].into());
2049    // the answer should be the same for lines with equivalent definitions
2050    let point : Point2 <f64> = [1.0, 3.0].into();
2051    let line_a = Line2::<f64>::new (
2052      [0.0, -1.0].into(), Unit2::normalize ([2.0, 1.0].into()));
2053    let line_b = Line2::<f64>::new (
2054      [2.0, 0.0].into(),  Unit2::normalize ([2.0, 1.0].into()));
2055    assert_relative_eq!(
2056      project_2d_point_on_line (&point, &line_a),
2057      project_2d_point_on_line (&point, &line_b));
2058    let line_a = Line2::<f64>::new (
2059      [0.0, 0.0].into(),   Unit2::normalize ([1.0, 1.0].into()));
2060    let line_b = Line2::<f64>::new (
2061      [-2.0, -2.0].into(), Unit2::normalize ([1.0, 1.0].into()));
2062    assert_ulps_eq!(
2063      project_2d_point_on_line (&point, &line_a),
2064      project_2d_point_on_line (&point, &line_b)
2065    );
2066    assert_relative_eq!(
2067      project_2d_point_on_line (&point, &line_a),
2068      [2.0, 2.0].into());
2069  }
2070
2071  #[test]
2072  fn test_project_3d_point_on_line() {
2073    use Unit3;
2074    // all the tests from 2d projection with 0.0 for the Z component
2075    let point : Point3 <f64> = [2.0, 2.0, 0.0].into();
2076    let line = Line3::<f64>::new ([0.0, 0.0, 0.0].into(), Unit3::axis_x());
2077    assert_eq!(project_3d_point_on_line (&point, &line), [2.0, 0.0, 0.0].into());
2078    let line = Line3::<f64>::new ([0.0, 0.0, 0.0].into(), Unit3::axis_y());
2079    assert_eq!(project_3d_point_on_line (&point, &line), [0.0, 2.0, 0.0].into());
2080    let point : Point3 <f64> = [0.0, 1.0, 0.0].into();
2081    let line = Line3::<f64>::new (
2082      [0.0, -1.0, 0.0].into(), Unit3::normalize ([1.0, 1.0, 0.0].into()));
2083    assert_relative_eq!(
2084      project_3d_point_on_line (&point, &line), [1.0, 0.0, 0.0].into());
2085    // the answer should be the same for lines with equivalent definitions
2086    let point : Point3 <f64> = [1.0, 3.0, 0.0].into();
2087    let line_a = Line3::<f64>::new (
2088      [0.0, -1.0, 0.0].into(), Unit3::normalize ([2.0, 1.0, 0.0].into()));
2089    let line_b = Line3::<f64>::new (
2090      [2.0, 0.0, 0.0].into(),  Unit3::normalize ([2.0, 1.0, 0.0].into()));
2091    assert_relative_eq!(
2092      project_3d_point_on_line (&point, &line_a),
2093      project_3d_point_on_line (&point, &line_b));
2094    let line_a = Line3::<f64>::new (
2095      [0.0, 0.0, 0.0].into(), Unit3::normalize ([1.0, 1.0, 0.0].into()));
2096    let line_b = Line3::<f64>::new (
2097      [2.0, 2.0, 0.0].into(), Unit3::normalize ([1.0, 1.0, 0.0].into()));
2098    assert_relative_eq!(
2099      project_3d_point_on_line (&point, &line_a),
2100      project_3d_point_on_line (&point, &line_b));
2101    assert_relative_eq!(
2102      project_3d_point_on_line (&point, &line_a),
2103      [2.0, 2.0, 0.0].into());
2104    // more tests
2105    let point : Point3 <f64> = [0.0, 0.0, 2.0].into();
2106    let line_a = Line3::<f64>::new (
2107      [-4.0, -4.0, -4.0].into(), Unit3::normalize ([1.0, 1.0, 1.0].into()));
2108    let line_b = Line3::<f64>::new (
2109      [4.0, 4.0, 4.0].into(),    Unit3::normalize ([1.0, 1.0, 1.0].into()));
2110    assert_relative_eq!(
2111      project_3d_point_on_line (&point, &line_a),
2112      project_3d_point_on_line (&point, &line_b),
2113      epsilon = 0.00000000000001
2114    );
2115    assert_relative_eq!(
2116      project_3d_point_on_line (&point, &line_a),
2117      [2.0/3.0, 2.0/3.0, 2.0/3.0].into(),
2118      epsilon = 0.00000000000001
2119    );
2120  }
2121
2122  #[test]
2123  fn test_smallest_rect() {
2124    // points are in counter-clockwise order
2125    let points : Vec <Point2 <f32>> = [
2126      [-3.0, -3.0],
2127      [ 3.0, -3.0],
2128      [ 3.0,  3.0],
2129      [ 0.0,  6.0],
2130      [-3.0,  3.0]
2131    ].map (Point2::from).into_iter().collect();
2132    let hull = Hull2::from_points (points.as_slice()).unwrap();
2133    assert_eq!(hull.points(), points);
2134    let rect = smallest_rect (0, 1, &hull);
2135    assert_eq!(rect.area, 54.0);
2136    // points are in counter-clockwise order
2137    let points : Vec <Point2 <f32>> = [
2138      [-1.0, -4.0],
2139      [ 2.0,  2.0],
2140      [-4.0, -1.0]
2141    ].map (Point2::from).into_iter().collect();
2142    let hull = Hull2::from_points (points.as_slice()).unwrap();
2143    assert_eq!(hull.points(), points);
2144    let rect = smallest_rect (0, 1, &hull);
2145    assert_eq!(rect.area, 27.0);
2146  }
2147
2148  #[test]
2149  fn test_obb2() {
2150    let points : Vec <Point2 <f32>> = [
2151      [-1.0, -4.0],
2152      [ 2.0,  2.0],
2153      [-4.0, -1.0]
2154    ].map (Point2::from).into_iter().collect();
2155    let obb = Obb2::containing_points (&points).unwrap();
2156    let corners = obb.corners();
2157    assert_eq!(obb.center, [-0.25, -0.25].into());
2158    approx::assert_relative_eq!(point2 (-4.0, -1.0), corners[0], max_relative = 0.00001);
2159    approx::assert_relative_eq!(point2 ( 0.5,  3.5), corners[1], max_relative = 0.00001);
2160    approx::assert_relative_eq!(point2 (-1.0, -4.0), corners[2], max_relative = 0.00001);
2161    approx::assert_relative_eq!(point2 ( 3.5,  0.5), corners[3], max_relative = 0.00001);
2162  }
2163} // end tests
math_utils/geometry/primitive/mod.rs

math_utils/geometry/primitive/
mod.rs
