geo/algorithm/convex_hull/
graham.rs

1use super::{swap_remove_to_first, trivial_hull};
2use crate::kernels::*;
3use crate::{Coord, GeoNum, LineString};
4
5/// The [Graham's scan] algorithm to compute the convex hull
6/// of a collection of points. This algorithm is less
7/// performant than the quick hull, but allows computing all
8/// the points on the convex hull, as opposed to a strict
9/// convex hull that does not include collinear points.
10///
11/// # References
12///
13/// Graham, R.L. (1972). ["An Efficient Algorithm for
14/// Determining the Convex Hull of a Finite Planar
15/// Set"](http://www.math.ucsd.edu/~ronspubs/72_10_convex_hull.pdf)
16/// (PDF). \
17/// Information Processing Letters. 1 (4): 132–133.
18/// [doi:10.1016/0020-0190(72)90045-2](https://doi.org/10.1016%2F0020-0190%2872%2990045-2)
19///
20/// [Graham's scan]: //en.wikipedia.org/wiki/Graham_scan
21pub fn graham_hull<T>(mut points: &mut [Coord<T>], include_on_hull: bool) -> LineString<T>
22where
23    T: GeoNum,
24{
25    if points.len() < 4 {
26        // Nothing to build with fewer than four points.
27        return trivial_hull(points, include_on_hull);
28    }
29
30    // Allocate output vector
31    let mut output = Vec::with_capacity(points.len());
32
33    // Find lexicographically least point and add to hull
34    use crate::utils::least_index;
35    use std::cmp::Ordering;
36    let min_idx = least_index(points);
37    let head = swap_remove_to_first(&mut points, min_idx);
38    output.push(*head);
39
40    // Sort rest of the points by angle it makes with head
41    // point. If two points are collinear with head, we sort
42    // by distance. We use kernel predicates here.
43    let cmp = |q: &Coord<T>, r: &Coord<T>| match T::Ker::orient2d(*q, *head, *r) {
44        Orientation::CounterClockwise => Ordering::Greater,
45        Orientation::Clockwise => Ordering::Less,
46        Orientation::Collinear => {
47            let dist1 = T::Ker::square_euclidean_distance(*head, *q);
48            let dist2 = T::Ker::square_euclidean_distance(*head, *r);
49            dist1.partial_cmp(&dist2).unwrap()
50        }
51    };
52    points.sort_unstable_by(cmp);
53
54    for pt in points.iter() {
55        while output.len() > 1 {
56            let len = output.len();
57            match T::Ker::orient2d(output[len - 2], output[len - 1], *pt) {
58                Orientation::CounterClockwise => {
59                    break;
60                }
61                Orientation::Clockwise => {
62                    output.pop();
63                }
64                Orientation::Collinear => {
65                    if include_on_hull {
66                        break;
67                    } else {
68                        output.pop();
69                    }
70                }
71            }
72        }
73        // Corner case: if the lex. least point added before
74        // this loop is repeated, then we should not end up
75        // adding it here (because output.len() == 1 in the
76        // first iteration)
77        if include_on_hull || pt != output.last().unwrap() {
78            output.push(*pt);
79        }
80    }
81
82    // Close and output the line string
83    let mut output = LineString::new(output);
84    output.close();
85    output
86}
87
88#[cfg(test)]
89mod test {
90    use super::*;
91    use crate::IsConvex;
92    fn test_convexity<T: GeoNum>(mut initial: Vec<Coord<T>>) {
93        let hull = graham_hull(&mut initial, false);
94        assert!(hull.is_strictly_ccw_convex());
95        let hull = graham_hull(&mut initial, true);
96        assert!(hull.is_ccw_convex());
97    }
98
99    #[test]
100    fn test_graham_hull_ccw() {
101        let initial = vec![
102            (1.0, 0.0),
103            (2.0, 1.0),
104            (1.75, 1.1),
105            (1.0, 2.0),
106            (0.0, 1.0),
107            (1.0, 0.0),
108        ];
109        let initial = initial.iter().map(|e| Coord::from((e.0, e.1))).collect();
110        test_convexity(initial);
111    }
112
113    #[test]
114    fn graham_hull_test1() {
115        let v: Vec<_> = vec![(0, 0), (4, 0), (4, 1), (1, 1), (1, 4), (0, 4), (0, 0)];
116        let initial = v.iter().map(|e| Coord::from((e.0, e.1))).collect();
117        test_convexity(initial);
118    }
119
120    #[test]
121    fn graham_hull_test2() {
122        let v = vec![
123            (0, 10),
124            (1, 1),
125            (10, 0),
126            (1, -1),
127            (0, -10),
128            (-1, -1),
129            (-10, 0),
130            (-1, 1),
131            (0, 10),
132        ];
133        let initial = v.iter().map(|e| Coord::from((e.0, e.1))).collect();
134        test_convexity(initial);
135    }
136
137    #[test]
138    fn graham_test_complex() {
139        test_convexity(geo_test_fixtures::poly1::<f64>().0);
140    }
141
142    #[test]
143    fn quick_hull_test_complex_2() {
144        test_convexity(geo_test_fixtures::poly2::<f64>().0);
145    }
146}
geo/algorithm/convex_hull/graham.rs

geo/algorithm/convex_hull/
graham.rs
