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

math_utils/geometry/primitive/
mod.rs
