rgeometry/
lib.rs

1// #![deny(warnings)]
2#![deny(clippy::cast_precision_loss)]
3#![deny(clippy::lossy_float_literal)]
4#![doc(test(no_crate_inject))]
5#![doc(html_playground_url = "https://rgeometry.org/playground.html")]
6#![doc(test(no_crate_inject))]
7use num_traits::*;
8use std::cmp::Ordering;
9use std::iter::Sum;
10use std::ops::*;
11
12pub mod algorithms;
13pub mod data;
14mod intersection;
15mod matrix;
16mod orientation;
17mod transformation;
18mod utils;
19
20pub use orientation::{Orientation, SoS};
21
22pub use intersection::Intersects;
23
24#[derive(Debug, Clone, Copy, PartialEq, Eq)]
25pub enum Error {
26  InsufficientVertices,
27  SelfIntersections,
28  DuplicatePoints,
29  /// Two consecutive line segments are either colinear or oriented clockwise.
30  ConvexViolation,
31  ClockWiseViolation,
32  CoLinearViolation,
33}
34
35impl std::fmt::Display for Error {
36  fn fmt(&self, f: &mut std::fmt::Formatter<'_>) -> Result<(), std::fmt::Error> {
37    match self {
38      Error::InsufficientVertices => write!(f, "Insufficient vertices"),
39      Error::SelfIntersections => write!(f, "Self intersections"),
40      Error::DuplicatePoints => write!(f, "Duplicate points"),
41      Error::ConvexViolation => write!(f, "Convex violation"),
42      Error::ClockWiseViolation => write!(f, "Clockwise violation"),
43      Error::CoLinearViolation => write!(
44        f,
45        "Two or more points are colinear and no valid solution exists"
46      ),
47    }
48  }
49}
50
51pub trait TotalOrd {
52  fn total_cmp(&self, other: &Self) -> Ordering;
53
54  fn total_min(self, other: Self) -> Self
55  where
56    Self: Sized,
57  {
58    std::cmp::min_by(self, other, TotalOrd::total_cmp)
59  }
60
61  fn total_max(self, other: Self) -> Self
62  where
63    Self: Sized,
64  {
65    std::cmp::max_by(self, other, TotalOrd::total_cmp)
66  }
67}
68
69impl TotalOrd for u32 {
70  fn total_cmp(&self, other: &Self) -> Ordering {
71    self.cmp(other)
72  }
73}
74
75impl<A: TotalOrd> TotalOrd for &A {
76  fn total_cmp(&self, other: &Self) -> Ordering {
77    (*self).total_cmp(*other)
78  }
79}
80
81impl<A: TotalOrd, B: TotalOrd> TotalOrd for (A, B) {
82  fn total_cmp(&self, other: &Self) -> Ordering {
83    self
84      .0
85      .total_cmp(&other.0)
86      .then_with(|| self.1.total_cmp(&other.1))
87  }
88}
89
90// FIXME: Should include ZHashable.
91pub trait PolygonScalar:
92  std::fmt::Debug
93  + Neg<Output = Self>
94  + NumAssignOps
95  + NumOps<Self, Self>
96  + TotalOrd
97  + PartialOrd
98  + Sum
99  + Clone
100{
101  fn from_constant(val: i8) -> Self;
102  fn cmp_dist(p: &[Self; 2], q: &[Self; 2], r: &[Self; 2]) -> std::cmp::Ordering;
103  fn cmp_slope(p: &[Self; 2], q: &[Self; 2], r: &[Self; 2]) -> std::cmp::Ordering;
104  fn cmp_vector_slope(p: &[Self; 2], q: &[Self; 2], r: &[Self; 2]) -> std::cmp::Ordering;
105  fn cmp_perp_vector_slope(p: &[Self; 2], q: &[Self; 2], r: &[Self; 2]) -> std::cmp::Ordering;
106}
107
108macro_rules! fixed_precision {
109  ( $ty:ty, $uty:ty, $long:ty, $ulong: ty ) => {
110    impl TotalOrd for $ty {
111      fn total_cmp(&self, other: &Self) -> Ordering {
112        self.cmp(other)
113      }
114    }
115
116    impl PolygonScalar for $ty {
117      fn from_constant(val: i8) -> Self {
118        val as $ty
119      }
120      fn cmp_dist(p: &[Self; 2], q: &[Self; 2], r: &[Self; 2]) -> std::cmp::Ordering {
121        fn diff(a: $ty, b: $ty) -> $ulong {
122          if b > a {
123            b.wrapping_sub(a) as $uty as $ulong
124          } else {
125            a.wrapping_sub(b) as $uty as $ulong
126          }
127        }
128        let pq_x = diff(p[0], q[0]);
129        let pq_y = diff(p[1], q[1]);
130        let (pq_dist_squared, pq_overflow) = (pq_x * pq_x).overflowing_add(pq_y * pq_y);
131        let pr_x = diff(p[0], r[0]);
132        let pr_y = diff(p[1], r[1]);
133        let (pr_dist_squared, pr_overflow) = (pr_x * pr_x).overflowing_add(pr_y * pr_y);
134        match (pq_overflow, pr_overflow) {
135          (true, false) => Ordering::Greater,
136          (false, true) => Ordering::Less,
137          _ => pq_dist_squared.cmp(&pr_dist_squared),
138        }
139      }
140
141      fn cmp_slope(p: &[Self; 2], q: &[Self; 2], r: &[Self; 2]) -> std::cmp::Ordering {
142        // Return the absolute difference along with its sign.
143        // diff(0, 10) => (10, true)
144        // diff(10, 0) => (10, false)
145        // diff(i8::MIN,i8:MAX) => (255_u16, true)
146        // diff(a,b) = (c, sign) where a = if sign { b-c } else { b+c }
147        fn diff(a: $ty, b: $ty) -> ($ulong, bool) {
148          if b > a {
149            (b.wrapping_sub(a) as $uty as $ulong, true)
150          } else {
151            (a.wrapping_sub(b) as $uty as $ulong, false)
152          }
153        }
154        let (ux, ux_neg) = diff(q[0], p[0]);
155        let (vy, vy_neg) = diff(r[1], p[1]);
156        let ux_vy_neg = ux_neg.bitxor(vy_neg) && ux != 0 && vy != 0;
157        let (uy, uy_neg) = diff(q[1], p[1]);
158        let (vx, vx_neg) = diff(r[0], p[0]);
159        let uy_vx_neg = uy_neg.bitxor(vx_neg) && uy != 0 && vx != 0;
160        match (ux_vy_neg, uy_vx_neg) {
161          (true, false) => Ordering::Less,
162          (false, true) => Ordering::Greater,
163          (true, true) => (uy * vx).cmp(&(ux * vy)),
164          (false, false) => (ux * vy).cmp(&(uy * vx)),
165        }
166      }
167
168      fn cmp_vector_slope(vector: &[Self; 2], p: &[Self; 2], q: &[Self; 2]) -> std::cmp::Ordering {
169        // Return the absolute difference along with its sign.
170        // diff(0, 10) => (10, true)
171        // diff(10, 0) => (10, false)
172        // diff(i8::MIN,i8:MAX) => (255_u16, true)
173        // diff(a,b) = (c, sign) where a = if sign { b-c } else { b+c }
174        fn diff(a: $ty, b: $ty) -> ($ulong, bool) {
175          if b > a {
176            (b.wrapping_sub(a) as $uty as $ulong, true)
177          } else {
178            (a.wrapping_sub(b) as $uty as $ulong, false)
179          }
180        }
181        let (ux, ux_neg) = (vector[0].unsigned_abs() as $ulong, vector[0] < 0);
182        let (vy, vy_neg) = diff(q[1], p[1]);
183        // neg xor neg = pos = 0
184        // neg xor pos = neg = 1
185        // pos xor neg = neg = 1
186        // pos xor pos = pos = 0
187        let ux_vy_neg = ux_neg.bitxor(vy_neg) && ux != 0 && vy != 0;
188        let (uy, uy_neg) = (vector[1].unsigned_abs() as $ulong, vector[1] < 0);
189        let (vx, vx_neg) = diff(q[0], p[0]);
190        let uy_vx_neg = uy_neg.bitxor(vx_neg) && uy != 0 && vx != 0;
191        match (ux_vy_neg, uy_vx_neg) {
192          (true, false) => Ordering::Less,
193          (false, true) => Ordering::Greater,
194          (true, true) => (uy * vx).cmp(&(ux * vy)),
195          (false, false) => (ux * vy).cmp(&(uy * vx)),
196        }
197      }
198
199      fn cmp_perp_vector_slope(
200        vector: &[Self; 2],
201        p: &[Self; 2],
202        q: &[Self; 2],
203      ) -> std::cmp::Ordering {
204        // Return the absolute difference along with its sign.
205        // diff(0, 10) => (10, true)
206        // diff(10, 0) => (10, false)
207        // diff(i8::MIN,i8:MAX) => (255_u16, true)
208        // diff(a,b) = (c, sign) where a = if sign { b-c } else { b+c }
209        fn diff(a: $ty, b: $ty) -> ($ulong, bool) {
210          if b > a {
211            (b.wrapping_sub(a) as $uty as $ulong, true)
212          } else {
213            (a.wrapping_sub(b) as $uty as $ulong, false)
214          }
215        }
216        let (ux, ux_neg) = (vector[1].unsigned_abs() as $ulong, vector[1] > 0);
217        let (vy, vy_neg) = diff(q[1], p[1]);
218        // neg xor neg = pos = 0
219        // neg xor pos = neg = 1
220        // pos xor neg = neg = 1
221        // pos xor pos = pos = 0
222        let ux_vy_neg = ux_neg.bitxor(vy_neg) && ux != 0 && vy != 0;
223        let (uy, uy_neg) = (vector[0].unsigned_abs() as $ulong, vector[0] < 0);
224        let (vx, vx_neg) = diff(q[0], p[0]);
225        let uy_vx_neg = uy_neg.bitxor(vx_neg) && uy != 0 && vx != 0;
226        match (ux_vy_neg, uy_vx_neg) {
227          (true, false) => Ordering::Less,
228          (false, true) => Ordering::Greater,
229          (true, true) => (uy * vx).cmp(&(ux * vy)),
230          (false, false) => (ux * vy).cmp(&(uy * vx)),
231        }
232      }
233    }
234  };
235}
236
237macro_rules! arbitrary_precision {
238  ( $( $ty:ty ),* ) => {
239    $(
240      impl TotalOrd for $ty {
241        fn total_cmp(&self, other: &Self) -> Ordering {
242          self.cmp(other)
243        }
244      }
245
246      impl PolygonScalar for $ty {
247      fn from_constant(val: i8) -> Self {
248        <$ty>::from_i8(val).unwrap()
249      }
250      fn cmp_dist(p: &[Self; 2], q: &[Self; 2], r: &[Self; 2]) -> std::cmp::Ordering {
251        let pq_x = &p[0] - &q[0];
252        let pq_y = &p[1] - &q[1];
253        let pq_dist_squared: Self = &pq_x*&pq_x + &pq_y*&pq_y;
254        let pr_x = &p[0] - &r[0];
255        let pr_y = &p[1] - &r[1];
256        let pr_dist_squared: Self = &pr_x*&pr_x + &pr_y*&pr_y;
257        pq_dist_squared.cmp(&pr_dist_squared)
258      }
259
260      fn cmp_slope(p: &[Self; 2], q: &[Self; 2], r: &[Self; 2]) -> std::cmp::Ordering {
261        let slope1 = (&r[1] - &q[1]) * (&q[0] - &p[0]);
262        let slope2 = (&q[1] - &p[1]) * (&r[0] - &q[0]);
263        slope1.cmp(&slope2)
264      }
265      fn cmp_vector_slope(vector: &[Self;2], p: &[Self; 2], q: &[Self; 2]) -> std::cmp::Ordering {
266        PolygonScalar::cmp_slope(
267          p,
268          &[&p[0] + &vector[0], &p[1] + &vector[1]],
269          q
270        )
271      }
272      fn cmp_perp_vector_slope(vector: &[Self;2], p: &[Self; 2], q: &[Self; 2]) -> std::cmp::Ordering {
273        PolygonScalar::cmp_slope(
274          p,
275          &[&p[0] - &vector[1], &p[1] + &vector[0]],
276          q
277        )
278      }
279    })*
280  };
281}
282
283macro_rules! wrapped_floating_precision {
284  ( $( $ty:ty ),* ) => {
285    $(
286      impl TotalOrd for $ty {
287        fn total_cmp(&self, other: &Self) -> Ordering {
288          self.cmp(other)
289        }
290      }
291
292      impl PolygonScalar for $ty {
293      fn from_constant(val: i8) -> Self {
294        <$ty>::from_i8(val).unwrap()
295      }
296      // FIXME: Use `geometry_predicates` to speed up calculation. Right now we're
297      // roughly 100x slower than necessary.
298      fn cmp_dist(p: &[Self; 2], q: &[Self; 2], r: &[Self; 2]) -> std::cmp::Ordering {
299        PolygonScalar::cmp_dist(
300          &[float_to_rational(p[0].into_inner()), float_to_rational(p[1].into_inner())],
301          &[float_to_rational(q[0].into_inner()), float_to_rational(q[1].into_inner())],
302          &[float_to_rational(r[0].into_inner()), float_to_rational(r[1].into_inner())],
303        )
304      }
305
306      // This function uses the arbitrary precision machinery of `geometry_predicates` to
307      // quickly compute the orientation of three 2D points. This is about 10x-50x slower
308      // than the inexact version.
309      fn cmp_slope(p: &[Self; 2], q: &[Self; 2], r: &[Self; 2]) -> std::cmp::Ordering {
310        let orient = geometry_predicates::predicates::orient2d(
311          [p[0].into_inner() as f64, p[1].into_inner() as f64],
312          [q[0].into_inner() as f64, q[1].into_inner() as f64],
313          [r[0].into_inner() as f64, r[1].into_inner() as f64],
314        );
315        if orient > 0.0 {
316          Ordering::Greater
317        } else if orient < 0.0 {
318          Ordering::Less
319        } else {
320          Ordering::Equal
321        }
322      }
323      // FIXME: Use `geometry_predicates` to speed up calculation. Right now we're
324      // roughly 100x slower than necessary.
325      fn cmp_vector_slope(vector: &[Self;2], p: &[Self; 2], q: &[Self; 2]) -> std::cmp::Ordering {
326        PolygonScalar::cmp_vector_slope(
327          &[float_to_rational(vector[0].into_inner()), float_to_rational(vector[1].into_inner())],
328          &[float_to_rational(p[0].into_inner()), float_to_rational(p[1].into_inner())],
329          &[float_to_rational(q[0].into_inner()), float_to_rational(q[1].into_inner())],
330        )
331      }
332      // FIXME: Use `geometry_predicates` to speed up calculation. Right now we're
333      // roughly 100x slower than necessary.
334      fn cmp_perp_vector_slope(vector: &[Self;2], p: &[Self; 2], q: &[Self; 2]) -> std::cmp::Ordering {
335        PolygonScalar::cmp_perp_vector_slope(
336          &[float_to_rational(vector[0].into_inner()), float_to_rational(vector[1].into_inner())],
337          &[float_to_rational(p[0].into_inner()), float_to_rational(p[1].into_inner())],
338          &[float_to_rational(q[0].into_inner()), float_to_rational(q[1].into_inner())],
339        )
340      }
341    })*
342  };
343}
344
345macro_rules! floating_precision {
346  ( $( $ty:ty ),* ) => {
347    $(
348      impl TotalOrd for $ty {
349        fn total_cmp(&self, other: &Self) -> Ordering {
350          <$ty>::total_cmp(self, other)
351        }
352      }
353
354      impl PolygonScalar for $ty {
355      fn from_constant(val: i8) -> Self {
356        <$ty>::from_i8(val).unwrap()
357      }
358      // FIXME: Use `geometry_predicates` to speed up calculation. Right now we're
359      // roughly 100x slower than necessary.
360      fn cmp_dist(p: &[Self; 2], q: &[Self; 2], r: &[Self; 2]) -> std::cmp::Ordering {
361        PolygonScalar::cmp_dist(
362          &[float_to_rational(p[0]), float_to_rational(p[1])],
363          &[float_to_rational(q[0]), float_to_rational(q[1])],
364          &[float_to_rational(r[0]), float_to_rational(r[1])],
365        )
366      }
367
368      // This function uses the arbitrary precision machinery of `geometry_predicates` to
369      // quickly compute the orientation of three 2D points. This is about 10x-50x slower
370      // than the inexact version.
371      fn cmp_slope(p: &[Self; 2], q: &[Self; 2], r: &[Self; 2]) -> std::cmp::Ordering {
372        let orient = geometry_predicates::predicates::orient2d(
373          [p[0] as f64, p[1] as f64],
374          [q[0] as f64, q[1] as f64],
375          [r[0] as f64, r[1] as f64],
376        );
377        if orient > 0.0 {
378          Ordering::Greater
379        } else if orient < 0.0 {
380          Ordering::Less
381        } else {
382          Ordering::Equal
383        }
384      }
385      // FIXME: Use `geometry_predicates` to speed up calculation. Right now we're
386      // roughly 100x slower than necessary.
387      fn cmp_vector_slope(vector: &[Self;2], p: &[Self; 2], q: &[Self; 2]) -> std::cmp::Ordering {
388        PolygonScalar::cmp_vector_slope(
389          &[float_to_rational(vector[0]), float_to_rational(vector[1])],
390          &[float_to_rational(p[0]), float_to_rational(p[1])],
391          &[float_to_rational(q[0]), float_to_rational(q[1])],
392        )
393      }
394      // FIXME: Use `geometry_predicates` to speed up calculation. Right now we're
395      // roughly 100x slower than necessary.
396      fn cmp_perp_vector_slope(vector: &[Self;2], p: &[Self; 2], q: &[Self; 2]) -> std::cmp::Ordering {
397        PolygonScalar::cmp_perp_vector_slope(
398          &[float_to_rational(vector[0]), float_to_rational(vector[1])],
399          &[float_to_rational(p[0]), float_to_rational(p[1])],
400          &[float_to_rational(q[0]), float_to_rational(q[1])],
401        )
402      }
403    })*
404  };
405}
406
407fixed_precision!(i8, u8, i16, u16);
408fixed_precision!(i16, u16, i32, u32);
409fixed_precision!(i32, u32, i64, u64);
410fixed_precision!(i64, u64, i128, u128);
411fixed_precision!(isize, usize, i128, u128);
412arbitrary_precision!(num_bigint::BigInt);
413arbitrary_precision!(num_rational::BigRational);
414wrapped_floating_precision!(ordered_float::OrderedFloat<f32>);
415wrapped_floating_precision!(ordered_float::OrderedFloat<f64>);
416wrapped_floating_precision!(ordered_float::NotNan<f32>);
417wrapped_floating_precision!(ordered_float::NotNan<f64>);
418floating_precision!(f32);
419floating_precision!(f64);
420
421#[cfg(feature = "rug")]
422impl TotalOrd for rug::Integer {
423  fn total_cmp(&self, other: &Self) -> Ordering {
424    self.cmp(other)
425  }
426}
427
428#[cfg(feature = "rug")]
429impl PolygonScalar for rug::Integer {
430  fn from_constant(val: i8) -> Self {
431    rug::Integer::from(val)
432  }
433  fn cmp_dist(p: &[Self; 2], q: &[Self; 2], r: &[Self; 2]) -> std::cmp::Ordering {
434    let [qx, qy] = q.clone();
435    let [px, py] = p.clone();
436    let pq_x = px - qx;
437    let pq_y = py - qy;
438    let pq_dist_squared: Self = pq_x.square() + pq_y.square();
439    let [rx, ry] = r.clone();
440    let [px, py] = p.clone();
441    let pr_x = px - rx;
442    let pr_y = py - ry;
443    let pr_dist_squared: Self = pr_x.square() + pr_y.square();
444    pq_dist_squared.cmp(&pr_dist_squared)
445  }
446
447  fn cmp_slope(p: &[Self; 2], q: &[Self; 2], r: &[Self; 2]) -> std::cmp::Ordering {
448    let [qx, qy] = q.clone();
449    let [rx, ry] = r.clone();
450    let ry_qy = ry - &q[1];
451    let qx_px = qx - &p[0];
452    let slope1 = ry_qy * qx_px;
453    let qy_py = qy - &p[1];
454    let rx_qx = rx - &q[0];
455    let slope2 = qy_py * rx_qx;
456    slope1.cmp(&slope2)
457  }
458  fn cmp_vector_slope(vector: &[Self; 2], p: &[Self; 2], q: &[Self; 2]) -> std::cmp::Ordering {
459    let new_x = rug::Integer::from(&p[0] + &vector[0]);
460    let new_y = rug::Integer::from(&p[1] + &vector[1]);
461    PolygonScalar::cmp_slope(p, &[new_x, new_y], q)
462  }
463  fn cmp_perp_vector_slope(vector: &[Self; 2], p: &[Self; 2], q: &[Self; 2]) -> std::cmp::Ordering {
464    let new_x = rug::Integer::from(&p[0] + &vector[1]);
465    let new_y = rug::Integer::from(&p[1] + &vector[0]);
466    PolygonScalar::cmp_slope(p, &[new_x, new_y], q)
467  }
468}
469
470fn float_to_rational(f: impl num::traits::float::FloatCore) -> num::BigRational {
471  num::BigRational::from_float(f).expect("cannot convert NaN or infinite to exact precision number")
472}
473
474#[cfg(test)]
475pub mod testing;
rgeometry/lib.rs

rgeometry/
lib.rs
