GridIntersection:: Add extra nodes to rectangular Polygon outputs so that coverage is valid

Author: dbaston

include/geos/operation/grid/Cell.h

@@ -40,6 +40,9 @@ class Cell
     {
     }
 
+    /// Get all points that fall on the specified side of this cell.
+    void getEdgePoints(Side s, std::vector<geom::CoordinateXY>& points) const;
+
     const geom::Envelope& box() const { return m_box; }
 
     double getWidth() const;

include/geos/operation/grid/Traversal.h

@@ -61,6 +61,8 @@ class GEOS_DLL Traversal
 
     Side getExitSide() const { return m_exit; }
 
+    const geom::CoordinateXY& getFirstCoordinate() const;
+
     const geom::CoordinateXY& getLastCoordinate() const;
 
     const geom::CoordinateXY& getExitCoordinate() const;

src/operation/grid/Cell.cpp

@@ -291,4 +291,20 @@ Cell::getCoveredPolygons(const GeometryFactory& gfact) const
     return TraversalAreas::getLeftHandRings(gfact, m_box, coord_lists);
 }
 
+void Cell::getEdgePoints(Side s, std::vector<CoordinateXY> &edgePoints) const
+{
+    for (const Traversal& t : m_traversals) {
+        const auto& coords = t.getCoordinates();
+
+        for (const auto& c : coords) {
+            if ((s == Side::LEFT && c.x == m_box.getMinX()) ||
+                (s == Side::RIGHT && c.x == m_box.getMaxX()) ||
+                (s == Side::BOTTOM && c.y == m_box.getMinY()) ||
+                (s == Side::TOP && c.y == m_box.getMaxY())) {
+                edgePoints.push_back(c);
+            }
+        }
+    }
+}
+
 }

src/operation/grid/GridIntersection.cpp

@@ -22,6 +22,7 @@
 #include <geos/operation/grid/Cell.h>
 #include <geos/operation/grid/FloodFill.h>
 #include <geos/operation/grid/GridIntersection.h>
+#include <geos/operation/grid/PerimeterDistance.h>
 #include <geos/operation/overlayng/CoverageUnion.h>
 #include <geos/util.h>
 
@@ -504,6 +505,58 @@ traverse_polygons(Matrix<std::unique_ptr<Cell>>& cells, const Grid<infinite_exte
     }
 }
 
+// Create a rectangular polygon from a set of points lying on the boundary
+// of a specified envelope. The provided points do not need to include the
+// corners of the envelope itself.
+static std::unique_ptr<Geometry>
+createRectangle(const geom::GeometryFactory& gfact, const Envelope& env, std::vector<CoordinateXY>& points)
+{
+    if (points.empty()) {
+        return gfact.toGeometry(&env);
+    } else {
+        auto perimeterDistanceCmp = [&env](const CoordinateXY& a, const CoordinateXY& b) {
+            return PerimeterDistance::getPerimeterDistance(env, a) <
+                PerimeterDistance::getPerimeterDistance(env, b);
+        };
+
+        points.emplace_back(env.getMinX(), env.getMinY());
+        points.emplace_back(env.getMinX(), env.getMaxY());
+        points.emplace_back(env.getMaxX(), env.getMinY());
+        points.emplace_back(env.getMaxX(), env.getMaxY());
+
+        std::sort(points.begin(), points.end(), perimeterDistanceCmp);
+        points.push_back(points.front());
+        auto seq = std::make_unique<geom::CoordinateSequence>(0, false, false);
+        seq->reserve(points.size());
+        for (const auto& c : points) {
+            seq->add(c, false);
+        }
+
+        auto ring = gfact.createLinearRing(std::move(seq));
+        return gfact.createPolygon(std::move(ring));
+    }
+}
+
+// Get any points from adjacent cells that are also on the boundary of this cell.
+static std::vector<CoordinateXY>
+getExtraNodes(const Matrix<std::unique_ptr<Cell> > &cells, size_t i, size_t j) {
+    std::vector<CoordinateXY> points;
+
+    if (const Cell *above = cells(i - 1, j).get()) {
+        above->getEdgePoints(Side::BOTTOM, points);
+    }
+    if (const Cell *below = cells(i + 1, j).get()) {
+        below->getEdgePoints(Side::TOP, points);
+    }
+    if (const Cell *left = cells(i, j - 1).get()) {
+        left->getEdgePoints(Side::RIGHT, points);
+    }
+    if (const Cell *right = cells(i, j + 1).get()) {
+        right->getEdgePoints(Side::LEFT, points);
+    }
+
+    return points;
+}
 
 std::unique_ptr<Geometry>
 GridIntersection::subdividePolygon(const Grid<bounded_extent>& p_grid, const Geometry& g, bool includeExterior)
