Port JTS locationtech/jts#655, Fix buffer to use largest enclosed are…

…a for invalid rings

NEWS

@@ -5,6 +5,7 @@ Changes in 3.9.1
 - Bug fixes / improvements:
   - Windows memory management quirk in createPolygon CAPI (#1050, Paul Ramsey)
   - Allow build on Apple ARM64 (Taras Zakharko)
+  - Fix buffer to use largest enclosed area for invalid rings (#732, Paul Ramsey)
 
 
 

include/geos/algorithm/Orientation.h

@@ -89,6 +89,32 @@ class GEOS_DLL Orientation {
     */
     static bool isCCW(const geom::CoordinateSequence* ring);
 
+    /**
+    * Tests if a ring defined by a CoordinateSequence is
+    * oriented counter-clockwise, using the signed area of the ring.
+    *
+    *  * The list of points is assumed to have the first and last points equal.
+    *  * This handles coordinate lists which contain repeated points.
+    *  * This handles rings which contain collapsed segments
+    *    (in particular, along the top of the ring).
+    *  * This handles rings which are invalid due to self-intersection
+    *
+    * This algorithm is guaranteed to work with valid rings.
+    * For invalid rings (containing self-intersections),
+    * the algorithm determines the orientation of
+    * the largest enclosed area (including overlaps).
+    * This provides a more useful result in some situations, such as buffering.
+    *
+    * However, this approach may be less accurate in the case of
+    * rings with almost zero area.
+    * (Note that the orientation of rings with zero area is essentially
+    * undefined, and hence non-deterministic.)
+    *
+    * @param ring a CoordinateSequence forming a ring (with first and last point identical)
+    * @return true if the ring is oriented counter-clockwise.
+    */
+    static bool isCCWArea(const geom::CoordinateSequence* ring);
+
 };
 
 

src/algorithm/Orientation.cpp

@@ -20,6 +20,7 @@
 #include <cmath>
 #include <vector>
 
+#include <geos/algorithm/Area.h>
 #include <geos/algorithm/Orientation.h>
 #include <geos/algorithm/CGAlgorithmsDD.h>
 #include <geos/geom/CoordinateSequence.h>
@@ -134,92 +135,12 @@ Orientation::isCCW(const geom::CoordinateSequence* ring)
     }
 }
 
-#if 0
 /* public static */
 bool
-Orientation::isCCW(const geom::CoordinateSequence* ring)
+Orientation::isCCWArea(const geom::CoordinateSequence* ring)
 {
-    // sanity check
-    if(ring->getSize() < 4) {
-        throw util::IllegalArgumentException("Ring has fewer than 4 points, so orientation cannot be determined");
-    }
-
-    // # of points without closing endpoint
-    const std::size_t nPts = ring->getSize() - 1;
-    assert(nPts >= 3); // This is here for scan-build
-
-    // find highest point
-    const geom::Coordinate* hiPt = &ring->getAt(0);
-    size_t hiIndex = 0;
-    for(std::size_t i = 1; i <= nPts; ++i) {
-        const geom::Coordinate* p = &ring->getAt(i);
-        if(p->y > hiPt->y) {
-            hiPt = p;
-            hiIndex = i;
-        }
-    }
-
-    // find distinct point before highest point
-    auto iPrev = hiIndex;
-    do {
-        if(iPrev == 0) {
-            iPrev = nPts;
-        }
-        iPrev = iPrev - 1;
-    }
-    while(ring->getAt(iPrev) == *hiPt && iPrev != hiIndex);
-
-    // find distinct point after highest point
-    auto iNext = hiIndex;
-    do {
-        iNext = (iNext + 1) % nPts;
-    }
-    while(ring->getAt(iNext) == *hiPt && iNext != hiIndex);
-
-    const geom::Coordinate* prev = &ring->getAt(iPrev);
-    const geom::Coordinate* next = &ring->getAt(iNext);
-
-    /*
-     * This check catches cases where the ring contains an A-B-A
-     * configuration of points.
-     * This can happen if the ring does not contain 3 distinct points
-     * (including the case where the input array has fewer than 4 elements),
-     * or it contains coincident line segments.
-     */
-    if(prev->equals2D(*hiPt) || next->equals2D(*hiPt) ||
-            prev->equals2D(*next)) {
-        return false;
-        // MD - don't bother throwing exception,
-        // since this isn't a complete check for ring validity
-        //throw  IllegalArgumentException("degenerate ring (does not contain 3 distinct points)");
-    }
-
-    int disc = Orientation::index(*prev, *hiPt, *next);
-
-    /*
-     *  If disc is exactly 0, lines are collinear.
-     * There are two possible cases:
-     *  (1) the lines lie along the x axis in opposite directions
-     *  (2) the lines lie on top of one another
-     *
-     *  (1) is handled by checking if next is left of prev ==> CCW
-     *  (2) should never happen, so we're going to ignore it!
-     *  (Might want to assert this)
-     */
-    bool isCCW = false;
-
-    if(disc == 0) {
-        // poly is CCW if prev x is right of next x
-        isCCW = (prev->x > next->x);
-    }
-    else {
-        // if area is positive, points are ordered CCW
-        isCCW = (disc > 0);
-    }
-
-    return isCCW;
+    return algorithm::Area::ofRingSigned(ring) < 0;
 }
