Add OverlayNG support for simple GeometryCollection inputs (#716)

include/geos/operation/overlayng/OverlayNG.h

@@ -57,9 +57,16 @@ namespace overlayng { // geos.operation.overlayng
  * * DIFFERENCE - all points which lie in the first geometry but not the second
  * * SYMDIFFERENCE - all points which lie in one geometry but not both
  *
- * Input geometries may have different dimension.
- * Input collections must be homogeneous
- * (all elements must have the same dimension).
+ * The requirements for overlay input are:
+ * * Input collections must be homogeneous
+ *   (all elements must have the same dimension).
+ * * Inputs may be simple link GeometryCollections.
+ *   A GeometryCollection is simple if it can be flattened into a valid Multi-geometry;
+ *   i.e. it is homogeneous and does not contain any overlapping Polygons.
+ * * In general, inputs must be valid geometries.
+ *   However, polygonal inputs may contain the following two kinds of "mild" invalid topology:
+ *   (i) rings which self-touch at discrete points (sometimes called inverted shells and exverted holes).
+ *   (ii) rings which touch along line segments (i.e. topology collapse).
  *
  * The precision model used for the computation can be supplied
  * independent of the precision model of the input geometry.
@@ -82,7 +89,7 @@ namespace overlayng { // geos.operation.overlayng
  * since the intersection clipping optimization can
  * interact with the snapping to alter the result.
  *
- * TOptionally the overlay computation can process using strict mode
+ * Optionally the overlay computation can process using strict mode
  * (via setStrictMode(boolean). In strict mode result semantics are:
  *
  *  - Lines and Points resulting from topology collapses are not included
@@ -101,7 +108,6 @@ namespace overlayng { // geos.operation.overlayng
  *
  * The original JTS overlay semantics correspond to non-strict mode.
  *
- *
  * If a robustness error occurs, a TopologyException is thrown.
  * These are usually caused by numerical rounding causing the noding
  * output to not be fully noded.
@@ -403,4 +409,3 @@ class GEOS_DLL OverlayNG {
 } // namespace geos.operation.overlayng
 } // namespace geos.operation
 } // namespace geos
-

include/geos/operation/overlayng/OverlayPoints.h

@@ -82,10 +82,6 @@ class GEOS_DLL OverlayPoints {
 
     std::map<Coordinate, std::unique_ptr<Point>> buildPointMap(const Geometry* geom);
 
-    Coordinate roundCoord(const Point* pt, const PrecisionModel* pm) const;
-
-
-
 public:
 
     /**
@@ -120,4 +116,3 @@ class GEOS_DLL OverlayPoints {
 } // namespace geos.operation.overlayng
 } // namespace geos.operation
 } // namespace geos
-

src/algorithm/locate/IndexedPointInAreaLocator.cpp

@@ -52,11 +52,19 @@ IndexedPointInAreaLocator::IntervalIndexedGeometry::init(const geom::Geometry& g
     // pre-compute size of segment vector
     std::size_t nsegs = 0;
     for(const geom::LineString* line : lines) {
+        //-- only include rings of Polygons or LinearRings
+        if (! line->isClosed())
+          continue;
+
         nsegs += line->getCoordinatesRO()->size() - 1;
     }
     index = decltype(index)(10, nsegs);
 
     for(const geom::LineString* line : lines) {
+        //-- only include rings of Polygons or LinearRings
+        if (! line->isClosed())
+          continue;
+
         addLine(line->getCoordinatesRO());
     }
 }
@@ -90,12 +98,6 @@ IndexedPointInAreaLocator::buildIndex(const geom::Geometry& g)
 IndexedPointInAreaLocator::IndexedPointInAreaLocator(const geom::Geometry& g)
     :	areaGeom(g)
 {
-    const std::type_info& areaGeomId = typeid(areaGeom);
-    if(areaGeomId != typeid(geom::Polygon)
-            &&	areaGeomId != typeid(geom::MultiPolygon)
-            &&	areaGeomId != typeid(geom::LinearRing)) {
-        throw util::IllegalArgumentException("Argument must be Polygonal or LinearRing");
-    }
 }
 
 geom::Location

src/operation/overlayng/OverlayMixedPoints.cpp

@@ -29,15 +29,48 @@
 #include <geos/operation/overlayng/OverlayUtil.h>
 #include <geos/util/Assert.h>
 
-
-
 namespace geos {      // geos
 namespace operation { // geos.operation
 namespace overlayng { // geos.operation.overlayng
 
-
 using namespace geos::geom;
 
+/**
+ * @brief Extracts and rounds coordinates from a geometry
+ *
+ */
+class CoordinateExtractingFilter: public geom::CoordinateFilter {
+public:
+    CoordinateExtractingFilter(CoordinateArraySequence& p_pts, const PrecisionModel& p_pm)
+        : pts(p_pts), pm(p_pm)
+    {}
+
+    /**
+     * Destructor.
+     * Virtual dctor promises appropriate behaviour when someone will
+     * delete a derived-class object via a base-class pointer.
+     * http://www.parashift.com/c++-faq-lite/virtual-functions.html#faq-20.7
+     */
+    ~CoordinateExtractingFilter() override {}
+
+    /**
+     * Performs a filtering operation with or on coord in "read-only" mode.
+     * @param coord The "read-only" Coordinate to which
+     * 				the filter is applied.
+     */
+    void
+    filter_ro(const geom::Coordinate* coord) override
+    {
+        Coordinate p(*coord);
+        pm.makePrecise(p);
+        pts.add(p);
+    }
+
+private:
+    CoordinateArraySequence& pts;
+    const PrecisionModel& pm;
+};
+
 /*public*/
 OverlayMixedPoints::OverlayMixedPoints(int p_opCode, const Geometry* geom0, const Geometry* geom1, const PrecisionModel* p_pm)
     : opCode(p_opCode)
@@ -234,16 +267,9 @@ std::unique_ptr<CoordinateArraySequence>
 OverlayMixedPoints::extractCoordinates(const Geometry* points, const PrecisionModel* p_pm) const
 {
     std::unique_ptr<CoordinateArraySequence> coords(new CoordinateArraySequence());
-    std::size_t n = points->getNumGeometries();
-    for (std::size_t i = 0; i < n; i++) {
-        const Point* point = static_cast<const Point*>(points->getGeometryN(i));
-        if (point->isEmpty()) {
-            continue;
-        }
-        Coordinate coord;
-        OverlayUtil::round(point, p_pm, coord);
-        coords->add(coord, true);
-    }
+
+    CoordinateExtractingFilter filter(*coords, *p_pm);
+    points->apply_ro(&filter);
     return coords;
 }
 
@@ -275,12 +301,6 @@ OverlayMixedPoints::extractLines(const Geometry* geom) const
     return list;
 }
 
-
-
-
-
-
-
 } // namespace geos.operation.overlayng
 } // namespace geos.operation
 } // namespace geos

