Fix ConcaveHullOfPolygons nested shell handling

include/geos/algorithm/hull/ConcaveHullOfPolygons.h

@@ -124,10 +124,6 @@ class GEOS_DLL ConcaveHullOfPolygons {
 
     std::unique_ptr<Geometry> createEmptyHull();
 
-    static void extractShellRings(
-        const Geometry* polygons,
-        std::vector<const LinearRing*>& rings);
-
     void buildHullTris();
 
     /**

include/geos/algorithm/hull/OuterShellsExtracter.h

@@ -0,0 +1,75 @@
+/**********************************************************************
+ *
+ * GEOS - Geometry Engine Open Source
+ * http://geos.osgeo.org
+ *
+ * Copyright (C) 2024 Martin Davis
+ *
+ * This is free software; you can redistribute and/or modify it under
+ * the terms of the GNU Lesser General Public Licence as published
+ * by the Free Software Foundation.
+ * See the COPYING file for more information.
+ *
+ **********************************************************************/
+
+#pragma once
+
+#include <vector>
+
+namespace geos {
+namespace geom {
+class Coordinate;
+class CoordinateSequence;
+class Envelope;
+class Geometry;
+class GeometryCollection;
+class GeometryFactory;
+class LinearRing;
+class Polygon;
+}
+}
+
+using geos::geom::Geometry;
+using geos::geom::LinearRing;
+
+namespace geos {
+namespace algorithm { // geos::algorithm
+namespace hull {      // geos::algorithm::hull
+
+/**
+ * Extracts the rings of outer shells from a polygonal geometry.
+ * Outer shells are the shells of polygon elements which
+ * are not nested inside holes of other polygons.
+ *
+ * \author Martin Davis
+ */
+class OuterShellsExtracter {
+private:
+
+    OuterShellsExtracter(const Geometry& g);
+
+    void extractOuterShells(std::vector<const LinearRing*>& outerShells);
+
+    bool isOuter(const LinearRing& shell, std::vector<const LinearRing*>& outerShells);
+
+    bool covers(const LinearRing& shellA, const LinearRing& shellB);
+
+    bool isPointInRing(const LinearRing& shell, const LinearRing& shellRing);
+
+    static void extractShellRings(const Geometry& polygons, std::vector<const LinearRing*>& shells);
+
+    static bool envelopeAreaComparator(
+        const LinearRing* g1,
+        const LinearRing* g2);
+
+    const Geometry& geom;
+
+public:
+    static void extractShells(const Geometry* polygons, std::vector<const LinearRing*>& shells);
+
+};
+
+} // geos::algorithm::hull
+} // geos::algorithm
+} // geos
+

src/algorithm/hull/ConcaveHullOfPolygons.cpp

@@ -13,6 +13,7 @@
  **********************************************************************/
 
 #include <geos/algorithm/hull/ConcaveHullOfPolygons.h>
+#include <geos/algorithm/hull/OuterShellsExtracter.h>
 #include <geos/geom/Coordinate.h>
 #include <geos/geom/Envelope.h>
 #include <geos/geom/Geometry.h>
@@ -182,26 +183,11 @@ ConcaveHullOfPolygons::createEmptyHull()
     return geomFactory->createPolygon();
 }
 
