shape_convex/convex2d/
solver.rs

1use super::*;
2
3impl<T> FastConvexHull<T> {
4    pub fn new(points: &[Point<T>]) -> Self where T: Real {
5        // TODO: Optimize using divide and conquer
6        graham_scan(points)
7    }
8}
9
10fn sort_by_min_angle<F: Real>(set: &[Point<F>], min: &Point<F>) -> Vec<Point<F>> {
11    let mut points: Vec<(F, F, Point<F>)> = set
12        .iter()
13        .map(|x| {
14            (
15                (x.y - min.y).atan2(x.x - min.x),
16                // angle
17                (x.y - min.y).hypot(x.x - min.x),
18                // distance (we want the closest to be first)
19                *x,
20            )
21        })
22        .collect();
23    points.sort_by(|(a1, d1, _), (a2, d2, _)| (a1, d1).partial_cmp(&(a2, d2)).unwrap_or(Equal));
24    points.into_iter().map(|x| x.2).collect()
25}
26
27fn z_vector_product<F: Real>(a: &Point<F>, b: &Point<F>, c: &Point<F>) -> F {
28    (b.x - a.x) * (c.y - a.y) - (c.x - a.x) * (b.y - a.y)
29}
30
31fn graham_scan<F: Real>(set: &[Point<F>]) -> FastConvexHull<F> {
32    if set.is_empty() {
33        return FastConvexHull {
34            bounds: vec![],
35            inners: vec![],
36        };
37    }
38
39    let mut inner = vec![];
40    let mut stack: Vec<Point<F>> = vec![];
41    let min = set
42        .iter()
43        .min_by(|a, b| {
44            let ord = a.y.partial_cmp(&b.y).unwrap_or(Equal);
45            match ord {
46                Equal => a.x.partial_cmp(&b.x).unwrap_or(Equal),
47                o => o,
48            }
49        })
50        .unwrap();
51    let points = sort_by_min_angle(set, min);
52
53    if points.len() <= 3 {
54        return FastConvexHull {
55            bounds: points,
56            inners: inner,
57        };
58    }
59
60    for point in points {
61        while stack.len() > 1
62            && z_vector_product(&stack[stack.len() - 2], &stack[stack.len() - 1], &point) <= F::zero()
63        {
64            // SAFETY: we know that stack.len() > 1
65            unsafe {
66                inner.push(stack.pop().unwrap_unchecked());
67            }
68        }
69        stack.push(point);
70    }
71
72    FastConvexHull {
73        bounds: stack,
74        inners: inner,
75    }
76}
77
78impl<T: Real> AddAssign<&[Point<T>]> for FastConvexHull<T> {
79    fn add_assign(&mut self, rhs: &[Point<T>]) {
80        self.bounds.extend_from_slice(rhs);
81        let Self { bounds, inners } = graham_scan(&self.bounds);
82        self.bounds = bounds;
83        self.inners.extend_from_slice(&inners);
84    }
85}
86
87impl<T: Real> AddAssign<&FastConvexHull<T>> for FastConvexHull<T> {
88    fn add_assign(&mut self, rhs: &FastConvexHull<T>) {
89        self.add_assign(rhs.bounds.as_slice());
90        self.inners.extend_from_slice(&rhs.inners);
91    }
92}
93
94
95#[test]
96// from https://codegolf.stackexchange.com/questions/11035/find-the-convex-hull-of-a-set-of-2d-points
97fn lots_of_points() {
98    let list = vec![
99        (4.4, 14.),
100        (6.7, 15.25),
101        (6.9, 12.8),
102        (2.1, 11.1),
103        (9.5, 14.9),
104        (13.2, 11.9),
105        (10.3, 12.3),
106        (6.8, 9.5),
107        (3.3, 7.7),
108        (0.6, 5.1),
109        (5.3, 2.4),
110        (8.45, 4.7),
111        (11.5, 9.6),
112        (13.8, 7.3),
113        (12.9, 3.1),
114        (11., 1.1),
115    ];
116    let ans = vec![
117        (11., 1.1),
118        (12.9, 3.1),
119        (13.8, 7.3),
120        (13.2, 11.9),
121        (9.5, 14.9),
122        (6.7, 15.25),
123        (4.4, 14.),
124        (2.1, 11.1),
125        (0.6, 5.1),
126        (5.3, 2.4),
127    ];
128
129    let points = PointSet::from_iter(list).points;
130    let ans = PointSet::from_iter(ans).points;
131    let ch = graham_scan(&points);
132    println!("{:#?}", ch);
133    assert_eq!(ch.bounds, ans);
134}
shape_convex/convex2d/solver.rs

shape_convex/convex2d/
solver.rs
