rgeometry/algorithms/convex_hull/
melkman.rs

1use crate::data::{Point, Polygon, PolygonConvex};
2
3use crate::PolygonScalar;
4use std::collections::VecDeque;
5
6// Wikipedia description of the convex hull problem: https://en.wikipedia.org/wiki/Convex_hull_of_a_simple_polygon
7// Melkman's algorithm: https://citeseerx.ist.psu.edu/viewdoc/download?doi=10.1.1.512.9681&rep=rep1&type=pdf
8
9// In the right region we pop from the left side, and vice versa.
10// In the composition region, which is the region shared by right and left, we pop the deque from both sides.
11
12pub fn convex_hull<T>(polygon: &Polygon<T>) -> PolygonConvex<T>
13where
14  T: PolygonScalar,
15{
16  let mut convex_hull: VecDeque<&Point<T, 2>> = VecDeque::new();
17  let mut last_idx = 0;
18  for p in polygon.iter() {
19    // Creat a deque with the first 3 points
20    if convex_hull.len() < 2 {
21      convex_hull.push_back(p);
22
23    // Check for colinear of the first 3 points and remove last verdix if so
24    } else if convex_hull.len() == 2 {
25      if Point::orient(convex_hull[0], convex_hull[1], p).is_colinear() {
26        convex_hull.pop_back();
27        convex_hull.push_back(p);
28        continue;
29      }
30      convex_hull.push_front(p);
31      convex_hull.push_back(p);
32
33      // Check and correct(if needed) the orientation if the first 3 verdices
34      if Point::orient(convex_hull[1], convex_hull[2], convex_hull[3]).is_cw() {
35        convex_hull.make_contiguous().reverse();
36      }
37      last_idx = convex_hull.len() - 1;
38
39    // If the new point is within the polygon, then don't add it to convex hull
40    } else if Point::orient(convex_hull[last_idx - 1], convex_hull[last_idx], p).is_ccw()
41      && Point::orient(p, convex_hull[1], convex_hull[0]).is_cw()
42    {
43      continue;
44    // Check for wrong rotations/colinear (fron and back) and remove verdices until correct
45    } else {
46      convex_hull.push_front(p);
47      convex_hull.push_back(p);
48      last_idx = convex_hull.len() - 1;
49
50      while Point::orient(
51        convex_hull[last_idx - 2],
52        convex_hull[last_idx - 1],
53        convex_hull[last_idx],
54      )
55      .is_cw()
56        || Point::orient(
57          convex_hull[last_idx - 2],
58          convex_hull[last_idx - 1],
59          convex_hull[last_idx],
60        )
61        .is_colinear()
62      {
63        convex_hull.remove(last_idx - 1);
64        last_idx = convex_hull.len() - 1;
65      }
66
67      while Point::orient(convex_hull[2], convex_hull[1], convex_hull[0]).is_ccw()
68        || Point::orient(convex_hull[2], convex_hull[1], convex_hull[0]).is_colinear()
69      {
70        convex_hull.remove(1);
71        last_idx = convex_hull.len() - 1;
72      }
73      if Point::orient(convex_hull[1], convex_hull[0], convex_hull[last_idx - 1]).is_colinear() {
74        convex_hull.remove(0);
75        last_idx = convex_hull.len() - 1;
76      }
77    }
78  }
79  // Pop last duplicated verdix
80  convex_hull.pop_back();
81
82  let polygon = Polygon::new_unchecked(convert_deque_to_vec(convex_hull));
83  PolygonConvex::new_unchecked(polygon)
84}
85
86fn convert_deque_to_vec<T: PolygonScalar>(dque: VecDeque<&Point<T, 2>>) -> Vec<Point<T, 2>> {
87  let mut vec: Vec<Point<T, 2>> = Vec::new();
88  for i in dque {
89    vec.push(i.clone());
90  }
91  vec
92}
93
94#[cfg(test)]
95mod tests {
96  use super::*;
97
98  use claims::assert_ok;
99
100  use test_strategy::proptest;
101
102  #[test]
103  fn unit_test_1() {
104    let input = Polygon::new(vec![
105      Point::new([0, 0]),
106      Point::new([2, 0]),
107      Point::new([2, 2]),
108      Point::new([1, 1]),
109      Point::new([0, 2]),
110    ])
111    .unwrap();
112    let output = Polygon::new(vec![
113      Point::new([0, 0]),
114      Point::new([2, 0]),
115      Point::new([2, 2]),
116      Point::new([0, 2]),
117    ])
118    .unwrap();
119    // println!("EXPECTED --- {:?}\n", output);
120    // println!("GOT --- {:?}\n", convex_hull(&input));
121    assert!(Polygon::equals(&convex_hull(&input), &output));
122  }
123
124  #[test]
125  fn unit_test_2() {
126    let input = Polygon::new(vec![
127      Point::new([0, 0]),
128      Point::new([1, 0]),
129      Point::new([2, 0]),
130      Point::new([2, 1]),
131      Point::new([2, 2]),
132      Point::new([1, 2]),
133      Point::new([0, 2]),
134      Point::new([0, 1]),
135    ])
136    .unwrap();
137    let output = Polygon::new(vec![
138      Point::new([0, 0]),
139      Point::new([2, 0]),
140      Point::new([2, 2]),
141      Point::new([0, 2]),
142    ])
143    .unwrap();
144    // println!("EXPECTED --- {:?}\n", output);
145    // println!("GOT --- {:?}\n", convex_hull(&input));
146    assert!(Polygon::equals(&convex_hull(&input), &output));
147  }
148
149  #[test]
150  fn unit_test_3() {
151    let input = Polygon::new(vec![
152      Point::new([0, 0]),
153      Point::new([2, 0]),
154      Point::new([3, 2]),
155      Point::new([2, 2]),
156      Point::new([-1, 3]),
157      Point::new([-1, 2]),
158      Point::new([-1, 1]),
159      Point::new([0, 2]),
160    ])
161    .unwrap();
162    let output = Polygon::new(vec![
163      Point::new([0, 0]),
164      Point::new([2, 0]),
165      Point::new([3, 2]),
166      Point::new([-1, 3]),
167      Point::new([-1, 1]),
168    ])
169    .unwrap();
170    // println!("EXPECTED --- {:?}\n", output);
171    // println!("GOT --- {:?}\n", convex_hull(&input));
172    assert!(Polygon::equals(&convex_hull(&input), &output));
173  }
174
175  #[proptest]
176  fn does_not_panic(poly: Polygon<i8>) {
177    convex_hull(&poly);
178  }
179
180  #[proptest]
181  fn is_valid_convex_polygon(poly: Polygon<i8>) {
182    assert_ok!(convex_hull(&poly).validate());
183  }
184
185  #[proptest]
186  fn is_idempotent(poly: PolygonConvex<i8>) {
187    assert!(Polygon::equals(&convex_hull(poly.polygon()), &poly));
188  }
189
190  #[proptest]
191  fn graham_scan_eq_prop(poly: Polygon<i8>) {
192    let points: Vec<Point<i8>> = poly
193      .iter_boundary()
194      .map(|cursor| cursor.point())
195      .cloned()
196      .collect();
197    let by_scan = crate::algorithms::convex_hull::graham_scan::convex_hull(points).unwrap();
198    assert!(Polygon::equals(&convex_hull(&poly), &by_scan));
199  }
200}
rgeometry/algorithms/convex_hull/melkman.rs

rgeometry/algorithms/convex_hull/
melkman.rs
