Add GEOSPolygonize_valid to CAPI

Closes #727

NEWS

@@ -4,9 +4,12 @@ Changes in 3.8.0
 - New things:
   - CAPI: GEOSBuildArea (#952, Even Rouault)
   - CAPI: GEOSMakeValid (#952, Even Rouault)
+  - CAPI: GEOSPolygonize_valid (#727, Dan Baston)
 
 - Improvements:
   - Improve performance and robustness of GEOSPointOnSurface (Martin Davis)
+  - Improve performance of GEOSPolygonize for cases with many potential
+    holes (#748, Dan Baston)
 
 
 Changes in 3.7.2

capi/geos_c.cpp

@@ -711,6 +711,12 @@ extern "C" {
         return GEOSPolygonize_r(handle, g, ngeoms);
     }
 
+    Geometry*
+    GEOSPolygonize_valid(const Geometry* const* g, unsigned int ngeoms)
+    {
+        return GEOSPolygonize_valid_r(handle, g, ngeoms);
+    }
+
     Geometry*
     GEOSPolygonizer_getCutEdges(const Geometry* const* g, unsigned int ngeoms)
     {

capi/geos_c.h.in

@@ -596,20 +596,20 @@ extern GEOSGeometry GEOS_DLL *GEOSClipByRect_r(GEOSContextHandle_t handle,
  * Polygonizes a set of Geometries which contain linework that
  * represents the edges of a planar graph.
  *
- * Any dimension of Geometry is handled - the constituent linework
- * is extracted to form the edges.
+ * All types of Geometry are accepted as input; the constituent
+ * linework is extracted as the edges to be polygonized.
  *
  * The edges must be correctly noded; that is, they must only meet
- * at their endpoints.
- * The Polygonizer will still run on incorrectly noded input
- * but will not form polygons from incorrectly noded edges.
+ * at their endpoints. Polygonization will accept incorrectly noded
+ * input but will not form polygons from non-noded edges, and reports
+ * them as errors.
  *
  * The Polygonizer reports the follow kinds of errors:
  *
  * - Dangles - edges which have one or both ends which are
  *   not incident on another edge endpoint
  * - Cut Edges - edges which are connected at both ends but
- *   which do not form part of polygon
+ *   which do not form part of a polygon
  * - Invalid Ring Lines - edges which form rings which are invalid
  *   (e.g. the component lines contain a self-intersection)
  *
@@ -618,11 +618,19 @@ extern GEOSGeometry GEOS_DLL *GEOSClipByRect_r(GEOSContextHandle_t handle,
  * collection. NULL is returned on exception. All returned
  * geometries must be destroyed by caller.
  *
+ * The GEOSPolygonize_valid_r variant allows extracting only polygons
+ * which form a valid polygonal result. The set of extracted polygons
+ * is guaranteed to be edge-disjoint. This is useful when it is known
+ * that the input lines form a valid polygonal geometry (which may
+ * include holes or nested polygons).
  */
 
 extern GEOSGeometry GEOS_DLL *GEOSPolygonize_r(GEOSContextHandle_t handle,
                               const GEOSGeometry *const geoms[],
                               unsigned int ngeoms);
+extern GEOSGeometry GEOS_DLL *GEOSPolygonize_valid_r(GEOSContextHandle_t handle,
+                                                     const GEOSGeometry *const geoms[],
+                                                     unsigned int ngems);
 extern GEOSGeometry GEOS_DLL *GEOSPolygonizer_getCutEdges_r(
                               GEOSContextHandle_t handle,
                               const GEOSGeometry * const geoms[],
@@ -1607,34 +1615,8 @@ extern GEOSGeometry GEOS_DLL *GEOSClipByRect(const GEOSGeometry* g, double xmin,
  * (both Geometries and pointers)
  */
 extern GEOSGeometry GEOS_DLL *GEOSPolygonize(const GEOSGeometry * const geoms[], unsigned int ngeoms);
+extern GEOSGeometry GEOS_DLL *GEOSPolygonize_valid(const GEOSGeometry * const geoms[], unsigned int ngeoms);
 extern GEOSGeometry GEOS_DLL *GEOSPolygonizer_getCutEdges(const GEOSGeometry * const geoms[], unsigned int ngeoms);
-/*
- * Polygonizes a set of Geometries which contain linework that
- * represents the edges of a planar graph.
- *
- * Any dimension of Geometry is handled - the constituent linework
- * is extracted to form the edges.
- *
- * The edges must be correctly noded; that is, they must only meet
- * at their endpoints.
- * The Polygonizer will still run on incorrectly noded input
- * but will not form polygons from incorrectly noded edges.
- *
- * The Polygonizer reports the follow kinds of errors:
- *
- * - Dangles - edges which have one or both ends which are
- *   not incident on another edge endpoint
- * - Cut Edges - edges which are connected at both ends but
- *   which do not form part of polygon
- * - Invalid Ring Lines - edges which form rings which are invalid
- *   (e.g. the component lines contain a self-intersection)
- *
- * Errors are reported to output parameters "cuts", "dangles" and
- * "invalid" (if not-null). Formed polygons are returned as a
- * collection. NULL is returned on exception. All returned
- * geometries must be destroyed by caller.
- *
- */
 extern GEOSGeometry GEOS_DLL *GEOSPolygonize_full(const GEOSGeometry* input,
     GEOSGeometry** cuts, GEOSGeometry** dangles, GEOSGeometry** invalid);
 

capi/geos_ts_c.cpp

@@ -3156,6 +3156,54 @@ extern "C" {
         return out;
     }
 
+    Geometry*
+    GEOSPolygonize_valid_r(GEOSContextHandle_t extHandle, const Geometry* const* g, unsigned int ngeoms)
+    {
+        if(0 == extHandle) {
+            return 0;
+        }
+
+        GEOSContextHandleInternal_t* handle = 0;
+        handle = reinterpret_cast<GEOSContextHandleInternal_t*>(extHandle);
+        if(0 == handle->initialized) {
+            return 0;
+        }
+
+        Geometry* out = 0;
+
+        try {
+            // Polygonize
+            using geos::operation::polygonize::Polygonizer;
+            Polygonizer plgnzr(true);
+            for(std::size_t i = 0; i < ngeoms; ++i) {
+                plgnzr.add(g[i]);
+            }
+
+            auto polys = plgnzr.getPolygons();
+            if (polys->empty()) {
+                out = handle->geomFactory->createGeometryCollection();
+            } else if (polys->size() == 1) {
+                out = (*polys)[0].release();
+            } else {
+                auto geoms = new std::vector<Geometry *>(polys->size());
+                for (size_t i = 0; i < polys->size(); i++) {
+                    (*geoms)[i] = (*polys)[i].release();
+                }
+
+                out = handle->geomFactory->createMultiPolygon(geoms);
+            }
+        }
+        catch(const std::exception& e) {
+            handle->ERROR_MESSAGE("%s", e.what());
+        }
+        catch(...) {
+            handle->ERROR_MESSAGE("Unknown exception thrown");
+        }
+
+        return out;
+    }
+
+
     Geometry*
     GEOSBuildArea_r(GEOSContextHandle_t extHandle, const Geometry* g)
     {

tests/unit/Makefile.am

@@ -160,7 +160,7 @@ geos_unit_SOURCES = \
 	capi/GEOSUserDataTest.cpp \
 	capi/GEOSPreparedGeometryTest.cpp \
 	capi/GEOSPointOnSurfaceTest.cpp \
-	capi/GEOSPolygonizer_getCutEdgesTest.cpp \
+	capi/GEOSPolygonizeTest.cpp \
 	capi/GEOSBufferTest.cpp \
 	capi/GEOSOffsetCurveTest.cpp \
 	capi/GEOSGeom_create.cpp \

tests/unit/capi/GEOSPolygonizeTest.cpp

@@ -0,0 +1,179 @@
+//
+// Test Suite for C-API GEOSPolygonize*
+
+#include <tut/tut.hpp>
+// geos
+#include <geos_c.h>
+// std
+#include <cstdarg>
+#include <cstdio>
+#include <cstdlib>
+#include <memory>
+
+namespace tut {
+//
+// Test Group
+//
+
+// Common data used in test cases.
+struct test_capigeospolygonize_data {
+    static void
+    notice(const char* fmt, ...)
+    {
+        std::fprintf(stdout, "NOTICE: ");
+
+        va_list ap;
+        va_start(ap, fmt);
+        std::vfprintf(stdout, fmt, ap);
+        va_end(ap);
+
+        std::fprintf(stdout, "\n");
+    }
+
+    test_capigeospolygonize_data()
+    {
+        initGEOS(notice, notice);
+    }
+
+    ~test_capigeospolygonize_data()
+    {
+        finishGEOS();
+    }
+
+};
+
+typedef test_group<test_capigeospolygonize_data> group;
+typedef group::object object;
+
+group test_capigeospolygonize_group("capi::GEOSPolygonize");
+
+//
+// Test Cases
+//
+
+template<>
+template<>
+void object::test<1>
+()
+{
+    constexpr int size = 2;
+    GEOSGeometry* geoms[size] = { nullptr };
+
+    geoms[0] = GEOSGeomFromWKT("LINESTRING(1 3, 3 3, 3 1, 1 1, 1 3)");
+    geoms[1] = GEOSGeomFromWKT("LINESTRING(1 3, 3 3, 3 1, 1 1, 1 3)");
+
+    GEOSGeometry* g = GEOSPolygonizer_getCutEdges(geoms, size);
+
+    ensure(nullptr != g);
+    ensure_equals(GEOSGetNumGeometries(g), size);
+
+    GEOSGeom_destroy(g);
+
+    for(auto& input : geoms) {
+        GEOSGeom_destroy(input);
+    }
+}
+
+template<>
+template<>
+void object::test<2>
+()
+{
+    constexpr int size = 6;
+    GEOSGeometry* geoms[size] = { nullptr };
+
+    // Example from JTS Developer's Guide, Chapter 6 - Polygonization
+    geoms[0] = GEOSGeomFromWKT("LINESTRING(0 0, 10 10)"); // isolated edge
+    geoms[1] = GEOSGeomFromWKT("LINESTRING(185 221, 100 100)"); // dangling edge
+    geoms[2] = GEOSGeomFromWKT("LINESTRING(185 221, 88 275, 180 316)");
+    geoms[3] = GEOSGeomFromWKT("LINESTRING(185 221, 292 281, 180 316)");
+    geoms[4] = GEOSGeomFromWKT("LINESTRING(189 98, 83 187, 185 221)");
+    geoms[5] = GEOSGeomFromWKT("LINESTRING(189 98, 325 168, 185 221)");
+
+    GEOSGeometry* g = GEOSPolygonizer_getCutEdges(geoms, size);
+
+    ensure(nullptr != g);
+    ensure_equals(GEOSGetNumGeometries(g), 0);
+
+    GEOSGeom_destroy(g);
+
+    for(auto& input : geoms) {
+        GEOSGeom_destroy(input);
+    }
+}
+
+template<>
+template<>
+void object::test<3>
+()
+{
+    constexpr int size = 2;
+    GEOSGeometry* geoms[size];
+    geoms[0] = GEOSGeomFromWKT("LINESTRING (100 100, 100 300, 300 300, 300 100, 100 100)");
+    geoms[1] = GEOSGeomFromWKT("LINESTRING (150 150, 150 250, 250 250, 250 150, 150 150)");
+
+    // GEOSPolygonize gives us a collection of two polygons
+    GEOSGeometry* g = GEOSPolygonize(geoms, size);
+    ensure(nullptr != g);
+    ensure_equals(GEOSGetNumGeometries(g), 2);
+    ensure_equals(GEOSGeomTypeId(g), GEOS_GEOMETRYCOLLECTION);
+    GEOSGeom_destroy(g);
+
+    // GEOSPolygonize_valid gives us a single polygon with a hole
+    g = GEOSPolygonize_valid(geoms, 2);
+
+    ensure(nullptr != g);
+    ensure_equals(GEOSGetNumGeometries(g), 1);
+    ensure_equals(GEOSGeomTypeId(g), GEOS_POLYGON);
+    GEOSGeom_destroy(g);
+
+    for(auto& input : geoms) {
+        GEOSGeom_destroy(input);
+    }
+}
+
+template<>
+template<>
+void object::test<4>
+()
+{
+    constexpr int size = 1;
+    GEOSGeometry* geoms[size];
+    geoms[0] = GEOSGeomFromWKT("LINESTRING (0 0, 1 1)");
+
+    GEOSGeometry* g = GEOSPolygonize_valid(geoms, size);
+
+    ensure(nullptr != g);
+    ensure_equals(GEOSGetNumGeometries(g), 0);
+    ensure_equals(GEOSGeomTypeId(g), GEOS_GEOMETRYCOLLECTION);
+    GEOSGeom_destroy(g);
+
+    for(auto& input : geoms) {
+        GEOSGeom_destroy(input);
+    }
+}
+
+template<>
+template<>
+void object::test<5>
+()
+{
+    constexpr int size = 2;
+    GEOSGeometry* geoms[size];
+    geoms[0] = GEOSGeomFromWKT("LINESTRING (0 0, 1 0, 1 1, 0 1, 0 0)");
+    geoms[1] = GEOSGeomFromWKT("LINESTRING (1 1, 2 1, 2 2, 1 2, 1 1)");
+
+    GEOSGeometry* g = GEOSPolygonize_valid(geoms, size);
+
+    ensure(nullptr != g);
+    ensure_equals(GEOSGetNumGeometries(g), 2);
+    ensure_equals(GEOSGeomTypeId(g), GEOS_MULTIPOLYGON);
+    GEOSGeom_destroy(g);
+
+    for(auto& input : geoms) {
+        GEOSGeom_destroy(input);
+    }
+}
+
+} // namespace tut
+