Use buffer hole erosion heuristic for rings

include/geos/operation/buffer/BufferCurveSetBuilder.h

@@ -71,6 +71,7 @@ namespace buffer { // geos.operation.buffer
  */
 class GEOS_DLL BufferCurveSetBuilder {
     using CoordinateSequence = geos::geom::CoordinateSequence;
+    using Envelope = geos::geom::Envelope;
 
 private:
 
@@ -193,7 +194,10 @@ class GEOS_DLL BufferCurveSetBuilder {
      * @param bufferDistance
      * @return
      */
-    bool isErodedCompletely(const geom::LinearRing* ringCoord,
+    bool isRingFullyEroded(const geom::LinearRing* ring, bool isHole,
+                            double bufferDistance);
+
+    bool isRingFullyEroded(const CoordinateSequence* ringCoord, const Envelope* env, bool isHole,
                             double bufferDistance);
 
     /**

src/operation/buffer/BufferCurveSetBuilder.cpp

@@ -237,7 +237,7 @@
 
     // optimization - don't bother computing buffer
     // if the polygon would be completely eroded
-    if(distance < 0.0 && isErodedCompletely(shell, distance)) {
+    if(distance < 0.0 && isRingFullyEroded(shell, false, distance)) {
 #if GEOS_DEBUG
         std::cerr << __FUNCTION__ << ": polygon is eroded completely " << std::endl;
 #endif
@@ -270,7 +270,7 @@
 
         // optimization - don't bother computing buffer for this hole
         // if the hole would be completely covered
-        if(distance > 0.0 && isErodedCompletely(hole, -distance)) {
+        if(distance > 0.0 && isRingFullyEroded(hole, true, distance)) {
             continue;
         }
 
@@ -292,14 +292,30 @@
 void
 BufferCurveSetBuilder::addRingBothSides(const CoordinateSequence* coord, double p_distance)
 {
-    addRingSide(coord, p_distance,
-                Position::LEFT,
-                Location::EXTERIOR, Location::INTERIOR);
-    /* Add the opposite side of the ring
-    */
-    addRingSide(coord, p_distance,
-                Position::RIGHT,
-                Location::INTERIOR, Location::EXTERIOR);
+    /*
+     * (f "hole" side will be eroded completely, avoid generating it.
+     * This prevents hole artifacts (e.g. https://github.com/libgeos/geos/issues/1223)
+     */
+    //-- distance is assumed positive, due to previous checks
+    Envelope env;
+    coord->expandEnvelope(env);
+    bool isHoleComputed = ! isRingFullyEroded(coord, &env, true, distance);
+
+    bool isCCW = isRingCCW(coord);
+
+    bool isShellLeft = ! isCCW;
+    if (isShellLeft || isHoleComputed) {
+        addRingSide(coord, p_distance,
+                    Position::LEFT,
+                    Location::EXTERIOR, Location::INTERIOR);
+    }
+
+    bool isShellRight = isCCW;
+    if (isShellRight || isHoleComputed) {
+        addRingSide(coord, p_distance,
+                    Position::RIGHT,
+                    Location::INTERIOR, Location::EXTERIOR);
+    }
 }
 
 
@@ -424,14 +440,22 @@
 
 /*private*/
 bool
-BufferCurveSetBuilder::isErodedCompletely(const LinearRing* ring,
+BufferCurveSetBuilder::isRingFullyEroded(const LinearRing* ring, bool isHole,
         double bufferDistance)
 {
     const CoordinateSequence* ringCoord = ring->getCoordinatesRO();
+    const Envelope* env = ring->getEnvelopeInternal();
+    return isRingFullyEroded(ringCoord, env, isHole, bufferDistance);
+}
 
+/*private*/
+bool
+BufferCurveSetBuilder::isRingFullyEroded(const CoordinateSequence* ringCoord, const Envelope* env, bool isHole,
+        double bufferDistance)
+{
     // degenerate ring has no area
     if(ringCoord->getSize() < 4) {
-        return bufferDistance < 0;
+        return true;
     }
 
     // important test to eliminate inverted triangle bug
@@ -440,10 +464,16 @@
         return isTriangleErodedCompletely(ringCoord, bufferDistance);
     }
 
-    const Envelope* env = ring->getEnvelopeInternal();
-    double envMinDimension = std::min(env->getHeight(), env->getWidth());
-    if(bufferDistance < 0.0 && 2 * std::abs(bufferDistance) > envMinDimension) {
-        return true;
+    bool isErodable = 
+        (  isHole && bufferDistance > 0) ||
+        (! isHole && bufferDistance < 0);
+
+    if (isErodable) {
+      //-- if envelope is narrower than twice the buffer distance, ring is eroded
+        double envMinDimension = std::min(env->getHeight(), env->getWidth());
+        if (2 * std::abs(bufferDistance) > envMinDimension) {
+            return true;
+        }
     }
     return false;
 }

tests/unit/operation/buffer/BufferOpTest.cpp

@@ -21,6 +21,8 @@
 #include <string>
 #include <vector>
 
+using namespace geos::operation::buffer;
+
 namespace tut {
 //
 // Test Group
@@ -39,11 +41,25 @@ struct test_bufferop_data {
     test_bufferop_data()
         : gf(*geos::geom::GeometryFactory::getDefaultInstance())
         , wktreader(&gf)
-        , default_quadrant_segments(geos::operation::buffer::BufferParameters::DEFAULT_QUADRANT_SEGMENTS)
+        , default_quadrant_segments(BufferParameters::DEFAULT_QUADRANT_SEGMENTS)
     {
         ensure_equals(default_quadrant_segments, int(8));
     }
 
+    std::unique_ptr<Geometry>
+    buffer(const std::string& wkt, double dist) {
+        std::unique_ptr<Geometry> geom = wktreader.read(wkt);
+        return BufferOp::bufferOp(geom.get(), dist);
+    }
+
+    void checkNumHoles(Geometry& geom, int nHoles) {
+        ensure_equals("Num Holes", dynamic_cast<const Polygon*>( &geom )->getNumInteriorRing(), std::size_t(nHoles));
+    }
+
+    void checkArea(Geometry& geom, double expectedArea, double tolerance) {
+        ensure_equals("Area", geom.getArea(), expectedArea, tolerance);
+    }
+
     void checkBufferEmpty(const std::string& wkt, double dist, bool isEmpty)
     {
         std::unique_ptr<Geometry> geom = wktreader.read(wkt);
@@ -596,4 +612,30 @@ void object::test<24>
     ensure_equals( (dynamic_cast<const Polygon*>( result.get() )->getNumInteriorRing()), 3u);
 }
 
+// testRingHoleEroded
+// See https://github.com/libgeos/geos/issues/1223
+template<>
+template<>
+void object::test<25>
+()
+{
+    std::string wkt("LINESTRING (25 44, 31 44, 32 38, 29 37, 25 37, 25 38, 24 40, 24 44, 25 44)");
+
+    std::unique_ptr<Geometry> result100 = buffer(wkt, 100);
+    checkValidPolygon(*result100);
+    checkNumHoles(*result100, 0);
+    checkArea(*result100, 34002, 100);
+
+    std::unique_ptr<Geometry> result10 = buffer(wkt, 10);
+    checkValidPolygon(*result10);
+    checkNumHoles(*result10, 0);
+    checkArea(*result10, 635.9, 1);
+
+    std::unique_ptr<Geometry> result2 = buffer(wkt, 2);
+    checkValidPolygon(*result2);
+    checkNumHoles(*result2, 1);
+    checkArea(*result2, 107, 1);
+
+}
+
 } // namespace tut