geo/algorithm/
centroid.rs

1use std::cmp::Ordering;
2
3use crate::area::{get_linestring_area, Area};
4use crate::dimensions::{Dimensions, Dimensions::*, HasDimensions};
5use crate::geometry::*;
6use crate::EuclideanLength;
7use crate::GeoFloat;
8
9/// Calculation of the centroid.
10/// The centroid is the arithmetic mean position of all points in the shape.
11/// Informally, it is the point at which a cutout of the shape could be perfectly
12/// balanced on the tip of a pin.
13/// The geometric centroid of a convex object always lies in the object.
14/// A non-convex object might have a centroid that _is outside the object itself_.
15///
16/// # Examples
17///
18/// ```
19/// use geo::Centroid;
20/// use geo::{point, polygon};
21///
22/// // rhombus shaped polygon
23/// let polygon = polygon![
24///     (x: -2., y: 1.),
25///     (x: 1., y: 3.),
26///     (x: 4., y: 1.),
27///     (x: 1., y: -1.),
28///     (x: -2., y: 1.),
29/// ];
30///
31/// assert_eq!(
32///     Some(point!(x: 1., y: 1.)),
33///     polygon.centroid(),
34/// );
35/// ```
36pub trait Centroid {
37    type Output;
38
39    /// See: <https://en.wikipedia.org/wiki/Centroid>
40    ///
41    /// # Examples
42    ///
43    /// ```
44    /// use geo::Centroid;
45    /// use geo::{line_string, point};
46    ///
47    /// let line_string = line_string![
48    ///     (x: 40.02f64, y: 116.34),
49    ///     (x: 40.02f64, y: 118.23),
50    /// ];
51    ///
52    /// assert_eq!(
53    ///     Some(point!(x: 40.02, y: 117.285)),
54    ///     line_string.centroid(),
55    /// );
56    /// ```
57    fn centroid(&self) -> Self::Output;
58}
59
60impl<T> Centroid for Line<T>
61where
62    T: GeoFloat,
63{
64    type Output = Point<T>;
65
66    /// The Centroid of a [`Line`] is its middle point
67    ///
68    /// # Examples
69    ///
70    /// ```
71    /// use geo::Centroid;
72    /// use geo::{Line, point};
73    ///
74    /// let line = Line::new(
75    ///     point!(x: 1.0f64, y: 3.0),
76    ///     point!(x: 2.0f64, y: 4.0),
77    /// );
78    ///
79    /// assert_eq!(
80    ///     point!(x: 1.5, y: 3.5),
81    ///     line.centroid(),
82    /// );
83    /// ```
84    fn centroid(&self) -> Self::Output {
85        let two = T::one() + T::one();
86        (self.start_point() + self.end_point()) / two
87    }
88}
89
90impl<T> Centroid for LineString<T>
91where
92    T: GeoFloat,
93{
94    type Output = Option<Point<T>>;
95
96    // The Centroid of a [`LineString`] is the mean of the middle of the segment
97    // weighted by the length of the segments.
98    ///
99    /// # Examples
100    ///
101    /// ```
102    /// use geo::Centroid;
103    /// use geo::{line_string, point};
104    ///
105    /// let line_string = line_string![
106    ///   (x: 1.0f32, y: 1.0),
107    ///   (x: 2.0, y: 2.0),
108    ///   (x: 4.0, y: 4.0)
109    ///   ];
110    ///
111    /// assert_eq!(
112    ///     // (1.0 * (1.5, 1.5) + 2.0 * (3.0, 3.0)) / 3.0
113    ///     Some(point!(x: 2.5, y: 2.5)),
114    ///     line_string.centroid(),
115    /// );
116    /// ```
117    fn centroid(&self) -> Self::Output {
118        let mut operation = CentroidOperation::new();
119        operation.add_line_string(self);
120        operation.centroid()
121    }
122}
123
124impl<T> Centroid for MultiLineString<T>
125where
126    T: GeoFloat,
127{
128    type Output = Option<Point<T>>;
129
130    /// The Centroid of a [`MultiLineString`] is the mean of the centroids of all the constituent linestrings,
131    /// weighted by the length of each linestring
132    ///
133    /// # Examples
134    ///
135    /// ```
136    /// use geo::Centroid;
137    /// use geo::{MultiLineString, line_string, point};
138    ///
139    /// let multi_line_string = MultiLineString::new(vec![
140    ///     // centroid: (2.5, 2.5)
141    ///     line_string![(x: 1.0f32, y: 1.0), (x: 2.0, y: 2.0), (x: 4.0, y: 4.0)],
142    ///     // centroid: (4.0, 4.0)
143    ///     line_string![(x: 1.0, y: 1.0), (x: 3.0, y: 3.0), (x: 7.0, y: 7.0)],
144    /// ]);
145    ///
146    /// assert_eq!(
147    ///     // ( 3.0 * (2.5, 2.5) + 6.0 * (4.0, 4.0) ) / 9.0
148    ///     Some(point!(x: 3.5, y: 3.5)),
149    ///     multi_line_string.centroid(),
150    /// );
151    /// ```
152    fn centroid(&self) -> Self::Output {
153        let mut operation = CentroidOperation::new();
154        operation.add_multi_line_string(self);
155        operation.centroid()
156    }
157}
158
159impl<T> Centroid for Polygon<T>
160where
161    T: GeoFloat,
162{
163    type Output = Option<Point<T>>;
164
165    /// The Centroid of a [`Polygon`] is the mean of its points
166    ///
167    /// # Examples
168    ///
169    /// ```
170    /// use geo::Centroid;
171    /// use geo::{polygon, point};
172    ///
173    /// let polygon = polygon![
174    ///     (x: 0.0f32, y: 0.0),
175    ///     (x: 2.0, y: 0.0),
176    ///     (x: 2.0, y: 1.0),
177    ///     (x: 0.0, y: 1.0),
178    /// ];
179    ///
180    /// assert_eq!(
181    ///     Some(point!(x: 1.0, y: 0.5)),
182    ///     polygon.centroid(),
183    /// );
184    /// ```
185    fn centroid(&self) -> Self::Output {
186        let mut operation = CentroidOperation::new();
187        operation.add_polygon(self);
188        operation.centroid()
189    }
190}
191
192impl<T> Centroid for MultiPolygon<T>
193where
194    T: GeoFloat,
195{
196    type Output = Option<Point<T>>;
197
198    /// The Centroid of a [`MultiPolygon`] is the mean of the centroids of its polygons, weighted
199    /// by the area of the polygons
200    ///
201    /// # Examples
202    ///
203    /// ```
204    /// use geo::Centroid;
205    /// use geo::{MultiPolygon, polygon, point};
206    ///
207    /// let multi_polygon = MultiPolygon::new(vec![
208    ///   // centroid (1.0, 0.5)
209    ///   polygon![
210    ///     (x: 0.0f32, y: 0.0),
211    ///     (x: 2.0, y: 0.0),
212    ///     (x: 2.0, y: 1.0),
213    ///     (x: 0.0, y: 1.0),
214    ///   ],
215    ///   // centroid (-0.5, 0.0)
216    ///   polygon![
217    ///     (x: 1.0, y: 1.0),
218    ///     (x: -2.0, y: 1.0),
219    ///     (x: -2.0, y: -1.0),
220    ///     (x: 1.0, y: -1.0),
221    ///   ]
222    /// ]);
223    ///
224    /// assert_eq!(
225    ///     // ( 2.0 * (1.0, 0.5) + 6.0 * (-0.5, 0.0) ) / 8.0
226    ///     Some(point!(x: -0.125, y: 0.125)),
227    ///     multi_polygon.centroid(),
228    /// );
229    /// ```
230    fn centroid(&self) -> Self::Output {
231        let mut operation = CentroidOperation::new();
232        operation.add_multi_polygon(self);
233        operation.centroid()
234    }
235}
236
237impl<T> Centroid for Rect<T>
238where
239    T: GeoFloat,
240{
241    type Output = Point<T>;
242
243    /// The Centroid of a [`Rect`] is the mean of its [`Point`]s
244    ///
245    /// # Examples
246    ///
247    /// ```
248    /// use geo::Centroid;
249    /// use geo::{Rect, point};
250    ///
251    /// let rect = Rect::new(
252    ///   point!(x: 0.0f32, y: 0.0),
253    ///   point!(x: 1.0, y: 1.0),
254    /// );
255    ///
256    /// assert_eq!(
257    ///     point!(x: 0.5, y: 0.5),
258    ///     rect.centroid(),
259    /// );
260    /// ```
261    fn centroid(&self) -> Self::Output {
262        self.center().into()
263    }
264}
265
266impl<T> Centroid for Triangle<T>
267where
268    T: GeoFloat,
269{
270    type Output = Point<T>;
271
272    /// The Centroid of a [`Triangle`] is the mean of its [`Point`]s
273    ///
274    /// # Examples
275    ///
276    /// ```
277    /// use geo::Centroid;
278    /// use geo::{Triangle, coord, point};
279    ///
280    /// let triangle = Triangle::new(
281    ///   coord!(x: 0.0f32, y: -1.0),
282    ///   coord!(x: 3.0, y: 0.0),
283    ///   coord!(x: 0.0, y: 1.0),
284    /// );
285    ///
286    /// assert_eq!(
287    ///     point!(x: 1.0, y: 0.0),
288    ///     triangle.centroid(),
289    /// );
290    /// ```
291    fn centroid(&self) -> Self::Output {
292        let mut operation = CentroidOperation::new();
293        operation.add_triangle(self);
294        operation
295            .centroid()
296            .expect("triangle cannot have an empty centroid")
297    }
298}
299
300impl<T> Centroid for Point<T>
301where
302    T: GeoFloat,
303{
304    type Output = Point<T>;
305
306    /// The Centroid of a [`Point`] is the point itself
307    ///
308    /// # Examples
309    ///
310    /// ```
311    /// use geo::Centroid;
312    /// use geo::point;
313    ///
314    /// let point = point!(x: 1.0f32, y: 2.0);
315    ///
316    /// assert_eq!(
317    ///     point!(x: 1.0f32, y: 2.0),
318    ///     point.centroid(),
319    /// );
320    /// ```
321    fn centroid(&self) -> Self::Output {
322        *self
323    }
324}
325
326impl<T> Centroid for MultiPoint<T>
327where
328    T: GeoFloat,
329{
330    type Output = Option<Point<T>>;
331
332    /// The Centroid of a [`MultiPoint`] is the mean of all [`Point`]s
333    ///
334    /// # Example
335    ///
336    /// ```
337    /// use geo::Centroid;
338    /// use geo::{MultiPoint, Point};
339    ///
340    /// let empty: Vec<Point> = Vec::new();
341    /// let empty_multi_points: MultiPoint<_> = empty.into();
342    /// assert_eq!(empty_multi_points.centroid(), None);
343    ///
344    /// let points: MultiPoint<_> = vec![(5., 1.), (1., 3.), (3., 2.)].into();
345    /// assert_eq!(points.centroid(), Some(Point::new(3., 2.)));
346    /// ```
347    fn centroid(&self) -> Self::Output {
348        let mut operation = CentroidOperation::new();
349        operation.add_multi_point(self);
350        operation.centroid()
351    }
352}
353
354impl<T> Centroid for Geometry<T>
355where
356    T: GeoFloat,
357{
358    type Output = Option<Point<T>>;
359
360    crate::geometry_delegate_impl! {
361        /// The Centroid of a [`Geometry`] is the centroid of its enum variant
362        ///
363        /// # Examples
364        ///
365        /// ```
366        /// use geo::Centroid;
367        /// use geo::{Geometry, Rect, point};
368        ///
369        /// let rect = Rect::new(
370        ///   point!(x: 0.0f32, y: 0.0),
371        ///   point!(x: 1.0, y: 1.0),
372        /// );
373        /// let geometry = Geometry::from(rect.clone());
374        ///
375        /// assert_eq!(
376        ///     Some(rect.centroid()),
377        ///     geometry.centroid(),
378        /// );
379        ///
380        /// assert_eq!(
381        ///     Some(point!(x: 0.5, y: 0.5)),
382        ///     geometry.centroid(),
383        /// );
384        /// ```
385        fn centroid(&self) -> Self::Output;
386    }
387}
388
389impl<T> Centroid for GeometryCollection<T>
390where
391    T: GeoFloat,
392{
393    type Output = Option<Point<T>>;
394
395    /// The Centroid of a [`GeometryCollection`] is the mean of the centroids of elements, weighted
396    /// by the area of its elements.
397    ///
398    /// Note that this means, that elements which have no area are not considered when calculating
399    /// the centroid.
400    ///
401    /// # Examples
402    ///
403    /// ```
404    /// use geo::Centroid;
405    /// use geo::{Geometry, GeometryCollection, Rect, Triangle, point, coord};
406    ///
407    /// let rect_geometry = Geometry::from(Rect::new(
408    ///   point!(x: 0.0f32, y: 0.0),
409    ///   point!(x: 1.0, y: 1.0),
410    /// ));
411    ///
412    /// let triangle_geometry = Geometry::from(Triangle::new(
413    ///     coord!(x: 0.0f32, y: -1.0),
414    ///     coord!(x: 3.0, y: 0.0),
415    ///     coord!(x: 0.0, y: 1.0),
416    /// ));
417    ///
418    /// let point_geometry = Geometry::from(
419    ///   point!(x: 12351.0, y: 129815.0)
420    /// );
421    ///
422    /// let geometry_collection = GeometryCollection::new_from(
423    ///   vec![
424    ///     rect_geometry,
425    ///     triangle_geometry,
426    ///     point_geometry
427    ///   ]
428    /// );
429    ///
430    /// assert_eq!(
431    ///     Some(point!(x: 0.875, y: 0.125)),
432    ///     geometry_collection.centroid(),
433    /// );
434    /// ```
435    fn centroid(&self) -> Self::Output {
436        let mut operation = CentroidOperation::new();
437        operation.add_geometry_collection(self);
438        operation.centroid()
439    }
440}
441
442struct CentroidOperation<T: GeoFloat>(Option<WeightedCentroid<T>>);
443impl<T: GeoFloat> CentroidOperation<T> {
444    fn new() -> Self {
445        CentroidOperation(None)
446    }
447
448    fn centroid(&self) -> Option<Point<T>> {
449        self.0.as_ref().map(|weighted_centroid| {
450            Point::from(weighted_centroid.accumulated / weighted_centroid.weight)
451        })
452    }
453
454    fn centroid_dimensions(&self) -> Dimensions {
455        self.0
456            .as_ref()
457            .map(|weighted_centroid| weighted_centroid.dimensions)
458            .unwrap_or(Empty)
459    }
460
461    fn add_coord(&mut self, coord: Coord<T>) {
462        self.add_centroid(ZeroDimensional, coord, T::one());
463    }
464
465    fn add_line(&mut self, line: &Line<T>) {
466        match line.dimensions() {
467            ZeroDimensional => self.add_coord(line.start),
468            OneDimensional => {
469                self.add_centroid(OneDimensional, line.centroid().0, line.euclidean_length())
470            }
471            _ => unreachable!("Line must be zero or one dimensional"),
472        }
473    }
474
475    fn add_line_string(&mut self, line_string: &LineString<T>) {
476        if self.centroid_dimensions() > OneDimensional {
477            return;
478        }
479
480        if line_string.0.len() == 1 {
481            self.add_coord(line_string.0[0]);
482            return;
483        }
484
485        for line in line_string.lines() {
486            self.add_line(&line);
487        }
488    }
489
490    fn add_multi_line_string(&mut self, multi_line_string: &MultiLineString<T>) {
491        if self.centroid_dimensions() > OneDimensional {
492            return;
493        }
494
495        for element in &multi_line_string.0 {
496            self.add_line_string(element);
497        }
498    }
499
500    fn add_polygon(&mut self, polygon: &Polygon<T>) {
501        // Polygons which are completely covered by their interior rings have zero area, and
502        // represent a unique degeneracy into a line_string which cannot be handled by accumulating
503        // directly into `self`. Instead, we perform a sub-operation, inspect the result, and only
504        // then incorporate the result into `self.
505
506        let mut exterior_operation = CentroidOperation::new();
507        exterior_operation.add_ring(polygon.exterior());
508
509        let mut interior_operation = CentroidOperation::new();
510        for interior in polygon.interiors() {
511            interior_operation.add_ring(interior);
512        }
513
514        if let Some(exterior_weighted_centroid) = exterior_operation.0 {
515            let mut poly_weighted_centroid = exterior_weighted_centroid;
516            if let Some(interior_weighted_centroid) = interior_operation.0 {
517                poly_weighted_centroid.sub_assign(interior_weighted_centroid);
518                if poly_weighted_centroid.weight.is_zero() {
519                    // A polygon with no area `interiors` completely covers `exterior`, degenerating to a linestring
520                    self.add_line_string(polygon.exterior());
521                    return;
522                }
523            }
524            self.add_weighted_centroid(poly_weighted_centroid);
525        }
526    }
527
528    fn add_multi_point(&mut self, multi_point: &MultiPoint<T>) {
529        if self.centroid_dimensions() > ZeroDimensional {
530            return;
531        }
532
533        for element in &multi_point.0 {
534            self.add_coord(element.0);
535        }
536    }
537
538    fn add_multi_polygon(&mut self, multi_polygon: &MultiPolygon<T>) {
539        for element in &multi_polygon.0 {
540            self.add_polygon(element);
541        }
542    }
543
544    fn add_geometry_collection(&mut self, geometry_collection: &GeometryCollection<T>) {
545        for element in &geometry_collection.0 {
546            self.add_geometry(element);
547        }
548    }
549
550    fn add_rect(&mut self, rect: &Rect<T>) {
551        match rect.dimensions() {
552            ZeroDimensional => self.add_coord(rect.min()),
553            OneDimensional => {
554                // Degenerate rect is a line, treat it the same way we treat flat polygons
555                self.add_line(&Line::new(rect.min(), rect.min()));
556                self.add_line(&Line::new(rect.min(), rect.max()));
557                self.add_line(&Line::new(rect.max(), rect.max()));
558                self.add_line(&Line::new(rect.max(), rect.min()));
559            }
560            TwoDimensional => {
561                self.add_centroid(TwoDimensional, rect.centroid().0, rect.unsigned_area())
562            }
563            Empty => unreachable!("Rect dimensions cannot be empty"),
564        }
565    }
566
567    fn add_triangle(&mut self, triangle: &Triangle<T>) {
568        match triangle.dimensions() {
569            ZeroDimensional => self.add_coord(triangle.0),
570            OneDimensional => {
571                // Degenerate triangle is a line, treat it the same way we treat flat
572                // polygons
573                let l0_1 = Line::new(triangle.0, triangle.1);
574                let l1_2 = Line::new(triangle.1, triangle.2);
575                let l2_0 = Line::new(triangle.2, triangle.0);
576                self.add_line(&l0_1);
577                self.add_line(&l1_2);
578                self.add_line(&l2_0);
579            }
580            TwoDimensional => {
581                let centroid = (triangle.0 + triangle.1 + triangle.2) / T::from(3).unwrap();
582                self.add_centroid(TwoDimensional, centroid, triangle.unsigned_area());
583            }
584            Empty => unreachable!("Rect dimensions cannot be empty"),
585        }
586    }
587
588    fn add_geometry(&mut self, geometry: &Geometry<T>) {
589        match geometry {
590            Geometry::Point(g) => self.add_coord(g.0),
591            Geometry::Line(g) => self.add_line(g),
592            Geometry::LineString(g) => self.add_line_string(g),
593            Geometry::Polygon(g) => self.add_polygon(g),
594            Geometry::MultiPoint(g) => self.add_multi_point(g),
595            Geometry::MultiLineString(g) => self.add_multi_line_string(g),
596            Geometry::MultiPolygon(g) => self.add_multi_polygon(g),
597            Geometry::GeometryCollection(g) => self.add_geometry_collection(g),
598            Geometry::Rect(g) => self.add_rect(g),
599            Geometry::Triangle(g) => self.add_triangle(g),
600        }
601    }
602
603    fn add_ring(&mut self, ring: &LineString<T>) {
604        debug_assert!(ring.is_closed());
605
606        let area = get_linestring_area(ring);
607        if area == T::zero() {
608            match ring.dimensions() {
609                // empty ring doesn't contribute to centroid
610                Empty => {}
611                // degenerate ring is a point
612                ZeroDimensional => self.add_coord(ring[0]),
613                // zero-area ring is a line string
614                _ => self.add_line_string(ring),
615            }
616            return;
617        }
618
619        // Since area is non-zero, we know the ring has at least one point
620        let shift = ring.0[0];
621
622        let accumulated_coord = ring.lines().fold(Coord::zero(), |accum, line| {
623            use crate::MapCoords;
624            let line = line.map_coords(|c| c - shift);
625            let tmp = line.determinant();
626            accum + (line.end + line.start) * tmp
627        });
628        let six = T::from(6).unwrap();
629        let centroid = accumulated_coord / (six * area) + shift;
630        let weight = area.abs();
631        self.add_centroid(TwoDimensional, centroid, weight);
632    }
633
634    fn add_centroid(&mut self, dimensions: Dimensions, centroid: Coord<T>, weight: T) {
635        let weighted_centroid = WeightedCentroid {
636            dimensions,
637            weight,
638            accumulated: centroid * weight,
639        };
640        self.add_weighted_centroid(weighted_centroid);
641    }
642
643    fn add_weighted_centroid(&mut self, other: WeightedCentroid<T>) {
644        match self.0.as_mut() {
645            Some(centroid) => centroid.add_assign(other),
646            None => self.0 = Some(other),
647        }
648    }
649}
650
651// Aggregated state for accumulating the centroid of a geometry or collection of geometries.
652struct WeightedCentroid<T: GeoFloat> {
653    weight: T,
654    accumulated: Coord<T>,
655    /// Collections of Geometries can have different dimensionality. Centroids must be considered
656    /// separately by dimensionality.
657    ///
658    /// e.g. If I have several Points, adding a new `Point` will affect their centroid.
659    ///
660    /// However, because a Point is zero dimensional, it is infinitely small when compared to
661    /// any 2-D Polygon. Thus a Point will not affect the centroid of any GeometryCollection
662    /// containing a 2-D Polygon.
663    ///
664    /// So, when accumulating a centroid, we must track the dimensionality of the centroid
665    dimensions: Dimensions,
666}
667
668impl<T: GeoFloat> WeightedCentroid<T> {
669    fn add_assign(&mut self, b: WeightedCentroid<T>) {
670        match self.dimensions.cmp(&b.dimensions) {
671            Ordering::Less => *self = b,
672            Ordering::Greater => {}
673            Ordering::Equal => {
674                self.accumulated = self.accumulated + b.accumulated;
675                self.weight = self.weight + b.weight;
676            }
677        }
678    }
679
680    fn sub_assign(&mut self, b: WeightedCentroid<T>) {
681        match self.dimensions.cmp(&b.dimensions) {
682            Ordering::Less => *self = b,
683            Ordering::Greater => {}
684            Ordering::Equal => {
685                self.accumulated = self.accumulated - b.accumulated;
686                self.weight = self.weight - b.weight;
687            }
688        }
689    }
690}
691
692#[cfg(test)]
693mod test {
694    use super::*;
695    use crate::{coord, line_string, point, polygon};
696
697    /// small helper to create a coordinate
698    fn c<T: GeoFloat>(x: T, y: T) -> Coord<T> {
699        coord! { x: x, y: y }
700    }
701
702    /// small helper to create a point
703    fn p<T: GeoFloat>(x: T, y: T) -> Point<T> {
704        point! { x: x, y: y }
705    }
706
707    // Tests: Centroid of LineString
708    #[test]
709    fn empty_linestring_test() {
710        let linestring: LineString<f32> = line_string![];
711        let centroid = linestring.centroid();
712        assert!(centroid.is_none());
713    }
714    #[test]
715    fn linestring_one_point_test() {
716        let coord = coord! {
717            x: 40.02f64,
718            y: 116.34,
719        };
720        let linestring = line_string![coord];
721        let centroid = linestring.centroid();
722        assert_eq!(centroid, Some(Point::from(coord)));
723    }
724    #[test]
725    fn linestring_test() {
726        let linestring = line_string![
727            (x: 1., y: 1.),
728            (x: 7., y: 1.),
729            (x: 8., y: 1.),
730            (x: 9., y: 1.),
731            (x: 10., y: 1.),
732            (x: 11., y: 1.)
733        ];
734        assert_eq!(linestring.centroid(), Some(point!(x: 6., y: 1. )));
735    }
736    #[test]
737    fn linestring_with_repeated_point_test() {
738        let l1 = LineString::from(vec![p(1., 1.), p(1., 1.), p(1., 1.)]);
739        assert_eq!(l1.centroid(), Some(p(1., 1.)));
740
741        let l2 = LineString::from(vec![p(2., 2.), p(2., 2.), p(2., 2.)]);
742        let mls = MultiLineString::new(vec![l1, l2]);
743        assert_eq!(mls.centroid(), Some(p(1.5, 1.5)));
744    }
745    // Tests: Centroid of MultiLineString
746    #[test]
747    fn empty_multilinestring_test() {
748        let mls: MultiLineString = MultiLineString::new(vec![]);
749        let centroid = mls.centroid();
750        assert!(centroid.is_none());
751    }
752    #[test]
753    fn multilinestring_with_empty_line_test() {
754        let mls: MultiLineString = MultiLineString::new(vec![line_string![]]);
755        let centroid = mls.centroid();
756        assert!(centroid.is_none());
757    }
758    #[test]
759    fn multilinestring_length_0_test() {
760        let coord = coord! {
761            x: 40.02f64,
762            y: 116.34,
763        };
764        let mls: MultiLineString = MultiLineString::new(vec![
765            line_string![coord],
766            line_string![coord],
767            line_string![coord],
768        ]);
769        assert_relative_eq!(mls.centroid().unwrap(), Point::from(coord));
770    }
771    #[test]
772    fn multilinestring_one_line_test() {
773        let linestring = line_string![
774            (x: 1., y: 1.),
775            (x: 7., y: 1.),
776            (x: 8., y: 1.),
777            (x: 9., y: 1.),
778            (x: 10., y: 1.),
779            (x: 11., y: 1.)
780        ];
781        let mls: MultiLineString = MultiLineString::new(vec![linestring]);
782        assert_relative_eq!(mls.centroid().unwrap(), point! { x: 6., y: 1. });
783    }
784    #[test]
785    fn multilinestring_test() {
786        let v1 = line_string![(x: 0.0, y: 0.0), (x: 1.0, y: 10.0)];
787        let v2 = line_string![(x: 1.0, y: 10.0), (x: 2.0, y: 0.0), (x: 3.0, y: 1.0)];
788        let v3 = line_string![(x: -12.0, y: -100.0), (x: 7.0, y: 8.0)];
789        let mls = MultiLineString::new(vec![v1, v2, v3]);
790        assert_relative_eq!(
791            mls.centroid().unwrap(),
792            point![x: -1.9097834383655845, y: -37.683866439745714]
793        );
794    }
795    // Tests: Centroid of Polygon
796    #[test]
797    fn empty_polygon_test() {
798        let poly: Polygon<f32> = polygon![];
799        assert!(poly.centroid().is_none());
800    }
801    #[test]
802    fn polygon_one_point_test() {
803        let p = point![ x: 2., y: 1. ];
804        let v = Vec::new();
805        let linestring = line_string![p.0];
806        let poly = Polygon::new(linestring, v);
807        assert_relative_eq!(poly.centroid().unwrap(), p);
808    }
809
810    #[test]
811    fn centroid_polygon_numerical_stability() {
812        let polygon = {
813            use std::f64::consts::PI;
814            const NUM_VERTICES: usize = 10;
815            const ANGLE_INC: f64 = 2. * PI / NUM_VERTICES as f64;
816
817            Polygon::new(
818                (0..NUM_VERTICES)
819                    .map(|i| {
820                        let angle = i as f64 * ANGLE_INC;
821                        coord! {
822                            x: angle.cos(),
823                            y: angle.sin(),
824                        }
825                    })
826                    .collect::<Vec<_>>()
827                    .into(),
828                vec![],
829            )
830        };
831
832        let centroid = polygon.centroid().unwrap();
833
834        let shift = coord! { x: 1.5e8, y: 1.5e8 };
835
836        use crate::map_coords::MapCoords;
837        let polygon = polygon.map_coords(|c| c + shift);
838
839        let new_centroid = polygon.centroid().unwrap().map_coords(|c| c - shift);
840        debug!("centroid {:?}", centroid.0);
841        debug!("new_centroid {:?}", new_centroid.0);
842        assert_relative_eq!(centroid.0.x, new_centroid.0.x, max_relative = 0.0001);
843        assert_relative_eq!(centroid.0.y, new_centroid.0.y, max_relative = 0.0001);
844    }
845
846    #[test]
847    fn polygon_test() {
848        let poly = polygon![
849            (x: 0., y: 0.),
850            (x: 2., y: 0.),
851            (x: 2., y: 2.),
852            (x: 0., y: 2.),
853            (x: 0., y: 0.)
854        ];
855        assert_relative_eq!(poly.centroid().unwrap(), point![x:1., y:1.]);
856    }
857    #[test]
858    fn polygon_hole_test() {
859        // hexagon
860        let ls1 = LineString::from(vec![
861            (5.0, 1.0),
862            (4.0, 2.0),
863            (4.0, 3.0),
864            (5.0, 4.0),
865            (6.0, 4.0),
866            (7.0, 3.0),
867            (7.0, 2.0),
868            (6.0, 1.0),
869            (5.0, 1.0),
870        ]);
871
872        let ls2 = LineString::from(vec![(5.0, 1.3), (5.5, 2.0), (6.0, 1.3), (5.0, 1.3)]);
873
874        let ls3 = LineString::from(vec![(5., 2.3), (5.5, 3.0), (6., 2.3), (5., 2.3)]);
875
876        let p1 = Polygon::new(ls1, vec![ls2, ls3]);
877        let centroid = p1.centroid().unwrap();
878        assert_relative_eq!(centroid, point!(x: 5.5, y: 2.5518518518518523));
879    }
880    #[test]
881    fn flat_polygon_test() {
882        let poly = Polygon::new(
883            LineString::from(vec![p(0., 1.), p(1., 1.), p(0., 1.)]),
884            vec![],
885        );
886        assert_eq!(poly.centroid(), Some(p(0.5, 1.)));
887    }
888    #[test]
889    fn multi_poly_with_flat_polygon_test() {
890        let poly = Polygon::new(
891            LineString::from(vec![p(0., 0.), p(1., 0.), p(0., 0.)]),
892            vec![],
893        );
894        let multipoly = MultiPolygon::new(vec![poly]);
895        assert_eq!(multipoly.centroid(), Some(p(0.5, 0.)));
896    }
897    #[test]
898    fn multi_poly_with_multiple_flat_polygon_test() {
899        let p1 = Polygon::new(
900            LineString::from(vec![p(1., 1.), p(1., 3.), p(1., 1.)]),
901            vec![],
902        );
903        let p2 = Polygon::new(
904            LineString::from(vec![p(2., 2.), p(6., 2.), p(2., 2.)]),
905            vec![],
906        );
907        let multipoly = MultiPolygon::new(vec![p1, p2]);
908        assert_eq!(multipoly.centroid(), Some(p(3., 2.)));
909    }
910    #[test]
911    fn multi_poly_with_only_points_test() {
912        let p1 = Polygon::new(
913            LineString::from(vec![p(1., 1.), p(1., 1.), p(1., 1.)]),
914            vec![],
915        );
916        assert_eq!(p1.centroid(), Some(p(1., 1.)));
917        let p2 = Polygon::new(
918            LineString::from(vec![p(2., 2.), p(2., 2.), p(2., 2.)]),
919            vec![],
920        );
921        let multipoly = MultiPolygon::new(vec![p1, p2]);
922        assert_eq!(multipoly.centroid(), Some(p(1.5, 1.5)));
923    }
924    #[test]
925    fn multi_poly_with_one_ring_and_one_real_poly() {
926        // if the multipolygon is composed of a 'normal' polygon (with an area not null)
927        // and a ring (a polygon with a null area)
928        // the centroid of the multipolygon is the centroid of the 'normal' polygon
929        let normal = Polygon::new(
930            LineString::from(vec![p(1., 1.), p(1., 3.), p(3., 1.), p(1., 1.)]),
931            vec![],
932        );
933        let flat = Polygon::new(
934            LineString::from(vec![p(2., 2.), p(6., 2.), p(2., 2.)]),
935            vec![],
936        );
937        let multipoly = MultiPolygon::new(vec![normal.clone(), flat]);
938        assert_eq!(multipoly.centroid(), normal.centroid());
939    }
940    #[test]
941    fn polygon_flat_interior_test() {
942        let poly = Polygon::new(
943            LineString::from(vec![p(0., 0.), p(0., 1.), p(1., 1.), p(1., 0.), p(0., 0.)]),
944            vec![LineString::from(vec![p(0., 0.), p(0., 1.), p(0., 0.)])],
945        );
946        assert_eq!(poly.centroid(), Some(p(0.5, 0.5)));
947    }
948    #[test]
949    fn empty_interior_polygon_test() {
950        let poly = Polygon::new(
951            LineString::from(vec![p(0., 0.), p(0., 1.), p(1., 1.), p(1., 0.), p(0., 0.)]),
952            vec![LineString::new(vec![])],
953        );
954        assert_eq!(poly.centroid(), Some(p(0.5, 0.5)));
955    }
956    #[test]
957    fn polygon_ring_test() {
958        let square = LineString::from(vec![p(0., 0.), p(0., 1.), p(1., 1.), p(1., 0.), p(0., 0.)]);
959        let poly = Polygon::new(square.clone(), vec![square]);
960        assert_eq!(poly.centroid(), Some(p(0.5, 0.5)));
961    }
962    #[test]
963    fn polygon_cell_test() {
964        // test the centroid of polygon with a null area
965        // this one a polygon with 2 interior polygon that makes a partition of the exterior
966        let square = LineString::from(vec![p(0., 0.), p(0., 2.), p(2., 2.), p(2., 0.), p(0., 0.)]);
967        let bottom = LineString::from(vec![p(0., 0.), p(2., 0.), p(2., 1.), p(0., 1.), p(0., 0.)]);
968        let top = LineString::from(vec![p(0., 1.), p(2., 1.), p(2., 2.), p(0., 2.), p(0., 1.)]);
969        let poly = Polygon::new(square, vec![top, bottom]);
970        assert_eq!(poly.centroid(), Some(p(1., 1.)));
971    }
972    // Tests: Centroid of MultiPolygon
973    #[test]
974    fn empty_multipolygon_polygon_test() {
975        assert!(MultiPolygon::<f64>::new(Vec::new()).centroid().is_none());
976    }
977
978    #[test]
979    fn multipolygon_one_polygon_test() {
980        let linestring =
981            LineString::from(vec![p(0., 0.), p(2., 0.), p(2., 2.), p(0., 2.), p(0., 0.)]);
982        let poly = Polygon::new(linestring, Vec::new());
983        assert_eq!(MultiPolygon::new(vec![poly]).centroid(), Some(p(1., 1.)));
984    }
985    #[test]
986    fn multipolygon_two_polygons_test() {
987        let linestring =
988            LineString::from(vec![p(2., 1.), p(5., 1.), p(5., 3.), p(2., 3.), p(2., 1.)]);
989        let poly1 = Polygon::new(linestring, Vec::new());
990        let linestring =
991            LineString::from(vec![p(7., 1.), p(8., 1.), p(8., 2.), p(7., 2.), p(7., 1.)]);
992        let poly2 = Polygon::new(linestring, Vec::new());
993        let centroid = MultiPolygon::new(vec![poly1, poly2]).centroid().unwrap();
994        assert_relative_eq!(
995            centroid,
996            point![x: 4.071428571428571, y: 1.9285714285714286]
997        );
998    }
999    #[test]
1000    fn multipolygon_two_polygons_of_opposite_clockwise_test() {
1001        let linestring = LineString::from(vec![(0., 0.), (2., 0.), (2., 2.), (0., 2.), (0., 0.)]);
1002        let poly1 = Polygon::new(linestring, Vec::new());
1003        let linestring = LineString::from(vec![(0., 0.), (-2., 0.), (-2., 2.), (0., 2.), (0., 0.)]);
1004        let poly2 = Polygon::new(linestring, Vec::new());
1005        assert_relative_eq!(
1006            MultiPolygon::new(vec![poly1, poly2]).centroid().unwrap(),
1007            point![x: 0., y: 1.]
1008        );
1009    }
1010    #[test]
1011    fn bounding_rect_test() {
1012        let bounding_rect = Rect::new(coord! { x: 0., y: 50. }, coord! { x: 4., y: 100. });
1013        let point = point![x: 2., y: 75.];
1014        assert_eq!(point, bounding_rect.centroid());
1015    }
1016    #[test]
1017    fn line_test() {
1018        let line1 = Line::new(c(0., 1.), c(1., 3.));
1019        assert_eq!(line1.centroid(), point![x: 0.5, y: 2.]);
1020    }
1021    #[test]
1022    fn collection_weighting() {
1023        let p0 = point!(x: 0.0, y: 0.0);
1024        let p1 = point!(x: 2.0, y: 0.0);
1025        let p2 = point!(x: 2.0, y: 2.0);
1026        let p3 = point!(x: 0.0, y: 2.0);
1027
1028        let multi_point = MultiPoint::new(vec![p0, p1, p2, p3]);
1029        assert_eq!(multi_point.centroid().unwrap(), point!(x: 1.0, y: 1.0));
1030
1031        let collection =
1032            GeometryCollection::new_from(vec![MultiPoint::new(vec![p1, p2, p3]).into(), p0.into()]);
1033
1034        assert_eq!(collection.centroid().unwrap(), point!(x: 1.0, y: 1.0));
1035    }
1036    #[test]
1037    fn triangles() {
1038        // boring triangle
1039        assert_eq!(
1040            Triangle::new(c(0., 0.), c(3., 0.), c(1.5, 3.)).centroid(),
1041            point!(x: 1.5, y: 1.0)
1042        );
1043
1044        // flat triangle
1045        assert_eq!(
1046            Triangle::new(c(0., 0.), c(3., 0.), c(1., 0.)).centroid(),
1047            point!(x: 1.5, y: 0.0)
1048        );
1049
1050        // flat triangle that's not axis-aligned
1051        assert_eq!(
1052            Triangle::new(c(0., 0.), c(3., 3.), c(1., 1.)).centroid(),
1053            point!(x: 1.5, y: 1.5)
1054        );
1055
1056        // triangle with some repeated points
1057        assert_eq!(
1058            Triangle::new(c(0., 0.), c(0., 0.), c(1., 0.)).centroid(),
1059            point!(x: 0.5, y: 0.0)
1060        );
1061
1062        // triangle with all repeated points
1063        assert_eq!(
1064            Triangle::new(c(0., 0.5), c(0., 0.5), c(0., 0.5)).centroid(),
1065            point!(x: 0., y: 0.5)
1066        )
1067    }
1068
1069    #[test]
1070    fn degenerate_triangle_like_ring() {
1071        let triangle = Triangle::new(c(0., 0.), c(1., 1.), c(2., 2.));
1072        let poly: Polygon<_> = triangle.into();
1073
1074        let line = Line::new(c(0., 1.), c(1., 3.));
1075
1076        let g1 = GeometryCollection::new_from(vec![triangle.into(), line.into()]);
1077        let g2 = GeometryCollection::new_from(vec![poly.into(), line.into()]);
1078        assert_eq!(g1.centroid(), g2.centroid());
1079    }
1080
1081    #[test]
1082    fn degenerate_rect_like_ring() {
1083        let rect = Rect::new(c(0., 0.), c(0., 4.));
1084        let poly: Polygon<_> = rect.into();
1085
1086        let line = Line::new(c(0., 1.), c(1., 3.));
1087
1088        let g1 = GeometryCollection::new_from(vec![rect.into(), line.into()]);
1089        let g2 = GeometryCollection::new_from(vec![poly.into(), line.into()]);
1090        assert_eq!(g1.centroid(), g2.centroid());
1091    }
1092
1093    #[test]
1094    fn rectangles() {
1095        // boring rect
1096        assert_eq!(
1097            Rect::new(c(0., 0.), c(4., 4.)).centroid(),
1098            point!(x: 2.0, y: 2.0)
1099        );
1100
1101        // flat rect
1102        assert_eq!(
1103            Rect::new(c(0., 0.), c(4., 0.)).centroid(),
1104            point!(x: 2.0, y: 0.0)
1105        );
1106
1107        // rect with all repeated points
1108        assert_eq!(
1109            Rect::new(c(4., 4.), c(4., 4.)).centroid(),
1110            point!(x: 4., y: 4.)
1111        );
1112
1113        // collection with rect
1114        let mut collection = GeometryCollection::new_from(vec![
1115            p(0., 0.).into(),
1116            p(6., 0.).into(),
1117            p(6., 6.).into(),
1118        ]);
1119        // sanity check
1120        assert_eq!(collection.centroid().unwrap(), point!(x: 4., y: 2.));
1121
1122        // 0-d rect treated like point
1123        collection.0.push(Rect::new(c(0., 6.), c(0., 6.)).into());
1124        assert_eq!(collection.centroid().unwrap(), point!(x: 3., y: 3.));
1125
1126        // 1-d rect treated like line. Since a line has higher dimensions than the rest of the
1127        // collection, its centroid clobbers everything else in the collection.
1128        collection.0.push(Rect::new(c(0., 0.), c(0., 2.)).into());
1129        assert_eq!(collection.centroid().unwrap(), point!(x: 0., y: 1.));
1130
1131        // 2-d has higher dimensions than the rest of the collection, so its centroid clobbers
1132        // everything else in the collection.
1133        collection
1134            .0
1135            .push(Rect::new(c(10., 10.), c(11., 11.)).into());
1136        assert_eq!(collection.centroid().unwrap(), point!(x: 10.5, y: 10.5));
1137    }
1138}
geo/algorithm/centroid.rs

geo/algorithm/
centroid.rs
