rgeometry/data/
point.rs

1use array_init::{array_init, try_array_init};
2use float::FloatCore;
3use num_bigint::BigInt;
4use num_rational::BigRational;
5use num_traits::*;
6use ordered_float::OrderedFloat;
7use ordered_float::{FloatIsNan, NotNan};
8use rand::distributions::{Distribution, Standard};
9use rand::Rng;
10use std::cmp::Ordering;
11use std::convert::TryFrom;
12use std::fmt::Debug;
13use std::iter::Sum;
14use std::ops::Deref;
15use std::ops::Index;
16
17use super::{Direction, Vector};
18use crate::{Orientation, PolygonScalar, TotalOrd};
19
20#[allow(clippy::derived_hash_with_manual_eq)]
21#[derive(Debug, Clone, Copy, Hash)]
22#[repr(transparent)] // Required for correctness!
23pub struct Point<T, const N: usize = 2> {
24  pub array: [T; N],
25}
26
27impl<T: TotalOrd, const N: usize> TotalOrd for Point<T, N> {
28  fn total_cmp(&self, other: &Self) -> Ordering {
29    for i in 0..N {
30      match self.array[i].total_cmp(&other.array[i]) {
31        Ordering::Equal => continue,
32        other => return other,
33      }
34    }
35    Ordering::Equal
36  }
37}
38
39impl<T: TotalOrd, const N: usize> PartialEq for Point<T, N> {
40  fn eq(&self, other: &Self) -> bool {
41    self.total_cmp(other).is_eq()
42  }
43}
44impl<T: TotalOrd, const N: usize> Eq for Point<T, N> {}
45
46impl<T: TotalOrd, const N: usize> PartialOrd for Point<T, N> {
47  fn partial_cmp(&self, other: &Self) -> Option<Ordering> {
48    Some(std::cmp::Ord::cmp(self, other))
49  }
50}
51
52impl<T: TotalOrd, const N: usize> Ord for Point<T, N> {
53  fn cmp(&self, other: &Self) -> Ordering {
54    self.total_cmp(other)
55  }
56}
57
58#[derive(Debug, Clone, Copy)]
59pub struct PointSoS<'a, T, const N: usize = 2> {
60  pub index: u32,
61  pub point: &'a Point<T, N>,
62}
63
64impl<T: TotalOrd, const N: usize> TotalOrd for PointSoS<'_, T, N> {
65  fn total_cmp(&self, other: &Self) -> Ordering {
66    (self.index, self.point).total_cmp(&(other.index, other.point))
67  }
68}
69
70impl<T: TotalOrd, const N: usize> PartialEq for PointSoS<'_, T, N> {
71  fn eq(&self, other: &Self) -> bool {
72    self.total_cmp(other).is_eq()
73  }
74}
75impl<T: TotalOrd, const N: usize> Eq for PointSoS<'_, T, N> {}
76
77impl<T: TotalOrd, const N: usize> PartialOrd for PointSoS<'_, T, N> {
78  fn partial_cmp(&self, other: &Self) -> Option<Ordering> {
79    Some(std::cmp::Ord::cmp(self, other))
80  }
81}
82
83impl<T: TotalOrd, const N: usize> Ord for PointSoS<'_, T, N> {
84  fn cmp(&self, other: &Self) -> Ordering {
85    self.total_cmp(other)
86  }
87}
88
89// Random sampling.
90impl<T, const N: usize> Distribution<Point<T, N>> for Standard
91where
92  Standard: Distribution<T>,
93{
94  // FIXME: Unify with code for Vector.
95  fn sample<R: Rng + ?Sized>(&self, rng: &mut R) -> Point<T, N> {
96    Point {
97      array: array_init(|_| rng.gen()),
98    }
99  }
100}
101
102// Methods on N-dimensional points.
103impl<T, const N: usize> Point<T, N> {
104  pub const fn new(array: [T; N]) -> Point<T, N> {
105    Point { array }
106  }
107
108  /// # Panics
109  ///
110  /// Panics if any of the inputs are NaN.
111  pub fn new_nn(array: [T; N]) -> Point<NotNan<T>, N>
112  where
113    T: FloatCore,
114  {
115    Point::new(array_init(|i| NotNan::new(array[i]).unwrap()))
116  }
117
118  pub fn as_vec(&self) -> &Vector<T, N> {
119    self.into()
120  }
121
122  // Warning: May cause arithmetic overflow.
123  // See `cmp_distance_to` for a safe way to compare distances in 2D.
124  pub fn squared_euclidean_distance<F>(&self, rhs: &Point<T, N>) -> F
125  where
126    T: Clone + Into<F>,
127    F: NumOps + Clone + Sum,
128  {
129    self
130      .array
131      .iter()
132      .zip(rhs.array.iter())
133      .map(|(a, b)| {
134        let diff = a.clone().into() - b.clone().into();
135        diff.clone() * diff
136      })
137      .sum()
138  }
139
140  // Similar to num_traits::identities::Zero but doesn't require an Add impl.
141  pub fn zero() -> Self
142  where
143    T: Sum,
144  {
145    Point {
146      array: array_init(|_| std::iter::empty().sum()),
147    }
148  }
149
150  pub fn map<U, F>(&self, f: F) -> Point<U, N>
151  where
152    T: Clone,
153    F: Fn(T) -> U,
154  {
155    Point {
156      array: array_init(|i| f(self.array[i].clone())),
157    }
158  }
159
160  pub fn cast<U>(&self) -> Point<U, N>
161  where
162    T: Clone + Into<U>,
163  {
164    Point {
165      array: array_init(|i| self.array[i].clone().into()),
166    }
167  }
168
169  pub fn to_float(&self) -> Point<OrderedFloat<f64>, N>
170  where
171    T: Clone + Into<f64>,
172  {
173    self.map(|v| OrderedFloat(v.into()))
174  }
175}
176
177impl<T, const N: usize> Index<usize> for Point<T, N> {
178  type Output = T;
179  fn index(&self, key: usize) -> &T {
180    self.array.index(key)
181  }
182}
183
184impl<const N: usize> TryFrom<Point<f64, N>> for Point<NotNan<f64>, N> {
185  type Error = FloatIsNan;
186  fn try_from(point: Point<f64, N>) -> Result<Point<NotNan<f64>, N>, FloatIsNan> {
187    Ok(Point {
188      array: try_array_init(|i| NotNan::try_from(point.array[i]))?,
189    })
190  }
191}
192
193impl<T> From<(T, T)> for Point<T, 2> {
194  fn from(point: (T, T)) -> Point<T, 2> {
195    Point {
196      array: [point.0, point.1],
197    }
198  }
199}
200
201impl From<Point<i64, 2>> for Point<BigInt, 2> {
202  fn from(point: Point<i64, 2>) -> Point<BigInt, 2> {
203    Point {
204      array: [point.array[0].into(), point.array[1].into()],
205    }
206  }
207}
208
209impl<const N: usize> From<&Point<BigRational, N>> for Point<f64, N> {
210  fn from(point: &Point<BigRational, N>) -> Point<f64, N> {
211    Point {
212      array: array_init(|i| point.array[i].to_f64().unwrap()),
213    }
214  }
215}
216
217impl<const N: usize> From<Point<BigRational, N>> for Point<f64, N> {
218  fn from(point: Point<BigRational, N>) -> Point<f64, N> {
219    Point {
220      array: array_init(|i| point.array[i].to_f64().unwrap()),
221    }
222  }
223}
224
225impl<const N: usize> From<&Point<f64, N>> for Point<BigRational, N> {
226  fn from(point: &Point<f64, N>) -> Point<BigRational, N> {
227    Point {
228      array: array_init(|i| BigRational::from_f64(point.array[i]).unwrap()),
229    }
230  }
231}
232
233impl<const N: usize> From<Point<f64, N>> for Point<BigRational, N> {
234  fn from(point: Point<f64, N>) -> Point<BigRational, N> {
235    Point {
236      array: array_init(|i| BigRational::from_f64(point.array[i]).unwrap()),
237    }
238  }
239}
240
241impl<T, const N: usize> From<Vector<T, N>> for Point<T, N> {
242  fn from(vector: Vector<T, N>) -> Point<T, N> {
243    Point { array: vector.0 }
244  }
245}
246
247// Methods on two-dimensional points.
248impl<T: PolygonScalar> Point<T> {
249  pub fn cmp_distance_to(&self, p: &Point<T, 2>, q: &Point<T, 2>) -> Ordering {
250    T::cmp_dist(self, p, q)
251  }
252
253  /// Determine the direction you have to turn if you walk from `p1`
254  /// to `p2` to `p3`.
255  ///
256  /// For fixed-precision types (i8,i16,i32,i64,etc), this function is
257  /// guaranteed to work for any input and never cause any arithmetic overflows.
258  ///
259  /// # Examples
260  ///
261  /// ```rust
262  /// # use rgeometry::data::Point;
263  /// # use rgeometry::Orientation;
264  /// let p1 = Point::new([ 0, 0 ]);
265  /// let p2 = Point::new([ 0, 1 ]); // One unit above p1.
266  /// // (0,0) -> (0,1) -> (0,2) == Orientation::CoLinear
267  /// assert!(Point::orient(&p1, &p2, &Point::new([ 0, 2 ])).is_colinear());
268  /// // (0,0) -> (0,1) -> (-1,2) == Orientation::CounterClockWise
269  /// assert!(Point::orient(&p1, &p2, &Point::new([ -1, 2 ])).is_ccw());
270  /// // (0,0) -> (0,1) -> (1,2) == Orientation::ClockWise
271  /// assert!(Point::orient(&p1, &p2, &Point::new([ 1, 2 ])).is_cw());
272  /// ```
273  ///
274  pub fn orient(p1: &Point<T, 2>, p2: &Point<T, 2>, p3: &Point<T, 2>) -> Orientation {
275    Orientation::new(p1, p2, p3)
276  }
277
278  pub fn orient_along_direction(
279    p1: &Point<T, 2>,
280    direction: Direction<'_, T, 2>,
281    p2: &Point<T, 2>,
282  ) -> Orientation {
283    match direction {
284      Direction::Vector(v) => Orientation::along_vector(p1, v, p2),
285      Direction::Through(p) => Orientation::new(p1, p, p2),
286    }
287  }
288
289  pub fn orient_along_vector(
290    p1: &Point<T, 2>,
291    vector: &Vector<T, 2>,
292    p2: &Point<T, 2>,
293  ) -> Orientation {
294    Orientation::along_vector(p1, vector, p2)
295  }
296
297  pub fn orient_along_perp_vector(
298    p1: &Point<T, 2>,
299    vector: &Vector<T, 2>,
300    p2: &Point<T, 2>,
301  ) -> Orientation {
302    Orientation::along_perp_vector(p1, vector, p2)
303  }
304
305  pub fn all_colinear(pts: &[Point<T>]) -> bool {
306    if pts.len() < 3 {
307      return true;
308    }
309    pts
310      .iter()
311      .all(|pt| Point::orient(&pts[0], &pts[1], pt).is_colinear())
312  }
313
314  /// Docs?
315  pub fn ccw_cmp_around(&self, p: &Point<T, 2>, q: &Point<T, 2>) -> Ordering {
316    self.ccw_cmp_around_with(&Vector([T::from_constant(1), T::from_constant(0)]), p, q)
317  }
318
319  pub fn ccw_cmp_around_with(
320    &self,
321    z: &Vector<T, 2>,
322    p: &Point<T, 2>,
323    q: &Point<T, 2>,
324  ) -> Ordering {
325    Orientation::ccw_cmp_around_with(z, self, p, q)
326  }
327}
328
329// FIXME: Use a macro
330impl<T> Point<T, 1> {
331  pub fn x_coord(&self) -> &T {
332    &self.array[0]
333  }
334}
335impl<T> Point<T, 2> {
336  pub fn x_coord(&self) -> &T {
337    &self.array[0]
338  }
339  pub fn y_coord(&self) -> &T {
340    &self.array[1]
341  }
342}
343impl<T> Point<T, 3> {
344  pub fn x_coord(&self) -> &T {
345    &self.array[0]
346  }
347  pub fn y_coord(&self) -> &T {
348    &self.array[1]
349  }
350  pub fn z_coord(&self) -> &T {
351    &self.array[2]
352  }
353}
354
355impl<T, const N: usize> Deref for Point<T, N> {
356  type Target = [T; N];
357  fn deref(&self) -> &[T; N] {
358    &self.array
359  }
360}
361
362mod add;
363mod sub;
364
365#[cfg(test)]
366#[cfg(not(tarpaulin_include))]
367pub mod tests {
368  use super::*;
369  use crate::testing::*;
370  use crate::Orientation::*;
371
372  use proptest::prelude::*;
373  use test_strategy::proptest;
374
375  #[proptest]
376  fn cmp_dist_i8_fuzz(pt1: Point<i8, 2>, pt2: Point<i8, 2>, pt3: Point<i8, 2>) {
377    let pt1_big: Point<BigInt, 2> = pt1.cast();
378    let pt2_big: Point<BigInt, 2> = pt2.cast();
379    let pt3_big: Point<BigInt, 2> = pt3.cast();
380    prop_assert_eq!(
381      pt1.cmp_distance_to(&pt2, &pt3),
382      pt1_big.cmp_distance_to(&pt2_big, &pt3_big)
383    );
384  }
385
386  #[proptest]
387  fn squared_euclidean_distance_fuzz(
388    #[strategy(any_nn::<2>())] pt1: Point<NotNan<f64>, 2>,
389    #[strategy(any_nn::<2>())] pt2: Point<NotNan<f64>, 2>,
390  ) {
391    let _out: f64 = pt1.squared_euclidean_distance(&pt2);
392  }
393
394  #[proptest]
395  fn squared_euclidean_distance_i8_fuzz(pt1: Point<i8, 2>, pt2: Point<i8, 2>) {
396    let _out: i64 = pt1.squared_euclidean_distance(&pt2);
397  }
398
399  #[proptest]
400  fn cmp_around_fuzz_nn(
401    #[strategy(any_nn::<2>())] pt1: Point<NotNan<f64>, 2>,
402    #[strategy(any_nn::<2>())] pt2: Point<NotNan<f64>, 2>,
403    #[strategy(any_nn::<2>())] pt3: Point<NotNan<f64>, 2>,
404  ) {
405    let _ = pt1.ccw_cmp_around(&pt2, &pt3);
406  }
407
408  #[proptest]
409  fn cmp_around_fuzz_i8(pt1: Point<i8, 2>, pt2: Point<i8, 2>, pt3: Point<i8, 2>) {
410    let _ = pt1.ccw_cmp_around(&pt2, &pt3);
411  }
412
413  #[proptest]
414  fn bigint_colinear(
415    #[strategy(any_r::<2>())] pt1: Point<BigInt, 2>,
416    #[strategy(any_r::<2>())] pt2: Point<BigInt, 2>,
417  ) {
418    let diff = &pt2 - &pt1;
419    let pt3 = &pt2 + &diff;
420    prop_assert!(Point::orient(&pt1, &pt2, &pt3).is_colinear())
421  }
422
423  #[proptest]
424  fn bigint_not_colinear(
425    #[strategy(any_r::<2>())] pt1: Point<BigInt, 2>,
426    #[strategy(any_r::<2>())] pt2: Point<BigInt, 2>,
427  ) {
428    let diff = &pt2 - &pt1;
429    let pt3 = &pt2 + &diff + &Vector([BigInt::from(1), BigInt::from(1)]);
430    prop_assert!(!Point::orient(&pt1, &pt2, &pt3).is_colinear())
431  }
432
433  // i8 and BigInt uses different algorithms for computing orientation but the results
434  // should be identical.
435  #[proptest]
436  fn orient_bigint_i8_prop(pt1: Point<i8, 2>, pt2: Point<i8, 2>, pt3: Point<i8, 2>) {
437    let pt1_big: Point<BigInt, 2> = pt1.cast();
438    let pt2_big: Point<BigInt, 2> = pt2.cast();
439    let pt3_big: Point<BigInt, 2> = pt3.cast();
440
441    prop_assert_eq!(
442      Point::orient(&pt1, &pt2, &pt3),
443      Point::orient(&pt1_big, &pt2_big, &pt3_big)
444    )
445  }
446
447  // i8 and f64 uses different algorithms for computing orientation but the results
448  // should be identical.
449  #[proptest]
450  fn orient_f64_i8_prop(pt1: Point<i8, 2>, pt2: Point<i8, 2>, pt3: Point<i8, 2>) {
451    let pt1_big: Point<OrderedFloat<f64>, 2> = pt1.to_float();
452    let pt2_big: Point<OrderedFloat<f64>, 2> = pt2.to_float();
453    let pt3_big: Point<OrderedFloat<f64>, 2> = pt3.to_float();
454
455    prop_assert_eq!(
456      Point::orient(&pt1, &pt2, &pt3),
457      Point::orient(&pt1_big, &pt2_big, &pt3_big)
458    )
459  }
460
461  #[proptest]
462  fn orientation_reverse(pt1: Point<i64, 2>, pt2: Point<i64, 2>, pt3: Point<i64, 2>) {
463    let abc = Point::orient(&pt1, &pt2, &pt3);
464    let cba = Point::orient(&pt3, &pt2, &pt1);
465    prop_assert_eq!(abc, cba.reverse())
466  }
467
468  #[test]
469  fn test_turns() {
470    assert_eq!(
471      Point::orient(
472        &Point::new([0, 0]),
473        &Point::new([1, 1]),
474        &Point::new([2, 2])
475      ),
476      CoLinear
477    );
478    assert_eq!(
479      Orientation::new(
480        &Point::new_nn([0.0, 0.0]),
481        &Point::new_nn([1.0, 1.0]),
482        &Point::new_nn([2.0, 2.0])
483      ),
484      CoLinear
485    );
486
487    assert_eq!(
488      Point::orient(
489        &Point::new([0, 0]),
490        &Point::new([0, 1]),
491        &Point::new([2, 2])
492      ),
493      ClockWise
494    );
495    assert_eq!(
496      Point::orient(
497        &Point::new_nn([0.0, 0.0]),
498        &Point::new_nn([0.0, 1.0]),
499        &Point::new_nn([2.0, 2.0])
500      ),
501      ClockWise
502    );
503
504    assert_eq!(
505      Point::orient(
506        &Point::new([0, 0]),
507        &Point::new([0, 1]),
508        &Point::new([-2, 2])
509      ),
510      CounterClockWise
511    );
512    assert_eq!(
513      Point::orient(
514        &Point::new([0, 0]),
515        &Point::new([0, 0]),
516        &Point::new([0, 0])
517      ),
518      CoLinear
519    );
520  }
521
522  #[test]
523  fn unit_1() {
524    assert_eq!(
525      Point::orient(
526        &Point::new([0, 0]),
527        &Point::new([1, 0]),
528        &Point::new([1, 0])
529      ),
530      CoLinear
531    );
532    assert_eq!(
533      Point::orient(
534        &Point::new([0, 0]),
535        &Point::new([1, 0]),
536        &Point::new([2, 0])
537      ),
538      CoLinear
539    );
540    assert_eq!(
541      Point::orient(
542        &Point::new([1, 0]),
543        &Point::new([2, 0]),
544        &Point::new([0, 0])
545      ),
546      CoLinear
547    );
548    assert_eq!(
549      Point::orient(
550        &Point::new([1, 0]),
551        &Point::new([2, 0]),
552        &Point::new([1, 0])
553      ),
554      CoLinear
555    );
556  }
557
558  #[test]
559  fn unit_2() {
560    assert_eq!(
561      Point::orient(
562        &Point::new([1, 0]),
563        &Point::new([0, 6]),
564        &Point::new([0, 8])
565      ),
566      ClockWise
567    );
568  }
569
570  #[test]
571  fn unit_3() {
572    assert_eq!(
573      Point::orient(
574        &Point::new([-12_i8, -126]),
575        &Point::new([-12, -126]),
576        &Point::new([0, -126])
577      ),
578      CoLinear
579    );
580  }
581}
rgeometry/data/point.rs

rgeometry/data/
point.rs
