rgeometry/algorithms/convex_hull/
graham_scan.rs

1use crate::data::{Point, Polygon, PolygonConvex};
2use crate::{Error, Orientation, PolygonScalar, TotalOrd};
3
4// https://en.wikipedia.org/wiki/Graham_scan
5
6// Doesn't allocate.
7// Properties:
8//    No panics.
9//    All Ok results are valid convex polygons.
10//    No points are outside the resulting convex polygon.
11/// Convex hull of a set of points.
12///
13/// [Graham scan][wiki] algorithm for finding the smallest convex polygon which
14/// contains all the given points.
15///
16/// # Errors
17/// Will return an error iff the input set contains less than three distinct points.
18///
19/// # Properties
20/// * No points from the input set will be outside the returned convex polygon.
21/// * All vertices in the convex polygon are from the input set.
22///
23/// # Time complexity
24/// $O(n \log n)$
25///
26/// # Examples
27///
28/// ```rust
29/// # pub fn main() {
30/// # use rgeometry::algorithms::convex_hull;
31/// # use rgeometry::data::Point;
32/// # use rgeometry::Error;
33/// let empty_set: Vec<Point<i32,2>> = vec![];
34/// assert_eq!(
35///   convex_hull(empty_set).err(),
36///   Some(Error::InsufficientVertices))
37/// # }
38/// ```
39///
40/// ```rust
41/// # pub fn main() {
42/// # use rgeometry::algorithms::convex_hull;
43/// # use rgeometry::data::Point;
44/// # use rgeometry::Error;
45/// let dups = vec![Point::new([0,0])].repeat(3);
46/// assert_eq!(
47///   convex_hull(dups).err(),
48///   Some(Error::InsufficientVertices))
49/// # }
50/// ```
51///
52/// [wiki]: https://en.wikipedia.org/wiki/Graham_scan
53// ```no_run
54// # pub fn main() {
55// #   use rgeometry::algorithms::convex_hull;
56// #   use rgeometry_wasm::playground::*;
57// #
58// #   static START: std::sync::Once = std::sync::Once::new();
59// #   START.call_once(|| on_mousemove(|_event| main()));
60// #
61// #   clear_screen();
62// #   set_viewport(2., 2.);
63// #
64// #   let pts = with_points(7);
65// #   let points = pts.clone();
66// if let Ok(convex) = convex_hull(points) {
67//   render_polygon(&convex);
68// #   context().set_fill_style(&"grey".into());
69// #   context().fill();
70// }
71// #   for pt in &pts {
72// #     render_point(&pt);
73// #   }
74// # }
75// ```
76//
77// <iframe src="https://web.rgeometry.org/wasm/gist/eac484cd855d001815d23a053919b5ca"></iframe>
78//
79pub fn convex_hull<T>(mut pts: Vec<Point<T>>) -> Result<PolygonConvex<T>, Error>
80where
81  T: PolygonScalar,
82{
83  let smallest = &smallest_point(&pts)?;
84
85  pts.sort_unstable_by(|a, b| {
86    smallest
87      .ccw_cmp_around(a, b)
88      .then_with(|| (a.y_coord(), a.x_coord()).total_cmp(&(b.y_coord(), b.x_coord())))
89    // .then_with(|| smallest.cmp_distance_to(a, b))
90  });
91  if pts.len() < 3 {
92    return Err(Error::InsufficientVertices);
93  }
94  debug_assert_eq!(&pts[0], smallest);
95  let mut write_idx = 1;
96  let mut read_idx = 2;
97  // Drop points that are co-linear with our origin.
98  {
99    let origin = pts[write_idx - 1].clone();
100    while read_idx < pts.len() {
101      if Point::orient(&origin, &pts[write_idx], &pts[read_idx]) == Orientation::CoLinear {
102        pts.swap(read_idx, write_idx);
103        read_idx += 1;
104      } else {
105        break;
106      }
107    }
108  }
109  // Filter out points until all consecutive points are oriented counter-clockwise.
110  while read_idx < pts.len() {
111    let p1 = &pts[read_idx];
112    let p2 = &pts[write_idx];
113    let p3 = &pts[write_idx - 1];
114    match Point::orient(p3, p2, p1) {
115      Orientation::CounterClockWise => {
116        pts.swap(read_idx, write_idx + 1);
117        read_idx += 1;
118        write_idx += 1;
119      }
120      Orientation::ClockWise | Orientation::CoLinear => {
121        write_idx -= 1;
122      }
123    }
124  }
125  pts.truncate(write_idx + 1);
126  if pts.len() < 3 {
127    return Err(Error::InsufficientVertices);
128  }
129  Ok(PolygonConvex::new_unchecked(Polygon::new_unchecked(pts)))
130}
131
132// Find the smallest point and remove it from the vector
133// O(n)
134fn smallest_point<T>(pts: &[Point<T>]) -> Result<Point<T>, Error>
135where
136  T: PolygonScalar,
137{
138  Ok(
139    pts
140      .iter()
141      .min_by(|a, b| TotalOrd::total_cmp(&(a.y_coord(), a.x_coord()), &(b.y_coord(), b.x_coord())))
142      .ok_or(Error::InsufficientVertices)?
143      .clone(),
144  )
145}
146
147#[cfg(test)]
148#[cfg(not(tarpaulin_include))]
149mod tests {
150  use super::*;
151  use crate::data::PointLocation;
152  use crate::testing::*;
153
154  use claims::assert_ok;
155  use num_bigint::BigInt;
156
157  use proptest::collection::*;
158  use proptest::prelude::*;
159  use test_strategy::proptest;
160
161  #[test]
162  fn convex_hull_colinear() {
163    let points = vec![
164      Point::new([0, 0]),
165      Point::new([1, 0]),
166      Point::new([2, 0]),
167      Point::new([3, 0]),
168      Point::new([4, 0]),
169      Point::new([1, 1]),
170    ];
171    let poly = convex_hull(points).unwrap();
172    assert_ok!(poly.validate());
173  }
174
175  #[test]
176  fn convex_hull_colinear_rev() {
177    let points = vec![
178      Point::new([0, 0]),
179      Point::new([1, 0]),
180      Point::new([0, 9]),
181      Point::new([0, 8]),
182      Point::new([0, 7]),
183      Point::new([0, 6]),
184    ];
185    let poly = convex_hull(points).unwrap();
186    assert_ok!(poly.validate());
187  }
188
189  #[test]
190  fn convex_hull_dups() {
191    let points = vec![
192      Point::new([0, 0]),
193      Point::new([1, 0]),
194      Point::new([0, 0]),
195      Point::new([1, 0]),
196      Point::new([2, 2]),
197      Point::new([2, 2]),
198      Point::new([5, 1]),
199      Point::new([5, 1]),
200    ];
201    let poly = convex_hull(points).unwrap();
202    assert_ok!(poly.validate());
203  }
204
205  #[test]
206  fn convex_hull_insufficient_dups() {
207    let points = vec![
208      Point::new([0, 0]),
209      Point::new([0, 0]),
210      Point::new([2, 2]),
211      Point::new([2, 2]),
212      Point::new([0, 0]),
213      Point::new([2, 2]),
214    ];
215    assert_eq!(convex_hull(points).err(), Some(Error::InsufficientVertices));
216  }
217
218  #[test]
219  fn convex_hull_invalid() {
220    let points: Vec<Point<i64>> = vec![
221      Point { array: [0, 0] },
222      Point { array: [100, 0] },
223      Point { array: [50, 1] },
224      Point { array: [40, 1] },
225      Point { array: [0, 100] },
226    ];
227    let points: Vec<Point<BigInt, 2>> = points.into_iter().map(|pt| pt.cast()).collect();
228    let poly = convex_hull(points).unwrap();
229    assert_ok!(poly.validate());
230  }
231
232  #[test]
233  fn unit_1() {
234    let points: Vec<Point<BigInt>> = vec![
235      Point::new([0, 0]).into(),
236      Point::new([-1, 1]).into(),
237      Point::new([0, 1]).into(),
238      Point::new([-717193444810564826, 1]).into(),
239    ];
240    let poly = convex_hull(points).unwrap();
241    assert_ok!(poly.validate());
242  }
243
244  #[test]
245  fn unit_2() {
246    let points: Vec<Point<i8>> = vec![
247      Point::new([0, 0]),
248      Point::new([0, -10]),
249      Point::new([-13, 0]),
250    ];
251    let poly = convex_hull(points).unwrap();
252    assert_ok!(poly.validate());
253  }
254
255  #[proptest]
256  fn convex_hull_prop(#[strategy(vec(any_r(), 0..100))] pts: Vec<Point<BigInt>>) {
257    if let Ok(poly) = convex_hull(pts.clone()) {
258      // Prop #1: Results are valid.
259      prop_assert_eq!(poly.validate().err(), None);
260      // Prop #2: No points from the input set are outside the polygon.
261      for pt in pts.iter() {
262        prop_assert_ne!(poly.locate(pt), PointLocation::Outside)
263      }
264      // Prop #3: All vertices are in the input set.
265      for pt in poly.iter() {
266        prop_assert!(pts.contains(pt))
267      }
268    }
269  }
270
271  #[proptest]
272  fn convex_hull_prop_i8(#[strategy(vec(any::<Point<i8>>(), 0..100))] pts: Vec<Point<i8>>) {
273    if let Ok(poly) = convex_hull(pts.clone()) {
274      // Prop #1: Results are valid.
275      prop_assert_eq!(poly.validate().err(), None);
276      // Prop #2: No points from the input set are outside the polygon.
277      for pt in pts.iter() {
278        prop_assert_ne!(poly.locate(pt), PointLocation::Outside)
279      }
280      // Prop #3: All vertices are in the input set.
281      for pt in poly.iter() {
282        prop_assert!(pts.contains(pt))
283      }
284    }
285  }
286}
rgeometry/algorithms/convex_hull/graham_scan.rs

rgeometry/algorithms/convex_hull/
graham_scan.rs
