jagua_rs/geometry/primitives/
simple_polygon.rs

1use std::borrow::Borrow;
2
3use itertools::Itertools;
4use ordered_float::{NotNan, OrderedFloat};
5
6use crate::geometry::Transformation;
7use crate::geometry::convex_hull::convex_hull_from_points;
8use crate::geometry::fail_fast::{SPSurrogate, SPSurrogateConfig, compute_pole};
9use crate::geometry::geo_enums::GeoPosition;
10use crate::geometry::geo_traits::{
11    CollidesWith, DistanceTo, SeparationDistance, Transformable, TransformableFrom,
12};
13use crate::geometry::primitives::Circle;
14use crate::geometry::primitives::Edge;
15use crate::geometry::primitives::Point;
16use crate::geometry::primitives::Rect;
17use crate::util::FPA;
18use anyhow::{Result, bail};
19
20/// A Simple Polygon is a polygon that does not intersect itself and contains no holes.
21/// It is a closed shape with a finite number of vertices and edges.
22/// [read more](https://en.wikipedia.org/wiki/Simple_polygon)
23#[derive(Clone, Debug)]
24pub struct SPolygon {
25    /// Set of points that form the polygon
26    pub vertices: Vec<Point>,
27    /// Bounding box
28    pub bbox: Rect,
29    /// Area of its interior
30    pub area: f32,
31    /// Maximum distance between any two points in the polygon
32    pub diameter: f32,
33    /// [Pole of inaccessibility](https://en.wikipedia.org/wiki/Pole_of_inaccessibility) represented as a circle
34    pub poi: Circle,
35    /// Optional surrogate representation of the polygon (subset of the original)
36    pub surrogate: Option<SPSurrogate>,
37}
38
39impl SPolygon {
40    /// Create a new simple polygon from a set of points, expensive operations are performed here! Use [Self::clone()] or [Self::transform()] to avoid recomputation.
41    pub fn new(mut points: Vec<Point>) -> Result<Self> {
42        if points.len() < 3 {
43            bail!("Simple polygon must have at least 3 points: {points:?}");
44        }
45        if points.iter().unique().count() != points.len() {
46            bail!("Simple polygon should not contain duplicate points: {points:?}");
47        }
48
49        let area = match SPolygon::calculate_area(&points) {
50            0.0 => bail!("Simple polygon has no area: {points:?}"),
51            area if area < 0.0 => {
52                //edges should always be ordered counterclockwise (positive area)
53                points.reverse();
54                -area
55            }
56            area => area,
57        };
58
59        let diameter = SPolygon::calculate_diameter(points.clone());
60        let bbox = SPolygon::generate_bounding_box(&points);
61        let poi = SPolygon::calculate_poi(&points, diameter)?;
62
63        Ok(SPolygon {
64            vertices: points,
65            bbox,
66            area,
67            diameter,
68            poi,
69            surrogate: None,
70        })
71    }
72
73    pub fn generate_surrogate(&mut self, config: SPSurrogateConfig) -> Result<()> {
74        //regenerate the surrogate if it is not present or if the config has changed
75        match &self.surrogate {
76            Some(surrogate) if surrogate.config == config => {}
77            _ => self.surrogate = Some(SPSurrogate::new(self, config)?),
78        }
79        Ok(())
80    }
81
82    pub fn vertex(&self, i: usize) -> Point {
83        self.vertices[i]
84    }
85
86    pub fn edge(&self, i: usize) -> Edge {
87        assert!(i < self.n_vertices(), "index out of bounds");
88        let j = if i == self.n_vertices() - 1 { 0 } else { i + 1 };
89        Edge {
90            start: self.vertices[i],
91            end: self.vertices[j],
92        }
93    }
94
95    pub fn edge_iter(&self) -> impl Iterator<Item = Edge> + '_ {
96        (0..self.n_vertices()).map(move |i| self.edge(i))
97    }
98
99    pub fn n_vertices(&self) -> usize {
100        self.vertices.len()
101    }
102
103    pub fn surrogate(&self) -> &SPSurrogate {
104        self.surrogate.as_ref().expect("surrogate not generated")
105    }
106
107    pub fn calculate_diameter(points: Vec<Point>) -> f32 {
108        //The two points furthest apart must be part of the convex hull
109        let ch = convex_hull_from_points(points);
110
111        //go through all pairs of points and find the pair with the largest distance
112        let sq_diam = ch
113            .iter()
114            .tuple_combinations()
115            .map(|(p1, p2)| p1.sq_distance_to(p2))
116            .max_by_key(|sq_d| NotNan::new(*sq_d).unwrap())
117            .expect("convex hull is empty");
118
119        sq_diam.sqrt()
120    }
121
122    pub fn generate_bounding_box(points: &[Point]) -> Rect {
123        let (mut x_min, mut y_min) = (f32::MAX, f32::MAX);
124        let (mut x_max, mut y_max) = (f32::MIN, f32::MIN);
125
126        for point in points.iter() {
127            x_min = x_min.min(point.0);
128            y_min = y_min.min(point.1);
129            x_max = x_max.max(point.0);
130            y_max = y_max.max(point.1);
131        }
132        Rect::try_new(x_min, y_min, x_max, y_max).unwrap()
133    }
134
135    //https://en.wikipedia.org/wiki/Shoelace_formula
136    //counterclockwise = positive area, clockwise = negative area
137    pub fn calculate_area(points: &[Point]) -> f32 {
138        let mut sigma: f32 = 0.0;
139        for i in 0..points.len() {
140            //next point
141            let j = (i + 1) % points.len();
142
143            let (x_i, y_i) = points[i].into();
144            let (x_j, y_j) = points[j].into();
145
146            sigma += (y_i + y_j) * (x_i - x_j)
147        }
148
149        0.5 * sigma
150    }
151
152    pub fn calculate_poi(points: &[Point], diameter: f32) -> Result<Circle> {
153        //need to make a dummy simple polygon, because the pole generation algorithm
154        //relies on many of the methods provided by the simple polygon struct
155        let dummy_sp = {
156            let bbox = SPolygon::generate_bounding_box(points);
157            let area = SPolygon::calculate_area(points);
158            let dummy_poi = Circle::try_new(Point(f32::MAX, f32::MAX), f32::MAX).unwrap();
159
160            SPolygon {
161                vertices: points.to_vec(),
162                bbox,
163                area,
164                diameter,
165                poi: dummy_poi,
166                surrogate: None,
167            }
168        };
169
170        compute_pole(&dummy_sp, &[])
171    }
172
173    pub fn centroid(&self) -> Point {
174        //based on: https://en.wikipedia.org/wiki/Centroid#Of_a_polygon
175
176        let area = self.area;
177        let mut c_x = 0.0;
178        let mut c_y = 0.0;
179
180        for i in 0..self.n_vertices() {
181            let j = if i == self.n_vertices() - 1 { 0 } else { i + 1 };
182            let Point(x_i, y_i) = self.vertex(i);
183            let Point(x_j, y_j) = self.vertex(j);
184            c_x += (x_i + x_j) * (x_i * y_j - x_j * y_i);
185            c_y += (y_i + y_j) * (x_i * y_j - x_j * y_i);
186        }
187
188        c_x /= 6.0 * area;
189        c_y /= 6.0 * area;
190
191        (c_x, c_y).into()
192    }
193}
194
195impl Transformable for SPolygon {
196    fn transform(&mut self, t: &Transformation) -> &mut Self {
197        //destructuring pattern to ensure that the code is updated when the struct changes
198        let SPolygon {
199            vertices: points,
200            bbox,
201            area: _,
202            diameter: _,
203            poi,
204            surrogate,
205        } = self;
206
207        //transform all points of the simple poly
208        points.iter_mut().for_each(|p| {
209            p.transform(t);
210        });
211
212        poi.transform(t);
213
214        //transform the surrogate
215        if let Some(surrogate) = surrogate.as_mut() {
216            surrogate.transform(t);
217        }
218
219        //regenerate bounding box
220        *bbox = SPolygon::generate_bounding_box(points);
221
222        self
223    }
224}
225
226impl TransformableFrom for SPolygon {
227    fn transform_from(&mut self, reference: &Self, t: &Transformation) -> &mut Self {
228        //destructuring pattern to ensure that the code is updated when the struct changes
229        let SPolygon {
230            vertices: points,
231            bbox,
232            area: _,
233            diameter: _,
234            poi,
235            surrogate,
236        } = self;
237
238        for (p, ref_p) in points.iter_mut().zip(&reference.vertices) {
239            p.transform_from(ref_p, t);
240        }
241
242        poi.transform_from(&reference.poi, t);
243
244        //transform the surrogate
245        if let Some(surrogate) = surrogate.as_mut() {
246            surrogate.transform_from(reference.surrogate(), t);
247        }
248        //regenerate bounding box
249        *bbox = SPolygon::generate_bounding_box(points);
250
251        self
252    }
253}
254
255impl CollidesWith<Point> for SPolygon {
256    fn collides_with(&self, point: &Point) -> bool {
257        //based on the ray casting algorithm: https://en.wikipedia.org/wiki/Point_in_polygon#Ray_casting_algorithm
258        match self.bbox.collides_with(point) {
259            false => false,
260            true => {
261                //horizontal ray shot to the right.
262                //Starting from the point to another point that is certainly outside the shape
263                let point_outside = Point(self.bbox.x_max + self.bbox.width(), point.1);
264                let ray = Edge {
265                    start: *point,
266                    end: point_outside,
267                };
268
269                let mut n_intersections = 0;
270                for edge in self.edge_iter() {
271                    //Check if the ray does not go through (or almost through) a vertex
272                    //This can result in funky behaviour, which could incorrect results
273                    //Therefore we handle this case
274                    let (s_x, s_y) = (FPA(edge.start.0), FPA(edge.start.1));
275                    let (e_x, e_y) = (FPA(edge.end.0), FPA(edge.end.1));
276                    let (p_x, p_y) = (FPA(point.0), FPA(point.1));
277
278                    if (s_y == p_y && s_x > p_x) || (e_y == p_y && e_x > p_x) {
279                        //in this case, the ray passes through (or dangerously close to) a vertex
280                        //We handle this case by only counting an intersection if the edge is below the ray
281                        if s_y < p_y || e_y < p_y {
282                            n_intersections += 1;
283                        }
284                    } else if ray.collides_with(&edge) {
285                        n_intersections += 1;
286                    }
287                }
288                n_intersections % 2 == 1
289            }
290        }
291    }
292}
293
294impl DistanceTo<Point> for SPolygon {
295    fn distance_to(&self, point: &Point) -> f32 {
296        self.sq_distance_to(point).sqrt()
297    }
298    fn sq_distance_to(&self, point: &Point) -> f32 {
299        match self.collides_with(point) {
300            true => 0.0,
301            false => self
302                .edge_iter()
303                .map(|edge| edge.sq_distance_to(point))
304                .min_by(|a, b| a.partial_cmp(b).unwrap())
305                .unwrap(),
306        }
307    }
308}
309
310impl SeparationDistance<Point> for SPolygon {
311    fn separation_distance(&self, point: &Point) -> (GeoPosition, f32) {
312        let (position, sq_distance) = self.sq_separation_distance(point);
313        (position, sq_distance.sqrt())
314    }
315
316    fn sq_separation_distance(&self, point: &Point) -> (GeoPosition, f32) {
317        let distance_to_closest_edge = self
318            .edge_iter()
319            .map(|edge| edge.sq_distance_to(point))
320            .min_by_key(|sq_d| OrderedFloat(*sq_d))
321            .unwrap();
322
323        match self.collides_with(point) {
324            true => (GeoPosition::Interior, distance_to_closest_edge),
325            false => (GeoPosition::Exterior, distance_to_closest_edge),
326        }
327    }
328}
329
330impl<T> From<T> for SPolygon
331where
332    T: Borrow<Rect>,
333{
334    fn from(r: T) -> Self {
335        let r = r.borrow();
336        SPolygon::new(vec![
337            (r.x_min, r.y_min).into(),
338            (r.x_max, r.y_min).into(),
339            (r.x_max, r.y_max).into(),
340            (r.x_min, r.y_max).into(),
341        ])
342        .unwrap()
343    }
344}
jagua_rs/geometry/primitives/simple_polygon.rs

jagua_rs/geometry/primitives/
simple_polygon.rs