src/operation/overlayng/OverlayPoints.cpp

@@ -27,6 +27,47 @@ namespace geos {      // geos
 namespace operation { // geos.operation
 namespace overlayng { // geos.operation.overlayng
 
+struct PointExtractingFilter: public GeometryComponentFilter {
+
+    PointExtractingFilter(std::map<Coordinate, std::unique_ptr<Point>>& p_ptMap, const PrecisionModel* p_pm)
+        : ptMap(p_ptMap), pm(p_pm)
+    {}
+
+    void
+    filter_ro(const Geometry* geom)
+    {
+        if (geom->getGeometryTypeId() != GEOS_POINT) return;
+
+        const Point* pt = static_cast<const Point*>(geom);
+        // don't add empty points
+        if (pt->isEmpty()) return;
+
+        Coordinate p = roundCoord(pt, pm);
+        /**
+        * Only add first occurrence of a point.
+        * This provides the merging semantics of overlay
+        */
+        if (ptMap.find(p) == ptMap.end()) {
+            std::unique_ptr<Point> newPt(pt->getFactory()->createPoint(p));
+            ptMap[p] = std::move(newPt);
+        }
+    }
+
+    static Coordinate
+    roundCoord(const Point* pt, const PrecisionModel* p_pm)
+    {
+        const Coordinate* p = pt->getCoordinate();
+        if (OverlayUtil::isFloating(p_pm))
+            return *p;
+        Coordinate p2 = *p;
+        p_pm->makePrecise(p2);
+        return p2;
+    }
+
+private:
+    std::map<Coordinate, std::unique_ptr<Point>>& ptMap;
+    const PrecisionModel* pm;
+};
 
 /*public static*/
 std::unique_ptr<Geometry>
@@ -126,54 +167,16 @@ OverlayPoints::computeUnion(std::map<Coordinate, std::unique_ptr<Point>>& map0,
     }
 }
 
-
 /*private*/
 std::map<Coordinate, std::unique_ptr<Point>>
 OverlayPoints::buildPointMap(const Geometry* geom)
 {
     std::map<Coordinate, std::unique_ptr<Point>> map;
-    for (std::size_t i = 0; i < geom->getNumGeometries(); i++) {
-        const Geometry* elt = geom->getGeometryN(i);
-        if (elt->getGeometryTypeId() != GEOS_POINT) {
-            throw util::IllegalArgumentException("Non-point geometry input to point overlay");
-        }
-        // don't add empty points
-        if (elt->isEmpty()) continue;
-
-        const Point* pt = static_cast<const Point*>(elt);
-        Coordinate p = roundCoord(pt, pm);
-        /**
-        * Only add first occurrence of a point.
-        * This provides the merging semantics of overlay
-        */
-        if (map.find(p) == map.end()) {
-            std::unique_ptr<Point> newPt(pt->getFactory()->createPoint(p));
-            map[p] = std::move(newPt);
-        }
-    }
+    PointExtractingFilter filter(map, pm);
+    geom->apply_ro(&filter);
     return map;
 }
 
-
-/*private*/
-Coordinate
-OverlayPoints::roundCoord(const Point* pt, const PrecisionModel* p_pm) const
-{
-    const Coordinate* p = pt->getCoordinate();
-    if (OverlayUtil::isFloating(p_pm))
-        return *p;
-    Coordinate p2 = *p;
-    p_pm->makePrecise(p2);
-    return p2;
-}
-
-
-
-
-
-
-
-
 } // namespace geos.operation.overlayng
 } // namespace geos.operation
 } // namespace geos