geo/algorithm/convex_hull/
qhull.rs

1use super::{swap_remove_to_first, trivial_hull};
2use crate::kernels::{HasKernel, Kernel, Orientation};
3use crate::utils::partition_slice;
4use crate::{coord, Coord, GeoNum, LineString};
5
6// Determines if `p_c` lies on the positive side of the
7// segment `p_a` to `p_b`. In other words, whether segment
8// `p_a` to `p_c` is a counter-clockwise rotation from the
9// segment. We use kernels to ensure this predicate is
10// exact.
11#[inline]
12fn is_ccw<T>(p_a: Coord<T>, p_b: Coord<T>, p_c: Coord<T>) -> bool
13where
14    T: GeoNum,
15{
16    let o = <T as HasKernel>::Ker::orient2d(p_a, p_b, p_c);
17    o == Orientation::CounterClockwise
18}
19
20// Adapted from https://web.archive.org/web/20180409175413/http://www.ahristov.com/tutorial/geometry-games/convex-hull.html
21pub fn quick_hull<T>(mut points: &mut [Coord<T>]) -> LineString<T>
22where
23    T: GeoNum,
24{
25    // can't build a hull from fewer than four points
26    if points.len() < 4 {
27        return trivial_hull(points, false);
28    }
29    let mut hull = vec![];
30
31    use crate::utils::least_and_greatest_index;
32    let (min, max) = {
33        let (min_idx, mut max_idx) = least_and_greatest_index(points);
34        let min = swap_remove_to_first(&mut points, min_idx);
35
36        // Two special cases to consider:
37        // (1) max_idx = 0, and got swapped
38        if max_idx == 0 {
39            max_idx = min_idx;
40        }
41
42        // (2) max_idx = min_idx: then any point could be
43        // chosen as max. But from case (1), it could now be
44        // 0, and we should not decrement it.
45        max_idx = max_idx.saturating_sub(1);
46
47        let max = swap_remove_to_first(&mut points, max_idx);
48        (min, max)
49    };
50
51    {
52        let (points, _) = partition_slice(points, |p| is_ccw(*max, *min, *p));
53        hull_set(*max, *min, points, &mut hull);
54    }
55    hull.push(*max);
56    let (points, _) = partition_slice(points, |p| is_ccw(*min, *max, *p));
57    hull_set(*min, *max, points, &mut hull);
58    hull.push(*min);
59    // close the polygon
60    let mut hull: LineString<_> = hull.into();
61    hull.close();
62    hull
63}
64
65// recursively calculate the convex hull of a subset of points
66fn hull_set<T>(p_a: Coord<T>, p_b: Coord<T>, mut set: &mut [Coord<T>], hull: &mut Vec<Coord<T>>)
67where
68    T: GeoNum,
69{
70    if set.is_empty() {
71        return;
72    }
73    if set.len() == 1 {
74        hull.push(set[0]);
75        return;
76    }
77
78    // Construct orthogonal vector to `p_b` - `p_a` We
79    // compute inner product of this with `v` - `p_a` to
80    // find the farthest point from the line segment a-b.
81    let p_orth = coord! {
82        x: p_a.y - p_b.y,
83        y: p_b.x - p_a.x,
84    };
85
86    let furthest_idx = set
87        .iter()
88        .map(|pt| {
89            let p_diff = coord! {
90                x: pt.x - p_a.x,
91                y: pt.y - p_a.y,
92            };
93            p_orth.x * p_diff.x + p_orth.y * p_diff.y
94        })
95        .enumerate()
96        .max_by(|(_, a), (_, b)| a.partial_cmp(b).unwrap())
97        .unwrap()
98        .0;
99
100    // move Coord at furthest_point from set into hull
101    let furthest_point = swap_remove_to_first(&mut set, furthest_idx);
102    // points over PB
103    {
104        let (points, _) = partition_slice(set, |p| is_ccw(*furthest_point, p_b, *p));
105        hull_set(*furthest_point, p_b, points, hull);
106    }
107    hull.push(*furthest_point);
108    // points over AP
109    let (points, _) = partition_slice(set, |p| is_ccw(p_a, *furthest_point, *p));
110    hull_set(p_a, *furthest_point, points, hull);
111}
112
113#[cfg(test)]
114mod test {
115    use super::*;
116    use crate::IsConvex;
117
118    #[test]
119    fn quick_hull_test1() {
120        let mut v = vec![
121            coord! { x: 0.0, y: 0.0 },
122            coord! { x: 4.0, y: 0.0 },
123            coord! { x: 4.0, y: 1.0 },
124            coord! { x: 1.0, y: 1.0 },
125            coord! { x: 1.0, y: 4.0 },
126            coord! { x: 0.0, y: 4.0 },
127            coord! { x: 0.0, y: 0.0 },
128        ];
129        let res = quick_hull(&mut v);
130        assert!(res.is_strictly_ccw_convex());
131    }
132
133    #[test]
134    fn quick_hull_test2() {
135        let mut v = vec![
136            coord! { x: 0, y: 10 },
137            coord! { x: 1, y: 1 },
138            coord! { x: 10, y: 0 },
139            coord! { x: 1, y: -1 },
140            coord! { x: 0, y: -10 },
141            coord! { x: -1, y: -1 },
142            coord! { x: -10, y: 0 },
143            coord! { x: -1, y: 1 },
144            coord! { x: 0, y: 10 },
145        ];
146        let correct = vec![
147            coord! { x: 0, y: -10 },
148            coord! { x: 10, y: 0 },
149            coord! { x: 0, y: 10 },
150            coord! { x: -10, y: 0 },
151            coord! { x: 0, y: -10 },
152        ];
153        let res = quick_hull(&mut v);
154        assert_eq!(res.0, correct);
155    }
156
157    #[test]
158    // test whether output is ccw
159    fn quick_hull_test_ccw() {
160        let initial = vec![
161            (1.0, 0.0),
162            (2.0, 1.0),
163            (1.75, 1.1),
164            (1.0, 2.0),
165            (0.0, 1.0),
166            (1.0, 0.0),
167        ];
168        let mut v: Vec<_> = initial.iter().map(|e| coord! { x: e.0, y: e.1 }).collect();
169        let correct = vec![(1.0, 0.0), (2.0, 1.0), (1.0, 2.0), (0.0, 1.0), (1.0, 0.0)];
170        let v_correct: Vec<_> = correct.iter().map(|e| coord! { x: e.0, y: e.1 }).collect();
171        let res = quick_hull(&mut v);
172        assert_eq!(res.0, v_correct);
173    }
174
175    #[test]
176    fn quick_hull_test_ccw_maintain() {
177        // initial input begins at min y, is oriented ccw
178        let initial = vec![
179            (0., 0.),
180            (2., 0.),
181            (2.5, 1.75),
182            (2.3, 1.7),
183            (1.75, 2.5),
184            (1.3, 2.),
185            (0., 2.),
186            (0., 0.),
187        ];
188        let mut v: Vec<_> = initial.iter().map(|e| coord! { x: e.0, y: e.1 }).collect();
189        let res = quick_hull(&mut v);
190        assert!(res.is_strictly_ccw_convex());
191    }
192
193    #[test]
194    fn quick_hull_test_complex() {
195        let mut coords = geo_test_fixtures::poly1::<f64>().0;
196        let correct = geo_test_fixtures::poly1_hull::<f64>().0;
197        let res = quick_hull(&mut coords);
198        assert_eq!(res.0, correct);
199    }
200
201    #[test]
202    fn quick_hull_test_complex_2() {
203        let mut coords = geo_test_fixtures::poly2::<f64>().0;
204        let correct = geo_test_fixtures::poly2_hull::<f64>().0;
205        let res = quick_hull(&mut coords);
206        assert_eq!(res.0, correct);
207    }
208
209    #[test]
210    fn quick_hull_test_collinear() {
211        // Initial input begins at min x, but not min y
212        // There are three points with same x.
213        // Output should not contain the middle point.
214        let initial = vec![
215            (-1., 0.),
216            (-1., -1.),
217            (-1., 1.),
218            (0., 0.),
219            (0., -1.),
220            (0., 1.),
221            (1., 0.),
222            (1., -1.),
223            (1., 1.),
224        ];
225        let mut v: Vec<_> = initial.iter().map(|e| coord! { x: e.0, y: e.1 }).collect();
226        let res = quick_hull(&mut v);
227        assert!(res.is_strictly_ccw_convex());
228    }
229}
geo/algorithm/convex_hull/qhull.rs

geo/algorithm/convex_hull/
qhull.rs
