imageproc/
geometry.rs

1//! Computational geometry functions, for example finding convex hulls.
2
3use crate::point::{distance, Line, Point, Rotation};
4use num::{cast, NumCast};
5use std::cmp::{Ord, Ordering};
6use std::f64::{self, consts::PI};
7
8/// Computes the length of an arc. If `closed` is set to `true` then the distance
9/// between the last and the first point is included in the total length.
10pub fn arc_length<T>(arc: &[Point<T>], closed: bool) -> f64
11where
12    T: NumCast + Copy,
13{
14    let mut length = arc.windows(2).map(|pts| distance(pts[0], pts[1])).sum();
15
16    if arc.len() > 2 && closed {
17        length += distance(arc[0], arc[arc.len() - 1]);
18    }
19
20    length
21}
22
23/// Approximates a polygon using the [Douglas–Peucker algorithm].
24///
25/// [Douglas–Peucker algorithm]: https://en.wikipedia.org/wiki/Ramer-Douglas-Peucker_algorithm
26pub fn approximate_polygon_dp<T>(curve: &[Point<T>], epsilon: f64, closed: bool) -> Vec<Point<T>>
27where
28    T: NumCast + Copy,
29{
30    if epsilon <= 0.0 {
31        panic!("epsilon must be greater than 0.0");
32    }
33
34    // Find the point with the maximum distance
35    let mut dmax = 0.0;
36    let mut index = 0;
37    let end = curve.len() - 1;
38    let line = Line::from_points(curve[0].to_f64(), curve[end].to_f64());
39    for (i, point) in curve.iter().enumerate().skip(1) {
40        let d = line.distance_from_point(point.to_f64());
41        if d > dmax {
42            index = i;
43            dmax = d;
44        }
45    }
46
47    let mut res = Vec::new();
48
49    // If max distance is greater than epsilon, recursively simplify
50    if dmax > epsilon {
51        // Recursive call
52        let mut partial1 = approximate_polygon_dp(&curve[0..=index], epsilon, false);
53        let mut partial2 = approximate_polygon_dp(&curve[index..=end], epsilon, false);
54
55        // Build the result list
56        partial1.pop();
57        res.append(&mut partial1);
58        res.append(&mut partial2);
59    } else {
60        res.push(curve[0]);
61        res.push(curve[end]);
62    }
63
64    if closed {
65        res.pop();
66    }
67
68    res
69}
70
71/// Finds the rectangle of least area that includes all input points. This rectangle need not be axis-aligned.
72///
73/// The returned points are the [top left, top right, bottom right, bottom left] points of this rectangle.
74pub fn min_area_rect<T>(points: &[Point<T>]) -> [Point<T>; 4]
75where
76    T: NumCast + Copy + Ord,
77{
78    let hull = convex_hull(points);
79    match hull.len() {
80        0 => panic!("no points are defined"),
81        1 => [hull[0]; 4],
82        2 => [hull[0], hull[1], hull[1], hull[0]],
83        _ => rotating_calipers(&hull),
84    }
85}
86
87/// An implementation of [rotating calipers] used for determining the
88/// bounding rectangle with the smallest area.
89///
90/// [rotating calipers]: https://en.wikipedia.org/wiki/Rotating_calipers
91fn rotating_calipers<T>(points: &[Point<T>]) -> [Point<T>; 4]
92where
93    T: NumCast + Copy,
94{
95    let mut edge_angles: Vec<f64> = points
96        .windows(2)
97        .map(|e| {
98            let edge = e[1].to_f64() - e[0].to_f64();
99            ((edge.y.atan2(edge.x) + PI) % (PI / 2.)).abs()
100        })
101        .collect();
102
103    edge_angles.dedup();
104
105    let mut min_area = f64::MAX;
106    let mut res = vec![Point::new(0.0, 0.0); 4];
107    for angle in edge_angles {
108        let rotation = Rotation::new(angle);
109        let rotated_points: Vec<Point<f64>> =
110            points.iter().map(|p| p.to_f64().rotate(rotation)).collect();
111
112        let (min_x, max_x, min_y, max_y) =
113            rotated_points
114                .iter()
115                .fold((f64::MAX, f64::MIN, f64::MAX, f64::MIN), |acc, p| {
116                    (
117                        acc.0.min(p.x),
118                        acc.1.max(p.x),
119                        acc.2.min(p.y),
120                        acc.3.max(p.y),
121                    )
122                });
123
124        let area = (max_x - min_x) * (max_y - min_y);
125        if area < min_area {
126            min_area = area;
127            res[0] = Point::new(max_x, min_y).invert_rotation(rotation);
128            res[1] = Point::new(min_x, min_y).invert_rotation(rotation);
129            res[2] = Point::new(min_x, max_y).invert_rotation(rotation);
130            res[3] = Point::new(max_x, max_y).invert_rotation(rotation);
131        }
132    }
133
134    res.sort_by(|a, b| a.x.partial_cmp(&b.x).unwrap());
135
136    let i1 = if res[1].y > res[0].y { 0 } else { 1 };
137    let i2 = if res[3].y > res[2].y { 2 } else { 3 };
138    let i3 = if res[3].y > res[2].y { 3 } else { 2 };
139    let i4 = if res[1].y > res[0].y { 1 } else { 0 };
140
141    [
142        Point::new(
143            cast(res[i1].x.floor()).unwrap(),
144            cast(res[i1].y.floor()).unwrap(),
145        ),
146        Point::new(
147            cast(res[i2].x.ceil()).unwrap(),
148            cast(res[i2].y.floor()).unwrap(),
149        ),
150        Point::new(
151            cast(res[i3].x.ceil()).unwrap(),
152            cast(res[i3].y.ceil()).unwrap(),
153        ),
154        Point::new(
155            cast(res[i4].x.floor()).unwrap(),
156            cast(res[i4].y.ceil()).unwrap(),
157        ),
158    ]
159}
160
161/// Finds the convex hull of a set of points, using the [Graham scan algorithm].
162///
163/// [Graham scan algorithm]: https://en.wikipedia.org/wiki/Graham_scan
164pub fn convex_hull<T>(points_slice: &[Point<T>]) -> Vec<Point<T>>
165where
166    T: NumCast + Copy + Ord,
167{
168    if points_slice.is_empty() {
169        return Vec::new();
170    }
171    let mut points: Vec<Point<T>> = points_slice.to_vec();
172    let mut start_point_pos = 0;
173    let mut start_point = points[0];
174    for (i, &point) in points.iter().enumerate().skip(1) {
175        if point.y < start_point.y || point.y == start_point.y && point.x < start_point.x {
176            start_point_pos = i;
177            start_point = point;
178        }
179    }
180    points.swap(0, start_point_pos);
181    points.remove(0);
182    points.sort_by(
183        |a, b| match orientation(start_point.to_i32(), a.to_i32(), b.to_i32()) {
184            Orientation::Collinear => {
185                if distance(start_point, *a) < distance(start_point, *b) {
186                    Ordering::Less
187                } else {
188                    Ordering::Greater
189                }
190            }
191            Orientation::Clockwise => Ordering::Greater,
192            Orientation::CounterClockwise => Ordering::Less,
193        },
194    );
195
196    let mut iter = points.iter().peekable();
197    let mut remaining_points = Vec::with_capacity(points.len());
198    while let Some(mut p) = iter.next() {
199        while iter.peek().is_some()
200            && orientation(
201                start_point.to_i32(),
202                p.to_i32(),
203                iter.peek().unwrap().to_i32(),
204            ) == Orientation::Collinear
205        {
206            p = iter.next().unwrap();
207        }
208        remaining_points.push(p);
209    }
210
211    let mut stack: Vec<Point<T>> = vec![Point::new(
212        cast(start_point.x).unwrap(),
213        cast(start_point.y).unwrap(),
214    )];
215
216    for p in points {
217        while stack.len() > 1
218            && orientation(
219                stack[stack.len() - 2].to_i32(),
220                stack[stack.len() - 1].to_i32(),
221                p.to_i32(),
222            ) != Orientation::CounterClockwise
223        {
224            stack.pop();
225        }
226        stack.push(p);
227    }
228    stack
229}
230
231#[derive(Debug, Copy, Clone, PartialEq, Eq)]
232enum Orientation {
233    Collinear,
234    Clockwise,
235    CounterClockwise,
236}
237
238/// Determines whether p -> q -> r is a left turn, a right turn, or the points are collinear.
239fn orientation(p: Point<i32>, q: Point<i32>, r: Point<i32>) -> Orientation {
240    let val = (q.y - p.y) * (r.x - q.x) - (q.x - p.x) * (r.y - q.y);
241    match val.cmp(&0) {
242        Ordering::Equal => Orientation::Collinear,
243        Ordering::Greater => Orientation::Clockwise,
244        Ordering::Less => Orientation::CounterClockwise,
245    }
246}
247
248#[cfg(test)]
249mod tests {
250    use super::*;
251    use crate::point::Point;
252
253    #[test]
254    fn test_arc_length() {
255        assert_eq!(arc_length::<f64>(&[], false), 0.0);
256        assert_eq!(arc_length(&[Point::new(1.0, 1.0)], false), 0.0);
257        assert_eq!(
258            arc_length(&[Point::new(1.0, 1.0), Point::new(4.0, 5.0)], false),
259            5.0
260        );
261        assert_eq!(
262            arc_length(
263                &[
264                    Point::new(1.0, 1.0),
265                    Point::new(4.0, 5.0),
266                    Point::new(9.0, 17.0)
267                ],
268                false
269            ),
270            18.0
271        );
272        assert_eq!(
273            arc_length(
274                &[
275                    Point::new(1.0, 1.0),
276                    Point::new(4.0, 5.0),
277                    Point::new(9.0, 17.0)
278                ],
279                true
280            ),
281            18.0 + (8f64.powf(2.0) + 16f64.powf(2.0)).sqrt()
282        );
283    }
284
285    #[test]
286    fn convex_hull_points() {
287        let star = vec![
288            Point::new(100, 20),
289            Point::new(90, 35),
290            Point::new(60, 25),
291            Point::new(90, 40),
292            Point::new(80, 55),
293            Point::new(101, 50),
294            Point::new(130, 60),
295            Point::new(115, 45),
296            Point::new(140, 30),
297            Point::new(120, 35),
298        ];
299        let points = convex_hull(&star);
300        assert_eq!(
301            points,
302            [
303                Point::new(100, 20),
304                Point::new(140, 30),
305                Point::new(130, 60),
306                Point::new(80, 55),
307                Point::new(60, 25)
308            ]
309        );
310    }
311
312    #[test]
313    fn convex_hull_points_empty_vec() {
314        let points = convex_hull::<i32>(&vec![]);
315        assert_eq!(points, []);
316    }
317
318    #[test]
319    fn convex_hull_points_with_negative_values() {
320        let star = vec![
321            Point::new(100, -20),
322            Point::new(90, 5),
323            Point::new(60, -15),
324            Point::new(90, 0),
325            Point::new(80, 15),
326            Point::new(101, 10),
327            Point::new(130, 20),
328            Point::new(115, 5),
329            Point::new(140, -10),
330            Point::new(120, -5),
331        ];
332        let points = convex_hull(&star);
333        assert_eq!(
334            points,
335            [
336                Point::new(100, -20),
337                Point::new(140, -10),
338                Point::new(130, 20),
339                Point::new(80, 15),
340                Point::new(60, -15)
341            ]
342        );
343    }
344
345    #[test]
346    fn test_min_area() {
347        assert_eq!(
348            min_area_rect(&[
349                Point::new(100, 20),
350                Point::new(140, 30),
351                Point::new(130, 60),
352                Point::new(80, 55),
353                Point::new(60, 25)
354            ]),
355            [
356                Point::new(60, 16),
357                Point::new(141, 24),
358                Point::new(137, 61),
359                Point::new(57, 53)
360            ]
361        )
362    }
363}
imageproc/geometry.rs

imageproc/
geometry.rs