@@ -547,9 +600,11 @@ GridIntersection::subdividePolygon(const Grid<bounded_extent>& p_grid, const Geo
                 }
             } else if (!edge && areas(i - 1, j - 1) == fill_values<float>::INTERIOR) {
                 // Cell is entirely covered by polygon
-                // TODO: Add nodes from adjacent polygons?
+                // In order to have the outputs forms a properly noded coverage,
+                // we need to add nodes from adjacent polygons.
                 Envelope env = cell_grid.getCellEnvelope(i, j);
-                geoms.push_back(gfact.toGeometry(&env));
+                std::vector<CoordinateXY> points = getExtraNodes(cells, i, j);
+                geoms.push_back(createRectangle(gfact, env, points));
             }
         }
     }

src/operation/grid/Traversal.cpp

@@ -96,10 +96,16 @@ Traversal::isTraversed() const
     return isEntered() && isExited();
 }
 
+const CoordinateXY&
+Traversal::getFirstCoordinate() const
+{
+    return m_coords.front();
+}
+
 const CoordinateXY&
 Traversal::getLastCoordinate() const
 {
-    return m_coords.at(m_coords.size() - 1);
+    return m_coords.back();
 }
 
 const CoordinateXY&

tests/unit/operation/grid/GridIntersectionTest.cpp

@@ -696,7 +696,7 @@ void object::test<37>()
 
     auto subdivided = GridIntersection::subdividePolygon(ext, *g, true);
 
-    ensure_equals("", subdivided->getArea(), g->getArea(), 1e-8);
+    check_subdivided_polygon(*g, *subdivided);
 }
 
 template<>
@@ -710,9 +710,8 @@ void object::test<38>()
 
     auto g = wkt_reader_.read("POINT (5 5)")->buffer(20);
 
-    auto subdivided = GridIntersection::subdividePolygon(ext, *g, true);
-
-    ensure_equals("", subdivided->getArea(), g->getArea(), 1e-8);
+    auto subd = GridIntersection::subdividePolygon(ext, *g, true);
+    check_subdivided_polygon(*g, *subd);
 }
 
 template<>
@@ -726,9 +725,8 @@ void object::test<39>()
 
     auto g = GeometryFactory::getDefaultInstance()->toGeometry(&e);
 
-    auto subdivided = GridIntersection::subdividePolygon(ext, *g, false);
-
-    ensure_equals("", subdivided->getArea(), e.getArea(), 1e-8);
+    auto subd = GridIntersection::subdividePolygon(ext, *g, false);
+    check_subdivided_polygon(*g, *subd);
 }
 
 template<>
@@ -742,9 +740,8 @@ void object::test<40>()
 
     auto g = wkt_reader_.read("MULTIPOLYGON (((1 1, 15 1, 15 25, 1 25, 1 1), (12 12, 12 14, 14 14, 14 12, 12 12)), ((16 1, 25 1, 25 25, 16 25, 16 1)))");
 
-    auto subdivided = GridIntersection::subdividePolygon(ext, *g, false);
-
-    ensure_equals("", subdivided->getArea(), g->getArea(), 1e-8);
+    auto subd = GridIntersection::subdividePolygon(ext, *g, false);
+    check_subdivided_polygon(*g, *subd);
 }
 
 template<>
@@ -827,10 +824,7 @@ void object::test<45>()
     check_area(*rci, ext, *g);
 
     auto subd = GridIntersection::subdividePolygon(ext, *g, false);
-    // TODO: Currently fails because polygons do not form a valid coverage.
-    // When we create a rectangle polygon for a fully interior cell, we do not add interior nodes
-    // along its edges if they are present in the input geometry.
-    //check_subdivided_polygon(*g, *subd);
+    check_subdivided_polygon(*g, *subd);
 }
 
 template<>
@@ -849,10 +843,7 @@ void object::test<46>()
 
     auto subd = GridIntersection::subdividePolygon(ext, *g, false);
 
-    // TODO: Currently fails because polygons do not form a valid coverage.
-    // When we create a rectangle polygon for a fully interior cell, we do not add interior nodes
-    // along its edges if they are present in the input geometry.
-    //check_subdivided_polygon(*g, *subd);
+    check_subdivided_polygon(*g, *subd);
 }
 
 template<>