-#endif
 
 
 } // namespace geos.algorithm

src/operation/buffer/OffsetCurveSetBuilder.cpp

@@ -320,8 +320,9 @@ OffsetCurveSetBuilder::addRingSide(const CoordinateSequence* coord,
 #if GEOS_DEBUG
     std::cerr << "OffsetCurveSetBuilder::addPolygonRing: CCW: " << Orientation::isCCW(coord) << std::endl;
 #endif
-    if(coord->size() >= LinearRing::MINIMUM_VALID_SIZE
-            && Orientation::isCCW(coord)) {
+    if(coord->size() >= LinearRing::MINIMUM_VALID_SIZE &&
+        Orientation::isCCWArea(coord))
+    {
         leftLoc = cwRightLoc;
         rightLoc = cwLeftLoc;
 #if GEOS_DEBUG

tests/unit/Makefile.am

@@ -44,7 +44,7 @@ geos_unit_SOURCES = \
 	algorithm/AngleTest.cpp \
 	algorithm/AreaTest.cpp \
 	algorithm/CGAlgorithms/computeOrientationTest.cpp \
-	algorithm/CGAlgorithms/isCCWTest.cpp \
+	algorithm/CGAlgorithms/OrientationIsCCWTest.cpp \
 	algorithm/CGAlgorithms/isPointInRingTest.cpp \
 	algorithm/CGAlgorithms/signedAreaTest.cpp \
 	algorithm/ConvexHullTest.cpp \

tests/unit/algorithm/CGAlgorithms/isCCWTest.cpp → tests/unit/algorithm/CGAlgorithms/OrientationIsCCWTest.cpp

@@ -40,7 +40,7 @@ struct test_isccw_data {
     }
 
     void
-    checkOrientationCCW(bool expectedCCW, const std::string& wkt)
+    checkCCW(bool expectedCCW, const std::string& wkt)
     {
         GeometryPtr geom(reader_.read(wkt));
         geos::geom::Polygon* poly = dynamic_cast<geos::geom::Polygon*>(geom.get());
@@ -50,6 +50,17 @@ struct test_isccw_data {
         ensure_equals("CoordinateSequence isCCW", expectedCCW, actualCCW);
     }
 
+    void
+    checkCCWArea(bool expectedCCWArea, const std::string& wkt)
+    {
+        GeometryPtr geom(reader_.read(wkt));
+        geos::geom::Polygon* poly = dynamic_cast<geos::geom::Polygon*>(geom.get());
+        ensure("WKT must be POLYGON)", poly != nullptr);
+        const geos::geom::CoordinateSequence* cs = poly->getExteriorRing()->getCoordinatesRO();
+        bool actualCCWArea = Orientation::isCCWArea(cs);
+        ensure_equals("CoordinateSequence isCCWArea", expectedCCWArea, actualCCWArea);
+    }
+
     void
     checkHexOrientationCCW(bool expectedCCW, std::istringstream& wkt)
     {
@@ -64,7 +75,7 @@ struct test_isccw_data {
 typedef test_group<test_isccw_data> group;
 typedef group::object object;
 
-group test_isccw_group("geos::algorithm::CGAlgorithms::isCCW");
+group test_isccw_group("geos::algorithm::CGAlgorithms::OrientationIsCCW");
 
 //
 // Test Cases
@@ -77,7 +88,7 @@ void object::test<1>
 ()
 {
     const std::string wkt("POLYGON ((60 180, 140 240, 140 240, 140 240, 200 180, 120 120, 60 180))");
-    checkOrientationCCW(false, wkt);
+    checkCCW(false, wkt);
 }
 
 // 2 - Test if coordinates of polygon are counter-clockwise oriented
@@ -87,7 +98,7 @@ void object::test<2>
 ()
 {
     const std::string wkt("POLYGON ((60 180, 140 120, 100 180, 140 240, 60 180))");
-    checkOrientationCCW(true, wkt);
+    checkCCW(true, wkt);
 }
 
 // 3 - Test the same polygon as in test No 2 but with duplicated top point
@@ -97,7 +108,7 @@ void object::test<3>
 ()
 {
     const std::string wkt("POLYGON ((60 180, 140 120, 100 180, 140 240, 140 240, 60 180))");
-    checkOrientationCCW(true, wkt);
+    checkCCW(true, wkt);
 }
 
 // 4 - Test orientation the narrow (almost collapsed) ring
@@ -133,7 +144,7 @@ void object::test<6>
 ()
 {
     const std::string wkt("POLYGON ((1 1, 9 1, 5 9, 1 1))");
-    checkOrientationCCW(true, wkt);
+    checkCCW(true, wkt);
 }
 
 // testFlatTopSegment
@@ -143,7 +154,7 @@ void object::test<7>
 ()
 {
     const std::string wkt("POLYGON ((100 200, 200 200, 200 100, 100 100, 100 200))");
-    checkOrientationCCW(false, wkt);
+    checkCCW(false, wkt);
 }
 
 // testFlatMultipleTopSegment
@@ -153,7 +164,7 @@ void object::test<8>
 ()
 {
     const std::string wkt("POLYGON ((100 200, 127 200, 151 200, 173 200, 200 200, 100 100, 100 200))");
-    checkOrientationCCW(false, wkt);
+    checkCCW(false, wkt);
 }
 
 // testDegenerateRingHorizontal
@@ -163,7 +174,7 @@ void object::test<9>
 ()
 {
     const std::string wkt("POLYGON ((100 200, 100 200, 200 200, 100 200))");
-    checkOrientationCCW(false, wkt);
+    checkCCW(false, wkt);
 }
 
 // testDegenerateRingAngled
@@ -173,7 +184,7 @@ void object::test<10>
 ()
 {
     const std::string wkt("POLYGON ((100 100, 100 100, 200 200, 100 100))");
-    checkOrientationCCW(false, wkt);
+    checkCCW(false, wkt);
 }
 
 // testDegenerateRingVertical
@@ -183,7 +194,7 @@ void object::test<11>
 ()
 {
     const std::string wkt("POLYGON ((200 100, 200 100, 200 200, 200 100))");
-    checkOrientationCCW(false, wkt);
+    checkCCW(false, wkt);
 }
 
 /**
@@ -196,7 +207,7 @@ void object::test<12>
 ()
 {
     const std::string wkt("POLYGON ((10 20, 61 20, 20 30, 50 60, 10 20))");
-    checkOrientationCCW(false, wkt);
+    checkCCW(false, wkt);
 }
 
 // testABATopFlatSegmentCollapse
@@ -206,7 +217,7 @@ void object::test<13>
 ()
 {
     const std::string wkt("POLYGON ((71 0, 40 40, 70 40, 40 40, 20 0, 71 0))");
-    checkOrientationCCW(true, wkt);
+    checkCCW(true, wkt);
 }
 
 // testABATopFlatSegmentCollapseMiddleStart
@@ -216,7 +227,7 @@ void object::test<14>
 ()
 {
     const std::string wkt("POLYGON ((90 90, 50 90, 10 10, 90 10, 50 90, 90 90))");
-    checkOrientationCCW(true, wkt);
+    checkCCW(true, wkt);
 }
 
 // testMultipleTopFlatSegmentCollapseSinglePoint
@@ -226,7 +237,7 @@ void object::test<15>
 ()
 {
     const std::string wkt("POLYGON ((100 100, 200 100, 150 200, 170 200, 200 200, 100 200, 150 200, 100 100))");
-    checkOrientationCCW(true, wkt);
+    checkCCW(true, wkt);
 }
 
 // testMultipleTopFlatSegmentCollapseFlatTop
@@ -236,7 +247,7 @@ void object::test<16>
 ()
 {
     const std::string wkt("POLYGON ((10 10, 90 10, 70 70, 90 70, 10 70, 30 70, 50 70, 10 10))");
-    checkOrientationCCW(true, wkt);
+    checkCCW(true, wkt);
 }
 
 

tests/unit/operation/buffer/BufferOpTest.cpp

@@ -3,6 +3,7 @@
 
 // tut
 #include <tut/tut.hpp>
+#include <utility.h>
 // geos
 #include <geos/operation/buffer/BufferOp.h>
 #include <geos/operation/buffer/BufferParameters.h>
@@ -447,5 +448,22 @@ void object::test<14>
     ensure_equals(gBuffer->getGeometryTypeId(), geos::geom::GEOS_MULTIPOLYGON);
 }
 
+
+// This now works since buffer ring orientation is changed to use signed-area test.
+// testBowtiePolygonLargestAreaRetained
+template<>
+template<>
+void object::test<15>
+()
+{
+    std::string wkt0("POLYGON ((10 10, 50 10, 25 35, 35 35, 10 10))");
+    GeomPtr g0(wktreader.read(wkt0));
+    GeomPtr gresult = g0->buffer(0.0);
+    std::string wkt1("POLYGON ((10 10, 30 30, 50 10, 10 10))");
+    GeomPtr gexpected(wktreader.read(wkt1));
+    ensure_equals_geometry(gresult.get(), gexpected.get());
+}
+
+
 } // namespace tut