diff --git a/Cargo.toml b/Cargo.toml index f791f902c..b1a72af0c 100644 --- a/Cargo.toml +++ b/Cargo.toml @@ -1,2 +1,8 @@ [workspace] members = ["geo", "geo-types", "geo-postgis"] + +[patch.crates-io] + +# Ensure any transitive dependencies also use the local geo/geo-types +geo = { path = "geo" } +geo-types = { path = "geo-types" } diff --git a/geo/CHANGES.md b/geo/CHANGES.md index 64db57d40..2c3abe63b 100644 --- a/geo/CHANGES.md +++ b/geo/CHANGES.md @@ -8,6 +8,14 @@ * Rewrite the crate documentation * +* Fix `Centroid` algorithm for `MultiLineString` when all members have only one + point. + * +* Implement `Centroid` algorithm on `Geometry` and its remaining variants. + * + +## 0.17.1 + * Add `GeodesicIntermediate` algorithm * diff --git a/geo/Cargo.toml b/geo/Cargo.toml index ad9e8e0b9..35804e4ab 100644 --- a/geo/Cargo.toml +++ b/geo/Cargo.toml @@ -22,7 +22,7 @@ geographiclib-rs = { version = "0.2" } proj = { version = "0.20.3", optional = true } -geo-types = { version = "0.7.0", path = "../geo-types", features = ["approx", "use-rstar"] } +geo-types = { version = "0.7.0", features = ["approx", "use-rstar"] } robust = { version = "0.2.2" } @@ -34,6 +34,9 @@ use-serde = ["serde", "geo-types/serde"] [dev-dependencies] approx = "0.4.0" criterion = { version = "0.3" } +# jts-test-runner is an internal crate which exists only to be part of the geo test suite. +# As such it's kept unpublished. It's in a separate repo primarily because it's kind of large. +jts-test-runner = { git = "https://github.com/georust/jts-test-runner", rev = "3294b9af1d3e64fcc9caf9646a93d0116f7d8321" } rand = "0.8.0" [[bench]] diff --git a/geo/fuzz/Cargo.toml b/geo/fuzz/Cargo.toml index 94472d1bc..2cf35a1d1 100644 --- a/geo/fuzz/Cargo.toml +++ b/geo/fuzz/Cargo.toml @@ -12,10 +12,8 @@ cargo-fuzz = true libfuzzer-sys = "0.4" [dependencies.geo] -path = ".." [dependencies.geo-types] -path = "../../geo-types" features = ["arbitrary"] # Prevent this from interfering with workspaces @@ -27,3 +25,7 @@ name = "simplify" path = "fuzz_targets/simplify.rs" test = false doc = false + +[patch.crates-io] +geo = { path = ".." } +geo-types = { path = "../../geo-types" } diff --git a/geo/src/algorithm/centroid.rs b/geo/src/algorithm/centroid.rs index 73a6c10fa..4d4392df3 100644 --- a/geo/src/algorithm/centroid.rs +++ b/geo/src/algorithm/centroid.rs @@ -1,10 +1,11 @@ -use num_traits::FromPrimitive; -use std::iter::Sum; +use std::cmp::Ordering; use crate::algorithm::area::{get_linestring_area, Area}; +use crate::algorithm::dimensions::{Dimensions, Dimensions::*, HasDimensions}; use crate::algorithm::euclidean_length::EuclideanLength; use crate::{ - CoordFloat, Line, LineString, MultiLineString, MultiPoint, MultiPolygon, Point, Polygon, Rect, + Coordinate, GeoFloat, Geometry, GeometryCollection, Line, LineString, MultiLineString, + MultiPoint, MultiPolygon, Point, Polygon, Rect, Triangle, }; /// Calculation of the centroid. @@ -58,250 +59,77 @@ pub trait Centroid { fn centroid(&self) -> Self::Output; } -// Calculation of a Polygon centroid without interior rings -fn simple_polygon_centroid(poly_ext: &LineString) -> Option> -where - T: CoordFloat + FromPrimitive + Sum, -{ - let area = get_linestring_area(poly_ext); - if area == T::zero() { - // if the polygon is flat (area = 0), it is considered as a linestring - return poly_ext.centroid(); - } - - // At this point, we know the exterior contains at least one point - let shift = poly_ext.0.first().unwrap(); - - let (sum_x, sum_y) = poly_ext - .lines() - .fold((T::zero(), T::zero()), |accum, line| { - use crate::algorithm::map_coords::MapCoords; - let line = line.map_coords(|&(x, y)| (x - shift.x, y - shift.y)); - let tmp = line.determinant(); - ( - accum.0 + ((line.end.x + line.start.x) * tmp), - accum.1 + ((line.end.y + line.start.y) * tmp), - ) - }); - let six = T::from_i32(6).unwrap(); - Some(Point::new( - sum_x / (six * area) + shift.x, - sum_y / (six * area) + shift.y, - )) -} - impl Centroid for Line where - T: CoordFloat, + T: GeoFloat, { type Output = Point; fn centroid(&self) -> Self::Output { let two = T::one() + T::one(); - let x = self.start.x + self.dx() / two; - let y = self.start.y + self.dy() / two; - Point::new(x, y) + (self.start_point() + self.end_point()) / two } } impl Centroid for LineString where - T: CoordFloat, + T: GeoFloat, { type Output = Option>; // The Centroid of a LineString is the mean of the middle of the segment // weighted by the length of the segments. fn centroid(&self) -> Self::Output { - if self.0.is_empty() { - return None; - } - if self.0.len() == 1 { - Some(Point(self.0[0])) - } else { - let (sum_x, sum_y, total_length) = - self.lines() - .fold((T::zero(), T::zero(), T::zero()), |accum, line| { - let segment_len = line.euclidean_length(); - let line_center = line.centroid(); - ( - accum.0 + segment_len * line_center.x(), - accum.1 + segment_len * line_center.y(), - accum.2 + segment_len, - ) - }); - if total_length == T::zero() { - // length == 0 means that all points were equal, we can just the first one - Some(Point(self.0[0])) - } else { - Some(Point::new(sum_x / total_length, sum_y / total_length)) - } - } + let mut operation = CentroidOperation::new(); + operation.add_line_string(self); + operation.centroid() } } impl Centroid for MultiLineString where - T: CoordFloat + FromPrimitive + Sum, + T: GeoFloat, { type Output = Option>; /// The Centroid of a MultiLineString is the mean of the centroids of all the constituent linestrings, /// weighted by the length of each linestring fn centroid(&self) -> Self::Output { - if self.0.is_empty() || self.iter().all(|ls| ls.0.is_empty()) { - return None; - } - if self.0.len() == 1 { - self.0[0].centroid() - } else { - let (sum_x, sum_y, total_length) = - self.0 - .iter() - .fold( - (T::zero(), T::zero(), T::zero()), - |accum, line| match line.centroid() { - Some(center) => { - let segment_len = line.euclidean_length(); - ( - accum.0 + segment_len * center.x(), - accum.1 + segment_len * center.y(), - accum.2 + segment_len, - ) - } - None => accum, - }, - ); - if total_length == T::zero() { - // length == 0 means that all points in all constituent linestrings were equal. - // we can just use the first defined point in this case. - for linestring in self.0.iter() { - if !linestring.0.is_empty() { - return Some(Point(linestring[0])); - } - } - None // this should never happen, since all linestrings being empty was previously checked - } else { - Some(Point::new(sum_x / total_length, sum_y / total_length)) - } - } + let mut operation = CentroidOperation::new(); + operation.add_multi_line_string(self); + operation.centroid() } } impl Centroid for Polygon where - T: CoordFloat + FromPrimitive + Sum, + T: GeoFloat, { type Output = Option>; - // Calculate the centroid of a Polygon. - // We distinguish between a simple polygon, which has no interior rings (holes), - // and a complex polygon, which has one or more interior rings. - // A complex polygon's centroid is the weighted average of its - // exterior shell centroid and the centroids of the interior ring(s). - // Both the shell and the ring(s) are considered simple polygons for the purposes of - // this calculation. - // See here for a formula: http://math.stackexchange.com/a/623849 - // See here for detail on alternative methods: https://fotino.me/calculating-centroids/ fn centroid(&self) -> Self::Output { - let linestring = &self.exterior(); - let vect = &linestring.0; - if vect.is_empty() { - return None; - } - if vect.len() == 1 { - Some(Point::new(vect[0].x, vect[0].y)) - } else { - let external_centroid = simple_polygon_centroid(self.exterior())?; - if self.interiors().is_empty() { - Some(external_centroid) - } else { - let external_area = get_linestring_area(self.exterior()).abs(); - // accumulate interior Polygons - let (totals_x, totals_y, internal_area) = self - .interiors() - .iter() - .filter_map(|ring| { - let area = get_linestring_area(ring).abs(); - let centroid = simple_polygon_centroid(ring)?; - Some((centroid.x() * area, centroid.y() * area, area)) - }) - .fold((T::zero(), T::zero(), T::zero()), |accum, val| { - (accum.0 + val.0, accum.1 + val.1, accum.2 + val.2) - }); - - let diff_area = external_area - internal_area; - if diff_area == T::zero() { - Some(external_centroid) - } else { - Some(Point::new( - ((external_centroid.x() * external_area) - totals_x) / diff_area, - ((external_centroid.y() * external_area) - totals_y) / diff_area, - )) - } - } - } + let mut operation = CentroidOperation::new(); + operation.add_polygon(self); + operation.centroid() } } impl Centroid for MultiPolygon where - T: CoordFloat + FromPrimitive + Sum, + T: GeoFloat, { type Output = Option>; fn centroid(&self) -> Self::Output { - let mut sum_area_x = T::zero(); - let mut sum_area_y = T::zero(); - let mut sum_seg_x = T::zero(); - let mut sum_seg_y = T::zero(); - let mut sum_x = T::zero(); - let mut sum_y = T::zero(); - let mut total_area = T::zero(); - let mut total_length = T::zero(); - let vect = &self.0; - if vect.is_empty() { - return None; - } - for poly in &self.0 { - let area = poly.unsigned_area(); - total_area = total_area + area; - if let Some(p) = poly.centroid() { - if area != T::zero() { - sum_area_x = sum_area_x + area * p.x(); - sum_area_y = sum_area_y + area * p.y(); - } else { - // the polygon is 'flat', we consider it as a linestring - let ls_len = poly.exterior().euclidean_length(); - if ls_len == T::zero() { - sum_x = sum_x + p.x(); - sum_y = sum_y + p.x(); - } else { - sum_seg_x = sum_seg_x + ls_len * p.x(); - sum_seg_y = sum_seg_y + ls_len * p.y(); - total_length = total_length + ls_len; - } - } - } - } - if total_area != T::zero() { - Some(Point::new(sum_area_x / total_area, sum_area_y / total_area)) - } else if total_length != T::zero() { - Some(Point::new( - sum_seg_x / total_length, - sum_seg_y / total_length, - )) - } else { - let nb_points = T::from_usize(self.0.len()).unwrap(); - // there was only "point" polygons, we do a simple centroid of all points - Some(Point::new(sum_x / nb_points, sum_y / nb_points)) - } + let mut operation = CentroidOperation::new(); + operation.add_multi_polygon(self); + operation.centroid() } } impl Centroid for Rect where - T: CoordFloat, + T: GeoFloat, { type Output = Point; @@ -312,12 +140,12 @@ where impl Centroid for Point where - T: CoordFloat, + T: GeoFloat, { type Output = Point; fn centroid(&self) -> Self::Output { - Point::new(self.x(), self.y()) + *self } } @@ -335,42 +163,318 @@ where /// ``` impl Centroid for MultiPoint where - T: CoordFloat, + T: GeoFloat, +{ + type Output = Option>; + + fn centroid(&self) -> Self::Output { + let mut operation = CentroidOperation::new(); + operation.add_multi_point(self); + operation.centroid() + } +} + +impl Centroid for Geometry +where + T: GeoFloat, +{ + type Output = Option>; + + crate::geometry_delegate_impl! { + fn centroid(&self) -> Self::Output; + } +} + +impl Centroid for GeometryCollection +where + T: GeoFloat, { type Output = Option>; fn centroid(&self) -> Self::Output { - if self.0.is_empty() { - return None; + let mut operation = CentroidOperation::new(); + operation.add_geometry_collection(self); + operation.centroid() + } +} + +impl Centroid for Triangle +where + T: GeoFloat, +{ + type Output = Point; + + fn centroid(&self) -> Self::Output { + let mut operation = CentroidOperation::new(); + operation.add_triangle(self); + operation + .centroid() + .expect("triangle cannot have an empty centroid") + } +} + +struct CentroidOperation(Option>); +impl CentroidOperation { + fn new() -> Self { + CentroidOperation(None) + } + + fn centroid(&self) -> Option> { + self.0.as_ref().map(|weighted_centroid| { + Point(weighted_centroid.accumulated / weighted_centroid.weight) + }) + } + + fn centroid_dimensions(&self) -> Dimensions { + self.0 + .as_ref() + .map(|weighted_centroid| weighted_centroid.dimensions) + .unwrap_or(Empty) + } + + fn add_coord(&mut self, coord: Coordinate) { + self.add_centroid(ZeroDimensional, coord, T::one()); + } + + fn add_line(&mut self, line: &Line) { + match line.dimensions() { + ZeroDimensional => self.add_coord(line.start), + OneDimensional => { + self.add_centroid(OneDimensional, line.centroid().0, line.euclidean_length()) + } + _ => unreachable!("Line must be zero or one dimensional"), + } + } + + fn add_line_string(&mut self, line_string: &LineString) { + if self.centroid_dimensions() > OneDimensional { + return; + } + + if line_string.0.len() == 1 { + self.add_coord(line_string.0[0]); + return; + } + + for line in line_string.lines() { + self.add_line(&line); + } + } + + fn add_multi_line_string(&mut self, multi_line_string: &MultiLineString) { + if self.centroid_dimensions() > OneDimensional { + return; + } + + for element in &multi_line_string.0 { + self.add_line_string(element); + } + } + + fn add_polygon(&mut self, polygon: &Polygon) { + // Polygons which are completely covered by their interior rings have zero area, and + // represent a unique degeneracy into a line_string which cannot be handled by accumulating + // directly into `self`. Instead, we perform a sub-operation, inspect the result, and only + // then incorporate the result into `self. + + let mut exterior_operation = CentroidOperation::new(); + exterior_operation.add_ring(polygon.exterior()); + + let mut interior_operation = CentroidOperation::new(); + for interior in polygon.interiors() { + interior_operation.add_ring(interior); + } + + if let Some(exterior_weighted_centroid) = exterior_operation.0 { + let mut poly_weighted_centroid = exterior_weighted_centroid; + if let Some(interior_weighted_centroid) = interior_operation.0 { + poly_weighted_centroid.sub_assign(interior_weighted_centroid); + if poly_weighted_centroid.weight.is_zero() { + // A polygon with no area `interiors` completely covers `exterior`, degenerating to a linestring + self.add_line_string(polygon.exterior()); + return; + } + } + self.add_weighted_centroid(poly_weighted_centroid); + } + } + + fn add_multi_point(&mut self, multi_point: &MultiPoint) { + if self.centroid_dimensions() > ZeroDimensional { + return; + } + + for element in &multi_point.0 { + self.add_coord(element.0); + } + } + + fn add_multi_polygon(&mut self, multi_polygon: &MultiPolygon) { + for element in &multi_polygon.0 { + self.add_polygon(element); + } + } + + fn add_geometry_collection(&mut self, geometry_collection: &GeometryCollection) { + for element in &geometry_collection.0 { + self.add_geometry(element); + } + } + + fn add_rect(&mut self, rect: &Rect) { + match rect.dimensions() { + ZeroDimensional => self.add_coord(rect.min()), + OneDimensional => { + // Degenerate rect is a line, treat it the same way we treat flat polygons + self.add_line(&Line::new(rect.min(), rect.min())); + self.add_line(&Line::new(rect.min(), rect.max())); + self.add_line(&Line::new(rect.max(), rect.max())); + self.add_line(&Line::new(rect.max(), rect.min())); + } + TwoDimensional => { + self.add_centroid(TwoDimensional, rect.centroid().0, rect.unsigned_area()) + } + Empty => unreachable!("Rect dimensions cannot be empty"), + } + } + + fn add_triangle(&mut self, triangle: &Triangle) { + match triangle.dimensions() { + ZeroDimensional => self.add_coord(triangle.0), + OneDimensional => { + // Degenerate triangle is a line, treat it the same way we treat flat + // polygons + let l0_1 = Line::new(triangle.0, triangle.1); + let l1_2 = Line::new(triangle.1, triangle.2); + let l2_0 = Line::new(triangle.2, triangle.0); + self.add_line(&l0_1); + self.add_line(&l1_2); + self.add_line(&l2_0); + } + TwoDimensional => { + let centroid = (triangle.0 + triangle.1 + triangle.2) / T::from(3).unwrap(); + self.add_centroid(TwoDimensional, centroid, triangle.unsigned_area()); + } + Empty => unreachable!("Rect dimensions cannot be empty"), + } + } + + fn add_geometry(&mut self, geometry: &Geometry) { + match geometry { + Geometry::Point(g) => self.add_coord(g.0), + Geometry::Line(g) => self.add_line(g), + Geometry::LineString(g) => self.add_line_string(g), + Geometry::Polygon(g) => self.add_polygon(g), + Geometry::MultiPoint(g) => self.add_multi_point(g), + Geometry::MultiLineString(g) => self.add_multi_line_string(g), + Geometry::MultiPolygon(g) => self.add_multi_polygon(g), + Geometry::GeometryCollection(g) => self.add_geometry_collection(g), + Geometry::Rect(g) => self.add_rect(g), + Geometry::Triangle(g) => self.add_triangle(g), + } + } + + fn add_ring(&mut self, ring: &LineString) { + debug_assert!(ring.is_closed()); + + let area = get_linestring_area(ring); + if area == T::zero() { + match ring.dimensions() { + // empty ring doesn't contribute to centroid + Empty => {} + // degenerate ring is a point + ZeroDimensional => self.add_coord(ring[0]), + // zero-area ring is a line string + _ => self.add_line_string(ring), + } + return; + } + + // Since area is non-zero, we know the ring has at least one point + let shift = ring.0[0]; + + let accumulated_coord = ring.lines().fold(Coordinate::zero(), |accum, line| { + use crate::algorithm::map_coords::MapCoords; + let line = line.map_coords(|&(x, y)| (x - shift.x, y - shift.y)); + let tmp = line.determinant(); + accum + (line.end + line.start) * tmp + }); + let six = T::from(6).unwrap(); + let centroid = accumulated_coord / (six * area) + shift; + let weight = area.abs(); + self.add_centroid(TwoDimensional, centroid, weight); + } + + fn add_centroid(&mut self, dimensions: Dimensions, centroid: Coordinate, weight: T) { + let weighted_centroid = WeightedCentroid { + dimensions, + weight: weight, + accumulated: centroid * weight, + }; + self.add_weighted_centroid(weighted_centroid); + } + + fn add_weighted_centroid(&mut self, other: WeightedCentroid) { + match self.0.as_mut() { + Some(centroid) => centroid.add_assign(other), + None => self.0 = Some(other), + } + } +} + +// Aggregated state for accumulating the centroid of a geometry or collection of geometries. +struct WeightedCentroid { + weight: T, + accumulated: Coordinate, + /// Collections of Geometries can have different dimensionality. Centroids must be considered + /// separately by dimensionality. + /// + /// e.g. If I have several Points, adding a new `Point` will affect their centroid. + /// + /// However, because a Point is zero dimensional, it is infinitely small when compared to + /// any 2-D Polygon. Thus a Point will not affect the centroid of any GeometryCollection + /// containing a 2-D Polygon. + /// + /// So, when accumulating a centroid, we must track the dimensionality of the centroid + dimensions: Dimensions, +} + +impl WeightedCentroid { + fn add_assign(&mut self, b: WeightedCentroid) { + match self.dimensions.cmp(&b.dimensions) { + Ordering::Less => *self = b, + Ordering::Greater => return, + Ordering::Equal => { + self.accumulated = self.accumulated + b.accumulated; + self.weight = self.weight + b.weight; + } + } + } + + fn sub_assign(&mut self, b: WeightedCentroid) { + match self.dimensions.cmp(&b.dimensions) { + Ordering::Less => *self = b, + Ordering::Greater => return, + Ordering::Equal => { + self.accumulated = self.accumulated - b.accumulated; + self.weight = self.weight - b.weight; + } } - let sum = self.iter().fold( - Point::new(T::zero(), T::zero()), - |a: Point, b: &Point| Point::new(a.x() + b.x(), a.y() + b.y()), - ); - Some(Point::new( - sum.x() / T::from(self.0.len()).unwrap(), - sum.y() / T::from(self.0.len()).unwrap(), - )) } } #[cfg(test)] mod test { - use crate::algorithm::centroid::Centroid; - use crate::algorithm::euclidean_distance::EuclideanDistance; - use crate::line_string; - use crate::{ - polygon, CoordFloat, Coordinate, Line, LineString, MultiLineString, MultiPolygon, Point, - Polygon, Rect, - }; + use super::*; + use crate::{line_string, point, polygon}; /// small helper to create a coordinate - fn c(x: T, y: T) -> Coordinate { + fn c(x: T, y: T) -> Coordinate { Coordinate { x, y } } /// small helper to create a point - fn p(x: T, y: T) -> Point { + fn p(x: T, y: T) -> Point { Point(c(x, y)) } @@ -401,10 +505,16 @@ mod test { (x: 10., y: 1.), (x: 11., y: 1.) ]; - assert_eq!( - linestring.centroid(), - Some(Point(Coordinate { x: 6., y: 1. })) - ); + assert_eq!(linestring.centroid(), Some(point!(x: 6., y: 1. ))); + } + #[test] + fn linestring_with_repeated_point_test() { + let l1 = LineString::from(vec![p(1., 1.), p(1., 1.), p(1., 1.)]); + assert_eq!(l1.centroid(), Some(p(1., 1.))); + + let l2 = LineString::from(vec![p(2., 2.), p(2., 2.), p(2., 2.)]); + let mls = MultiLineString(vec![l1, l2]); + assert_eq!(mls.centroid(), Some(p(1.5, 1.5))); } // Tests: Centroid of MultiLineString #[test] @@ -430,7 +540,7 @@ mod test { line_string![coord], line_string![coord], ]); - assert_eq!(mls.centroid(), Some(Point(coord))); + assert_relative_eq!(mls.centroid().unwrap(), Point(coord)); } #[test] fn multilinestring_one_line_test() { @@ -443,7 +553,7 @@ mod test { (x: 11., y: 1.) ]; let mls: MultiLineString = MultiLineString(vec![linestring]); - assert_eq!(mls.centroid(), Some(Point(Coordinate { x: 6., y: 1. }))); + assert_relative_eq!(mls.centroid().unwrap(), Point(Coordinate { x: 6., y: 1. })); } #[test] fn multilinestring_test() { @@ -451,12 +561,9 @@ mod test { let v2 = line_string![(x: 1.0, y: 10.0), (x: 2.0, y: 0.0), (x: 3.0, y: 1.0)]; let v3 = line_string![(x: -12.0, y: -100.0), (x: 7.0, y: 8.0)]; let mls = MultiLineString(vec![v1, v2, v3]); - assert_eq!( - mls.centroid(), - Some(Point(Coordinate { - x: -1.9097834383655845, - y: -37.683866439745714 - })) + assert_relative_eq!( + mls.centroid().unwrap(), + point![x: -1.9097834383655845, y: -37.683866439745714] ); } // Tests: Centroid of Polygon @@ -467,11 +574,11 @@ mod test { } #[test] fn polygon_one_point_test() { - let p = Point(Coordinate { x: 2., y: 1. }); + let p = point![ x: 2., y: 1. ]; let v = Vec::new(); let linestring = line_string![p.0]; let poly = Polygon::new(linestring, v); - assert_eq!(poly.centroid(), Some(p)); + assert_relative_eq!(poly.centroid().unwrap(), p); } #[test] @@ -523,10 +630,11 @@ mod test { (x: 0., y: 2.), (x: 0., y: 0.) ]; - assert_eq!(poly.centroid(), Some(Point::new(1., 1.))); + assert_relative_eq!(poly.centroid().unwrap(), point![x:1., y:1.]); } #[test] fn polygon_hole_test() { + // hexagon let ls1 = LineString::from(vec![ (5.0, 1.0), (4.0, 2.0), @@ -545,8 +653,7 @@ mod test { let p1 = Polygon::new(ls1, vec![ls2, ls3]); let centroid = p1.centroid().unwrap(); - assert_relative_eq!(centroid.x(), 5.5, max_relative = 1e-6); - assert_relative_eq!(centroid.y(), 2.5518518518518514, max_relative = 1e-6); + assert_relative_eq!(centroid, point!(x: 5.5, y: 2.5518518518518523)); } #[test] fn flat_polygon_test() { @@ -660,11 +767,11 @@ mod test { let linestring = LineString::from(vec![p(7., 1.), p(8., 1.), p(8., 2.), p(7., 2.), p(7., 1.)]); let poly2 = Polygon::new(linestring, Vec::new()); - let dist = MultiPolygon(vec![poly1, poly2]) - .centroid() - .unwrap() - .euclidean_distance(&p(4.07142857142857, 1.92857142857143)); - assert_relative_eq!(dist, 0.0, epsilon = 1e-14); // the results is larger than f64::EPSILON + let centroid = MultiPolygon(vec![poly1, poly2]).centroid().unwrap(); + assert_relative_eq!( + centroid, + point![x: 4.071428571428571, y: 1.9285714285714286] + ); } #[test] fn multipolygon_two_polygons_of_opposite_clockwise_test() { @@ -672,20 +779,133 @@ mod test { let poly1 = Polygon::new(linestring, Vec::new()); let linestring = LineString::from(vec![(0., 0.), (-2., 0.), (-2., 2.), (0., 2.), (0., 0.)]); let poly2 = Polygon::new(linestring, Vec::new()); - assert_eq!( - MultiPolygon(vec![poly1, poly2]).centroid(), - Some(Point::new(0., 1.)) + assert_relative_eq!( + MultiPolygon(vec![poly1, poly2]).centroid().unwrap(), + point![x: 0., y: 1.] ); } #[test] fn bounding_rect_test() { let bounding_rect = Rect::new(Coordinate { x: 0., y: 50. }, Coordinate { x: 4., y: 100. }); - let point = Point(Coordinate { x: 2., y: 75. }); + let point = point![x: 2., y: 75.]; assert_eq!(point, bounding_rect.centroid()); } #[test] fn line_test() { let line1 = Line::new(c(0., 1.), c(1., 3.)); - assert_eq!(line1.centroid(), Point::new(0.5, 2.)); + assert_eq!(line1.centroid(), point![x: 0.5, y: 2.]); + } + #[test] + fn collection_weighting() { + let p0 = point!(x: 0.0, y: 0.0); + let p1 = point!(x: 2.0, y: 0.0); + let p2 = point!(x: 2.0, y: 2.0); + let p3 = point!(x: 0.0, y: 2.0); + + let multi_point = MultiPoint(vec![p0, p1, p2, p3]); + assert_eq!(multi_point.centroid().unwrap(), point!(x: 1.0, y: 1.0)); + + let collection = GeometryCollection(vec![MultiPoint(vec![p1, p2, p3]).into(), p0.into()]); + + assert_eq!(collection.centroid().unwrap(), point!(x: 1.0, y: 1.0)); + } + #[test] + fn triangles() { + // boring triangle + assert_eq!( + Triangle(c(0., 0.), c(3., 0.), c(1.5, 3.)).centroid(), + point!(x: 1.5, y: 1.0) + ); + + // flat triangle + assert_eq!( + Triangle(c(0., 0.), c(3., 0.), c(1., 0.)).centroid(), + point!(x: 1.5, y: 0.0) + ); + + // flat triangle that's not axis-aligned + assert_eq!( + Triangle(c(0., 0.), c(3., 3.), c(1., 1.)).centroid(), + point!(x: 1.5, y: 1.5) + ); + + // triangle with some repeated points + assert_eq!( + Triangle(c(0., 0.), c(0., 0.), c(1., 0.)).centroid(), + point!(x: 0.5, y: 0.0) + ); + + // triangle with all repeated points + assert_eq!( + Triangle(c(0., 0.5), c(0., 0.5), c(0., 0.5)).centroid(), + point!(x: 0., y: 0.5) + ) + } + + #[test] + fn degenerate_triangle_like_ring() { + let triangle = Triangle(c(0., 0.), c(1., 1.), c(2., 2.)); + let poly: Polygon<_> = triangle.clone().into(); + + let line = Line::new(c(0., 1.), c(1., 3.)); + + let g1 = GeometryCollection(vec![triangle.into(), line.into()]); + let g2 = GeometryCollection(vec![poly.into(), line.into()]); + assert_eq!(g1.centroid(), g2.centroid()); + } + + #[test] + fn degenerate_rect_like_ring() { + let rect = Rect::new(c(0., 0.), c(0., 4.)); + let poly: Polygon<_> = rect.clone().into(); + + let line = Line::new(c(0., 1.), c(1., 3.)); + + let g1 = GeometryCollection(vec![rect.into(), line.into()]); + let g2 = GeometryCollection(vec![poly.into(), line.into()]); + assert_eq!(g1.centroid(), g2.centroid()); + } + + #[test] + fn rectangles() { + // boring rect + assert_eq!( + Rect::new(c(0., 0.), c(4., 4.)).centroid(), + point!(x: 2.0, y: 2.0) + ); + + // flat rect + assert_eq!( + Rect::new(c(0., 0.), c(4., 0.)).centroid(), + point!(x: 2.0, y: 0.0) + ); + + // rect with all repeated points + assert_eq!( + Rect::new(c(4., 4.), c(4., 4.)).centroid(), + point!(x: 4., y: 4.) + ); + + // collection with rect + let mut collection = + GeometryCollection(vec![p(0., 0.).into(), p(6., 0.).into(), p(6., 6.).into()]); + // sanity check + assert_eq!(collection.centroid().unwrap(), point!(x: 4., y: 2.)); + + // 0-d rect treated like point + collection.0.push(Rect::new(c(0., 6.), c(0., 6.)).into()); + assert_eq!(collection.centroid().unwrap(), point!(x: 3., y: 3.)); + + // 1-d rect treated like line. Since a line has higher dimensions than the rest of the + // collection, it's centroid clobbers everything else in the collection. + collection.0.push(Rect::new(c(0., 0.), c(0., 2.)).into()); + assert_eq!(collection.centroid().unwrap(), point!(x: 0., y: 1.)); + + // 2-d has higher dimensions than the rest of the collection, so it's centroid clobbers + // everything else in the collection. + collection + .0 + .push(Rect::new(c(10., 10.), c(11., 11.)).into()); + assert_eq!(collection.centroid().unwrap(), point!(x: 10.5, y: 10.5)); } } diff --git a/geo/src/algorithm/dimensions.rs b/geo/src/algorithm/dimensions.rs index cfd34013b..08f2fc969 100644 --- a/geo/src/algorithm/dimensions.rs +++ b/geo/src/algorithm/dimensions.rs @@ -1,5 +1,6 @@ +use crate::algorithm::kernels::Orientation::Collinear; use crate::{ - CoordNum, Geometry, GeometryCollection, Line, LineString, MultiLineString, MultiPoint, + CoordNum, GeoNum, Geometry, GeometryCollection, Line, LineString, MultiLineString, MultiPoint, MultiPolygon, Point, Polygon, Rect, Triangle, }; @@ -131,7 +132,7 @@ pub trait HasDimensions { fn boundary_dimensions(&self) -> Dimensions; } -impl HasDimensions for Geometry { +impl HasDimensions for Geometry { crate::geometry_delegate_impl! { fn is_empty(&self) -> bool; fn dimensions(&self) -> Dimensions; @@ -187,7 +188,6 @@ impl HasDimensions for LineString { return Dimensions::Empty; } - debug_assert!(self.0.len() > 1, "invalid line_string with 1 coord"); let first = self.0[0]; if self.0.iter().any(|&coord| first != coord) { Dimensions::OneDimensional @@ -226,7 +226,19 @@ impl HasDimensions for Polygon { } fn dimensions(&self) -> Dimensions { - Dimensions::TwoDimensional + use crate::algorithm::coords_iter::CoordsIter; + let mut coords = self.exterior_coords_iter(); + match coords.next() { + None => Dimensions::Empty, + Some(coord_0) => { + if coords.all(|coord_n| coord_0 == coord_n) { + // all coords are a single point + Dimensions::ZeroDimensional + } else { + Dimensions::TwoDimensional + } + } + } } fn boundary_dimensions(&self) -> Dimensions { @@ -309,7 +321,7 @@ impl HasDimensions for MultiPolygon { } } -impl HasDimensions for GeometryCollection { +impl HasDimensions for GeometryCollection { fn is_empty(&self) -> bool { if self.0.is_empty() { true @@ -375,18 +387,21 @@ impl HasDimensions for Rect { } } -impl HasDimensions for Triangle { +impl HasDimensions for Triangle { fn is_empty(&self) -> bool { false } fn dimensions(&self) -> Dimensions { - if self.0 == self.1 && self.1 == self.2 { - // degenerate triangle is a point - Dimensions::ZeroDimensional - } else if self.0 == self.1 || self.1 == self.2 || self.2 == self.0 { - // degenerate triangle is a line - Dimensions::OneDimensional + use crate::algorithm::kernels::Kernel; + if Collinear == C::Ker::orient2d(self.0, self.1, self.2) { + if self.0 == self.1 && self.1 == self.2 { + // degenerate triangle is a point + Dimensions::ZeroDimensional + } else { + // degenerate triangle is a line + Dimensions::OneDimensional + } } else { Dimensions::TwoDimensional } diff --git a/geo/src/algorithm/rotate.rs b/geo/src/algorithm/rotate.rs index 33fc1f681..c676c7c4e 100644 --- a/geo/src/algorithm/rotate.rs +++ b/geo/src/algorithm/rotate.rs @@ -1,10 +1,9 @@ use crate::algorithm::centroid::Centroid; use crate::algorithm::map_coords::MapCoords; use crate::{ - CoordFloat, Line, LineString, MultiLineString, MultiPoint, MultiPolygon, Point, Polygon, + CoordFloat, GeoFloat, Line, LineString, MultiLineString, MultiPoint, MultiPolygon, Point, + Polygon, }; -use num_traits::FromPrimitive; -use std::iter::Sum; #[inline] fn rotate_inner(x: T, y: T, x0: T, y0: T, sin_theta: T, cos_theta: T) -> Point @@ -146,7 +145,7 @@ where impl Rotate for Line where - T: CoordFloat, + T: GeoFloat, { fn rotate(&self, angle: T) -> Self { let centroid = self.centroid(); @@ -159,7 +158,7 @@ where impl Rotate for LineString where - T: CoordFloat, + T: GeoFloat, { /// Rotate the LineString about its centroid by the given number of degrees fn rotate(&self, angle: T) -> Self { @@ -169,7 +168,7 @@ where impl Rotate for Polygon where - T: CoordFloat + FromPrimitive + Sum, + T: GeoFloat, { /// Rotate the Polygon about its centroid by the given number of degrees fn rotate(&self, angle: T) -> Self { @@ -191,7 +190,7 @@ where impl Rotate for MultiPolygon where - T: CoordFloat + FromPrimitive + Sum, + T: GeoFloat, { /// Rotate the contained Polygons about their centroids by the given number of degrees fn rotate(&self, angle: T) -> Self { @@ -201,7 +200,7 @@ where impl Rotate for MultiLineString where - T: CoordFloat + FromPrimitive, + T: GeoFloat, { /// Rotate the contained LineStrings about their centroids by the given number of degrees fn rotate(&self, angle: T) -> Self { @@ -211,7 +210,7 @@ where impl Rotate for MultiPoint where - T: CoordFloat + FromPrimitive, + T: CoordFloat, { /// Rotate the contained Points about their centroids by the given number of degrees fn rotate(&self, angle: T) -> Self { diff --git a/geo/tests/jts_tests.rs b/geo/tests/jts_tests.rs new file mode 100644 index 000000000..14441936f --- /dev/null +++ b/geo/tests/jts_tests.rs @@ -0,0 +1,57 @@ +use jts_test_runner::TestRunner; + +#[test] +fn test_centroid() { + let mut runner = TestRunner::new().matching_filename_glob("*Centroid.xml"); + runner.run().expect("test cases failed"); + + // sanity check that *something* was run + assert!( + runner.failures().len() + runner.successes().len() > 0, + "No tests were run." + ); + + if !runner.failures().is_empty() { + let failure_text = runner + .failures() + .iter() + .map(|failure| format!("{}", failure)) + .collect::>() + .join("\n"); + panic!( + "{} failures / {} successes in JTS test suite:\n{}", + runner.failures().len(), + runner.successes().len(), + failure_text + ); + } +} + +#[test] +// several of the ConvexHull tests are currently failing +#[ignore] +fn test_all_general() { + let mut runner = TestRunner::new(); + runner.run().expect("test cases failed"); + + // sanity check that *something* was run + assert!( + runner.failures().len() + runner.successes().len() > 0, + "No tests were run." + ); + + if !runner.failures().is_empty() { + let failure_text = runner + .failures() + .iter() + .map(|failure| format!("{}", failure)) + .collect::>() + .join("\n"); + panic!( + "{} failures / {} successes in JTS test suite:\n{}", + runner.failures().len(), + runner.successes().len(), + failure_text + ); + } +}