gistools/geometry/tools/polys/
inside.rs

1use crate::geometry::orient2d;
2use alloc::vec::Vec;
3use s2json::{
4    Feature, Geometry, GetXY, MultiLineString, MultiLineString3D, MultiLineString3DGeometry,
5    MultiLineStringGeometry, MultiPoint, MultiPoint3D, MultiPoint3DGeometry, MultiPointGeometry,
6    MultiPolygon, MultiPolygon3D, MultiPolygon3DGeometry, MultiPolygonGeometry, Point, Point3D,
7    Point3DGeometry, PointGeometry, VectorFeature, VectorGeometry, VectorMultiLineString,
8    VectorMultiLineStringGeometry, VectorMultiPoint, VectorMultiPointGeometry, VectorMultiPolygon,
9    VectorMultiPolygonGeometry, VectorPoint, VectorPointGeometry,
10};
11
12/// The result of a point being inside, outside, or on the boundary
13#[derive(Debug, Clone, Copy, PartialEq)]
14pub enum InsideResult {
15    /// The point is inside
16    Inside,
17    /// The point is outside
18    Outside,
19    /// The point is on the boundary
20    Boundary,
21}
22impl InsideResult {
23    /// Get the dominant result. Inside > Boundary > Outside
24    ///
25    /// Inside takes precedence over Boundary and Outside
26    ///
27    /// Boundary takes precedence over Outside
28    pub fn get_dominant(&self, other: InsideResult) -> InsideResult {
29        if *self == other {
30            *self
31        } else {
32            match (self, other) {
33                (InsideResult::Inside, _) => InsideResult::Inside,
34                (_, InsideResult::Inside) => InsideResult::Inside,
35                (InsideResult::Boundary, _) => InsideResult::Boundary,
36                (_, InsideResult::Boundary) => InsideResult::Boundary,
37                _ => InsideResult::Outside,
38            }
39        }
40    }
41}
42
43/// Check if a point is inside a geometry.
44///
45/// If the geoemtry we check against is a Point, check if the points are equal
46///
47/// If the geometry we check against is a MultiPoint or LineString, treat as a polygon with no holes
48///
49/// If the geometry we check against is a Polygon or MultiPolygon, check if the point is inside
50///
51/// This trait is implemented for:
52/// - [`Feature`]
53/// - [`Geometry`]
54/// - [`PointGeometry`]
55/// - [`MultiPointGeometry`]
56/// - [`MultiLineStringGeometry`]
57/// - [`MultiPolygonGeometry`]
58/// - [`Point3DGeometry`]
59/// - [`MultiPoint3DGeometry`]
60/// - [`MultiLineString3DGeometry`]
61/// - [`MultiPolygon3DGeometry`]
62/// - [`VectorFeature`]
63/// - [`VectorGeometry`]
64/// - [`VectorPointGeometry`]
65/// - [`VectorMultiPointGeometry`]
66/// - [`VectorMultiLineStringGeometry`]
67/// - [`VectorMultiPolygonGeometry`]
68/// - [`VectorMultiPoint`]
69/// - [`VectorMultiLineString`]
70/// - [`VectorMultiPolygon`]
71/// - `&Vec<Vec<P>>` or `&[Vec<P>]` where P implements [`GetXY`]
72///
73/// And all specific geometries of the above enums
74pub trait Inside {
75    /// Check if a point is inside a geometry.
76    ///
77    /// If the geoemtry we check against is a Point, check if the points are equal
78    ///
79    /// If the geometry we check against is a MultiPoint or LineString, treat as a polygon with no holes
80    ///
81    /// If the geometry we check against is a Polygon or MultiPolygon, check if the point is inside
82    fn inside<P: GetXY>(&self, b: &P) -> InsideResult;
83}
84
85// Feature and below
86
87impl<P: GetXY> Inside for &Vec<Vec<P>> {
88    fn inside<B: GetXY>(&self, b: &B) -> InsideResult {
89        point_in_polygon(b, self)
90    }
91}
92impl<P: GetXY> Inside for &[Vec<P>] {
93    fn inside<B: GetXY>(&self, b: &B) -> InsideResult {
94        point_in_polygon(b, self)
95    }
96}
97
98impl<M, P: Clone + Default, D: Clone + Default> Inside for Feature<M, P, D> {
99    fn inside<B: GetXY>(&self, b: &B) -> InsideResult {
100        self.geometry.inside(b)
101    }
102}
103impl<M: Clone + Default> Inside for Geometry<M> {
104    fn inside<P: GetXY>(&self, b: &P) -> InsideResult {
105        match self {
106            Geometry::Point(g) => g.inside(b),
107            Geometry::MultiPoint(g) => g.inside(b),
108            Geometry::LineString(g) => g.inside(b),
109            Geometry::MultiLineString(g) => g.inside(b),
110            Geometry::Polygon(g) => g.inside(b),
111            Geometry::MultiPolygon(g) => g.inside(b),
112            Geometry::Point3D(g) => g.inside(b),
113            Geometry::MultiPoint3D(g) => g.inside(b),
114            Geometry::LineString3D(g) => g.inside(b),
115            Geometry::MultiLineString3D(g) => g.inside(b),
116            Geometry::Polygon3D(g) => g.inside(b),
117            Geometry::MultiPolygon3D(g) => g.inside(b),
118        }
119    }
120}
121impl<M: Clone + Default> Inside for PointGeometry<M> {
122    fn inside<P: GetXY>(&self, b: &P) -> InsideResult {
123        self.coordinates.inside(b)
124    }
125}
126impl<M: Clone + Default> Inside for MultiPointGeometry<M> {
127    fn inside<P: GetXY>(&self, b: &P) -> InsideResult {
128        self.coordinates.inside(b)
129    }
130}
131impl<M: Clone + Default> Inside for MultiLineStringGeometry<M> {
132    fn inside<P: GetXY>(&self, b: &P) -> InsideResult {
133        self.coordinates.inside(b)
134    }
135}
136impl<M: Clone + Default> Inside for MultiPolygonGeometry<M> {
137    fn inside<P: GetXY>(&self, b: &P) -> InsideResult {
138        self.coordinates.inside(b)
139    }
140}
141impl<M: Clone + Default> Inside for Point3DGeometry<M> {
142    fn inside<P: GetXY>(&self, b: &P) -> InsideResult {
143        self.coordinates.inside(b)
144    }
145}
146impl<M: Clone + Default> Inside for MultiPoint3DGeometry<M> {
147    fn inside<P: GetXY>(&self, b: &P) -> InsideResult {
148        self.coordinates.inside(b)
149    }
150}
151impl<M: Clone + Default> Inside for MultiLineString3DGeometry<M> {
152    fn inside<P: GetXY>(&self, b: &P) -> InsideResult {
153        self.coordinates.inside(b)
154    }
155}
156impl<M: Clone + Default> Inside for MultiPolygon3DGeometry<M> {
157    fn inside<P: GetXY>(&self, b: &P) -> InsideResult {
158        self.coordinates.inside(b)
159    }
160}
161
162// Feature Point types
163
164impl Inside for Point {
165    fn inside<P: GetXY>(&self, b: &P) -> InsideResult {
166        if self.x() == b.x() && self.y() == b.y() {
167            InsideResult::Inside
168        } else {
169            InsideResult::Outside
170        }
171    }
172}
173impl Inside for MultiPoint {
174    fn inside<P: GetXY>(&self, b: &P) -> InsideResult {
175        point_in_polyline(b, self)
176    }
177}
178impl Inside for MultiLineString {
179    fn inside<P: GetXY>(&self, b: &P) -> InsideResult {
180        point_in_polygon(b, self)
181    }
182}
183impl Inside for MultiPolygon {
184    fn inside<P: GetXY>(&self, b: &P) -> InsideResult {
185        let mut res = InsideResult::Outside;
186        for poly in self {
187            res = res.get_dominant(poly.inside(b));
188        }
189        res
190    }
191}
192impl Inside for Point3D {
193    fn inside<P: GetXY>(&self, b: &P) -> InsideResult {
194        if self.x() == b.x() && self.y() == b.y() {
195            InsideResult::Inside
196        } else {
197            InsideResult::Outside
198        }
199    }
200}
201impl Inside for MultiPoint3D {
202    fn inside<P: GetXY>(&self, b: &P) -> InsideResult {
203        point_in_polyline(b, self)
204    }
205}
206impl Inside for MultiLineString3D {
207    fn inside<P: GetXY>(&self, b: &P) -> InsideResult {
208        point_in_polygon(b, self)
209    }
210}
211impl Inside for MultiPolygon3D {
212    fn inside<P: GetXY>(&self, b: &P) -> InsideResult {
213        let mut res = InsideResult::Outside;
214        for poly in self {
215            res = res.get_dominant(poly.inside(b));
216        }
217        res
218    }
219}
220
221// Vector Feature and below
222
223impl<M, P: Clone + Default, D: Clone + Default> Inside for VectorFeature<M, P, D> {
224    fn inside<B: GetXY>(&self, b: &B) -> InsideResult {
225        self.geometry.inside(b)
226    }
227}
228impl<M: Clone + Default> Inside for VectorGeometry<M> {
229    fn inside<P: GetXY>(&self, b: &P) -> InsideResult {
230        match self {
231            VectorGeometry::Point(g) => g.inside(b),
232            VectorGeometry::MultiPoint(g) => g.inside(b),
233            VectorGeometry::LineString(g) => g.inside(b),
234            VectorGeometry::MultiLineString(g) => g.inside(b),
235            VectorGeometry::Polygon(g) => g.inside(b),
236            VectorGeometry::MultiPolygon(g) => g.inside(b),
237        }
238    }
239}
240impl<M: Clone + Default> Inside for VectorPointGeometry<M> {
241    fn inside<P: GetXY>(&self, b: &P) -> InsideResult {
242        self.coordinates.inside(b)
243    }
244}
245impl<M: Clone + Default> Inside for VectorMultiPointGeometry<M> {
246    fn inside<P: GetXY>(&self, b: &P) -> InsideResult {
247        self.coordinates.inside(b)
248    }
249}
250impl<M: Clone + Default> Inside for VectorMultiLineStringGeometry<M> {
251    fn inside<P: GetXY>(&self, b: &P) -> InsideResult {
252        self.coordinates.inside(b)
253    }
254}
255impl<M: Clone + Default> Inside for VectorMultiPolygonGeometry<M> {
256    fn inside<P: GetXY>(&self, b: &P) -> InsideResult {
257        self.coordinates.inside(b)
258    }
259}
260
261// Vector Point Types
262
263impl<M: Clone + Default> Inside for VectorPoint<M> {
264    fn inside<P: GetXY>(&self, b: &P) -> InsideResult {
265        if self.x() == b.x() && self.y() == b.y() {
266            InsideResult::Inside
267        } else {
268            InsideResult::Outside
269        }
270    }
271}
272impl<M: Clone + Default> Inside for VectorMultiPoint<M> {
273    fn inside<P: GetXY>(&self, b: &P) -> InsideResult {
274        point_in_polyline(b, self)
275    }
276}
277impl<M: Clone + Default> Inside for VectorMultiLineString<M> {
278    fn inside<P: GetXY>(&self, b: &P) -> InsideResult {
279        point_in_polygon(b, self)
280    }
281}
282impl<M: Clone + Default> Inside for VectorMultiPolygon<M> {
283    fn inside<P: GetXY>(&self, b: &P) -> InsideResult {
284        let mut res = InsideResult::Outside;
285        for poly in self {
286            res = res.get_dominant(poly.inside(b));
287        }
288        res
289    }
290}
291
292/// A Robust point in polygon test
293///
294/// ## Parameters
295/// - `point`: the point
296/// - `polygon`: the polygon
297///
298/// ## Returns
299/// true if the point is in the polygon, 0 if on the boundary, false otherwise
300pub fn point_in_polygon<P1: GetXY, P2: GetXY>(point: &P1, polygon: &[Vec<P2>]) -> InsideResult {
301    let mut k = 0;
302    let mut f;
303    let mut u1;
304    let mut v1;
305    let mut u2;
306    let mut v2;
307    let mut current_p: &P2;
308    let mut next_p: &P2;
309
310    let x = point.x();
311    let y = point.y();
312    for contour in polygon {
313        let contour_len = contour.len() - 1;
314
315        current_p = &contour[0];
316        if current_p.x() != contour[contour_len].x() && current_p.y() != contour[contour_len].y() {
317            // since the first and last coordinates in a ring are not the same, assume it's not a polygon and return false
318            return InsideResult::Outside;
319        }
320
321        u1 = current_p.x() - x;
322        v1 = current_p.y() - y;
323
324        for ii in 0..contour_len {
325            next_p = &contour[ii + 1];
326
327            u2 = next_p.x() - x;
328            v2 = next_p.y() - y;
329
330            if v1 == 0. && v2 == 0. {
331                if (u2 <= 0. && u1 >= 0.) || (u1 <= 0. && u2 >= 0.) {
332                    return InsideResult::Boundary;
333                }
334            } else if (v2 >= 0. && v1 <= 0.) || (v2 <= 0. && v1 >= 0.) {
335                f = orient2d(u1, u2, v1, v2, 0., 0.);
336                if f == 0. {
337                    return InsideResult::Boundary;
338                }
339                if (f > 0. && v2 > 0. && v1 <= 0.) || (f < 0. && v2 <= 0. && v1 > 0.) {
340                    k += 1;
341                }
342            }
343            // current_p = next_p;
344            v1 = v2;
345            u1 = u2;
346        }
347    }
348
349    if k % 2 == 0 { InsideResult::Outside } else { InsideResult::Inside }
350}
351
352/// A Robust point in polygon test
353///
354/// ## Parameters
355/// - `point`: the point
356/// - `polygon`: the polygon
357///
358/// ## Returns
359/// true if the point is in the polygon, 0 if on the boundary, false otherwise
360pub fn point_in_polyline<P1: GetXY, P2: GetXY>(point: &P1, contour: &[P2]) -> InsideResult {
361    let mut k = 0;
362    let mut f;
363    let mut u1;
364    let mut v1;
365    let mut u2;
366    let mut v2;
367    let mut next_p: &P2;
368
369    let x = point.x();
370    let y = point.y();
371    let contour_len = contour.len() - 1;
372
373    let current_p = &contour[0];
374    if current_p.x() != contour[contour_len].x() && current_p.y() != contour[contour_len].y() {
375        // since the first and last coordinates in a ring are not the same, assume it's not a polygon and return false
376        return InsideResult::Outside;
377    }
378
379    u1 = current_p.x() - x;
380    v1 = current_p.y() - y;
381
382    for ii in 0..contour_len {
383        next_p = &contour[ii + 1];
384
385        u2 = next_p.x() - x;
386        v2 = next_p.y() - y;
387
388        if v1 == 0. && v2 == 0. {
389            if (u2 <= 0. && u1 >= 0.) || (u1 <= 0. && u2 >= 0.) {
390                return InsideResult::Boundary;
391            }
392        } else if (v2 >= 0. && v1 <= 0.) || (v2 <= 0. && v1 >= 0.) {
393            f = orient2d(u1, u2, v1, v2, 0., 0.);
394            if f == 0. {
395                return InsideResult::Boundary;
396            }
397            if (f > 0. && v2 > 0. && v1 <= 0.) || (f < 0. && v2 <= 0. && v1 > 0.) {
398                k += 1;
399            }
400        }
401        // current_p = next_p;
402        v1 = v2;
403        u1 = u2;
404    }
405
406    if k % 2 == 0 { InsideResult::Outside } else { InsideResult::Inside }
407}
408
409/// Check if a polyline/hole is inside another polyline/outer ring
410///
411/// ## Parameters
412/// - `hole`: the hole to test if inside the outer
413/// - `outer`: the outer to test against
414///
415/// ## Returns
416/// true if the hole is inside the outer
417pub fn polyline_in_polyline<P: GetXY>(hole: &[P], outer: &[P]) -> bool {
418    for point in hole {
419        match point_in_polyline(point, outer) {
420            InsideResult::Inside => return true,
421            InsideResult::Outside => return false,
422            InsideResult::Boundary => continue,
423        }
424    }
425    // If we make it means all points of the hole were on the boundary therefore its inside the outer
426    true
427}
gistools/geometry/tools/polys/inside.rs

gistools/geometry/tools/polys/
inside.rs
