Fix overlay for GCs with empty

include/geos/geom/HeuristicOverlay.h

@@ -63,6 +63,8 @@ class StructuredCollection {
         , poly_union(nullptr)
     {};
 
+    static std::unique_ptr<Geometry> overlay(const Geometry* g0, const Geometry* g1, int opCode);
+
     void readCollection(const Geometry* g);
     const Geometry* getPolyUnion()  const { return poly_union.get(); }
     const Geometry* getLineUnion()  const { return line_union.get(); }

src/geom/Geometry.cpp

@@ -619,43 +619,6 @@ Geometry::intersection(const Geometry* other) const
 std::unique_ptr<Geometry>
 Geometry::Union(const Geometry* other) const
 {
-#ifdef SHORTCIRCUIT_PREDICATES
-    // if envelopes are disjoint return a MULTI geom or
-    // a geometrycollection
-    if(! getEnvelopeInternal()->intersects(other->getEnvelopeInternal())) {
-//cerr<<"SHORTCIRCUITED-UNION engaged"<<endl;
-        const GeometryCollection* coll;
-
-        std::size_t ngeomsThis = getNumGeometries();
-        std::size_t ngeomsOther = other->getNumGeometries();
-
-        // Allocated for ownership transfer
-        std::vector<std::unique_ptr<Geometry>> v;
-        v.reserve(ngeomsThis + ngeomsOther);
-
-
-        if(nullptr != (coll = dynamic_cast<const GeometryCollection*>(this))) {
-            for(std::size_t i = 0; i < ngeomsThis; ++i) {
-                v.push_back(coll->getGeometryN(i)->clone());
-            }
-        }
-        else {
-            v.push_back(this->clone());
-        }
-
-        if(nullptr != (coll = dynamic_cast<const GeometryCollection*>(other))) {
-            for(std::size_t i = 0; i < ngeomsOther; ++i) {
-                v.push_back(coll->getGeometryN(i)->clone());
-            }
-        }
-        else {
-            v.push_back(other->clone());
-        }
-
-        return _factory->buildGeometry(std::move(v));
-    }
-#endif
-
     return HeuristicOverlay(this, other, OverlayNG::UNION);
 }
 
@@ -676,40 +639,6 @@ Geometry::difference(const Geometry* other) const
 std::unique_ptr<Geometry>
 Geometry::symDifference(const Geometry* other) const
 {
-    // if envelopes are disjoint return a MULTI geom or
-    // a geometrycollection
-    if(! getEnvelopeInternal()->intersects(other->getEnvelopeInternal()) && !(isEmpty() && other->isEmpty())) {
-        const GeometryCollection* coll;
-
-        std::size_t ngeomsThis = getNumGeometries();
-        std::size_t ngeomsOther = other->getNumGeometries();
-
-        // Allocated for ownership transfer
-        std::vector<std::unique_ptr<Geometry>> v;
-        v.reserve(ngeomsThis + ngeomsOther);
-
-
-        if(nullptr != (coll = dynamic_cast<const GeometryCollection*>(this))) {
-            for(std::size_t i = 0; i < ngeomsThis; ++i) {
-                v.push_back(coll->getGeometryN(i)->clone());
-            }
-        }
-        else {
-            v.push_back(this->clone());
-        }
-
-        if(nullptr != (coll = dynamic_cast<const GeometryCollection*>(other))) {
-            for(std::size_t i = 0; i < ngeomsOther; ++i) {
-                v.push_back(coll->getGeometryN(i)->clone());
-            }
-        }
-        else {
-            v.push_back(other->clone());
-        }
-
-        return _factory->buildGeometry(std::move(v));
-    }
-
     return HeuristicOverlay(this, other, OverlayNG::SYMDIFFERENCE);
 }
 

src/geom/HeuristicOverlay.cpp

@@ -40,31 +40,63 @@ namespace geom { // geos::geom
 using operation::overlayng::OverlayNG;
 using operation::overlayng::OverlayNGRobust;
 
+bool isDisjointNonEmpty(const Geometry* g0, const Geometry* g1)
+{
+    if (g0->isEmpty() && g1->isEmpty())
+        return false;
+    return ! g0->getEnvelopeInternal()->intersects(g1->getEnvelopeInternal());
+}
+
+void 
+extractELements(const Geometry* g, std::vector<std::unique_ptr<Geometry>>& v)
+{
+    const GeometryCollection* coll;
+    if (nullptr != (coll = dynamic_cast<const GeometryCollection*>(g))) {
+        //-- buffer num geoms, since can be expensive for nested collections
+        std::size_t numGeoms = g->getNumGeometries();
+        for(std::size_t i = 0; i < numGeoms; ++i) {
+            //-- recurse to handle nested GCs
+            extractELements(coll->getGeometryN(i), v);
+        }
+    }
+    else if (g->isEmpty()) {
+        return;
+    }
+    else {
+        v.push_back(g->clone());
+    }
+}
+
+std::unique_ptr<Geometry>
+reduceCombine(const Geometry* g0, const Geometry* g1)
+{
+    // Allocated for ownership transfer
+    std::vector<std::unique_ptr<Geometry>> v;
+    v.reserve(g0->getNumGeometries() + g1->getNumGeometries());
+    extractELements(g0, v);
+    extractELements(g1, v);
+    return g0->getFactory()->buildGeometry(std::move(v));
+}
+
 std::unique_ptr<Geometry>
 HeuristicOverlay(const Geometry* g0, const Geometry* g1, int opCode)
 {
-    std::unique_ptr<Geometry> ret;
+    //-- if both are empty use OverlayNG return-type logic
+    if ( ( opCode == OverlayNG::UNION 
+        || opCode == OverlayNG::SYMDIFFERENCE)
+            && isDisjointNonEmpty(g0, g1)) {
+        return reduceCombine(g0, g1);
+    }
 
     /*
     * overlayng::OverlayNGRobust does not currently handle
     * GeometryCollection (collections of mixed dimension)
     * so we handle that case here.
     */
-    if ((g0->isMixedDimension() && !g0->isEmpty()) ||
-        (g1->isMixedDimension() && !g1->isEmpty()))
+    if ((g0->isMixedDimension() && ! g0->isEmpty()) ||
+        (g1->isMixedDimension() && ! g1->isEmpty()))
     {
-        StructuredCollection s0(g0);
-        StructuredCollection s1(g1);
-        switch (opCode) {
-        case OverlayNG::UNION:
-            return s0.doUnion(s1);
-        case OverlayNG::DIFFERENCE:
-            return s0.doDifference(s1);
-        case OverlayNG::SYMDIFFERENCE:
-            return s0.doSymDifference(s1);
-        case OverlayNG::INTERSECTION:
-            return s0.doIntersection(s1);
-        }
+        return StructuredCollection::overlay(g0, g1, opCode);
     }
 
     /*
@@ -86,6 +118,7 @@ HeuristicOverlay(const Geometry* g0, const Geometry* g1, int opCode)
     * Running overlayng::OverlayNGRobust at this stage should guarantee
     * that none of the other heuristics are ever needed.
     */
+    std::unique_ptr<Geometry> ret;
     if (g0 == nullptr && g1 == nullptr) {
         return std::unique_ptr<Geometry>(nullptr);
     }
@@ -108,6 +141,24 @@ HeuristicOverlay(const Geometry* g0, const Geometry* g1, int opCode)
     return ret;
 }
 
+/* public static */
+std::unique_ptr<Geometry>
+StructuredCollection::overlay(const Geometry* g0, const Geometry* g1, int opCode)
+{
+    StructuredCollection s0(g0);
+    StructuredCollection s1(g1);
+    switch (opCode) {
+    case OverlayNG::UNION:
+        return s0.doUnion(s1);
+    case OverlayNG::DIFFERENCE:
+        return s0.doDifference(s1);
+    case OverlayNG::SYMDIFFERENCE:
+        return s0.doSymDifference(s1);
+    default: // only OverlayNG::INTERSECTION
+        return s0.doIntersection(s1);
+    }
+}
+
 /* public */
 void
 StructuredCollection::readCollection(const Geometry* g)

src/operation/union/CascadedPolygonUnion.cpp

@@ -20,11 +20,11 @@
 
 #include <geos/geom/Geometry.h>
 #include <geos/geom/GeometryFactory.h>
-#include <geos/geom/HeuristicOverlay.h>
 #include <geos/geom/MultiPolygon.h>
 #include <geos/geom/Polygon.h>
 #include <geos/index/strtree/TemplateSTRtree.h>
 #include <geos/operation/overlayng/OverlayNG.h>
+#include <geos/operation/overlayng/OverlayNGRobust.h>
 #include <geos/operation/union/CascadedPolygonUnion.h>
 #include <geos/operation/valid/IsValidOp.h>
 #include <geos/operation/valid/IsSimpleOp.h>
@@ -201,10 +201,8 @@ CascadedPolygonUnion::restrictToPolygons(std::unique_ptr<geom::Geometry> g)
 std::unique_ptr<geom::Geometry>
 ClassicUnionStrategy::Union(const geom::Geometry* g0, const geom::Geometry* g1)
 {
-    // TODO make an rvalue overload for this that can consume its inputs.
-    // At a minimum, a copy in the buffer fallback can be eliminated.
     try {
-        return geom::HeuristicOverlay(g0, g1, operation::overlayng::OverlayNG::UNION);
+        return operation::overlayng::OverlayNGRobust::Overlay(g0, g1, operation::overlayng::OverlayNG::UNION);
     }
     catch (const util::TopologyException &ex) {
         ::geos::ignore_unused_variable_warning(ex);

src/operation/valid/MakeValid.cpp

@@ -23,9 +23,9 @@
 #include <geos/operation/valid/IsValidOp.h>
 
 #include <geos/operation/overlayng/OverlayNG.h>
+#include <geos/operation/overlayng/OverlayNGRobust.h>
 #include <geos/operation/polygonize/BuildArea.h>
 #include <geos/operation/union/UnaryUnionOp.h>
-#include <geos/geom/HeuristicOverlay.h>
 #include <geos/geom/Geometry.h>
 #include <geos/geom/GeometryCollection.h>
 #include <geos/geom/GeometryFactory.h>
@@ -51,6 +51,7 @@
 
 using namespace geos::geom;
 using geos::operation::overlayng::OverlayNG;
+using geos::operation::overlayng::OverlayNGRobust;
 
 namespace geos {
 namespace operation { // geos.operation
@@ -60,19 +61,19 @@ namespace valid { // geos.operation.valid
 static std::unique_ptr<geom::Geometry>
 makeValidSymDifference(const geom::Geometry* g0, const geom::Geometry* g1)
 {
-    return HeuristicOverlay(g0, g1, OverlayNG::SYMDIFFERENCE);
+    return OverlayNGRobust::Overlay(g0, g1, OverlayNG::SYMDIFFERENCE);
 }
 
 static std::unique_ptr<geom::Geometry>
 makeValidDifference(const geom::Geometry* g0, const geom::Geometry* g1)
 {
-    return HeuristicOverlay(g0, g1, OverlayNG::DIFFERENCE);
+    return OverlayNGRobust::Overlay(g0, g1, OverlayNG::DIFFERENCE);
 }
 
 static std::unique_ptr<geom::Geometry>
 makeValidUnion(const geom::Geometry* g0, const geom::Geometry* g1)
 {
-    return HeuristicOverlay(g0, g1, OverlayNG::UNION);
+    return OverlayNGRobust::Overlay(g0, g1, OverlayNG::UNION);
 }
 
 /*

tests/xmltester/BufferResultMatcher.cpp

@@ -19,7 +19,6 @@
 #include "BufferResultMatcher.h"
 
 #include <geos/geom/Geometry.h>
-#include <geos/geom/HeuristicOverlay.h>
 #include <geos/operation/overlayng/OverlayNG.h>
 #include <geos/algorithm/distance/DiscreteHausdorffDistance.h>
 
@@ -65,15 +64,8 @@ BufferResultMatcher::isSymDiffAreaInTolerance(
     const geom::Geometry& actualBuffer,
     const geom::Geometry& expectedBuffer)
 {
-    typedef std::unique_ptr<geom::Geometry> GeomPtr;
-
-    using geos::geom::HeuristicOverlay;
-
     double area = expectedBuffer.getArea();
-    GeomPtr diff = HeuristicOverlay(&actualBuffer, &expectedBuffer,
-                            operation::overlayng::OverlayNG::SYMDIFFERENCE);
-
-    double areaDiff = diff->getArea();
+    double areaDiff = actualBuffer.symDifference(&expectedBuffer)->getArea();
 
     // can't get closer than difference area = 0 !
     // This also handles case when symDiff is empty