Add GEOSDensify to CAPI

NEWS

@@ -6,6 +6,7 @@ Changes in 3.10.0
     library (#868, Paul Ramsey)
   - geosop CLI for GEOS (Martin Davis)
   - Full doxygen of the C-API (Paul Ramsey)
+  - CAPI: GEOSDensify (Brendan Ward)
 
 - Improvements:
   - Preserve ordering of lines in overlay results (Martin Davis)

capi/geos_c.cpp

@@ -427,6 +427,13 @@ extern "C" {
                                      joinStyle, mitreLimit);
     }
 
+    Geometry*
+    GEOSDensify(const Geometry* g, double tolerance)
+    {
+        return GEOSDensify_r(handle, g, tolerance);
+    }
+
+
     Geometry*
     GEOSSingleSidedBuffer(const Geometry* g, double width, int quadsegs,
                           int joinStyle, double mitreLimit, int leftSide)

capi/geos_c.h.in

@@ -582,6 +582,12 @@ extern GEOSGeometry GEOS_DLL *GEOSBufferWithStyle_r(
     double width, int quadsegs, int endCapStyle,
 	int joinStyle, double mitreLimit);
 
+/** \see GEOSDensify */
+extern GEOSGeometry GEOS_DLL *GEOSDensify_r(
+    GEOSContextHandle_t handle,
+	const GEOSGeometry* g,
+    double tolerance);
+
 /** \see GEOSOffsetCurve */
 extern GEOSGeometry GEOS_DLL *GEOSOffsetCurve_r(
     GEOSContextHandle_t handle,
@@ -2011,6 +2017,21 @@ extern GEOSGeometry GEOS_DLL *GEOSBufferWithStyle(
     int joinStyle,
     double mitreLimit);
 
+/**
+* Densifies a geometry using a given distance tolerance.
+* Additional vertices will be added to every line segment
+* that is greater this tolerance; these vertices will
+* evenly subdivide that segment.
+* Only linear components of input geometry are densified.
+* \param g The geometry to densify
+* \param tolerance the distance tolerance to densify
+* \return The densified geometry, or NULL on exception.
+* Caller is responsible for freeing with GEOSGeom_destroy().
+*/
+extern GEOSGeometry GEOS_DLL *GEOSDensify(
+    const GEOSGeometry* g,
+    double tolerance);
+
 /**
 * Generates offset curve for linear geometry.
 * Only LineStrings accepted as input.

capi/geos_ts_c.cpp

@@ -37,6 +37,7 @@
 #include <geos/geom/Coordinate.h>
 #include <geos/geom/IntersectionMatrix.h>
 #include <geos/geom/Envelope.h>
+#include <geos/geom/util/Densifier.h>
 #include <geos/index/strtree/SimpleSTRtree.h>
 #include <geos/index/strtree/GeometryItemDistance.h>
 #include <geos/index/ItemVisitor.h>
@@ -1110,6 +1111,20 @@ extern "C" {
         });
     }
 
+    Geometry*
+    GEOSDensify_r(GEOSContextHandle_t extHandle, const Geometry* g, double tolerance)
+    {
+        using geos::geom::util::Densifier;
+
+        return execute(extHandle, [&]() {
+            Densifier densifier(g);
+            densifier.setDistanceTolerance(tolerance);
+            auto g3 = densifier.getResultGeometry();
+            g3->setSRID(g->getSRID());
+            return g3.release();
+        });
+    }
+
     Geometry*
     GEOSOffsetCurve_r(GEOSContextHandle_t extHandle, const Geometry* g1, double width, int quadsegs, int joinStyle,
                       double mitreLimit)

src/geom/util/Densifier.cpp

@@ -18,6 +18,8 @@
  *
  **********************************************************************/
 
+#include <math.h>
+
 #include <geos/geom/util/Densifier.h>
 #include <geos/geom/CoordinateSequenceFactory.h>
 #include <geos/geom/CoordinateList.h>
@@ -113,7 +115,8 @@ Densifier::densifyPoints(const Coordinate::Vect pts, double distanceTolerance, c
         seg.p1 = *(it + 1);
         coordList.insert(coordList.end(), seg.p0, false);
         double len = seg.getLength();
-        int densifiedSegCount = (int)(len / distanceTolerance) + 1;
+
+        int densifiedSegCount = (int)(ceil(len / distanceTolerance));
         if(densifiedSegCount > 1) {
             double densifiedSegLen = len / densifiedSegCount;
             for(int j = 1; j < densifiedSegCount; j++) {
@@ -124,6 +127,10 @@ Densifier::densifyPoints(const Coordinate::Vect pts, double distanceTolerance, c
                 coordList.insert(coordList.end(), p, false);
             }
         }
+        else {
+            // no densification required; insert the last coordinate and continue
+            coordList.insert(coordList.end(), seg.p1, false);
+        }
     }
     coordList.insert(coordList.end(), pts[pts.size() - 1], false);
     return coordList.toCoordinateArray();
@@ -157,6 +164,10 @@ Densifier::setDistanceTolerance(double tol)
 Geometry::Ptr
 Densifier::getResultGeometry() const
 {
+    if (inputGeom->isEmpty()) {
+        return inputGeom->clone();
+    }
+
     DensifyTransformer dt(distanceTolerance);
     return dt.transform(inputGeom);
 }

tests/unit/capi/GEOSDensifyTest.cpp

@@ -0,0 +1,214 @@
+//
+// Test Suite for C-API GEOSBuffer and GEOSBufferWithStyle
+
+#include <tut/tut.hpp>
+// geos
+#include <geos_c.h>
+
+#include "capi_test_utils.h"
+
+
+namespace tut {
+//
+// Test Group
+//
+
+// Common data used in test cases.
+struct test_capigeosdensify_data : public capitest::utility {
+    GEOSGeometry* geom1_;
+    GEOSGeometry* geom2_;
+    GEOSWKTWriter* w_;
+    char* wkt_;
+
+    test_capigeosdensify_data()
+        : geom1_(nullptr), geom2_(nullptr), w_(nullptr)
+    {
+        initGEOS(notice, notice);
+        w_ = GEOSWKTWriter_create();
+        GEOSWKTWriter_setTrim(w_, 1);
+    }
+
+    ~test_capigeosdensify_data()
+    {
+        GEOSGeom_destroy(geom1_);
+        GEOSGeom_destroy(geom2_);
+        GEOSWKTWriter_destroy(w_);
+        geom1_ = nullptr;
+        geom2_ = nullptr;
+        wkt_ = nullptr;
+        finishGEOS();
+    }
+
+};
+
+typedef test_group<test_capigeosdensify_data> group;
+typedef group::object object;
+
+group test_capigeosdensify_group("capi::GEOSDensify");
+
+//
+// Test Cases
+//
+
+// Densify with a tolerance slightly larger than length of all edges.
+// Result should match inputs.
+template<>
+template<>
+void object::test<1>()
+{
+    geom1_ = GEOSGeomFromWKT("POLYGON ((0 0, 10 0, 10 10, 0 10, 0 0), (1 1, 1 2, 2 2, 2 1, 1 1))");
+    ensure(geom1_ != nullptr);
+    GEOSSetSRID(geom1_, 3857);
+
+    geom2_ = GEOSDensify(geom1_, 10.0);
+
+    ensure("result not null", geom2_ != nullptr);
+    ensure_geometry_equals(geom2_, geom1_);
+    ensure_equals("result SRID == expected SRID", GEOSGetSRID(geom2_), 3857);
+}
+
+
+// Densify with a tolerance that evenly subdivides all outer and inner edges.
+template<>
+template<>
+void object::test<2>()
+{
+    geom1_ = GEOSGeomFromWKT("POLYGON ((0 0, 10 0, 10 10, 0 10, 0 0), (1 1, 1 7, 7 7, 7 1, 1 1))");
+    ensure(geom1_ != nullptr);
+    GEOSSetSRID(geom1_, 3857);
+
+    geom2_ = GEOSDensify(geom1_, 5.0);
+
+    ensure("result not null", geom2_ != nullptr);
+    ensure_geometry_equals(
+        geom2_,
+        "POLYGON ((0 0, 5 0, 10 0, 10 5, 10 10, 5 10, 0 10, 0 5, 0 0), (1 1, 1 4, 1 7, 4 7, 7 7, 7 4, 7 1, 4 1, 1 1))");
+
+    ensure_equals("result SRID == expected SRID", GEOSGetSRID(geom2_), 3857);
+}
+
+// Densify a LINESTRING
+template<>
+template<>
+void object::test<3>()
+{
+    geom1_ = GEOSGeomFromWKT("LINESTRING (0 0, 0 6 )");
+    ensure(geom1_ != nullptr);
+    GEOSSetSRID(geom1_, 3857);
+
+    geom2_ = GEOSDensify(geom1_, 3);
+
+    ensure("result not null", geom2_ != nullptr);
+    ensure_geometry_equals(
+        geom2_,
+        "LINESTRING (0 0, 0 3, 0 6)");
+
+    ensure_equals("result SRID == expected SRID", GEOSGetSRID(geom2_), 3857);
+}
+
+// Ensure that tolerance results in the right number of subdivisions
+// ceil(6 / 2.9999999) = 3 new segments; 2 new vertices
+template<>
+template<>
+void object::test<4>()
+{
+    geom1_ = GEOSGeomFromWKT("LINESTRING (0 0, 0 6 )");
+    ensure(geom1_ != nullptr);
+    GEOSSetSRID(geom1_, 3857);
+
+    geom2_ = GEOSDensify(geom1_, 2.9999999);
+
+    ensure("result not null", geom2_ != nullptr);
+    ensure_geometry_equals(
+        geom2_,
+        "LINESTRING (0 0, 0 2, 0 4, 0 6)");
+
+    ensure_equals("result SRID == expected SRID", GEOSGetSRID(geom2_), 3857);
+}
+
+
+// Densify a LINEARRING
+template<>
+template<>
+void object::test<5>()
+{
+    geom1_ = GEOSGeomFromWKT("LINEARRING (0 0, 0 6, 6 6, 0 0)");
+    ensure(geom1_ != nullptr);
+    GEOSSetSRID(geom1_, 3857);
+
+    geom2_ = GEOSDensify(geom1_, 3.0);
+
+    ensure("result not null", geom2_ != nullptr);
+    ensure_geometry_equals(
+        geom2_,
+        "LINEARRING (0 0, 0 3, 0 6, 3 6, 6 6, 4 4, 2 2, 0 0)");
+    ensure_equals("result SRID == expected SRID", GEOSGetSRID(geom2_), 3857);
+}
+
+// Densify a POINT
+// Results should match inputs
+template<>
+template<>
+void object::test<6>()
+{
+    geom1_ = GEOSGeomFromWKT("POINT (0 0)");
+    ensure(geom1_ != nullptr);
+    GEOSSetSRID(geom1_, 3857);
+
+    geom2_ = GEOSDensify(geom1_, 3.0);
+
+    ensure("result not null", geom2_ != nullptr);
+    ensure_geometry_equals(geom2_, geom1_);
+    ensure_equals("result SRID == expected SRID", GEOSGetSRID(geom2_), 3857);
+}
+
+// Densify a MULTIPOINT
+// Results should match inputs
+template<>
+template<>
+void object::test<7>()
+{
+    geom1_ = GEOSGeomFromWKT("MULTIPOINT ((0 0), (10 10))");
+    ensure(geom1_ != nullptr);
+    GEOSSetSRID(geom1_, 3857);
+
+    geom2_ = GEOSDensify(geom1_, 3.0);
+
+    ensure("result not null", geom2_ != nullptr);
+    ensure_geometry_equals(geom2_, geom1_);
+    ensure_equals("result SRID == expected SRID", GEOSGetSRID(geom2_), 3857);
+}
+
+// Densify an empty polygon
+// Results should match inputs
+template<>
+template<>
+void object::test<8>()
+{
+    geom1_ = GEOSGeomFromWKT("POLYGON EMPTY");
+    ensure(geom1_ != nullptr);
+
+    geom2_ = GEOSDensify(geom1_, 3.0);
+
+    ensure("result not null", geom2_ != nullptr);
+    ensure_geometry_equals(geom2_, geom1_);
+}
+
+// Densify with an invalid tolerances should fail
+// Note: this raises "IllegalArgumentException: Tolerance must be positive:
+template<>
+template<>
+void object::test<9>()
+{
+    geom1_ = GEOSGeomFromWKT("POLYGON ((0 0, 10 0, 10 10, 0 10, 0 0))");
+    ensure(geom1_ != nullptr);
+
+    geom2_ = GEOSDensify(geom1_, 0.0);
+    ensure(geom2_ == nullptr);
+
+    geom2_ = GEOSDensify(geom1_, -1.0);
+    ensure(geom2_ == nullptr);
+}
+
+} // namespace tut
+