GridIntersection: Avoid passing self-touching rings to Polygonizer

Author: dbaston

include/geos/operation/grid/Cell.h

@@ -76,13 +76,14 @@ class Cell
      * @param c             Coordinate to process
      * @param prev_original The last *uninterpolated* coordinate preceding `c` in the
      *                      boundary being processed
+     * @param parentage an optional pointer indicating the source of the coordinate.
      *
      * @return `true` if the Coordinate was inside this cell, `false` otherwise
      */
-    bool take(const geom::CoordinateXY& c, const geom::CoordinateXY* prev_original = nullptr);
+    bool take(const geom::CoordinateXY& c, const geom::CoordinateXY* prev_original, const void* parentage);
 
   private:
-    std::vector<const std::vector<geom::CoordinateXY>*> getCoordLists() const;
+    std::vector<const Traversal*> getTraversals() const;
 
     enum class Location
     {

include/geos/operation/grid/Traversal.h

@@ -14,24 +14,26 @@
 
 #pragma once
 
-#include <vector>
-
+#include <geos/export.h>
 #include <geos/geom/Coordinate.h>
 #include <geos/operation/grid/Side.h>
 
+#include <vector>
+
 namespace geos::operation::grid {
 
 /**
  * @brief The Traversal class records the coordinates of a line that are
  *        within a grid cell, as well as the `Side` from which the line
  *        entered and exited the cell.
  */
-class Traversal
+class GEOS_DLL Traversal
 {
   public:
     Traversal()
       : m_entry{ Side::NONE }
       , m_exit{ Side::NONE }
+      , m_parentage{ nullptr }
     {
     }
 
@@ -48,7 +50,7 @@ class Traversal
     bool hasMultipleUniqueCoordinates() const;
 
     /// Begin a Traversal on the specified `Side`
-    void enter(const geom::CoordinateXY& c, Side s);
+    void enter(const geom::CoordinateXY& c, Side s, const void* parentage);
 
     /// Complete a Traversal on the specified `Side`
     void exit(const geom::CoordinateXY& c, Side s);
@@ -67,10 +69,13 @@ class Traversal
 
     const std::vector<geom::CoordinateXY>& getCoordinates() const { return m_coords; }
 
+    const void* getParentage() const { return m_parentage; }
+
   private:
     std::vector<geom::CoordinateXY> m_coords;
     Side m_entry;
     Side m_exit;
+    const void* m_parentage;
 };
 
 }

include/geos/operation/grid/TraversalAreas.h

@@ -19,6 +19,8 @@
 
 namespace geos::operation::grid {
 
+class Traversal;
+
 class GEOS_DLL TraversalAreas {
 public:
     /**
@@ -32,7 +34,7 @@ class GEOS_DLL TraversalAreas {
      * @return total area
      */
     static double
-    getLeftHandArea(const geom::Envelope& box, const std::vector<const std::vector<geom::CoordinateXY>*>& coord_lists);
+    getLeftHandArea(const geom::Envelope& box, const std::vector<const Traversal*>& coord_lists);
 
     /**
      * @brief Return an areal geometry representing the closed rings formed by this box and the provided Coordinate sequences
@@ -46,7 +48,7 @@ class GEOS_DLL TraversalAreas {
      * @return a Polygon or MultiPolygon geometry
      */
     static std::unique_ptr<geom::Geometry>
-    getLeftHandRings(const geom::GeometryFactory& gfact, const geom::Envelope& box, const std::vector<const std::vector<geom::CoordinateXY>*>& coord_lists);
+    getLeftHandRings(const geom::GeometryFactory& gfact, const geom::Envelope& box, const std::vector<const Traversal*>& coord_lists);
 
 };
 

src/operation/grid/Cell.cpp

@@ -191,7 +191,7 @@ Cell::getLastTraversal()
 }
 
 bool
-Cell::take(const CoordinateXY& c, const CoordinateXY* prev_original)
+Cell::take(const CoordinateXY& c, const CoordinateXY* prev_original, const void* parentage)
 {
     Traversal& t = traversal_in_progress();
 
@@ -201,7 +201,7 @@ Cell::take(const CoordinateXY& c, const CoordinateXY* prev_original)
         std::cout << "Entering " << m_box << " from " << getSide(c) << " at " << c << std::endl;
 #endif
 
-        t.enter(c, getSide(c));
+        t.enter(c, getSide(c), parentage);
         return true;
     }
 
@@ -258,32 +258,32 @@ Cell::isDetermined() const
     return false;
 }
 
