From 63013dc2949586512d946d92b4e68dabefee8bda Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Stephan=20H=C3=BCgel?= Date: Mon, 27 Mar 2017 21:48:54 +0100 Subject: [PATCH 1/6] Correctly calculate centroids for complex polygons (i.e. with interior holes) --- src/algorithm/centroid.rs | 114 ++++++++++++++++++++++++++++++-------- 1 file changed, 91 insertions(+), 23 deletions(-) diff --git a/src/algorithm/centroid.rs b/src/algorithm/centroid.rs index ed30afc17..22855751a 100644 --- a/src/algorithm/centroid.rs +++ b/src/algorithm/centroid.rs @@ -23,21 +23,47 @@ pub trait Centroid { fn centroid(&self) -> Option>; } +// Calculation of simple (no interior holes) Polygon area +fn simple_polygon_area(linestring: &LineString) -> T where T: Float { + if linestring.0.is_empty() || linestring.0.len() == 1 { + return T::zero(); + } + let mut tmp = T::zero(); + for ps in linestring.0.windows(2) { + tmp = tmp + (ps[0].x() * ps[1].y() - ps[1].x() * ps[0].y()); + } + tmp / (T::one() + T::one()) +} + +// Calculation of a Polygon centroid without interior rings +fn simple_polygon_centroid(poly_ext: &LineString) -> Option> + where T: Float + FromPrimitive +{ + let vect = &poly_ext.0; + let area = simple_polygon_area(poly_ext); + let mut sum_x = T::zero(); + let mut sum_y = T::zero(); + for ps in vect.windows(2) { + let tmp = ps[0].x() * ps[1].y() - ps[1].x() * ps[0].y(); + sum_x = sum_x + ((ps[1].x() + ps[0].x()) * tmp); + sum_y = sum_y + ((ps[1].y() + ps[0].y()) * tmp); + } + let six = T::from_i32(6).unwrap(); + Some(Point::new(sum_x / (six * area), sum_y / (six * area))) +} + impl Centroid for LineString where T: Float { - /// - /// The Centroid of a LineString is the mean of the middle of the segment - /// weighted by the length of the segments. - /// + // The Centroid of a LineString is the mean of the middle of the segment + // weighted by the length of the segments. fn centroid(&self) -> Option> { let vect = &self.0; if vect.is_empty() { return None; } if vect.len() == 1 { - Some(Point::new(vect[0].x(), - vect[0].y())) + Some(Point::new(vect[0].x(), vect[0].y())) } else { let mut sum_x = T::zero(); let mut sum_y = T::zero(); @@ -57,12 +83,15 @@ impl Centroid for LineString impl Centroid for Polygon where T: Float + FromPrimitive { - /// - /// Centroid on a Polygon. - /// See: https://en.wikipedia.org/wiki/Centroid - /// + // Centroid of a Polygon. + // See: https://en.wikipedia.org/wiki/Centroid + // We distinguish between a simple polygon, which has no interior holes, + // and a complex polygon, which has one or more interior holes. + // A complex polygon's centroid is the weighted average of its + // exterior shell centroid and the centroids of the interior ring(s), + // which are treated as *polygons* for the purposes of this calculation. + // See here for a formula: http://math.stackexchange.com/a/623849 fn centroid(&self) -> Option> { - // TODO: consideration of inner polygons; let linestring = &self.exterior; let vect = &linestring.0; if vect.is_empty() { @@ -71,16 +100,27 @@ impl Centroid for Polygon if vect.len() == 1 { Some(Point::new(vect[0].x(), vect[0].y())) } else { - let area = self.area(); - let mut sum_x = T::zero(); - let mut sum_y = T::zero(); - for ps in vect.windows(2) { - let tmp = ps[0].x() * ps[1].y() - ps[1].x() * ps[0].y(); - sum_x = sum_x + ((ps[1].x() + ps[0].x()) * tmp); - sum_y = sum_y + ((ps[1].y() + ps[0].y()) * tmp); + let external_centroid = simple_polygon_centroid(&self.exterior).unwrap(); + if !self.interiors.is_empty() { + let external_area = simple_polygon_area(&self.exterior).abs(); + // accumulate interior Polygons + let (totals_x, totals_y, internal_area) = + self.interiors + .iter() + .map(|ring| { + let area = simple_polygon_area(ring).abs(); + let centroid = simple_polygon_centroid(ring).unwrap(); + ((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 centroid_x = ((external_centroid.x() * external_area) - totals_x) / + (external_area - internal_area); + let centroid_y = ((external_centroid.y() * external_area) - totals_y) / + (external_area - internal_area); + return Some(Point::new(centroid_x, centroid_y)); } - let six = T::from_i32(6).unwrap(); - Some(Point::new(sum_x / (six * area), sum_y / (six * area))) + Some(external_centroid) } } } @@ -118,7 +158,7 @@ impl Centroid for Bbox /// fn centroid(&self) -> Option> { let two = T::one() + T::one(); - Some(Point::new((self.xmax + self.xmin)/two, (self.ymax + self.ymin)/two)) + Some(Point::new((self.xmax + self.xmin) / two, (self.ymax + self.ymin) / two)) } } @@ -176,6 +216,26 @@ mod test { let poly = Polygon::new(linestring, v); assert_eq!(poly.centroid(), Some(p(1., 1.))); } + #[test] + fn polygon_hole_test() { + let ls1 = LineString(vec![Point::new(5.0, 1.0), + Point::new(4.0, 2.0), + Point::new(4.0, 3.0), + Point::new(5.0, 4.0), + Point::new(6.0, 4.0), + Point::new(7.0, 3.0), + Point::new(7.0, 2.0), + Point::new(6.0, 1.0), + Point::new(5.0, 1.0)]); + + let ls2 = LineString(vec![Point::new(5.0, 1.3), + Point::new(5.5, 2.0), + Point::new(6.0, 1.3), + Point::new(5.0, 1.3)]); + let p1 = Polygon::new(ls1, vec![ls2]); + let centroid = p1.centroid().unwrap(); + assert_eq!(centroid, Point::new(5.5, 2.550877192982456)); + } /// Tests: Centroid of MultiPolygon #[test] fn empty_multipolygon_polygon_test() { @@ -195,7 +255,10 @@ mod test { let poly1 = Polygon::new(linestring, Vec::new()); let linestring = LineString(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().distance(&p(4.07142857142857, 1.92857142857143)); + let dist = MultiPolygon(vec![poly1, poly2]) + .centroid() + .unwrap() + .distance(&p(4.07142857142857, 1.92857142857143)); assert!(dist < COORD_PRECISION); } #[test] @@ -209,7 +272,12 @@ mod test { } #[test] fn bbox_test() { - let bbox = Bbox{ xmax: 4., xmin: 0., ymax: 100., ymin: 50.}; + let bbox = Bbox { + xmax: 4., + xmin: 0., + ymax: 100., + ymin: 50., + }; let point = Point(Coordinate { x: 2., y: 75. }); assert_eq!(point, bbox.centroid().unwrap()); } From b688345c6d46c31a356d27e71b11a111549c0764 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Stephan=20H=C3=BCgel?= Date: Wed, 29 Mar 2017 23:14:05 +0100 Subject: [PATCH 2/6] Add second interior ring --- src/algorithm/centroid.rs | 14 +++++++++++--- 1 file changed, 11 insertions(+), 3 deletions(-) diff --git a/src/algorithm/centroid.rs b/src/algorithm/centroid.rs index 22855751a..46f7ad102 100644 --- a/src/algorithm/centroid.rs +++ b/src/algorithm/centroid.rs @@ -24,7 +24,9 @@ pub trait Centroid { } // Calculation of simple (no interior holes) Polygon area -fn simple_polygon_area(linestring: &LineString) -> T where T: Float { +fn simple_polygon_area(linestring: &LineString) -> T + where T: Float +{ if linestring.0.is_empty() || linestring.0.len() == 1 { return T::zero(); } @@ -232,9 +234,15 @@ mod test { Point::new(5.5, 2.0), Point::new(6.0, 1.3), Point::new(5.0, 1.3)]); - let p1 = Polygon::new(ls1, vec![ls2]); + + let ls3 = LineString(vec![Point::new(5., 2.3), + Point::new(5.5, 3.0), + Point::new(6., 2.3), + Point::new(5., 2.3)]); + + let p1 = Polygon::new(ls1, vec![ls2, ls3]); let centroid = p1.centroid().unwrap(); - assert_eq!(centroid, Point::new(5.5, 2.550877192982456)); + assert_eq!(centroid, Point::new(5.5, 2.5518518518518514)); } /// Tests: Centroid of MultiPolygon #[test] From 0ce5cf2b9ec38fd3abd5380e6c789b527a9e59c1 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Stephan=20H=C3=BCgel?= Date: Thu, 30 Mar 2017 15:25:15 +0100 Subject: [PATCH 3/6] No need to assign to intermediate variables --- src/algorithm/centroid.rs | 9 ++++----- 1 file changed, 4 insertions(+), 5 deletions(-) diff --git a/src/algorithm/centroid.rs b/src/algorithm/centroid.rs index 46f7ad102..f751871cc 100644 --- a/src/algorithm/centroid.rs +++ b/src/algorithm/centroid.rs @@ -116,11 +116,10 @@ impl Centroid for Polygon }) .fold((T::zero(), T::zero(), T::zero()), |accum, val| (accum.0 + val.0, accum.1 + val.1, accum.2 + val.2)); - let centroid_x = ((external_centroid.x() * external_area) - totals_x) / - (external_area - internal_area); - let centroid_y = ((external_centroid.y() * external_area) - totals_y) / - (external_area - internal_area); - return Some(Point::new(centroid_x, centroid_y)); + return Some(Point::new(((external_centroid.x() * external_area) - totals_x) / + (external_area - internal_area), + ((external_centroid.y() * external_area) - totals_y) / + (external_area - internal_area))); } Some(external_centroid) } From ca25bd1cfd51436c3a63fc680ed50ae242aa66e4 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Stephan=20H=C3=BCgel?= Date: Fri, 31 Mar 2017 15:37:53 +0100 Subject: [PATCH 4/6] Remove unused import, and fix example Showing a numerical result is preferable to (NaN, NaN) --- src/algorithm/centroid.rs | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/src/algorithm/centroid.rs b/src/algorithm/centroid.rs index f751871cc..3520b51ca 100644 --- a/src/algorithm/centroid.rs +++ b/src/algorithm/centroid.rs @@ -9,15 +9,15 @@ pub trait Centroid { /// Calculation of the centroid, see: https://en.wikipedia.org/wiki/Centroid /// /// ``` - /// use geo::{Point, LineString, Coordinate}; + /// use geo::{Point, LineString}; /// use geo::algorithm::centroid::Centroid; /// /// let mut vec = Vec::new(); /// vec.push(Point::new(40.02f64, 116.34)); - /// vec.push(Point::new(40.02f64, 116.34)); + /// vec.push(Point::new(40.02f64, 118.23)); /// let linestring = LineString(vec); /// - /// println!("Centroid {:?}", linestring.centroid()); + /// assert_eq!(linestring.centroid().unwrap(), Point::new(40.02, 117.285)); /// ``` /// fn centroid(&self) -> Option>; From dbfaddfbec98b87b090a6226e67f43fbd11a13cf Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Stephan=20H=C3=BCgel?= Date: Fri, 31 Mar 2017 15:45:45 +0100 Subject: [PATCH 5/6] Tidy up docs and comments --- src/algorithm/centroid.rs | 19 ++++++++----------- 1 file changed, 8 insertions(+), 11 deletions(-) diff --git a/src/algorithm/centroid.rs b/src/algorithm/centroid.rs index 3520b51ca..dc8e8fab0 100644 --- a/src/algorithm/centroid.rs +++ b/src/algorithm/centroid.rs @@ -6,7 +6,7 @@ use algorithm::distance::Distance; /// Calculation of the centroid. pub trait Centroid { - /// Calculation of the centroid, see: https://en.wikipedia.org/wiki/Centroid + /// See: https://en.wikipedia.org/wiki/Centroid /// /// ``` /// use geo::{Point, LineString}; @@ -85,14 +85,15 @@ impl Centroid for LineString impl Centroid for Polygon where T: Float + FromPrimitive { - // Centroid of a Polygon. - // See: https://en.wikipedia.org/wiki/Centroid + // Calculate the centroid of a Polygon. // We distinguish between a simple polygon, which has no interior holes, // and a complex polygon, which has one or more interior holes. // A complex polygon's centroid is the weighted average of its // exterior shell centroid and the centroids of the interior ring(s), - // which are treated as *polygons* for the purposes of this calculation. + // which are both considered as as * 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: See: https://fotino.me/calculating-centroids/ fn centroid(&self) -> Option> { let linestring = &self.exterior; let vect = &linestring.0; @@ -129,7 +130,6 @@ impl Centroid for Polygon impl Centroid for MultiPolygon where T: Float + FromPrimitive { - // See: https://fotino.me/calculating-centroids/ fn centroid(&self) -> Option> { let mut sum_x = T::zero(); let mut sum_y = T::zero(); @@ -154,9 +154,6 @@ impl Centroid for MultiPolygon impl Centroid for Bbox where T: Float { - /// - /// Calculate the Centroid of a Bbox. - /// fn centroid(&self) -> Option> { let two = T::one() + T::one(); Some(Point::new((self.xmax + self.xmin) / two, (self.ymax + self.ymin) / two)) @@ -168,7 +165,7 @@ mod test { use types::{COORD_PRECISION, Coordinate, Point, LineString, Polygon, MultiPolygon, Bbox}; use algorithm::centroid::Centroid; use algorithm::distance::Distance; - /// Tests: Centroid of LineString + // Tests: Centroid of LineString #[test] fn empty_linestring_test() { let vec = Vec::>::new(); @@ -192,7 +189,7 @@ mod test { assert_eq!(linestring.centroid(), Some(Point(Coordinate { x: 6., y: 1. }))); } - /// Tests: Centroid of Polygon + // Tests: Centroid of Polygon #[test] fn empty_polygon_test() { let v1 = Vec::new(); @@ -243,7 +240,7 @@ mod test { let centroid = p1.centroid().unwrap(); assert_eq!(centroid, Point::new(5.5, 2.5518518518518514)); } - /// Tests: Centroid of MultiPolygon + // Tests: Centroid of MultiPolygon #[test] fn empty_multipolygon_polygon_test() { assert!(MultiPolygon::(Vec::new()).centroid().is_none()); From ad6f0c555ad1154fdcdbbbe874d34ff04e77c913 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Stephan=20H=C3=BCgel?= Date: Sun, 2 Apr 2017 16:59:37 +0100 Subject: [PATCH 6/6] Typo fix --- src/algorithm/centroid.rs | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/algorithm/centroid.rs b/src/algorithm/centroid.rs index dc8e8fab0..864d5a711 100644 --- a/src/algorithm/centroid.rs +++ b/src/algorithm/centroid.rs @@ -90,10 +90,10 @@ impl Centroid for Polygon // and a complex polygon, which has one or more interior holes. // A complex polygon's centroid is the weighted average of its // exterior shell centroid and the centroids of the interior ring(s), - // which are both considered as as * simple polygons* for the purposes of + // which are both 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: See: https://fotino.me/calculating-centroids/ + // See here for detail on alternative methods: https://fotino.me/calculating-centroids/ fn centroid(&self) -> Option> { let linestring = &self.exterior; let vect = &linestring.0;