Write islands as separate polys in WKT (#2549)

filters/private/hexer/HexGrid.cpp

@@ -368,27 +368,23 @@ void HexGrid::cleanPossibleRoot(Segment s, Path *p)
     }
 }
 
-void HexGrid::toWKT(std::ostream& output) const
+void HexGrid::toWKT(std::ostream& out) const
 {
-    auto outputPath = [this,&output](size_t pathNum)
+    auto writePath = [this, &out](size_t pathNum)
     {
-        Path *p = rootPaths()[pathNum];
-
-        output << "(";
-        p->toWKT(output);
-        output << ")";
+       rootPaths()[pathNum]->toWKT(out);
     };
 
-    output << "MULTIPOLYGON (";
+    out << "MULTIPOLYGON (";
 
     if (rootPaths().size())
-        outputPath(0);
+        writePath(0);
     for (size_t pi = 1; pi < rootPaths().size(); ++pi)
     {
-        output << ",";
-        outputPath(pi);
+        out << ",";
+        writePath(pi);
     }
-    output << ")";
+    out << ")";
 }
 
 size_t HexGrid::densePointCount() const

filters/private/hexer/HexGrid.hpp

@@ -37,6 +37,8 @@
 #include <unordered_map>
 #include <set>
 
+#include <pdal/pdal_export.hpp>
+
 #include "exception.hpp"
 #include "Hexagon.hpp"
 #include "Mathpair.hpp"
@@ -54,7 +56,7 @@ class HexGrid
 {
     friend class HexIter;
 public:
-    HexGrid(int dense_limit);
+    PDAL_DLL HexGrid(int dense_limit);
     HexGrid(double height, int dense_limit) : m_dense_limit(dense_limit)
         { initialize(height); }
 
@@ -64,20 +66,24 @@ class HexGrid
             delete m_paths[i];
     }
 
+    // Exported for testing.
+    PDAL_DLL void findShapes();
+    PDAL_DLL void findParentPaths();
+    PDAL_DLL void toWKT(std::ostream&) const;
+    PDAL_DLL void addDenseHexagon(int x, int y);
+
     bool dense(Hexagon *h);
     void addPoint(double x, double y)
         { addPoint(Point(x, y)); }
     void addPoint(Point p);
     void processSample();
-    void findShapes();
-    void findParentPaths();
+
     void extractShapes();
     void dumpInfo();
     void drawHexagons();
     Hexagon *getHexagon(int x, int y);
     Hexagon *getHexagon(const Coord& c)
         { return getHexagon(c.m_x, c.m_y); }
-    void addDenseHexagon(int x, int y);
     HexIter hexBegin();
     HexIter hexEnd();
     double width() const
@@ -96,7 +102,6 @@ class HexGrid
         { return m_paths; }
     void setSampleSize(unsigned sampleSize)
         { m_maxSample = sampleSize; }
-    void toWKT(std::ostream&) const;
     size_t densePointCount() const;
 
 private:

filters/private/hexer/Path.cpp

@@ -34,6 +34,8 @@
 
 #include "Path.hpp"
 
+#include <cassert>
+
 using namespace std;
 
 namespace hexer
@@ -66,34 +68,63 @@ vector<Point> Path::points() const
     return points;
 }
 
-void Path::toWKT( std::ostream& output) const
+void Path::writeRing(std::ostream& out) const
 {
-    vector<Point> pts = points();
-
-    auto outputPoint = [&output](Point& p)
+    auto outputPoint = [&out](const Point& p)
     {
-        output << p.m_x << " " << p.m_y;
+        out << p.m_x << " " << p.m_y;
     };
 
-    output << "(";
+    const vector<Point>& pts = points();
+    assert(pts.size() > 2);
+    out << "(";
+    outputPoint(pts.front());
+    for (auto it = pts.begin() + 1; it != pts.end(); ++it)
+    {
+        out << ", ";
+        outputPoint(*it);
+    }
+    out << ")";
+}
+
+// WKT (or GeoJSON) doesn't allow nesting of polygons.  You can just have
+// polygons and holes.  Islands within the holes need to be described as
+// separate polygons.  To that end, we gather the islands from all holes
+// and return them to be processed as separate polygons.
+PathPtrList Path::writePolygon(std::ostream& out) const
+{
+    PathPtrList islands;
 
-    auto pi = pts.begin();
-    if (pi != pts.end())
-        outputPoint(*pi++);
-    for (; pi != pts.end(); ++pi)
+    out << "(";
+    writeRing(out);
+    const PathPtrList& paths = subPaths();
+    for (auto& p : paths)
     {
-        output << ", ";
-        outputPoint(*pi);
+        out << ", ";
+        p->writeRing(out);
+        const PathPtrList& subs(p->subPaths());
+        islands.insert(islands.end(), subs.begin(), subs.end());
     }
+    out << ")";
+    return islands;
+}
 
-    output << ")";
 
-    vector<Path *> paths = subPaths();
-    for (size_t pi = 0; pi != paths.size(); ++pi)
+void Path::toWKT(std::ostream& out) const
+{
+    PathPtrList islands = writePolygon(out);
+
+    // See the note on writePolygon()
+    while (islands.size())
     {
-        Path* p = paths[pi];
-        output <<",";
-        p->toWKT(output);
+        PathPtrList paths;
+        paths.swap(islands);
+        for (Path *p : paths)
+        {
+            out << ", ";
+            PathPtrList subIslands = p->writePolygon(out);
+            islands.insert(islands.end(), subIslands.begin(), subIslands.end());
+        }
     }
 }
 

filters/private/hexer/Path.hpp

@@ -49,6 +49,8 @@ enum Orientation
 };
 
 class HexGrid;
