Add support for centroid() on geometry collections

The implementation isn't really efficient, because it needs to copy
(part of) the geometry collection, but it can be improved later on
if needed.

Author: joto

src/geom-functions.cpp

@@ -700,16 +700,55 @@ geometry_t line_merge(geometry_t const &input)
     return output;
 }
 
+/**
+ * This helper function is used to calculate centroids of geometry collections.
+ * It first creates a multi geometry that only contains the geometries of
+ * dimension N from the input collection. This is done by copying the geometry,
+ * which isn't very efficient, but hopefully the centroid of a geometry
+ * collection isn't used very often. This can be optimized if needed.
+ *
+ * Then the centroid of this new collection is calculated.
+ *
+ * Nested geometry collections are not allowed.
+ */
+template <std::size_t N, typename T>
+static void filtered_centroid(collection_t const &collection, point_t *center)
+{
+    multigeometry_t<T> multi;
+    for (auto const &geom : collection) {
+        assert(!geom.is_collection());
+        if (!geom.is_null() && dimension(geom) == N) {
+            if (geom.is_multi()) {
+                for (auto const &sgeom : geom.get<multigeometry_t<T>>()) {
+                    multi.add_geometry() = sgeom;
+                }
+            } else {
+                multi.add_geometry() = geom.get<T>();
+            }
+        }
+    }
+    boost::geometry::centroid(multi, *center);
+}
+
 geometry_t centroid(geometry_t const &geom)
 {
     geom::geometry_t output{point_t{}, geom.srid()};
     auto &center = output.get<point_t>();
 
     geom.visit(overloaded{
         [&](geom::nullgeom_t const & /*input*/) { output.reset(); },
-        [&](geom::collection_t const & /*input*/) {
-            throw std::runtime_error{
-                "Centroid of geometry collection not implemented yet"};
+        [&](geom::collection_t const &input) {
+            switch (dimension(input)) {
+            case 0:
+                filtered_centroid<0, point_t>(input, &center);
+                break;
+            case 1:
+                filtered_centroid<1, linestring_t>(input, &center);
+                break;
+            default: // 2
+                filtered_centroid<2, polygon_t>(input, &center);
+                break;
+            }
         },
         [&](auto const &input) { boost::geometry::centroid(input, center); }});
 

tests/test-geom-collections.cpp

@@ -27,10 +27,30 @@ TEST_CASE("geometry collection with point", "[NoDB]")
     REQUIRE(num_geometries(geom) == 1);
     REQUIRE(area(geom) == Approx(0.0));
     REQUIRE(length(geom) == Approx(0.0));
-    REQUIRE_THROWS(centroid(geom));
+    REQUIRE(centroid(geom) == geom::geometry_t{geom::point_t{1, 1}});
     REQUIRE(geometry_n(geom, 1) == geom::geometry_t{geom::point_t{1, 1}});
 }
 
+TEST_CASE("geometry collection with multipoint", "[NoDB]")
+{
+    geom::geometry_t mpgeom{geom::multipoint_t{}};
+    auto &mp = mpgeom.get<geom::multipoint_t>();
+    mp.add_geometry(geom::point_t{1, 1});
+    mp.add_geometry(geom::point_t{1, 2});
+    mp.add_geometry(geom::point_t{2, 1});
+    mp.add_geometry(geom::point_t{2, 2});
+
+    geom::geometry_t geom{geom::collection_t{}};
+    auto &c = geom.get<geom::collection_t>();
+    c.add_geometry(std::move(mpgeom));
+
+    REQUIRE(geometry_type(geom) == "GEOMETRYCOLLECTION");
+    REQUIRE(dimension(geom) == 0);
+    REQUIRE(num_geometries(geom) == 1);
+    REQUIRE(area(geom) == Approx(0.0));
+    REQUIRE(centroid(geom) == geom::geometry_t{geom::point_t{1.5, 1.5}});
+}
+
 TEST_CASE("geometry collection with several geometries", "[NoDB]")
 {
     geom::geometry_t geom{geom::collection_t{}};
@@ -45,7 +65,7 @@ TEST_CASE("geometry collection with several geometries", "[NoDB]")
     REQUIRE(num_geometries(geom) == 3);
     REQUIRE(area(geom) == Approx(0.0));
     REQUIRE(length(geom) == Approx(1.41421));
-    REQUIRE_THROWS(centroid(geom));
+    REQUIRE(centroid(geom) == geom::geometry_t{geom::point_t{1.5, 1.5}});
     REQUIRE(geometry_n(geom, 1) == geom::geometry_t{geom::point_t{1, 1}});
     REQUIRE(geometry_n(geom, 2) ==
             geom::geometry_t{geom::linestring_t{{1, 1}, {2, 2}}});
@@ -57,7 +77,7 @@ TEST_CASE("create_collection from OSM data", "[NoDB]")
     test_buffer_t buffer;
     buffer.add_node("n1 x1 y1");
     buffer.add_way("w20 Nn1x1y1,n2x2y1,n3x2y2,n4x1y2,n1x1y1");
-    buffer.add_way("w21 Nn5x10y10,n6x10y20");
+    buffer.add_way("w21 Nn5x10y10,n6x10y11");
     buffer.add_relation("r30 Mw20@");
 
     auto const geom = geom::create_collection(buffer.buffer());
@@ -70,11 +90,11 @@ TEST_CASE("create_collection from OSM data", "[NoDB]")
     REQUIRE(c[0] == geom::geometry_t{geom::point_t{1, 1}});
     REQUIRE(c[1] == geom::geometry_t{geom::linestring_t{
                         {1, 1}, {2, 1}, {2, 2}, {1, 2}, {1, 1}}});
-    REQUIRE(c[2] == geom::geometry_t{geom::linestring_t{{10, 10}, {10, 20}}});
+    REQUIRE(c[2] == geom::geometry_t{geom::linestring_t{{10, 10}, {10, 11}}});
 
     REQUIRE(area(geom) == Approx(0.0));
-    REQUIRE(length(geom) == Approx(14.0));
-    REQUIRE_THROWS(centroid(geom));
+    REQUIRE(length(geom) == Approx(5.0));
+    REQUIRE(centroid(geom) == geom::geometry_t{geom::point_t{3.2, 3.3}});
 }
 
 TEST_CASE("create_collection from no OSM data returns null geometry", "[NoDB]")