diff --git a/geo-types/CHANGES.md b/geo-types/CHANGES.md index abd44c63f..18270844c 100644 --- a/geo-types/CHANGES.md +++ b/geo-types/CHANGES.md @@ -1,5 +1,10 @@ # Changes +## Master (not released) + +* Add vector-space operations to `Coordinate` and `Point` + * + ## 0.6.0 * Remove `COORD_PRECISION` which was an arbitrary constant of 0.1m @@ -114,4 +119,4 @@ ## 0.1.0 * New crate with core types from `geo` - * \ No newline at end of file + * diff --git a/geo-types/src/line.rs b/geo-types/src/line.rs index 0395f1b14..4b4d91efe 100644 --- a/geo-types/src/line.rs +++ b/geo-types/src/line.rs @@ -1,6 +1,9 @@ use crate::{Coordinate, CoordinateType, Point}; -/// A line segment made up of exactly two [`Point`s](struct.Point.html). +/// A line segment made up of exactly two +/// [`Coordinate`s](struct.Coordinate.html). The interior and +/// boundaries are defined as with a `LineString` with the +/// two end points. #[derive(Eq, PartialEq, Clone, Copy, Debug, Hash)] #[cfg_attr(feature = "serde", derive(Serialize, Deserialize))] pub struct Line @@ -37,6 +40,11 @@ where } } + /// Calculate the difference in coordinates (Δx, Δy). + pub fn delta(&self) -> Coordinate { + self.end - self.start + } + /// Calculate the difference in ‘x’ components (Δx). /// /// Equivalent to: @@ -53,7 +61,7 @@ where /// # ); /// ``` pub fn dx(&self) -> T { - self.end.x - self.start.x + self.delta().x } /// Calculate the difference in ‘y’ components (Δy). @@ -72,7 +80,7 @@ where /// # ); /// ``` pub fn dy(&self) -> T { - self.end.y - self.start.y + self.delta().y } /// Calculate the slope (Δy/Δx). diff --git a/geo-types/src/line_string.rs b/geo-types/src/line_string.rs index bfef99934..6d3ee2786 100644 --- a/geo-types/src/line_string.rs +++ b/geo-types/src/line_string.rs @@ -2,9 +2,25 @@ use crate::{Coordinate, CoordinateType, Line, Point, Triangle}; use std::iter::FromIterator; use std::ops::{Index, IndexMut}; -/// An ordered collection of two or more [`Coordinate`s](struct.Coordinate.html), representing a +/// An ordered collection of two or more +/// [`Coordinate`s](struct.Coordinate.html), representing a /// path between locations. /// +/// A `LineString` is _closed_ if it is empty, or if the +/// first and last coordinates are the same. The _boundary_ +/// of a `LineString` is empty if closed, and otherwise the +/// end points. The interior is the (infinite) set of all +/// points along the linestring _not including_ the +/// boundary. +/// +/// # Validity +/// +/// A `LineString` is valid if it is either empty or +/// contains 2 or more coordinates. Further, a closed +/// `LineString` must not self intersect. Note that the +/// validity is not enforced, and the operations and +/// predicates are undefined on invalid linestrings. +/// /// # Examples /// /// Create a `LineString` by calling it directly: @@ -180,11 +196,9 @@ impl LineString { self.0.len() } - /// Checks if the linestring is closed; i.e. the first - /// and last points have the same coords. Note that a - /// single point is considered closed, but the output on - /// an empty linestring is _unspecified_ and must not be - /// relied upon. + /// Checks if the linestring is closed; i.e. it is + /// either empty or, the first and last points are the + /// same. /// /// # Examples /// diff --git a/geo-types/src/multi_line_string.rs b/geo-types/src/multi_line_string.rs index c323cf340..fe93b4970 100644 --- a/geo-types/src/multi_line_string.rs +++ b/geo-types/src/multi_line_string.rs @@ -1,11 +1,16 @@ use crate::{CoordinateType, LineString}; use std::iter::FromIterator; -/// A collection of [`LineString`s](line_string/struct.LineString.html). +/// A collection of +/// [`LineString`s](line_string/struct.LineString.html). The +/// interior and the boundary are the union of the interior +/// or the boundary of the constituent line strings. /// -/// Can be created from a `Vec` of `LineString`s, or from an Iterator which yields `LineString`s. +/// Can be created from a `Vec` of `LineString`s, or from an +/// Iterator which yields `LineString`s. /// -/// Iterating over this objects, yields the component `LineString`s. +/// Iterating over this objects, yields the component +/// `LineString`s. #[derive(Eq, PartialEq, Clone, Debug, Hash)] #[cfg_attr(feature = "serde", derive(Serialize, Deserialize))] pub struct MultiLineString(pub Vec>) diff --git a/geo-types/src/multi_point.rs b/geo-types/src/multi_point.rs index ea54fc34d..3db979b9b 100644 --- a/geo-types/src/multi_point.rs +++ b/geo-types/src/multi_point.rs @@ -1,7 +1,9 @@ use crate::{CoordinateType, Point}; use std::iter::FromIterator; -/// A collection of [`Point`s](struct.Point.html). +/// A collection of [`Point`s](struct.Point.html). The +/// interior and the boundary are the union of the interior +/// or the boundary of the constituent points. /// /// # Examples /// diff --git a/geo-types/src/multi_polygon.rs b/geo-types/src/multi_polygon.rs index b0b14fad1..db816ae90 100644 --- a/geo-types/src/multi_polygon.rs +++ b/geo-types/src/multi_polygon.rs @@ -1,7 +1,9 @@ use crate::{CoordinateType, Polygon}; use std::iter::FromIterator; -/// A collection of [`Polygon`s](struct.Polygon.html). +/// A collection of [`Polygon`s](struct.Polygon.html). The +/// interior and the boundary are the union of the interior +/// or the boundary of the constituent polygons. /// /// Can be created from a `Vec` of `Polygon`s, or `collect`ed from an Iterator which yields `Polygon`s. /// diff --git a/geo-types/src/point.rs b/geo-types/src/point.rs index 7e4151fc4..4a4890fe8 100644 --- a/geo-types/src/point.rs +++ b/geo-types/src/point.rs @@ -9,6 +9,10 @@ use std::ops::{Add, Div, Mul, Neg, Sub}; /// tuples, or arrays – see the `From` impl section for a /// complete list. /// +/// A point is _valid_ iff neither coordinate is `NaN`. The +/// _interior_ of the point is itself (a singleton set), and +/// its _boundary_ is empty. +/// /// # Examples /// /// ``` diff --git a/geo-types/src/polygon.rs b/geo-types/src/polygon.rs index 0114b7beb..d8da9235f 100644 --- a/geo-types/src/polygon.rs +++ b/geo-types/src/polygon.rs @@ -7,6 +7,11 @@ use num_traits::{Float, Signed}; /// [`LineString`]. It may contain zero or more holes (_interior rings_), also /// represented by `LineString`s. /// +/// The _boundary_ of the polygon is the union of the +/// boundaries of the exterior and interiors. The interior +/// is all the points inside the polygon (not on the +/// boundary). +/// /// The `Polygon` structure guarantees that all exterior and interior rings will /// be _closed_, such that the first and last `Coordinate` of each ring has /// the same value. diff --git a/geo/CHANGES.md b/geo/CHANGES.md index 814544db4..bbe703763 100644 --- a/geo/CHANGES.md +++ b/geo/CHANGES.md @@ -1,5 +1,15 @@ # Changes +## Master (not released) + +* Add robust predicates + * + * + * + * +* Improve numerical stability in centroid computation + * + ## 0.14.2 * Bump proj version to 0.20.3 diff --git a/geo/src/algorithm/closest_point.rs b/geo/src/algorithm/closest_point.rs index 10600bc97..f629ba213 100644 --- a/geo/src/algorithm/closest_point.rs +++ b/geo/src/algorithm/closest_point.rs @@ -1,3 +1,4 @@ +use crate::kernels::*; use crate::prelude::*; use crate::{Closest, Line, LineString, MultiLineString, MultiPoint, MultiPolygon, Point, Polygon}; use num_traits::Float; @@ -47,7 +48,7 @@ impl ClosestPoint for Point { } } -impl ClosestPoint for Line { +impl ClosestPoint for Line { fn closest_point(&self, p: &Point) -> Closest { let line_length = self.euclidean_length(); if line_length == F::zero() { @@ -77,7 +78,7 @@ impl ClosestPoint for Line { let y = direction_vector.y(); let c = Point(self.start + (t * x, t * y).into()); - if self.contains(p) { + if self.intersects(p) { Closest::Intersection(c) } else { Closest::SinglePoint(c) @@ -106,20 +107,20 @@ where best } -impl ClosestPoint for LineString { +impl ClosestPoint for LineString { fn closest_point(&self, p: &Point) -> Closest { closest_of(self.lines(), *p) } } -impl ClosestPoint for Polygon { +impl ClosestPoint for Polygon { fn closest_point(&self, p: &Point) -> Closest { let prospectives = self.interiors().iter().chain(iter::once(self.exterior())); closest_of(prospectives, *p) } } -impl ClosestPoint for MultiPolygon { +impl ClosestPoint for MultiPolygon { fn closest_point(&self, p: &Point) -> Closest { closest_of(self.0.iter(), *p) } @@ -131,7 +132,7 @@ impl ClosestPoint for MultiPoint { } } -impl ClosestPoint for MultiLineString { +impl ClosestPoint for MultiLineString { fn closest_point(&self, p: &Point) -> Closest { closest_of(self.0.iter(), *p) } diff --git a/geo/src/algorithm/contains/geometry.rs b/geo/src/algorithm/contains/geometry.rs new file mode 100644 index 000000000..da1121206 --- /dev/null +++ b/geo/src/algorithm/contains/geometry.rs @@ -0,0 +1,58 @@ +use super::Contains; +use crate::kernels::*; +use crate::*; + +// ┌──────────────────────────────┐ +// │ Implementations for Geometry │ +// └──────────────────────────────┘ + +impl Contains> for Geometry +where + T: HasKernel, +{ + fn contains(&self, coord: &Coordinate) -> bool { + match self { + Geometry::Point(g) => g.contains(coord), + Geometry::Line(g) => g.contains(coord), + Geometry::LineString(g) => g.contains(coord), + Geometry::Polygon(g) => g.contains(coord), + Geometry::MultiPoint(g) => g.contains(coord), + Geometry::MultiLineString(g) => g.contains(coord), + Geometry::MultiPolygon(g) => g.contains(coord), + Geometry::GeometryCollection(g) => g.contains(coord), + Geometry::Rect(g) => g.contains(coord), + Geometry::Triangle(g) => g.contains(coord), + } + } +} + +impl Contains> for Geometry +where + T: HasKernel, +{ + fn contains(&self, point: &Point) -> bool { + self.contains(&point.0) + } +} + +// ┌────────────────────────────────────────┐ +// │ Implementations for GeometryCollection │ +// └────────────────────────────────────────┘ + +impl Contains> for GeometryCollection +where + T: HasKernel, +{ + fn contains(&self, coord: &Coordinate) -> bool { + self.0.iter().any(|geometry| geometry.contains(coord)) + } +} + +impl Contains> for GeometryCollection +where + T: HasKernel, +{ + fn contains(&self, point: &Point) -> bool { + self.contains(&point.0) + } +} diff --git a/geo/src/algorithm/contains/line.rs b/geo/src/algorithm/contains/line.rs new file mode 100644 index 000000000..9b60b2e8e --- /dev/null +++ b/geo/src/algorithm/contains/line.rs @@ -0,0 +1,84 @@ +use super::Contains; +use crate::intersects::Intersects; +use crate::kernels::*; +use crate::*; + +// ┌──────────────────────────┐ +// │ Implementations for Line │ +// └──────────────────────────┘ + +impl Contains> for Line +where + T: HasKernel, +{ + fn contains(&self, coord: &Coordinate) -> bool { + if self.start == self.end { + &self.start == coord + } else { + coord != &self.start && coord != &self.end && self.intersects(coord) + } + } +} + +impl Contains> for Line +where + T: HasKernel, +{ + fn contains(&self, p: &Point) -> bool { + self.contains(&p.0) + } +} + +impl Contains> for Line +where + T: HasKernel, +{ + fn contains(&self, line: &Line) -> bool { + if line.start == line.end { + self.contains(&line.start) + } else { + self.intersects(&line.start) && self.intersects(&line.end) + } + } +} + +impl Contains> for Line +where + T: HasKernel, +{ + fn contains(&self, linestring: &LineString) -> bool { + // Empty linestring has no interior, and not + // contained in anything. + if linestring.0.is_empty() { + return false; + } + + // The interior of the linestring should have some + // intersection with the interior of self. Two cases + // arise: + // + // 1. There are at least two distinct points in the + // linestring. Then, if both intersect, the interior + // between these two must have non-empty intersection. + // + // 2. Otherwise, all the points on the linestring + // are the same. In this case, the interior is this + // specific point, and it should be contained in the + // line. + let first = linestring.0.first().unwrap(); + let mut all_equal = true; + + // If all the vertices of the linestring intersect + // self, then the interior or boundary of the + // linestring cannot have non-empty intersection + // with the exterior. + let all_intersects = linestring.0.iter().all(|c| { + if c != first { + all_equal = false; + } + self.intersects(c) + }); + + all_intersects && (!all_equal || self.contains(first)) + } +} diff --git a/geo/src/algorithm/contains/line_string.rs b/geo/src/algorithm/contains/line_string.rs new file mode 100644 index 000000000..9175110f1 --- /dev/null +++ b/geo/src/algorithm/contains/line_string.rs @@ -0,0 +1,128 @@ +use super::Contains; +use crate::kernels::*; +use crate::*; +use intersects::Intersects; + +// ┌────────────────────────────────┐ +// │ Implementations for LineString │ +// └────────────────────────────────┘ + +impl Contains> for LineString +where + T: HasKernel, +{ + fn contains(&self, coord: &Coordinate) -> bool { + if self.0.is_empty() { + return false; + } + + if self.is_closed() && coord == &self.0[0] { + return true; + } + + self.lines() + .enumerate() + .any(|(i, line)| line.contains(coord) || (i > 0 && coord == &line.start)) + } +} + +impl Contains> for LineString +where + T: HasKernel, +{ + fn contains(&self, p: &Point) -> bool { + self.contains(&p.0) + } +} + +impl Contains> for LineString +where + T: HasKernel, +{ + fn contains(&self, line: &Line) -> bool { + if line.start == line.end { + return self.contains(&line.start); + } + + // We copy the line as we may truncate the line as + // we find partial matches. + let mut line = *line; + let mut first_cut = None; + + let lines_iter = self.lines(); + let num_lines = lines_iter.len(); + + // We need to repeat the logic twice to handle cases + // where the linestring starts at the middle of the line. + for (i, segment) in self.lines().chain(lines_iter).enumerate() { + if i >= num_lines { + // The first loop was done. If we never cut + // the line, or at the cut segment again, we + // can exit now. + if let Some(upto_i) = first_cut { + if i >= num_lines + upto_i { + break; + } + } else { + break; + } + } + // Look for a segment that intersects at least + // one of the end points. + let other = if segment.intersects(&line.start) { + line.end + } else if segment.intersects(&line.end) { + line.start + } else { + continue; + }; + + // If the other end point also intersects this + // segment, then we are done. + let new_inside = if segment.intersects(&other) { + return true; + } + // otoh, if the line contains one of the ends of + // the segments, then we truncate the line to + // the part outside. + else if line.contains(&segment.start) { + segment.start + } else if line.contains(&segment.end) { + segment.end + } else { + continue; + }; + + first_cut = first_cut.or(Some(i)); + if other == line.start { + line.end = new_inside; + } else { + line.start = new_inside; + } + } + + false + } +} + +impl Contains> for LineString +where + T: HasKernel, +{ + fn contains(&self, rhs: &LineString) -> bool { + rhs.lines().all(|l| self.contains(&l)) + } +} + +// ┌─────────────────────────────────────┐ +// │ Implementations for MultiLineString │ +// └─────────────────────────────────────┘ +impl Contains for MultiLineString +where + T: CoordinateType, + LineString: Contains, +{ + fn contains(&self, rhs: &G) -> bool { + self.0.iter().any(|p| p.contains(rhs)) + } +} diff --git a/geo/src/algorithm/contains.rs b/geo/src/algorithm/contains/mod.rs similarity index 51% rename from geo/src/algorithm/contains.rs rename to geo/src/algorithm/contains/mod.rs index 55b30c976..737996187 100644 --- a/geo/src/algorithm/contains.rs +++ b/geo/src/algorithm/contains/mod.rs @@ -1,434 +1,48 @@ -use num_traits::Float; - -use crate::algorithm::intersects::Intersects; -use crate::utils; -use crate::{ - Coordinate, CoordinateType, Geometry, GeometryCollection, Line, LineString, MultiLineString, - MultiPoint, MultiPolygon, Point, Polygon, Rect, Triangle, -}; - -/// Checks if the geometry A is completely inside the B geometry +/// Checks if `rhs` is completely contained within `self`. +/// More formally, the interior of `rhs` has non-empty +/// (set-theoretic) intersection but neither the interior, +/// nor the boundary of `rhs` intersects the exterior of +/// `self`. In other words, the [DE-9IM] intersection matrix +/// of `(rhs, self)` is `T*F**F***`. +/// +/// [DE-9IM]: https://en.wikipedia.org/wiki/DE-9IM +/// +/// # Examples +/// +/// ``` +/// use geo::algorithm::contains::Contains; +/// use geo::{line_string, point, Polygon}; +/// +/// let line_string = line_string![ +/// (x: 0., y: 0.), +/// (x: 2., y: 0.), +/// (x: 2., y: 2.), +/// (x: 0., y: 2.), +/// (x: 0., y: 0.), +/// ]; +/// +/// let polygon = Polygon::new(line_string.clone(), vec![]); +/// +/// // Point in Point +/// assert!(point!(x: 2., y: 0.).contains(&point!(x: 2., y: 0.))); +/// +/// // Point in Linestring +/// assert!(line_string.contains(&point!(x: 2., y: 0.))); +/// +/// // Point in Polygon +/// assert!(polygon.contains(&point!(x: 1., y: 1.))); +/// ``` pub trait Contains { - /// Checks if `rhs` is completely contained within `self`. - /// - /// # Examples - /// - /// ``` - /// use geo::algorithm::contains::Contains; - /// use geo::{line_string, point, Polygon}; - /// - /// let line_string = line_string![ - /// (x: 0., y: 0.), - /// (x: 2., y: 0.), - /// (x: 2., y: 2.), - /// (x: 0., y: 2.), - /// (x: 0., y: 0.), - /// ]; - /// - /// let polygon = Polygon::new(line_string.clone(), vec![]); - /// - /// // Point in Point - /// assert!(point!(x: 2., y: 0.).contains(&point!(x: 2., y: 0.))); - /// - /// // Point in Linestring - /// assert!(line_string.contains(&point!(x: 2., y: 0.))); - /// - /// // Point in Polygon - /// assert!(polygon.contains(&point!(x: 1., y: 1.))); - /// ``` fn contains(&self, rhs: &Rhs) -> bool; } -// ┌───────────────────────────┐ -// │ Implementations for Point │ -// └───────────────────────────┘ - -impl Contains> for Point -where - T: Float, -{ - fn contains(&self, coord: &Coordinate) -> bool { - self.contains(&Point(*coord)) - } -} - -impl Contains> for Point -where - T: Float, -{ - fn contains(&self, p: &Point) -> bool { - ::geo_types::private_utils::point_contains_point(*self, *p) - } -} - -// ┌────────────────────────────────┐ -// │ Implementations for MultiPoint │ -// └────────────────────────────────┘ - -impl Contains> for MultiPoint -where - T: Float, -{ - fn contains(&self, coord: &Coordinate) -> bool { - self.0.iter().any(|point| point.contains(coord)) - } -} - -impl Contains> for MultiPoint -where - T: Float, -{ - fn contains(&self, point: &Point) -> bool { - self.contains(&point.0) - } -} - -// ┌──────────────────────────┐ -// │ Implementations for Line │ -// └──────────────────────────┘ - -impl Contains> for Line -where - T: Float, -{ - fn contains(&self, coord: &Coordinate) -> bool { - self.contains(&Point(*coord)) - } -} - -impl Contains> for Line -where - T: Float, -{ - fn contains(&self, p: &Point) -> bool { - self.intersects(p) - } -} - -impl Contains> for Line -where - T: Float, -{ - fn contains(&self, line: &Line) -> bool { - self.contains(&line.start_point()) && self.contains(&line.end_point()) - } -} - -impl Contains> for Line -where - T: Float, -{ - fn contains(&self, linestring: &LineString) -> bool { - linestring.points_iter().all(|pt| self.contains(&pt)) - } -} - -// ┌────────────────────────────────┐ -// │ Implementations for LineString │ -// └────────────────────────────────┘ - -impl Contains> for LineString -where - T: Float, -{ - fn contains(&self, coord: &Coordinate) -> bool { - self.contains(&Point(*coord)) - } -} - -impl Contains> for LineString -where - T: Float, -{ - fn contains(&self, p: &Point) -> bool { - self.lines().any(|line| line.contains(p)) - } -} - -impl Contains> for LineString -where - T: Float, -{ - fn contains(&self, line: &Line) -> bool { - let (p0, p1) = line.points(); - let mut look_for: Option> = None; - for segment in self.lines() { - if look_for.is_none() { - // If segment contains an endpoint of line, we mark the other endpoint as the - // one we are looking for. - if segment.contains(&p0) { - look_for = Some(p1); - } else if segment.contains(&p1) { - look_for = Some(p0); - } - } - if let Some(p) = look_for { - // If we are looking for an endpoint, we need to either find it, or show that we - // should continue to look for it - if segment.contains(&p) { - // If the segment contains the endpoint we are looking for we are done - return true; - } else if !line.contains(&segment.end_point()) { - // If not, and the end of the segment is not on the line, we should stop - // looking - look_for = None - } - } - } - false - } -} - -// ┌─────────────────────────────────────┐ -// │ Implementations for MultiLineString │ -// └─────────────────────────────────────┘ - -impl Contains> for MultiLineString -where - T: Float, -{ - fn contains(&self, coord: &Coordinate) -> bool { - self.0.iter().any(|line_string| line_string.contains(coord)) - } -} - -impl Contains> for MultiLineString -where - T: Float, -{ - fn contains(&self, point: &Point) -> bool { - self.contains(&point.0) - } -} - -// ┌─────────────────────────────┐ -// │ Implementations for Polygon │ -// └─────────────────────────────┘ - -impl Contains> for Polygon -where - T: Float, -{ - fn contains(&self, coord: &Coordinate) -> bool { - match utils::coord_pos_relative_to_line_string(*coord, &self.exterior()) { - utils::CoordPos::OnBoundary | utils::CoordPos::Outside => false, - _ => self.interiors().iter().all(|ls| { - utils::coord_pos_relative_to_line_string(*coord, ls) == utils::CoordPos::Outside - }), - } - } -} - -impl Contains> for Polygon -where - T: Float, -{ - fn contains(&self, p: &Point) -> bool { - self.contains(&p.0) - } -} - -impl Contains> for Polygon -where - T: Float, -{ - fn contains(&self, line: &Line) -> bool { - // both endpoints are contained in the polygon and the line - // does NOT intersect the exterior or any of the interior boundaries - self.contains(&line.start_point()) - && self.contains(&line.end_point()) - && !self.exterior().intersects(line) - && !self.interiors().iter().any(|inner| inner.intersects(line)) - } -} - -impl Contains> for Polygon -where - T: Float, -{ - fn contains(&self, poly: &Polygon) -> bool { - // decompose poly's exterior ring into Lines, and check each for containment - poly.exterior().lines().all(|line| self.contains(&line)) - } -} - -impl Contains> for Polygon -where - T: Float, -{ - fn contains(&self, linestring: &LineString) -> bool { - // All LineString points must be inside the Polygon - if linestring.points_iter().all(|point| self.contains(&point)) { - // The Polygon interior is allowed to intersect with the LineString - // but the Polygon's rings are not - !self - .interiors() - .iter() - .any(|ring| ring.intersects(linestring)) - } else { - false - } - } -} - -// ┌──────────────────────────────────┐ -// │ Implementations for MultiPolygon │ -// └──────────────────────────────────┘ - -impl Contains> for MultiPolygon -where - T: Float, -{ - fn contains(&self, coord: &Coordinate) -> bool { - self.0.iter().any(|poly| poly.contains(coord)) - } -} - -impl Contains> for MultiPolygon -where - T: Float, -{ - fn contains(&self, p: &Point) -> bool { - self.contains(&p.0) - } -} - -// ┌──────────────────────────┐ -// │ Implementations for Rect │ -// └──────────────────────────┘ - -impl Contains> for Rect -where - T: CoordinateType, -{ - fn contains(&self, coord: &Coordinate) -> bool { - coord.x >= self.min().x - && coord.x <= self.max().x - && coord.y >= self.min().y - && coord.y <= self.max().y - } -} - -impl Contains> for Rect -where - T: CoordinateType, -{ - fn contains(&self, p: &Point) -> bool { - self.contains(&p.0) - } -} - -impl Contains> for Rect -where - T: CoordinateType, -{ - fn contains(&self, bounding_rect: &Rect) -> bool { - // All points of LineString must be in the polygon ? - self.min().x <= bounding_rect.min().x - && self.max().x >= bounding_rect.max().x - && self.min().y <= bounding_rect.min().y - && self.max().y >= bounding_rect.max().y - } -} - -// ┌──────────────────────────────┐ -// │ Implementations for Triangle │ -// └──────────────────────────────┘ - -impl Contains> for Triangle -where - T: CoordinateType, -{ - fn contains(&self, coord: &Coordinate) -> bool { - let s = self.0.y * self.2.x - self.0.x * self.2.y - + (self.2.y - self.0.y) * coord.x - + (self.0.x - self.2.x) * coord.y; - let t = self.0.x * self.1.y - self.0.y * self.1.x - + (self.0.y - self.1.y) * coord.x - + (self.1.x - self.0.x) * coord.y; - - if (s < T::zero()) != (t < T::zero()) { - return false; - } - - let a = (T::zero() - self.1.y) * self.2.x - + self.0.y * (self.2.x - self.1.x) - + self.0.x * (self.1.y - self.2.y) - + self.1.x * self.2.y; - - if a == T::zero() { - false - } else if a < T::zero() { - s <= T::zero() && s + t >= a - } else { - s >= T::zero() && s + t <= a - } - } -} - -impl Contains> for Triangle -where - T: CoordinateType, -{ - fn contains(&self, point: &Point) -> bool { - self.contains(&point.0) - } -} - -// ┌──────────────────────────────┐ -// │ Implementations for Geometry │ -// └──────────────────────────────┘ - -impl Contains> for Geometry -where - T: Float, -{ - fn contains(&self, coord: &Coordinate) -> bool { - match self { - Geometry::Point(g) => g.contains(coord), - Geometry::Line(g) => g.contains(coord), - Geometry::LineString(g) => g.contains(coord), - Geometry::Polygon(g) => g.contains(coord), - Geometry::MultiPoint(g) => g.contains(coord), - Geometry::MultiLineString(g) => g.contains(coord), - Geometry::MultiPolygon(g) => g.contains(coord), - Geometry::GeometryCollection(g) => g.contains(coord), - Geometry::Rect(g) => g.contains(coord), - Geometry::Triangle(g) => g.contains(coord), - } - } -} - -impl Contains> for Geometry -where - T: Float, -{ - fn contains(&self, point: &Point) -> bool { - self.contains(&point.0) - } -} - -// ┌────────────────────────────────────────┐ -// │ Implementations for GeometryCollection │ -// └────────────────────────────────────────┘ - -impl Contains> for GeometryCollection -where - T: Float, -{ - fn contains(&self, coord: &Coordinate) -> bool { - self.0.iter().any(|geometry| geometry.contains(coord)) - } -} - -impl Contains> for GeometryCollection -where - T: Float, -{ - fn contains(&self, point: &Point) -> bool { - self.contains(&point.0) - } -} +mod geometry; +mod line; +mod line_string; +mod point; +mod polygon; +mod rect; +mod triangle; // ┌───────┐ // │ Tests │ @@ -524,7 +138,11 @@ mod test { #[test] fn linestring_point_is_vertex_test() { let linestring = LineString::from(vec![(0., 0.), (2., 0.), (2., 2.)]); - assert!(linestring.contains(&Point::new(2., 2.))); + // Note: the end points of a linestring are not + // considered to be "contained" + assert!(linestring.contains(&Point::new(2., 0.))); + assert!(!linestring.contains(&Point::new(0., 0.))); + assert!(!linestring.contains(&Point::new(2., 2.))); } #[test] fn linestring_test() { @@ -709,32 +327,64 @@ mod test { } #[test] fn linestring_in_line_test() { - let line = Line::from([(0., 1.), (3., 4.)]); + let line = Line::from([(0, 10), (30, 40)]); // linestring0 in line - let linestring0 = LineString::from(vec![(0.1, 1.1), (1., 2.), (1.5, 2.5)]); + let linestring0 = LineString::from(vec![(01, 11), (10, 20), (15, 25)]); // linestring1 starts and ends in line, but wanders in the middle - let linestring1 = LineString::from(vec![(0.1, 1.1), (2., 2.), (1.5, 2.5)]); + let linestring1 = LineString::from(vec![(01, 11), (20, 20), (15, 25)]); // linestring2 is co-linear, but extends beyond line - let linestring2 = LineString::from(vec![(0.1, 1.1), (1., 2.), (4., 5.)]); + let linestring2 = LineString::from(vec![(01, 11), (10, 20), (40, 50)]); // no part of linestring3 is contained in line - let linestring3 = LineString::from(vec![(1.1, 1.1), (2., 2.), (2.5, 2.5)]); + let linestring3 = LineString::from(vec![(11, 11), (20, 20), (25, 25)]); + // a linestring with singleton interior on the boundary of the line + let linestring4 = LineString::from(vec![(0, 10), (0, 10), (0, 10)]); + // a linestring with singleton interior that is contained in the line + let linestring5 = LineString::from(vec![(1, 11), (1, 11), (1, 11)]); assert!(line.contains(&linestring0)); assert!(!line.contains(&linestring1)); assert!(!line.contains(&linestring2)); assert!(!line.contains(&linestring3)); + assert!(!line.contains(&linestring4)); + assert!(line.contains(&linestring5)); } #[test] fn line_in_polygon_test() { let c = |x, y| Coordinate { x, y }; - let line = Line::new(c(0., 1.), c(3., 4.)); - let linestring0 = line_string![c(-1., 0.), c(5., 0.), c(5., 5.), c(0., 5.), c(-1., 0.)]; + let line = Line::new(c(0, 10), c(30, 40)); + let linestring0 = line_string![c(-10, 0), c(50, 0), c(50, 50), c(0, 50), c(-10, 0)]; let poly0 = Polygon::new(linestring0, Vec::new()); - let linestring1 = line_string![c(0., 0.), c(0., 2.), c(2., 2.), c(2., 0.), c(0., 0.)]; + let linestring1 = line_string![c(0, 0), c(0, 20), c(20, 20), c(20, 0), c(0, 0)]; let poly1 = Polygon::new(linestring1, Vec::new()); assert!(poly0.contains(&line)); assert!(!poly1.contains(&line)); } #[test] + #[ignore] + fn line_in_polygon_edgecases_test() { + // Some DE-9IM edge cases for checking line is + // inside polygon The end points of the line can be + // on the boundary of the polygon but we don't allow + // that yet. + let c = |x, y| Coordinate { x, y }; + // A non-convex polygon + let linestring0 = line_string![c(0, 0), c(1, 1), c(1, -1), c(-1, -1), c(-1, 1)]; + let poly = Polygon::new(linestring0, Vec::new()); + + assert!(poly.contains(&Line::new(c(0, 0), c(1, -1)))); + assert!(poly.contains(&Line::new(c(-1, 1), c(1, -1)))); + assert!(!poly.contains(&Line::new(c(-1, 1), c(1, 1)))); + } + #[test] + fn line_in_linestring_edgecases() { + let c = |x, y| Coordinate { x, y }; + use crate::line_string; + let mut ls = line_string![c(0, 0), c(1, 0), c(0, 1), c(-1, 0)]; + assert!(!ls.contains(&Line::from([(0, 0), (0, 0)]))); + ls.close(); + assert!(ls.contains(&Line::from([(0, 0), (0, 0)]))); + assert!(ls.contains(&Line::from([(-1, 0), (1, 0)]))); + } + #[test] fn line_in_linestring_test() { let line0 = Line::from([(1., 1.), (2., 2.)]); // line0 is completely contained in the second segment @@ -775,17 +425,17 @@ mod test { } #[test] - fn triangle_contains_point_on_edge() { + fn triangle_not_contains_point_on_edge() { let t = Triangle::from([(0.0, 0.0), (2.0, 0.0), (2.0, 2.0)]); let p = Point::new(1.0, 0.0); - assert!(t.contains(&p)); + assert!(!t.contains(&p)); } #[test] - fn triangle_contains_point_on_vertex() { + fn triangle_not_contains_point_on_vertex() { let t = Triangle::from([(0.0, 0.0), (2.0, 0.0), (2.0, 2.0)]); let p = Point::new(2.0, 0.0); - assert!(t.contains(&p)); + assert!(!t.contains(&p)); } #[test] diff --git a/geo/src/algorithm/contains/point.rs b/geo/src/algorithm/contains/point.rs new file mode 100644 index 000000000..6e6210e30 --- /dev/null +++ b/geo/src/algorithm/contains/point.rs @@ -0,0 +1,37 @@ +use super::Contains; +use crate::*; + +// ┌────────────────────────────────┐ +// │ Implementations for Point │ +// └────────────────────────────────┘ + +impl Contains> for Point +where + T: CoordinateType, +{ + fn contains(&self, coord: &Coordinate) -> bool { + &self.0 == coord + } +} + +impl Contains> for Point +where + T: CoordinateType, +{ + fn contains(&self, p: &Point) -> bool { + self.contains(&p.0) + } +} + +// ┌────────────────────────────────┐ +// │ Implementations for MultiPoint │ +// └────────────────────────────────┘ +impl Contains for MultiPoint +where + T: CoordinateType, + Point: Contains, +{ + fn contains(&self, rhs: &G) -> bool { + self.0.iter().any(|p| p.contains(rhs)) + } +} diff --git a/geo/src/algorithm/contains/polygon.rs b/geo/src/algorithm/contains/polygon.rs new file mode 100644 index 000000000..b0abd6ab6 --- /dev/null +++ b/geo/src/algorithm/contains/polygon.rs @@ -0,0 +1,92 @@ +use super::Contains; +use crate::intersects::Intersects; +use crate::kernels::*; +use crate::*; + +// ┌─────────────────────────────┐ +// │ Implementations for Polygon │ +// └─────────────────────────────┘ + +impl Contains> for Polygon +where + T: HasKernel, +{ + fn contains(&self, coord: &Coordinate) -> bool { + match utils::coord_pos_relative_to_ring(*coord, &self.exterior()) { + utils::CoordPos::OnBoundary | utils::CoordPos::Outside => false, + _ => self.interiors().iter().all(|ls| { + utils::coord_pos_relative_to_ring(*coord, ls) == utils::CoordPos::Outside + }), + } + } +} + +impl Contains> for Polygon +where + T: HasKernel, +{ + fn contains(&self, p: &Point) -> bool { + self.contains(&p.0) + } +} + +// TODO: ensure DE-9IM compliance: esp., when +// line.start and line.end is on the boundaries +impl Contains> for Polygon +where + T: HasKernel, +{ + fn contains(&self, line: &Line) -> bool { + // both endpoints are contained in the polygon and the line + // does NOT intersect the exterior or any of the interior boundaries + self.contains(&line.start) + && self.contains(&line.end) + && !self.exterior().intersects(line) + && !self.interiors().iter().any(|inner| inner.intersects(line)) + } +} + +// TODO: also check interiors +impl Contains> for Polygon +where + T: HasKernel, +{ + fn contains(&self, poly: &Polygon) -> bool { + // decompose poly's exterior ring into Lines, and check each for containment + poly.exterior().lines().all(|line| self.contains(&line)) + } +} + +// TODO: ensure DE-9IM compliance +impl Contains> for Polygon +where + T: HasKernel, +{ + fn contains(&self, linestring: &LineString) -> bool { + // All LineString points must be inside the Polygon + if linestring.points_iter().all(|point| self.contains(&point)) { + // The Polygon interior is allowed to intersect with the LineString + // but the Polygon's rings are not + !self + .interiors() + .iter() + .any(|ring| ring.intersects(linestring)) + } else { + false + } + } +} + +// ┌──────────────────────────────────┐ +// │ Implementations for MultiPolygon │ +// └──────────────────────────────────┘ +// TODO: ensure DE-9IM compliance +impl Contains for MultiPolygon +where + T: CoordinateType, + Polygon: Contains, +{ + fn contains(&self, rhs: &G) -> bool { + self.0.iter().any(|p| p.contains(rhs)) + } +} diff --git a/geo/src/algorithm/contains/rect.rs b/geo/src/algorithm/contains/rect.rs new file mode 100644 index 000000000..d32a0eab4 --- /dev/null +++ b/geo/src/algorithm/contains/rect.rs @@ -0,0 +1,41 @@ +use super::Contains; +use crate::*; + +// ┌──────────────────────────┐ +// │ Implementations for Rect │ +// └──────────────────────────┘ + +impl Contains> for Rect +where + T: CoordinateType, +{ + fn contains(&self, coord: &Coordinate) -> bool { + coord.x > self.min().x + && coord.x < self.max().x + && coord.y > self.min().y + && coord.y < self.max().y + } +} + +impl Contains> for Rect +where + T: CoordinateType, +{ + fn contains(&self, p: &Point) -> bool { + self.contains(&p.0) + } +} + +impl Contains> for Rect +where + T: CoordinateType, +{ + fn contains(&self, other: &Rect) -> bool { + // TODO: check for degenerate rectangle (which is a line or a point) + // All points of LineString must be in the polygon ? + self.min().x <= other.min().x + && self.max().x >= other.max().x + && self.min().y <= other.min().y + && self.max().y >= other.max().y + } +} diff --git a/geo/src/algorithm/contains/triangle.rs b/geo/src/algorithm/contains/triangle.rs new file mode 100644 index 000000000..a8cf84183 --- /dev/null +++ b/geo/src/algorithm/contains/triangle.rs @@ -0,0 +1,27 @@ +use super::Contains; +use crate::kernels::*; +use crate::*; + +// ┌──────────────────────────────┐ +// │ Implementations for Triangle │ +// └──────────────────────────────┘ + +impl Contains> for Triangle +where + T: HasKernel, +{ + fn contains(&self, coord: &Coordinate) -> bool { + let ls = LineString(vec![self.0, self.1, self.2, self.0]); + use utils::*; + coord_pos_relative_to_ring(*coord, &ls) == CoordPos::Inside + } +} + +impl Contains> for Triangle +where + T: HasKernel, +{ + fn contains(&self, point: &Point) -> bool { + self.contains(&point.0) + } +} diff --git a/geo/src/algorithm/convex_hull/graham.rs b/geo/src/algorithm/convex_hull/graham.rs index 1fd503920..186d8bd74 100644 --- a/geo/src/algorithm/convex_hull/graham.rs +++ b/geo/src/algorithm/convex_hull/graham.rs @@ -89,20 +89,13 @@ where mod test { use super::*; use crate::algorithm::is_convex::IsConvex; - use geo_types::CoordinateType; use std::fmt::Debug; - fn test_convexity(initial: &[(T, T)]) { + fn test_convexity(initial: &[(T, T)]) { let mut v: Vec<_> = initial .iter() .map(|e| Coordinate::from((e.0, e.1))) .collect(); let hull = graham_hull(&mut v, false); - eprintln!("Strict hull"); - for v in hull.0.iter() { - eprintln!("{:?}", v); - } - eprintln!("{}", hull.is_closed()); - eprintln!("{:?}", hull.convex_orientation(true, None)); assert!(hull.is_strictly_ccw_convex()); let hull = graham_hull(&mut v, true); assert!(hull.is_ccw_convex()); diff --git a/geo/src/algorithm/convex_hull/mod.rs b/geo/src/algorithm/convex_hull/mod.rs index 373f3c4fe..8ee8c7c22 100644 --- a/geo/src/algorithm/convex_hull/mod.rs +++ b/geo/src/algorithm/convex_hull/mod.rs @@ -117,10 +117,17 @@ where // Remove repeated points unless collinear points // are to be included. - let mut ls: LineString = points.iter().copied().collect(); + let mut ls: Vec> = points.iter().copied().collect(); if !include_on_hull { - ls.0.dedup(); + ls.dedup(); } + + // A linestring with a single point is invalid. + if ls.len() == 1 { + ls.push(ls[0]); + } + + let mut ls = LineString(ls); ls.close(); // Maintain the CCW invariance diff --git a/geo/src/algorithm/euclidean_distance.rs b/geo/src/algorithm/euclidean_distance.rs index 7a43e4c5a..81c314025 100644 --- a/geo/src/algorithm/euclidean_distance.rs +++ b/geo/src/algorithm/euclidean_distance.rs @@ -3,7 +3,7 @@ use crate::algorithm::euclidean_length::EuclideanLength; use crate::algorithm::intersects::Intersects; use crate::algorithm::polygon_distance_fast_path::*; use crate::kernels::*; -use crate::utils::{coord_pos_relative_to_line_string, CoordPos}; +use crate::utils::{coord_pos_relative_to_ring, CoordPos}; use crate::{ Line, LineString, MultiLineString, MultiPoint, MultiPolygon, Point, Polygon, Triangle, }; @@ -132,7 +132,7 @@ where impl EuclideanDistance> for Point where - T: Float, + T: Float + HasKernel, { /// Minimum distance from a Point to a Polygon fn euclidean_distance(&self, polygon: &Polygon) -> T { @@ -165,7 +165,7 @@ where impl EuclideanDistance> for Polygon where - T: Float, + T: Float + HasKernel, { /// Minimum distance from a Polygon to a Point fn euclidean_distance(&self, point: &Point) -> T { @@ -175,7 +175,7 @@ where impl EuclideanDistance> for Point where - T: Float, + T: Float + HasKernel, { /// Minimum distance from a Point to a MultiPolygon fn euclidean_distance(&self, mpolygon: &MultiPolygon) -> T { @@ -189,7 +189,7 @@ where impl EuclideanDistance> for MultiPolygon where - T: Float, + T: Float + HasKernel, { /// Minimum distance from a MultiPolygon to a Point fn euclidean_distance(&self, point: &Point) -> T { @@ -263,7 +263,7 @@ where /// LineString-LineString distance impl EuclideanDistance> for LineString where - T: Float + Signed + RTreeNum, + T: Float + HasKernel + Signed + RTreeNum, { fn euclidean_distance(&self, other: &LineString) -> T { if self.intersects(other) { @@ -280,9 +280,9 @@ where /// contain a point from the candidate Polygon's outer shell in their simple representations fn ring_contains_point(poly: &Polygon, p: Point) -> bool where - T: Float, + T: HasKernel, { - match coord_pos_relative_to_line_string(p.0, &poly.exterior()) { + match coord_pos_relative_to_ring(p.0, &poly.exterior()) { CoordPos::Inside => true, CoordPos::OnBoundary | CoordPos::Outside => false, } @@ -291,7 +291,7 @@ where /// LineString to Line impl EuclideanDistance> for LineString where - T: Float + FloatConst + Signed + RTreeNum, + T: Float + FloatConst + Signed + RTreeNum + HasKernel, { fn euclidean_distance(&self, other: &Line) -> T { self.lines().fold(Bounded::max_value(), |acc, line| { @@ -303,7 +303,7 @@ where /// Line to LineString impl EuclideanDistance> for Line where - T: Float + FloatConst + Signed + RTreeNum, + T: Float + FloatConst + Signed + RTreeNum + HasKernel, { fn euclidean_distance(&self, other: &LineString) -> T { other.euclidean_distance(self) @@ -313,7 +313,7 @@ where /// LineString to Polygon impl EuclideanDistance> for LineString where - T: Float + FloatConst + Signed + RTreeNum, + T: Float + FloatConst + Signed + RTreeNum + HasKernel, { fn euclidean_distance(&self, other: &Polygon) -> T { if self.intersects(other) || other.contains(self) { @@ -334,7 +334,7 @@ where /// Polygon to LineString distance impl EuclideanDistance> for Polygon where - T: Float + FloatConst + Signed + RTreeNum, + T: Float + FloatConst + Signed + RTreeNum + HasKernel, { fn euclidean_distance(&self, other: &LineString) -> T { other.euclidean_distance(self) @@ -344,7 +344,7 @@ where /// Line to MultiPolygon distance impl EuclideanDistance> for Line where - T: Float + FloatConst + Signed + RTreeNum, + T: Float + FloatConst + Signed + RTreeNum + HasKernel, { fn euclidean_distance(&self, mpolygon: &MultiPolygon) -> T { mpolygon @@ -358,7 +358,7 @@ where /// MultiPolygon to Line distance impl EuclideanDistance> for MultiPolygon where - T: Float + FloatConst + Signed + RTreeNum, + T: Float + FloatConst + Signed + RTreeNum + HasKernel, { fn euclidean_distance(&self, other: &Line) -> T { other.euclidean_distance(self) @@ -368,7 +368,7 @@ where /// Line to Line distance impl EuclideanDistance> for Line where - T: Float + FloatConst + Signed + RTreeNum, + T: Float + FloatConst + Signed + RTreeNum + HasKernel, { fn euclidean_distance(&self, other: &Line) -> T { if self.intersects(other) || self.contains(other) { @@ -386,7 +386,7 @@ where // Line to Polygon distance impl EuclideanDistance> for Line where - T: Float + Signed + RTreeNum + FloatConst, + T: Float + Signed + RTreeNum + FloatConst + HasKernel, { fn euclidean_distance(&self, other: &Polygon) -> T { if other.contains(self) || self.intersects(other) { @@ -420,7 +420,7 @@ where // Polygon to Line distance impl EuclideanDistance> for Polygon where - T: Float + FloatConst + Signed + RTreeNum, + T: Float + FloatConst + Signed + RTreeNum + HasKernel, { fn euclidean_distance(&self, other: &Line) -> T { other.euclidean_distance(self) @@ -468,7 +468,7 @@ where impl EuclideanDistance> for Triangle where - T: Float, + T: Float + HasKernel, { fn euclidean_distance(&self, point: &Point) -> T { if self.contains(point) { @@ -659,6 +659,7 @@ mod test { // A point inside the cutout triangle let p = Point::new(3.5, 2.5); let dist = p.euclidean_distance(&poly); + // 0.41036467732879783 <-- Shapely assert_relative_eq!(dist, 0.41036467732879767); } diff --git a/geo/src/algorithm/extremes.rs b/geo/src/algorithm/extremes.rs index 98c27993b..fef8c126c 100644 --- a/geo/src/algorithm/extremes.rs +++ b/geo/src/algorithm/extremes.rs @@ -1,20 +1,7 @@ use super::kernels::*; use crate::prelude::*; use crate::*; -use num_traits::{Signed, Zero}; - -/// Compute the sign of the dot product of `u` and `v` using -/// robust predicates. The output is `CounterClockwise` if -/// the sign is positive, `Clockwise` if negative, and -/// `Collinear` if zero. -fn dot_product_sign(u: Coordinate, v: Coordinate) -> Orientation -where - T: CoordinateType + Signed + HasKernel, -{ - let zero = Coordinate::zero(); - let vdash = Coordinate { x: -v.y, y: v.x }; - T::Ker::orient2d(zero, u, vdash) -} +use num_traits::Signed; // Useful direction vectors, aligned with x and y axes: // 1., 0. = largest x @@ -28,7 +15,7 @@ fn above(u: Coordinate, vi: Coordinate, vj: Coordinate) -> bool where T: CoordinateType + Signed + HasKernel, { - dot_product_sign(u, vi - vj) == Orientation::CounterClockwise + T::Ker::dot_product_sign(u, vi - vj) == Orientation::CounterClockwise } /// Predicate that returns `true` if `vi` is (strictly) below `vj` @@ -38,16 +25,15 @@ fn below(u: Coordinate, vi: Coordinate, vj: Coordinate) -> bool where T: CoordinateType + Signed + HasKernel, { - dot_product_sign(u, vi - vj) == Orientation::Clockwise + T::Ker::dot_product_sign(u, vi - vj) == Orientation::Clockwise } // wrapper for extreme-finding function fn find_extreme_indices(func: F, polygon: &Polygon) -> Result where - T: CoordinateType + HasKernel + Signed, + T: HasKernel + Signed, F: Fn(Coordinate, &Polygon) -> Result, { - use crate::is_convex::IsConvex; if !polygon.exterior().is_convex() { return Err(()); } @@ -68,7 +54,7 @@ where // u: a direction vector. We're using a point to represent this, which is a hack but works fine fn polymax_naive_indices(u: Coordinate, poly: &Polygon) -> Result where - T: CoordinateType + HasKernel + Signed, + T: HasKernel + Signed, { let vertices = &poly.exterior().0; let mut max: usize = 0; diff --git a/geo/src/algorithm/intersects/coordinate.rs b/geo/src/algorithm/intersects/coordinate.rs new file mode 100644 index 000000000..037832cc6 --- /dev/null +++ b/geo/src/algorithm/intersects/coordinate.rs @@ -0,0 +1,11 @@ +use super::Intersects; +use crate::*; + +impl Intersects> for Coordinate +where + T: CoordinateType, +{ + fn intersects(&self, rhs: &Coordinate) -> bool { + self == rhs + } +} diff --git a/geo/src/algorithm/intersects/line.rs b/geo/src/algorithm/intersects/line.rs new file mode 100644 index 000000000..b08a2c663 --- /dev/null +++ b/geo/src/algorithm/intersects/line.rs @@ -0,0 +1,56 @@ +use super::{point_in_rect, Intersects}; +use crate::kernels::*; +use crate::*; + +impl Intersects> for Line +where + T: HasKernel, +{ + fn intersects(&self, rhs: &Coordinate) -> bool { + T::Ker::orient2d(self.start, self.end, *rhs) == Orientation::Collinear + && point_in_rect(*rhs, self.start, self.end) + } +} +symmetric_intersects_impl!(Coordinate, Line, HasKernel); + +impl Intersects> for Line +where + T: HasKernel, +{ + fn intersects(&self, p: &Point) -> bool { + self.intersects(&p.0) + } +} +symmetric_intersects_impl!(Point, Line, HasKernel); + +impl Intersects> for Line +where + T: HasKernel, +{ + fn intersects(&self, line: &Line) -> bool { + if self.start == self.end { + return line.intersects(&self.start); + } + let check_1_1 = T::Ker::orient2d(self.start, self.end, line.start); + let check_1_2 = T::Ker::orient2d(self.start, self.end, line.end); + + if check_1_1 != check_1_2 { + let check_2_1 = T::Ker::orient2d(line.start, line.end, self.start); + let check_2_2 = T::Ker::orient2d(line.start, line.end, self.end); + + check_2_1 != check_2_2 + } else if check_1_1 == Orientation::Collinear { + // Special case: collinear line segments. + + // Equivalent to the point-line intersection + // impl., but removes the calls to the kernel + // predicates. + point_in_rect(line.start, self.start, self.end) + || point_in_rect(line.end, self.start, self.end) + || point_in_rect(self.end, line.start, line.end) + || point_in_rect(self.end, line.start, line.end) + } else { + false + } + } +} diff --git a/geo/src/algorithm/intersects/line_string.rs b/geo/src/algorithm/intersects/line_string.rs new file mode 100644 index 000000000..2cdfb2546 --- /dev/null +++ b/geo/src/algorithm/intersects/line_string.rs @@ -0,0 +1,53 @@ +use super::Intersects; +use crate::kernels::*; +use crate::*; + +impl Intersects> for LineString +where + T: HasKernel, +{ + fn intersects(&self, coord: &Coordinate) -> bool { + self.lines().any(|l| coord.intersects(&l)) + } +} +symmetric_intersects_impl!(Coordinate, LineString, HasKernel); + +impl Intersects> for LineString +where + T: HasKernel, +{ + fn intersects(&self, point: &Point) -> bool { + self.intersects(&point.0) + } +} +symmetric_intersects_impl!(Point, LineString, HasKernel); + +impl Intersects> for LineString +where + T: HasKernel, +{ + fn intersects(&self, line: &Line) -> bool { + self.lines().any(|l| line.intersects(&l)) + } +} +symmetric_intersects_impl!(Line, LineString, HasKernel); + +impl Intersects> for LineString +where + T: HasKernel, +{ + fn intersects(&self, linestring: &LineString) -> bool { + if self.0.is_empty() || linestring.0.is_empty() { + return false; + } + for a in self.lines() { + for b in linestring.lines() { + if a.intersects(&b) { + return true; + } + } + } + + false + } +} diff --git a/geo/src/algorithm/intersects.rs b/geo/src/algorithm/intersects/mod.rs similarity index 69% rename from geo/src/algorithm/intersects.rs rename to geo/src/algorithm/intersects/mod.rs index 55532634a..7c145f03e 100644 --- a/geo/src/algorithm/intersects.rs +++ b/geo/src/algorithm/intersects/mod.rs @@ -1,278 +1,104 @@ -use crate::algorithm::contains::Contains; -use crate::{Line, LineString, Point, Polygon, Rect}; -use num_traits::Float; - -/// Checks if the geometry A intersects the geometry B. - +use crate::*; + +/// Checks if the geometry Self intersects the geometry Rhs. +/// More formally, either boundary or interior of Self has +/// non-empty (set-theoretic) intersection with the boundary +/// or interior of Rhs. In other words, the [DE-9IM] +/// intersection matrix for (Self, Rhs) is _not_ `FF*FF****`. +/// +/// This predicate is symmetric: `a.intersects(b)` iff +/// `b.intersects(a)`. +/// +/// [DE-9IM]: https://en.wikipedia.org/wiki/DE-9IM +/// +/// # Examples +/// +/// ``` +/// use geo::algorithm::intersects::Intersects; +/// use geo::line_string; +/// +/// let line_string_a = line_string![ +/// (x: 3., y: 2.), +/// (x: 7., y: 6.), +/// ]; +/// +/// let line_string_b = line_string![ +/// (x: 3., y: 4.), +/// (x: 8., y: 4.), +/// ]; +/// +/// let line_string_c = line_string![ +/// (x: 9., y: 2.), +/// (x: 11., y: 5.), +/// ]; +/// +/// assert!(line_string_a.intersects(&line_string_b)); +/// assert!(!line_string_a.intersects(&line_string_c)); +/// ``` pub trait Intersects { - /// Checks if the geometry A intersects the geometry B. - /// - /// # Examples - /// - /// ``` - /// use geo::algorithm::intersects::Intersects; - /// use geo::line_string; - /// - /// let line_string_a = line_string![ - /// (x: 3., y: 2.), - /// (x: 7., y: 6.), - /// ]; - /// - /// let line_string_b = line_string![ - /// (x: 3., y: 4.), - /// (x: 8., y: 4.), - /// ]; - /// - /// let line_string_c = line_string![ - /// (x: 9., y: 2.), - /// (x: 11., y: 5.), - /// ]; - /// - /// assert!(line_string_a.intersects(&line_string_b)); - /// assert!(!line_string_a.intersects(&line_string_c)); - /// ``` fn intersects(&self, rhs: &Rhs) -> bool; } -impl Intersects> for Line -where - T: Float, -{ - fn intersects(&self, p: &Point) -> bool { - let tx = if self.dx() == T::zero() { - None - } else { - Some((p.x() - self.start.x) / self.dx()) - }; - let ty = if self.dy() == T::zero() { - None - } else { - Some((p.y() - self.start.y) / self.dy()) - }; - match (tx, ty) { - (None, None) => { - // Degenerate line - p.0 == self.start - } - (Some(t), None) => { - // Horizontal line - p.y() == self.start.y && T::zero() <= t && t <= T::one() - } - (None, Some(t)) => { - // Vertical line - p.x() == self.start.x && T::zero() <= t && t <= T::one() - } - (Some(t_x), Some(t_y)) => { - // All other lines - (t_x - t_y).abs() <= T::epsilon() && T::zero() <= t_x && t_x <= T::one() - } - } - } -} - -impl Intersects> for Point -where - T: Float, -{ - fn intersects(&self, line: &Line) -> bool { - line.intersects(self) - } -} - -impl Intersects> for Line -where - T: Float, -{ - fn intersects(&self, line: &Line) -> bool { - // Using Cramer's Rule: - // https://en.wikipedia.org/wiki/Intersection_%28Euclidean_geometry%29#Two_line_segments - let a1 = self.dx(); - let a2 = self.dy(); - let b1 = -line.dx(); - let b2 = -line.dy(); - let c1 = line.start.x - self.start.x; - let c2 = line.start.y - self.start.y; - - let d = a1 * b2 - a2 * b1; - if d == T::zero() { - let (self_start, self_end) = self.points(); - let (other_start, other_end) = line.points(); - // lines are parallel - // return true iff at least one endpoint intersects the other line - self_start.intersects(line) - || self_end.intersects(line) - || other_start.intersects(self) - || other_end.intersects(self) - } else { - let s = (c1 * b2 - c2 * b1) / d; - let t = (a1 * c2 - a2 * c1) / d; - (T::zero() <= s) && (s <= T::one()) && (T::zero() <= t) && (t <= T::one()) - } - } -} - -impl Intersects> for Line -where - T: Float, -{ - fn intersects(&self, linestring: &LineString) -> bool { - linestring.lines().any(|line| self.intersects(&line)) - } -} - -impl Intersects> for LineString -where - T: Float, -{ - fn intersects(&self, line: &Line) -> bool { - line.intersects(self) - } -} - -impl Intersects> for Line -where - T: Float, -{ - fn intersects(&self, p: &Polygon) -> bool { - p.exterior().intersects(self) - || p.interiors().iter().any(|inner| inner.intersects(self)) - || p.contains(&self.start_point()) - || p.contains(&self.end_point()) - } -} - -impl Intersects> for Polygon -where - T: Float, -{ - fn intersects(&self, line: &Line) -> bool { - line.intersects(self) - } -} - -impl Intersects> for LineString -where - T: Float, -{ - // See: https://github.com/brandonxiang/geojson-python-utils/blob/33b4c00c6cf27921fb296052d0c0341bd6ca1af2/geojson_utils.py - fn intersects(&self, linestring: &LineString) -> bool { - if self.0.is_empty() || linestring.0.is_empty() { - return false; - } - for a in self.lines() { - for b in linestring.lines() { - let u_b = b.dy() * a.dx() - b.dx() * a.dy(); - if u_b == T::zero() { - continue; - } - let ua_t = b.dx() * (a.start.y - b.start.y) - b.dy() * (a.start.x - b.start.x); - let ub_t = a.dx() * (a.start.y - b.start.y) - a.dy() * (a.start.x - b.start.x); - let u_a = ua_t / u_b; - let u_b = ub_t / u_b; - if (T::zero() <= u_a) - && (u_a <= T::one()) - && (T::zero() <= u_b) - && (u_b <= T::one()) - { - return true; - } +// Since `Intersects` is symmetric, we use a macro to +// implement `T: Intersects` if `S: Intersects` is +// available. As a convention, we provide explicit impl. +// whenever the Rhs is a "simpler geometry" than the target +// type, and use the macro for the reverse impl. +#[macro_use] +macro_rules! symmetric_intersects_impl { + ($t:ty, $k:ty, $bounds:tt $(+ $more_bounds:tt )*) => { + impl $crate::algorithm::intersects::Intersects<$k> for $t + where T: $bounds $(+ $more_bounds)* + { + fn intersects(&self, rhs: &$k) -> bool { + rhs.intersects(self) } - } - false - } -} -impl Intersects> for Polygon -where - T: Float, -{ - fn intersects(&self, linestring: &LineString) -> bool { - // line intersects inner or outer polygon edge - if self.exterior().intersects(linestring) - || self - .interiors() - .iter() - .any(|inner| inner.intersects(linestring)) - { - true - } else { - // or if it's contained in the polygon - linestring.points_iter().any(|point| self.contains(&point)) } - } -} - -impl Intersects> for LineString -where - T: Float, -{ - fn intersects(&self, polygon: &Polygon) -> bool { - polygon.intersects(self) - } + }; } -// helper function for intersection check +mod coordinate; +mod line; +mod line_string; +mod point; +mod polygon; +mod rect; +mod triangle; + +// Helper function to check value lies between min and max. +// Only makes sense if min <= max (or always false) +#[inline] fn value_in_range(value: T, min: T, max: T) -> bool where - T: Float, -{ - (value >= min) && (value <= max) -} - -impl Intersects> for Rect -where - T: Float, + T: std::cmp::PartialOrd, { - fn intersects(&self, bounding_rect: &Rect) -> bool { - let x_overlap = value_in_range(self.min().x, bounding_rect.min().x, bounding_rect.max().x) - || value_in_range(bounding_rect.min().x, self.min().x, self.max().x); - - let y_overlap = value_in_range(self.min().y, bounding_rect.min().y, bounding_rect.max().y) - || value_in_range(bounding_rect.min().y, self.min().y, self.max().y); - - x_overlap && y_overlap - } + value >= min && value <= max } -impl Intersects> for Rect +// Helper function to check value lies between two bounds, +// where the ordering of the bounds is not known +#[inline] +fn value_in_between(value: T, bound_1: T, bound_2: T) -> bool where - T: Float, + T: std::cmp::PartialOrd, { - fn intersects(&self, polygon: &Polygon) -> bool { - polygon.intersects(self) + if bound_1 < bound_2 { + value_in_range(value, bound_1, bound_2) + } else { + value_in_range(value, bound_2, bound_1) } } -impl Intersects> for Polygon +// Helper function to check point lies inside rect given by +// bounds. The first bound need not be min. +#[inline] +fn point_in_rect(value: Coordinate, bound_1: Coordinate, bound_2: Coordinate) -> bool where - T: Float, + T: CoordinateType, { - fn intersects(&self, bounding_rect: &Rect) -> bool { - let p = Polygon::new( - LineString::from(vec![ - (bounding_rect.min().x, bounding_rect.min().y), - (bounding_rect.min().x, bounding_rect.max().y), - (bounding_rect.max().x, bounding_rect.max().y), - (bounding_rect.max().x, bounding_rect.min().y), - (bounding_rect.min().x, bounding_rect.min().y), - ]), - vec![], - ); - self.intersects(&p) - } -} - -impl Intersects> for Polygon -where - T: Float, -{ - fn intersects(&self, polygon: &Polygon) -> bool { - // self intersects (or contains) any line in polygon - self.intersects(polygon.exterior()) || - polygon.interiors().iter().any(|inner_line_string| self.intersects(inner_line_string)) || - // self is contained inside polygon - polygon.intersects(self.exterior()) - } + value_in_between(value.x, bound_1.x, bound_2.x) + && value_in_between(value.y, bound_1.y, bound_2.y) } #[cfg(test)] diff --git a/geo/src/algorithm/intersects/point.rs b/geo/src/algorithm/intersects/point.rs new file mode 100644 index 000000000..2d0a2f952 --- /dev/null +++ b/geo/src/algorithm/intersects/point.rs @@ -0,0 +1,37 @@ +use super::Intersects; +use crate::kernels::*; +use crate::*; + +impl Intersects> for Point +where + T: CoordinateType, +{ + fn intersects(&self, rhs: &Coordinate) -> bool { + self.0.intersects(rhs) + } +} +symmetric_intersects_impl!(Coordinate, Point, CoordinateType); + +impl Intersects> for Point +where + T: CoordinateType, +{ + fn intersects(&self, rhs: &Point) -> bool { + self.intersects(&rhs.0) + } +} + +impl Intersects for MultiPoint +where + T: CoordinateType, + Point: Intersects, +{ + fn intersects(&self, rhs: &G) -> bool { + self.0.iter().any(|p| p.intersects(rhs)) + } +} +symmetric_intersects_impl!(Coordinate, MultiPoint, CoordinateType); +symmetric_intersects_impl!(Point, MultiPoint, CoordinateType); +symmetric_intersects_impl!(Line, MultiPoint, HasKernel); +symmetric_intersects_impl!(LineString, MultiPoint, HasKernel); +symmetric_intersects_impl!(Polygon, MultiPoint, HasKernel); diff --git a/geo/src/algorithm/intersects/polygon.rs b/geo/src/algorithm/intersects/polygon.rs new file mode 100644 index 000000000..088620af3 --- /dev/null +++ b/geo/src/algorithm/intersects/polygon.rs @@ -0,0 +1,86 @@ +use super::Intersects; +use crate::contains::Contains; +use crate::kernels::*; +use crate::utils::{coord_pos_relative_to_ring, CoordPos}; +use crate::*; + +impl Intersects> for Polygon +where + T: HasKernel, +{ + fn intersects(&self, p: &Coordinate) -> bool { + coord_pos_relative_to_ring(*p, &self.exterior()) != CoordPos::Outside + && self + .interiors() + .iter() + .all(|int| coord_pos_relative_to_ring(*p, int) != CoordPos::Inside) + } +} +symmetric_intersects_impl!(Coordinate, Polygon, HasKernel); + +impl Intersects> for Polygon +where + T: HasKernel, +{ + fn intersects(&self, p: &Point) -> bool { + self.intersects(&p.0) + } +} +symmetric_intersects_impl!(Point, Polygon, HasKernel); + +impl Intersects> for Polygon +where + T: HasKernel, +{ + fn intersects(&self, line: &Line) -> bool { + self.exterior().intersects(line) + || self.interiors().iter().any(|inner| inner.intersects(line)) + || self.contains(&line.start) + || self.contains(&line.end) + } +} +symmetric_intersects_impl!(Line, Polygon, HasKernel); + +impl Intersects> for Polygon +where + T: HasKernel, +{ + fn intersects(&self, linestring: &LineString) -> bool { + // line intersects inner or outer polygon edge + if self.exterior().intersects(linestring) + || self + .interiors() + .iter() + .any(|inner| inner.intersects(linestring)) + { + true + } else { + // or if it's contained in the polygon + linestring.0.iter().any(|c| self.contains(c)) + } + } +} +symmetric_intersects_impl!(LineString, Polygon, HasKernel); + +impl Intersects> for Polygon +where + T: HasKernel, +{ + fn intersects(&self, rect: &Rect) -> bool { + self.intersects(&rect.clone().to_polygon()) + } +} +symmetric_intersects_impl!(Rect, Polygon, HasKernel); + +impl Intersects> for Polygon +where + T: HasKernel, +{ + fn intersects(&self, polygon: &Polygon) -> bool { + // self intersects (or contains) any line in polygon + self.intersects(polygon.exterior()) || + polygon.interiors().iter().any(|inner_line_string| self.intersects(inner_line_string)) || + // self is contained inside polygon + polygon.intersects(self.exterior()) + } +} diff --git a/geo/src/algorithm/intersects/rect.rs b/geo/src/algorithm/intersects/rect.rs new file mode 100644 index 000000000..8078be842 --- /dev/null +++ b/geo/src/algorithm/intersects/rect.rs @@ -0,0 +1,17 @@ +use super::{value_in_range, Intersects}; +use crate::*; + +impl Intersects> for Rect +where + T: CoordinateType, +{ + fn intersects(&self, other: &Rect) -> bool { + let x_overlap = value_in_range(self.min().x, other.min().x, other.max().x) + || value_in_range(other.min().x, self.min().x, self.max().x); + + let y_overlap = value_in_range(self.min().y, other.min().y, other.max().y) + || value_in_range(other.min().y, self.min().y, self.max().y); + + x_overlap && y_overlap + } +} diff --git a/geo/src/algorithm/intersects/triangle.rs b/geo/src/algorithm/intersects/triangle.rs new file mode 100644 index 000000000..2f2d7fe60 --- /dev/null +++ b/geo/src/algorithm/intersects/triangle.rs @@ -0,0 +1,18 @@ +use super::Intersects; +use crate::kernels::*; +use crate::*; + +impl Intersects for Triangle +where + T: CoordinateType, + Polygon: Intersects, +{ + fn intersects(&self, rhs: &G) -> bool { + self.clone().to_polygon().intersects(rhs) + } +} +symmetric_intersects_impl!(Coordinate, Triangle, HasKernel); +symmetric_intersects_impl!(Point, Triangle, HasKernel); +symmetric_intersects_impl!(Line, Triangle, HasKernel); +symmetric_intersects_impl!(LineString, Triangle, HasKernel); +symmetric_intersects_impl!(Polygon, Triangle, HasKernel); diff --git a/geo/src/algorithm/is_convex.rs b/geo/src/algorithm/is_convex.rs index c02c998f4..5159e16fa 100644 --- a/geo/src/algorithm/is_convex.rs +++ b/geo/src/algorithm/is_convex.rs @@ -1,5 +1,5 @@ use crate::kernels::*; -use crate::{Coordinate, CoordinateType, LineString}; +use crate::{Coordinate, LineString}; /// Predicates to test the convexity of a [ `LineString` ]. /// A closed `LineString` is said to be _convex_ if it @@ -108,7 +108,7 @@ pub trait IsConvex { fn is_collinear(&self) -> bool; } -impl IsConvex for LineString { +impl IsConvex for LineString { fn convex_orientation( &self, allow_collinear: bool, @@ -140,7 +140,7 @@ fn is_convex_shaped( specific_orientation: Option, ) -> Option where - T: CoordinateType + HasKernel, + T: HasKernel, { let n = coords.len(); diff --git a/geo/src/algorithm/kernels/mod.rs b/geo/src/algorithm/kernels/mod.rs index fb56baf80..7d8a95c2f 100644 --- a/geo/src/algorithm/kernels/mod.rs +++ b/geo/src/algorithm/kernels/mod.rs @@ -1,4 +1,5 @@ use crate::{Coordinate, CoordinateType}; +use num_traits::Zero; #[derive(Debug, PartialEq, Eq, PartialOrd, Ord, Clone, Copy)] pub enum Orientation { @@ -9,18 +10,11 @@ pub enum Orientation { /// Kernel trait to provide predicates to operate on /// different scalar types. -pub trait Kernel { - type Scalar: CoordinateType; - +pub trait Kernel { /// Gives the orientation of 3 2-dimensional points: /// ccw, cw or collinear (None) - fn orient2d( - p: Coordinate, - q: Coordinate, - r: Coordinate, - ) -> Orientation { + fn orient2d(p: Coordinate, q: Coordinate, r: Coordinate) -> Orientation { let res = (q.x - p.x) * (r.y - q.y) - (q.y - p.y) * (r.x - q.x); - use num_traits::Zero; if res > Zero::zero() { Orientation::CounterClockwise } else if res < Zero::zero() { @@ -30,17 +24,27 @@ pub trait Kernel { } } - fn square_euclidean_distance( - p: Coordinate, - q: Coordinate, - ) -> Self::Scalar { + fn square_euclidean_distance(p: Coordinate, q: Coordinate) -> T { (p.x - q.x) * (p.x - q.x) + (p.y - q.y) * (p.y - q.y) } + + /// Compute the sign of the dot product of `u` and `v` using + /// robust predicates. The output is `CounterClockwise` if + /// the sign is positive, `Clockwise` if negative, and + /// `Collinear` if zero. + fn dot_product_sign(u: Coordinate, v: Coordinate) -> Orientation { + let zero = Coordinate::zero(); + let vdash = Coordinate { + x: T::zero() - v.y, + y: v.x, + }; + Self::orient2d(zero, u, vdash) + } } /// Marker trait to assign Kernel for scalars pub trait HasKernel: CoordinateType { - type Ker: Kernel; + type Ker: Kernel; } // Helper macro to implement `HasKernel` on a a scalar type @@ -51,7 +55,7 @@ pub trait HasKernel: CoordinateType { macro_rules! has_kernel { ($t:ident, $k:ident) => { impl $crate::algorithm::kernels::HasKernel for $t { - type Ker = $k<$t>; + type Ker = $k; } }; } diff --git a/geo/src/algorithm/kernels/robust.rs b/geo/src/algorithm/kernels/robust.rs index b622a44dd..a2b9f9830 100644 --- a/geo/src/algorithm/kernels/robust.rs +++ b/geo/src/algorithm/kernels/robust.rs @@ -1,6 +1,5 @@ use super::{Kernel, Orientation}; use crate::Coordinate; -use std::marker::PhantomData; /// Robust kernel that uses [fast robust /// predicates](//www.cs.cmu.edu/~quake/robust.html) to @@ -8,17 +7,11 @@ use std::marker::PhantomData; /// used with types that can _always_ be casted to `f64` /// _without loss in precision_. #[derive(Default)] -pub struct RobustKernel(PhantomData); +pub struct RobustKernel; use num_traits::{Float, NumCast}; -impl Kernel for RobustKernel { - type Scalar = T; - - fn orient2d( - p: Coordinate, - q: Coordinate, - r: Coordinate, - ) -> Orientation { +impl Kernel for RobustKernel { + fn orient2d(p: Coordinate, q: Coordinate, r: Coordinate) -> Orientation { use robust::{orient2d, Coord}; let orientation = orient2d( diff --git a/geo/src/algorithm/kernels/simple.rs b/geo/src/algorithm/kernels/simple.rs index a50c08a59..340097eae 100644 --- a/geo/src/algorithm/kernels/simple.rs +++ b/geo/src/algorithm/kernels/simple.rs @@ -1,13 +1,10 @@ use super::Kernel; use crate::CoordinateType; -use std::marker::PhantomData; /// Simple kernel provides the direct implementation of the /// predicates. These are meant to be used with exact /// arithmetic signed tpyes (eg. i32, i64). #[derive(Default)] -pub struct SimpleKernel(PhantomData); +pub struct SimpleKernel; -impl Kernel for SimpleKernel { - type Scalar = T; -} +impl Kernel for SimpleKernel {} diff --git a/geo/src/algorithm/orient.rs b/geo/src/algorithm/orient.rs index b41a12472..a4fb8ef50 100644 --- a/geo/src/algorithm/orient.rs +++ b/geo/src/algorithm/orient.rs @@ -1,5 +1,5 @@ use super::kernels::*; -use crate::{CoordinateType, MultiPolygon, Polygon}; +use crate::{MultiPolygon, Polygon}; use crate::algorithm::winding_order::{Winding, WindingOrder}; @@ -68,7 +68,7 @@ pub trait Orient { impl Orient for Polygon where - T: CoordinateType + HasKernel, + T: HasKernel, { fn orient(&self, direction: Direction) -> Polygon { orient(self, direction) @@ -77,7 +77,7 @@ where impl Orient for MultiPolygon where - T: CoordinateType + HasKernel, + T: HasKernel, { fn orient(&self, direction: Direction) -> MultiPolygon { MultiPolygon(self.0.iter().map(|poly| poly.orient(direction)).collect()) @@ -100,7 +100,7 @@ pub enum Direction { // and the interior ring(s) will be oriented clockwise fn orient(poly: &Polygon, direction: Direction) -> Polygon where - T: CoordinateType + HasKernel, + T: HasKernel, { let interiors = poly .interiors() diff --git a/geo/src/algorithm/winding_order.rs b/geo/src/algorithm/winding_order.rs index a4c66a29a..3aa37d590 100644 --- a/geo/src/algorithm/winding_order.rs +++ b/geo/src/algorithm/winding_order.rs @@ -91,8 +91,8 @@ pub trait Winding { impl Winding for LineString where - T: CoordinateType + HasKernel, - K: Kernel, + T: HasKernel, + K: Kernel, { type Scalar = T; diff --git a/geo/src/lib.rs b/geo/src/lib.rs index 44e454e60..ef459400e 100644 --- a/geo/src/lib.rs +++ b/geo/src/lib.rs @@ -105,6 +105,7 @@ pub mod prelude { pub use crate::algorithm::haversine_intermediate::HaversineIntermediate; pub use crate::algorithm::haversine_length::HaversineLength; pub use crate::algorithm::intersects::Intersects; + pub use crate::algorithm::is_convex::IsConvex; pub use crate::algorithm::map_coords::MapCoords; pub use crate::algorithm::orient::Orient; #[cfg(feature = "use-proj")] diff --git a/geo/src/utils.rs b/geo/src/utils.rs index 96bf415a6..62bca87e4 100644 --- a/geo/src/utils.rs +++ b/geo/src/utils.rs @@ -1,7 +1,8 @@ //! Internal utility functions, types, and data structures. -use crate::contains::Contains; -use geo_types::{Coordinate, CoordinateType}; +use crate::intersects::Intersects; +use crate::kernels::*; +use geo_types::{Coordinate, CoordinateType, Line}; /// Partition a mutable slice in-place so that it contains all elements for /// which `predicate(e)` is `true`, followed by all elements for which @@ -83,42 +84,81 @@ pub enum CoordPos { Outside, } -/// Calculate the position of a `Coordinate` relative to a `LineString` -pub fn coord_pos_relative_to_line_string( +/// Calculate the position of a `Coordinate` relative to a +/// closed `LineString`. +pub fn coord_pos_relative_to_ring( coord: crate::Coordinate, linestring: &crate::LineString, ) -> CoordPos where - T: num_traits::Float, + T: HasKernel, { - // See: http://www.ecse.rpi.edu/Homepages/wrf/Research/Short_Notes/pnpoly.html - // http://geospatialpython.com/search - // ?updated-min=2011-01-01T00:00:00-06:00&updated-max=2012-01-01T00:00:00-06:00&max-results=19 + // Use the ray-tracing algorithm: count #times a + // horizontal ray from point (to positive infinity). + // + // See: https://en.wikipedia.org/wiki/Point_in_polygon + + assert!(linestring.is_closed()); // LineString without points if linestring.0.is_empty() { return CoordPos::Outside; } - // Point is on linestring - if linestring.contains(&coord) { - return CoordPos::OnBoundary; + if linestring.0.len() == 1 { + // If LineString has one point, it will not generate + // any lines. So, we handle this edge case separately. + return if coord == linestring.0[0] { + CoordPos::OnBoundary + } else { + CoordPos::Outside + }; } - let mut xints = T::zero(); let mut crossings = 0; for line in linestring.lines() { - if coord.y > line.start.y.min(line.end.y) - && coord.y <= line.start.y.max(line.end.y) - && coord.x <= line.start.x.max(line.end.x) - { - if line.start.y != line.end.y { - xints = (coord.y - line.start.y) * (line.end.x - line.start.x) - / (line.end.y - line.start.y) - + line.start.x; - } - if (line.start.x == line.end.x) || (coord.x <= xints) { - crossings += 1; - } + // Check if coord lies on the line + if line.intersects(&coord) { + return CoordPos::OnBoundary; + } + + // Ignore if the line is strictly to the left of the coord. + let max_x = if line.start.x < line.end.x { + line.end.x + } else { + line.start.x + }; + if max_x < coord.x { + continue; + } + + // Ignore if line is horizontal. This includes an + // edge case where the ray would intersect a + // horizontal segment of the ring infinitely many + // times, and is irrelevant for the calculation. + if line.start.y == line.end.y { + continue; + } + + // Ignore if the intersection of the line is + // possibly at the beginning of the line. This is to + // prevent a double counting when the ray passes + // through a vertex of the plygon + if line.start.y == coord.y { + continue; + } + + // Otherwise, check if ray intersects the line + // segment. Enough to consider ray upto the max_x + // coordinate of the current segment. + let ray = Line::new( + coord, + Coordinate { + x: max_x, + y: coord.y, + }, + ); + if ray.intersects(&line) { + crossings += 1; } } if crossings % 2 == 1 {