rgeometry/algorithms/triangulation/
earclip.rs

1use crate::algorithms::zhash::{ZHashBox, ZHashable};
2use crate::data::{Point, PointId, PointLocation, Polygon, TriangleView};
3use crate::Orientation;
4use crate::PolygonScalar;
5
6// use rand::rngs::mock::StepRng;
7use rand::rngs::SmallRng;
8use rand::Rng;
9use rand::SeedableRng;
10
11/// $O(n^2)$ Polygon triangulation. Ears are selected in a pseudo-random manner.
12pub fn earclip<T>(poly: &Polygon<T>) -> impl Iterator<Item = (PointId, PointId, PointId)> + '_
13where
14  T: PolygonScalar,
15{
16  // FIXME: Support holes.
17  assert!(poly.rings.len() == 1);
18  // let rng = StepRng::new(0, 0);
19  let rng = SmallRng::seed_from_u64(0xDEADBEEF);
20  triangulate_list(&poly.points, &poly.rings[0], rng)
21}
22
23/// $O(n)$ Polygon triangulation. Ears are selected in a pseudo-random manner.
24///
25/// $O(n^2)$ worst case complexity. Expected time is linear.
26pub fn earclip_hashed<T>(
27  poly: &Polygon<T>,
28) -> impl Iterator<Item = (PointId, PointId, PointId)> + '_
29where
30  T: PolygonScalar + ZHashable,
31{
32  // FIXME: Support holes.
33  assert!(poly.rings.len() == 1);
34  // let rng = StepRng::new(0, 0);
35  let rng = SmallRng::seed_from_u64(0xDEADBEEF);
36  triangulate_list_hashed(&poly.points, &poly.rings[0], rng)
37}
38
39// triangulate: Vec<Point<T,2>> + Vec<PointId> -> Vec<(PointId,PointId,PointId)>
40// triangulate: Polygon -> Vec<Polygon>
41// triangulate: Polygon -> Vec<(PointId, PointId, PointId)>
42
43// Create a linked list of points. O(n)
44// Clone it as a list of possible ears. O(n)
45// For each vertex in the list of possible ears:
46//   Check if the vertex is an ear. O(n)
47//   Delete it from the list of possible-ears.
48//   If it is:
49//     Delete it from the vertex list.
50//     Insert vertex->prev and vertex->next into possible-ears if they aren't there already.
51//     Emit an edge: (vertex.prev, ear, vertex.next)
52pub fn triangulate_list<'a, T, R>(
53  points: &'a [Point<T>],
54  order: &'a [PointId],
55  mut rng: R,
56) -> impl Iterator<Item = (PointId, PointId, PointId)> + 'a
57where
58  T: PolygonScalar,
59  R: Rng + 'static,
60{
61  let mut len = order.len();
62  let mut vertices = List::new(points, order);
63  let mut possible_ears = EarStore::new(order.len());
64  std::iter::from_fn(move || match len {
65    0..=2 => None,
66    _ => loop {
67      let focus = vertices.cursor(possible_ears.pop(&mut rng).unwrap());
68      let prev = focus.prev();
69      let next = focus.next();
70      if is_ear(prev, focus, next) {
71        possible_ears.new_possible_ear(prev.position);
72        possible_ears.new_possible_ear(next.position);
73        let out = (prev.point_id(), focus.point_id(), next.point_id());
74        let focus = focus.position;
75        vertices.delete(focus);
76        len -= 1;
77        return Some(out);
78      }
79    },
80  })
81}
82
83fn is_ear<T>(a: Cursor<'_, T>, b: Cursor<'_, T>, c: Cursor<'_, T>) -> bool
84where
85  T: PolygonScalar,
86{
87  let trig = TriangleView::new_unchecked([a.point(), b.point(), c.point()]);
88  if trig.orientation() == Orientation::CounterClockWise {
89    let mut focus = c.next();
90    while focus != a {
91      if trig.locate(focus.point()) != PointLocation::Outside {
92        return false;
93      }
94      focus = focus.next();
95    }
96    true
97  } else {
98    false
99  }
100}
101
102///////////////////////////////////////////////////////////////////////////////
103// Z-order hash ear clipping
104
105pub fn triangulate_list_hashed<'a, T, R>(
106  points: &'a [Point<T>],
107  order: &'a [PointId],
108  mut rng: R,
109) -> impl Iterator<Item = (PointId, PointId, PointId)> + 'a
110where
111  T: PolygonScalar + ZHashable,
112  R: Rng + 'static,
113{
114  let mut len = order.len();
115  let mut vertices = List::new(points, order);
116  let mut possible_ears = EarStore::new(order.len());
117
118  let zbox = zbox_slice(points, order);
119  let key = ZHashable::zhash_key(zbox);
120  let zhashes: Vec<u64> = order
121    .iter()
122    .map(|&pid| ZHashable::zhash_fn(key, &points[pid.usize()]))
123    .collect();
124  let mut zorder = List::new_sorted(points, order, zhashes);
125
126  std::iter::from_fn(move || match len {
127    0..=2 => None,
128    _ => loop {
129      let focus = possible_ears.pop(&mut rng).unwrap();
130      let prev = vertices.prev(focus);
131      let next = vertices.next(focus);
132      if is_ear_hashed(
133        key,
134        zorder.cursor(prev),
135        zorder.cursor(focus),
136        zorder.cursor(next),
137      ) {
138        possible_ears.new_possible_ear(prev);
139        possible_ears.new_possible_ear(next);
140        vertices.delete(focus);
141        zorder.delete(focus);
142        len -= 1;
143        let out = (order[prev], order[focus], order[next]);
144        return Some(out);
145      }
146    },
147  })
148}
149
150fn is_ear_hashed<T: PolygonScalar + ZHashable>(
151  key: <T as ZHashable>::ZHashKey,
152  a: Cursor<'_, T>,
153  b: Cursor<'_, T>,
154  c: Cursor<'_, T>,
155) -> bool {
156  let trig = TriangleView::new_unchecked([a.point(), b.point(), c.point()]);
157  if trig.orientation() == Orientation::CounterClockWise {
158    // Points inside the triangle are guaranteed to have a zhash
159    // between the hashes of the bounding box.
160    let (min, max) = trig.bounding_box();
161    let min_hash = ZHashable::zhash_fn(key, &min);
162    let max_hash = ZHashable::zhash_fn(key, &max);
163
164    // Points inside the triangle are likely to be nearby so start
165    // by searching in both directions.
166    let mut up_focus = b.next();
167    let mut down_focus = b.prev();
168
169    let cond = |cursor: Cursor<'_, T>| cursor.valid() && cursor.between(min_hash, max_hash);
170    let check = |cursor: Cursor<'_, T>| {
171      cursor != a
172        && cursor != b
173        && cursor != c
174        && trig.locate(cursor.point()) != PointLocation::Outside
175    };
176
177    while cond(up_focus) && cond(down_focus) {
178      if check(up_focus) || check(down_focus) {
179        return false;
180      }
181      up_focus = up_focus.next();
182      down_focus = down_focus.prev();
183    }
184
185    // Look upwards
186    while cond(up_focus) {
187      if check(up_focus) {
188        return false;
189      }
190      up_focus = up_focus.next();
191    }
192
193    // Look downwards
194    while cond(down_focus) {
195      if check(down_focus) {
196        return false;
197      }
198      down_focus = down_focus.prev();
199    }
200
201    true
202  } else {
203    false
204  }
205}
206
207fn zbox_slice<'a, T>(slice: &'a [Point<T>], order: &'a [PointId]) -> ZHashBox<'a, T>
208where
209  T: PolygonScalar + ZHashable,
210{
211  if slice.is_empty() {
212    panic!("Cannot zhash an empty slice");
213    // return ZHashBox {
214    //   min_x: &T::zero(),
215    //   max_x: &T::zero(),
216    //   min_y: &T::zero(),
217    //   max_y: &T::zero(),
218    // };
219  }
220  let mut zbox = ZHashBox {
221    min_x: slice[0].x_coord(),
222    max_x: slice[0].x_coord(),
223    min_y: slice[0].y_coord(),
224    max_y: slice[0].y_coord(),
225  };
226  for &pid in order {
227    let pt = &slice[pid.usize()];
228    if zbox.min_x > pt.x_coord() {
229      zbox.min_x = pt.x_coord();
230    }
231    if zbox.max_x < pt.x_coord() {
232      zbox.max_x = pt.x_coord();
233    }
234    if zbox.min_y > pt.y_coord() {
235      zbox.min_y = pt.y_coord();
236    }
237    if zbox.max_y < pt.y_coord() {
238      zbox.max_y = pt.y_coord();
239    }
240  }
241  zbox
242}
243
244///////////////////////////////////////////////////////////////////////////////
245// Linked List that supports deletions and re-insertions (of deleted items)
246
247struct Cursor<'a, T> {
248  list: &'a List<'a, T>,
249  position: usize,
250}
251
252impl<T> Eq for Cursor<'_, T> {}
253
254impl<'a, T> PartialEq for Cursor<'a, T> {
255  fn eq(&self, other: &Cursor<'a, T>) -> bool {
256    self.position == other.position
257  }
258}
259
260impl<T> Copy for Cursor<'_, T> {}
261
262impl<'a, T> Clone for Cursor<'a, T> {
263  fn clone(&self) -> Cursor<'a, T> {
264    *self
265  }
266}
267
268impl<'a, T> Cursor<'a, T> {
269  fn valid(&self) -> bool {
270    self.position != usize::MAX
271  }
272
273  fn next(mut self) -> Cursor<'a, T> {
274    self.position = self.list.next(self.position);
275    self
276  }
277
278  fn prev(mut self) -> Cursor<'a, T> {
279    self.position = self.list.prev(self.position);
280    self
281  }
282
283  fn point(&self) -> &Point<T> {
284    self.list.point(self.position)
285  }
286
287  fn point_id(&self) -> PointId {
288    self.list.point_id(self.position)
289  }
290
291  fn hash(&self) -> u64 {
292    self.list.hash(self.position)
293  }
294
295  fn between(&self, min: u64, max: u64) -> bool {
296    self.hash() >= min && self.hash() <= max
297  }
298}
299
300struct List<'a, T> {
301  points: &'a [Point<T>],
302  order: &'a [PointId],
303  hashes: Vec<u64>,
304  prev: Vec<usize>,
305  next: Vec<usize>,
306}
307
308impl<'a, T> List<'a, T> {
309  fn new(points: &'a [Point<T>], order: &'a [PointId]) -> List<'a, T> {
310    let size = order.len();
311    let mut prev = Vec::with_capacity(size);
312    let mut next = Vec::with_capacity(size);
313    prev.resize(size, 0);
314    next.resize(size, 0);
315    for i in 0..size {
316      prev[(i + 1) % size] = i;
317      next[i] = (i + 1) % size;
318    }
319    List {
320      points,
321      order,
322      hashes: vec![],
323      prev,
324      next,
325    }
326  }
327
328  fn prev(&self, vertex: usize) -> usize {
329    self.prev[vertex]
330  }
331
332  fn next(&self, vertex: usize) -> usize {
333    self.next[vertex]
334  }
335
336  fn delete(&mut self, vertex: usize) {
337    let prev = self.prev[vertex];
338    let next = self.next[vertex];
339    if prev != usize::MAX {
340      self.next[prev] = next;
341    }
342    if next != usize::MAX {
343      self.prev[next] = prev;
344    }
345  }
346
347  fn cursor(&self, vertex: usize) -> Cursor<'_, T> {
348    Cursor {
349      list: self,
350      position: vertex,
351    }
352  }
353
354  fn point(&self, vertex: usize) -> &Point<T> {
355    &self.points[self.point_id(vertex).usize()]
356  }
357
358  fn point_id(&self, vertex: usize) -> PointId {
359    self.order[vertex]
360  }
361
362  fn hash(&self, vertex: usize) -> u64 {
363    self.hashes[vertex]
364  }
365}
366
367impl<'a, T> List<'a, T> {
368  fn new_sorted(points: &'a [Point<T>], order: &'a [PointId], keys: Vec<u64>) -> List<'a, T> {
369    let size = keys.len();
370    let mut v: Vec<usize> = (0..size).collect();
371    v.sort_unstable_by_key(|&idx| keys[idx]);
372    let mut prev = Vec::with_capacity(size);
373    let mut next = Vec::with_capacity(size);
374    prev.resize(size, 0);
375    next.resize(size, 0);
376    for i in 0..size {
377      prev[v[(i + 1) % size]] = v[i];
378      next[v[i]] = v[(i + 1) % size];
379    }
380    next[v[size - 1]] = usize::MAX;
381    prev[v[0]] = usize::MAX;
382    List {
383      points,
384      order,
385      hashes: keys,
386      prev,
387      next,
388    }
389  }
390}
391
392///////////////////////////////////////////////////////////////////////////////
393// Collection of possible ears.
394
395struct EarStore {
396  possible_ears_vec: Vec<usize>,
397  possible_ears_set: IntSet,
398}
399
400impl EarStore {
401  fn new(size: usize) -> EarStore {
402    EarStore {
403      possible_ears_vec: (0..size).collect(),
404      possible_ears_set: IntSet::with_capacity(size),
405    }
406  }
407
408  fn new_possible_ear(&mut self, possible_ear: usize) {
409    if !self.possible_ears_set.contains(possible_ear) {
410      self.possible_ears_set.insert(possible_ear);
411      self.possible_ears_vec.push(possible_ear);
412    }
413  }
414
415  fn pop<R>(&mut self, rng: &mut R) -> Option<usize>
416  where
417    R: Rng + ?Sized,
418  {
419    let n = rng.gen_range(0..self.possible_ears_vec.len());
420    let next = self.possible_ears_vec.swap_remove(n);
421    self.possible_ears_set.delete(next);
422    Some(next)
423  }
424}
425
426///////////////////////////////////////////////////////////////////////////////
427// IntSet
428
429// FIXME: Use a bitset?
430struct IntSet {
431  set: Vec<bool>,
432}
433
434impl IntSet {
435  fn with_capacity(capacity: usize) -> IntSet {
436    let mut set = Vec::with_capacity(capacity);
437    set.resize(capacity, true);
438    IntSet { set }
439  }
440
441  fn contains(&self, value: usize) -> bool {
442    self.set[value]
443  }
444
445  fn insert(&mut self, value: usize) {
446    self.set[value] = true
447  }
448
449  fn delete(&mut self, value: usize) {
450    self.set[value] = false
451  }
452}
453
454///////////////////////////////////////////////////////////////////////////////
455// Tests
456
457#[cfg(test)]
458#[cfg(not(tarpaulin_include))]
459mod tests {
460  use super::*;
461  use crate::data::*;
462  use num_bigint::BigInt;
463  use rand::rngs::SmallRng;
464  use rand::SeedableRng;
465
466  fn trig_area_2x<F: PolygonScalar + Into<BigInt>>(p: &Polygon<F>) -> BigInt {
467    let mut trig_area_2x = BigInt::from(0);
468    // let mut rng = StepRng::new(0,0);
469    let rng = SmallRng::seed_from_u64(0);
470    for (a, b, c) in triangulate_list(&p.points, &p.rings[0], rng) {
471      let trig = TriangleView::new_unchecked([p.point(a), p.point(b), p.point(c)]);
472      trig_area_2x += trig.signed_area_2x::<BigInt>();
473    }
474    trig_area_2x
475  }
476
477  #[test]
478  fn basic_1() {
479    let p = Polygon::new(vec![
480      Point::new([0, 0]),
481      Point::new([1, 0]),
482      Point::new([1, 1]),
483    ])
484    .unwrap();
485
486    assert_eq!(p.signed_area_2x::<BigInt>(), trig_area_2x(&p));
487  }
488
489  #[test]
490  fn basic_2() {
491    let p = Polygon::new(vec![
492      Point::new([0, 0]),
493      Point::new([1, 0]),
494      Point::new([2, 0]),
495      Point::new([3, 0]),
496      Point::new([4, 0]),
497      Point::new([1, 1]),
498    ])
499    .unwrap();
500
501    assert_eq!(p.signed_area_2x::<BigInt>(), trig_area_2x(&p));
502  }
503
504  #[test]
505  fn basic_3() {
506    let rng = SmallRng::seed_from_u64(0);
507    let p: Polygon<i8> = Polygon::new(vec![
508      Point { array: [-44, -11] },
509      Point { array: [-43, 23] },
510      Point { array: [-64, 44] },
511      Point { array: [-52, 114] },
512      Point { array: [-82, 69] },
513    ])
514    .unwrap();
515
516    triangulate_list_hashed(&p.points, &p.rings[0], rng).count();
517  }
518
519  #[test]
520  fn basic_4() {
521    let rng = SmallRng::seed_from_u64(0);
522    let p: Polygon<i8> = Polygon::new(vec![
523      Point { array: [12, 5] },   // 0
524      Point { array: [0, 8] },    // 1
525      Point { array: [-10, -6] }, // 2 cut
526      Point { array: [-3, 3] },   // 3
527      Point { array: [-2, 4] },   // 4
528    ])
529    .unwrap();
530
531    triangulate_list_hashed(&p.points, &p.rings[0], rng).count();
532  }
533
534  use proptest::prelude::*;
535  use test_strategy::proptest;
536
537  #[proptest]
538  fn equal_area_prop(poly: Polygon<i64>) {
539    prop_assert_eq!(poly.signed_area_2x::<BigInt>(), trig_area_2x(&poly));
540  }
541
542  #[proptest]
543  fn hashed_identity_prop(poly: Polygon<i8>) {
544    let rng = SmallRng::seed_from_u64(0);
545    let not_hashed: Vec<(PointId, PointId, PointId)> =
546      triangulate_list(&poly.points, &poly.rings[0], rng.clone()).collect();
547    let hashed: Vec<(PointId, PointId, PointId)> =
548      triangulate_list_hashed(&poly.points, &poly.rings[0], rng).collect();
549    prop_assert_eq!(not_hashed, hashed);
550  }
551}
rgeometry/algorithms/triangulation/earclip.rs

rgeometry/algorithms/triangulation/
earclip.rs