+class Path;
+using PathPtrList = std::vector<Path *>;
 
 class Path
 {
@@ -59,8 +61,8 @@ class Path
 
     ~Path()
     {
-        for (std::vector<Path*>::size_type i = 0; i < m_children.size(); ++i)
-            delete m_children[i];
+        for (auto p : m_children)
+            delete p;
     }
 
     void push_back(const Segment& s)
@@ -85,7 +87,7 @@ class Path
     Orientation orientation() const
         { return m_orientation; }
     std::vector<Point> points() const;
-    std::vector<Path *> subPaths() const
+    PathPtrList subPaths() const
         { return m_children; }
     void toWKT( std::ostream& output) const;
 
@@ -95,12 +97,15 @@ class Path
     /// Parent path (NULL if root path)
     Path *m_parent;
     /// Children
-    std::vector<Path *> m_children;
+    PathPtrList m_children;
     /// Orientation of path AT EXTRACTION - segments are ALWAYS ordered
     /// clockwise.
     Orientation m_orientation;
     /// List of segments that make up the path.
     std::vector<Segment> m_segs;
+
+    void writeRing(std::ostream& out) const;
+    PathPtrList writePolygon(std::ostream& out) const;
 };
 
 } //namespace hexer

test/unit/filters/HexbinFilterTest.cpp

@@ -35,6 +35,7 @@
 #include <pdal/pdal_test_main.hpp>
 
 #include <io/LasReader.hpp>
+#include <filters/private/hexer/HexGrid.hpp>
 
 #include <pdal/SpatialReference.hpp>
 #include <pdal/StageFactory.hpp>
@@ -94,3 +95,40 @@ TEST(HexbinFilterTest, HexbinFilterTest_test_1)
     out.close();
     FileUtils::deleteFile(filename);
 }
+
+// Test that we create proper WKT for geometry with islands.
+TEST(HexbinFilterTest, issue_2507)
+{
+    hexer::HexGrid grid(1);
+
+    // This is an arrangement with two holes.  One of the holes has two
+    // islands and one of those islands has a hole.
+    // Here's a link to the picture.  The green numbers indicate the hexes
+    // listed here:
+    // https://photos.app.goo.gl/P3B3mU4Zre6zADEQ6
+    std::vector<std::pair<int, int>> hexes {
+        {0, 3}, {0, 4}, {0,5}, {0, 6},
+        {1, 2}, {1, 6},
+        {2, 2}, {2, 4}, {2, 5}, {2, 7},
+        {3, 1}, {3, 3}, {3, 5}, {3, 7},
+        {4, 1}, {4, 2}, {4, 4}, {4, 5}, {4, 8},
+        {5, 0}, {5, 2}, {5, 6}, {5, 8},
+        {6, 1}, {6, 3}, {6, 4}, {6, 8},
+        {7, 1}, {7, 3}, {7, 4}, {7, 5}, {7, 7},
+        {8, 2}, {8, 3}, {8, 4}, {8, 5}, {8, 6}, {8, 7}
+    };
+
+    for (auto p : hexes)
+        grid.addDenseHexagon(p.first, p.second);
+    grid.findShapes();
+    grid.findParentPaths();
+
+    std::ostringstream oss;
+    grid.toWKT(oss);
+    std::string s(oss.str());
+
+    std::string test =
+    R"delim(MULTIPOLYGON (((-5 -0.5, -5 -0.5, -6 -1, -6 -1, -7 -1.5, -7 -1.5, -8 -2, -8 -2, -8 -2, -8 -3, -8 -3, -8 -4, -8 -4, -8 -5, -8 -5, -8 -6, -8 -6, -8 -7, -8 -7, -8 -7, -7 -7.5, -7 -7.5, -6 -8, -6 -8, -5 -8.5, -5 -8.5, -5 -8.5, -4 -8, -4 -8, -3 -7.5, -3 -7.5, -2 -7, -2 -7, -1 -6.5, -1 -6.5, 0 -6, 0 -6, 0 -6, 0 -5, 0 -5, 0 -4, 0 -4, 0 -3, 0 -3, 0 -3, -1 -2.5, -1 -2.5, -2 -2, -2 -2, -3 -1.5, -3 -1.5, -4 -1, -4 -1, -5 -0.5, -5 -0.5), (-4 -2, -4 -1, -5 -0.5, -6 -1, -6 -1, -7 -1.5, -7 -1.5, -8 -2, -8 -3, -7 -3.5, -6 -3, -6 -3, -5 -2.5, -5 -2.5, -4 -2), (0 -6, 0 -5, 0 -5, 0 -4, 0 -4, 0 -3, -1 -2.5, -1 -2.5, -2 -2, -2 -2, -3 -1.5, -4 -2, -4 -2, -5 -2.5, -5 -2.5, -6 -3, -6 -4, -6 -4, -6 -4, -7 -4.5, -7 -5.5, -7 -5.5, -7 -5.5, -8 -6, -8 -7, -7 -7.5, -7 -7.5, -6 -8, -6 -8, -5 -8.5, -4 -8, -4 -8, -3 -7.5, -3 -7.5, -2 -7, -2 -7, -1 -6.5, -1 -6.5, 0 -6)), ((-3 -3.5, -3 -3.5, -4 -4, -4 -4, -4 -4, -4 -5, -4 -5, -4 -5, -3 -5.5, -3 -5.5, -3 -5.5, -2 -5, -2 -5, -2 -5, -2 -4, -2 -4, -2 -4, -3 -3.5, -3 -3.5), (-2 -5, -2 -4, -3 -3.5, -4 -4, -4 -5, -3 -5.5, -2 -5)), ((-5 -6.5, -5 -6.5, -5 -6.5, -5 -6.5, -5 -6.5, -5 -6.5, -5 -6.5))))delim";
+    EXPECT_EQ(s, test);
+}
+