-std::vector<const std::vector<CoordinateXY>*>
-Cell::getCoordLists() const
+std::vector<const Traversal*>
+Cell::getTraversals() const
 {
-    std::vector<const std::vector<CoordinateXY>*> coord_lists;
-    coord_lists.reserve(m_traversals.size());
+    std::vector<const Traversal*> traversals;
+    traversals.reserve(m_traversals.size());
 
     for (const auto& t : m_traversals) {
         if (t.isTraversed() || t.isClosedRing()) {
-            coord_lists.push_back(&t.getCoordinates());
+            traversals.push_back(&t);
         }
     }
 
-    return coord_lists;
+    return traversals;
 }
 
 double
 Cell::getCoveredFraction() const
 {
-    auto coord_lists = getCoordLists();
+    auto coord_lists = getTraversals();
     return TraversalAreas::getLeftHandArea(m_box, coord_lists) / getArea();
 }
 
 std::unique_ptr<Geometry>
 Cell::getCoveredPolygons(const GeometryFactory& gfact) const
 {
-    auto coord_lists = getCoordLists();
+    auto coord_lists = getTraversals();
     return TraversalAreas::getLeftHandRings(gfact, m_box, coord_lists);
 }
 

src/operation/grid/GridIntersection.cpp

@@ -308,8 +308,8 @@ collect_lengths(const Matrix<std::unique_ptr<Cell>>& cells)
     return lengths;
 }
 
-void
-traverse_cells(Matrix<std::unique_ptr<Cell>>& cells, std::vector<CoordinateXY>& coords, const Grid<infinite_extent>& ring_grid, bool areal)
+static void
+traverse_cells(Matrix<std::unique_ptr<Cell>>& cells, std::vector<CoordinateXY>& coords, const Grid<infinite_extent>& ring_grid, bool areal, const void* parentage)
 {
     size_t pos = 0;
     size_t row = ring_grid.getRow(coords.front().y);
@@ -323,7 +323,7 @@ traverse_cells(Matrix<std::unique_ptr<Cell>>& cells, std::vector<CoordinateXY>&
             const CoordinateXY* next_coord = last_exit ? last_exit : &coords[pos];
             const CoordinateXY* prev_coord = pos > 0 ? &coords[pos - 1] : nullptr;
 
-            cell.take(*next_coord, prev_coord);
+            cell.take(*next_coord, prev_coord, parentage);
 
             if (cell.getLastTraversal().isExited()) {
                 // Only push our exit coordinate if it's not same as the
@@ -437,7 +437,7 @@ GridIntersection::processLine(const LineString& ls, bool exterior_ring)
     }
 
     Matrix<std::unique_ptr<Cell>> cells(ring_grid.getNumRows(), ring_grid.getNumCols());
-    traverse_cells(cells, coordsVec, ring_grid, m_areal);
+    traverse_cells(cells, coordsVec, ring_grid, m_areal, &ls);
 
     // Compute the fraction covered for all cells and assign it to
     // the area matrix
@@ -473,8 +473,8 @@ traverse_ring(Matrix<std::unique_ptr<Cell>>& cells, const Grid<infinite_extent>&
 
     if (want_ccw != algorithm::Orientation::isCCW(&seq)) {
       std::reverse(coords.begin(), coords.end());
-   }
-    traverse_cells(cells, coords, grid, true);
+    }
+    traverse_cells(cells, coords, grid, true, &g);
 }
 
 void

src/operation/grid/Traversal.cpp

@@ -34,14 +34,15 @@ Traversal::isEmpty() const
 }
 
 void
-Traversal::enter(const CoordinateXY& c, Side s)
+Traversal::enter(const CoordinateXY& c, Side s, const void* parentage)
 {
     if (!m_coords.empty()) {
         throw std::runtime_error("Traversal already started");
     }
 
     add(c);
     m_entry = s;
+    m_parentage = parentage;
 }
 
 void