-/* private static */
-void
-ConcaveHullOfPolygons::extractShellRings(const Geometry* polygons, std::vector<const LinearRing*>& rings)
-{
-    rings.clear();
-    for (std::size_t i = 0; i < polygons->getNumGeometries(); i++) {
-        const Geometry* consGeom = polygons->getGeometryN(i);
-        const Polygon* consPoly = static_cast<const Polygon*>(consGeom);
-        const LinearRing* lr = consPoly->getExteriorRing();
-        rings.push_back(lr);
-    }
-    return;
-}
-
-
 /* private */
 void
 ConcaveHullOfPolygons::buildHullTris()
 {
-    extractShellRings(inputPolygons, polygonRings);
+    OuterShellsExtracter::extractShells(inputPolygons, polygonRings);
     std::unique_ptr<Polygon> frame = createFrame(inputPolygons->getEnvelopeInternal());
     ConstrainedDelaunayTriangulator::triangulatePolygon(frame.get(), triList);
     //System.out.println(tris);

src/algorithm/hull/OuterShellsExtracter.cpp

@@ -0,0 +1,123 @@
+
+/**********************************************************************
+ *
+ * GEOS - Geometry Engine Open Source
+ * http://geos.osgeo.org
+ *
+ * Copyright (C) 2022 Paul Ramsey <[email protected]>
+ *
+ * This is free software; you can redistribute and/or modify it under
+ * the terms of the GNU Lesser General Public Licence as published
+ * by the Free Software Foundation.
+ * See the COPYING file for more information.
+ *
+ **********************************************************************/
+
+#include <geos/algorithm/hull/OuterShellsExtracter.h>
+#include <geos/algorithm/PointLocation.h>
+#include <geos/geom/Polygon.h>
+#include <geos/geom/Coordinate.h>
+
+using geos::algorithm::PointLocation;
+using geos::geom::Polygon;
+using geos::geom::CoordinateXY;
+
+namespace geos {
+namespace algorithm {
+namespace hull {
+
+/* public static */
+void 
+OuterShellsExtracter::extractShells(const Geometry* polygons, std::vector<const LinearRing*>& shells)
+{
+    OuterShellsExtracter extracter(*polygons);
+    extracter.extractOuterShells(shells);
+}
+
+/* private */
+OuterShellsExtracter::OuterShellsExtracter(const geom::Geometry& g)
+    : geom(g)
+{
+}
+
+/* private */
+void 
+OuterShellsExtracter::extractOuterShells(std::vector<const LinearRing*>& outerShells) 
+{
+    std::vector<const LinearRing*> shells;
+    extractShellRings(geom, shells);
+    //-- sort shells in order of increasing envelope area
+    std::sort(shells.begin(), shells.end(), envelopeAreaComparator);
+
+    //-- Scan shells by decreasing area to ensure that shells are added before any nested shells
+    for (auto i = shells.rbegin(); i != shells.rend(); ++i) { 
+      const LinearRing* shell = *i;
+      if (outerShells.size() == 0 
+          || isOuter(*shell, outerShells)) {
+        outerShells.push_back(shell);
+      }
+    }
+}
+
+bool
+OuterShellsExtracter::envelopeAreaComparator(
+        const LinearRing* g1,
+        const LinearRing* g2)
+{
+    double area1 = g1->getEnvelopeInternal()->getArea();
+    double area2 = g2->getEnvelopeInternal()->getArea();
+    if (area1 < area2)
+        return true;
+    else
+        return false;
+}
+
+/* private */
+bool 
+OuterShellsExtracter::isOuter(const LinearRing& shell, std::vector<const LinearRing*>& outerShells)
+{
+    for (const LinearRing* outShell : outerShells) {
+      if (covers(*outShell, shell)) {
+        return false;
+      }
+    }
+    return true;
+}
+
+/* private */
+bool 
+OuterShellsExtracter::covers(const LinearRing& shellA, const LinearRing& shellB)
+{
+    //-- if shellB envelope is not covered then shell is not covered
+    if (! shellA.getEnvelopeInternal()->covers(shellB.getEnvelopeInternal()))
+      return false;
+    //-- if a shellB point lies inside shellA, shell is covered (since shells do not overlap)
+    if (isPointInRing(shellB, shellA))
+      return true;
+    return false;
+}
+
+bool 
+OuterShellsExtracter::isPointInRing(const LinearRing& shell, const LinearRing& shellRing)
+{
+    //TODO: optimize this with cached index
+    const CoordinateXY* pt = shell.getCoordinate();
+    return PointLocation::isInRing(*pt, shellRing.getCoordinatesRO());
+}
+
+void 
+OuterShellsExtracter::extractShellRings(const Geometry& polygons, std::vector<const LinearRing*>& shells)
+{
+    shells.clear();
+    for (std::size_t i = 0; i < polygons.getNumGeometries(); i++) {
+        const Geometry* consGeom = polygons.getGeometryN(i);
+        const Polygon* consPoly = static_cast<const Polygon*>(consGeom);
+        const LinearRing* lr = consPoly->getExteriorRing();
+        shells.push_back(lr);
+    }
+}
+
+
+} // geos::algorithm::hull
+} // geos::algorithm
+} // geos

tests/unit/algorithm/hull/ConcaveHullOfPolygonsTest.cpp

@@ -191,5 +191,27 @@ void object::test<7>()
         "POLYGON ((6 9, 8 9, 9 5, 8 0, 6 0, 5 0, 1 0, 1 4, 1 5, 1 9, 5 9, 6 9))");
 }
 
+// testPolygonHole
+template<>
+template<>
+void object::test<8>()
+{
+    checkHullByLenRatio(
+        "MULTIPOLYGON (((1 1, 10 3, 19 1, 16 8, 19 7, 19 19, 10 20, 8 17, 1 19, 1 1), (3 4, 5 10, 3 16, 9 14, 14 15, 15 9, 13 5, 3 4)))", 
+        0.9, 
+        "POLYGON ((10 20, 19 19, 19 7, 19 1, 10 3, 1 1, 1 19, 10 20), (13 5, 15 9, 14 15, 9 14, 3 16, 5 10, 3 4, 13 5))" );    
+}
+
+
+// testPolygonNestedPoly
+template<>
+template<>
+void object::test<9>()
+{
+    checkHullByLenRatio(
+        "MULTIPOLYGON (((1 1, 10 3, 19 1, 16 8, 19 7, 19 19, 10 20, 8 17, 1 19, 1 1), (3 4, 5 10, 3 16, 9 14, 14 15, 15 9, 13 5, 3 4)), ((6 10, 7 13, 10 12, 12 13, 13 11, 11 9, 13 8, 9 6, 6 6, 7 8, 6 10)))", 
+        0.9, 
+        "MULTIPOLYGON (((10 20, 19 19, 19 7, 19 1, 10 3, 1 1, 1 19, 10 20), (13 5, 15 9, 14 15, 9 14, 3 16, 5 10, 3 4, 13 5)), ((7 13, 10 12, 12 13, 13 11, 11 9, 13 8, 9 6, 6 6, 7 8, 6 10, 7 13)))" );    
+}
 
 } // namespace tut