feat: curve A/P PIP for collections (rebase of #1322)

include/geos/algorithm/locate/SimplePointInAreaLocator.h

@@ -87,6 +87,17 @@ class GEOS_DLL SimplePointInAreaLocator : public PointOnGeometryLocator {
     static bool isContained(const geom::CoordinateXY& p,
                             const geom::Geometry* geom);
 
+    static bool isAnyPointContained(const geom::Geometry& pt, const geom::Geometry& areaGeom);
+
+    static bool isEveryPointContained(const geom::Geometry& pt, const geom::Geometry& areaGeom);
+
+    /**
+     * DE-9IM Contains for areal vs point geometry:
+     * non-empty, no exterior points, and at least one interior point.
+     * Boundary-only point sets return false (unlike Covers).
+     */
+    static bool isAreaContainingPoints(const geom::Geometry& pt, const geom::Geometry& areaGeom);
+
     SimplePointInAreaLocator(const geom::Geometry* p_g)
         : g(p_g)
     { }

src/algorithm/locate/SimplePointInAreaLocator.cpp

@@ -46,6 +46,98 @@ SimplePointInAreaLocator::isContained(const CoordinateXY& p, const Geometry* geo
     return Location::EXTERIOR != locate(p, geom);
 }
 
+bool
+SimplePointInAreaLocator::isAnyPointContained(const geom::Geometry& pt, const geom::Geometry& area)
+{
+    // Empty point sets do not intersect any area.
+    if (pt.isEmpty()) {
+        return false;
+    }
+    if (pt.getNumGeometries() > 1) {
+        for (size_t i = 0; i < pt.getNumGeometries(); i++ ) {
+            if (isAnyPointContained(*pt.getGeometryN(i), area)) {
+                return true;
+            }
+        }
+        return false;
+    }
+
+    const CoordinateXY* c = pt.getCoordinate();
+    if (c == nullptr) {
+        return false;
+    }
+    return isContained(*c, &area);
+}
+
+bool
+SimplePointInAreaLocator::isEveryPointContained(const geom::Geometry &pt, const geom::Geometry &area)
+{
+    // Vacuous truth: every point of an empty set is contained.
+    if (pt.isEmpty()) {
+        return true;
+    }
+    if (pt.getNumGeometries() > 1) {
+        for (size_t i = 0; i < pt.getNumGeometries(); i++ ) {
+            if (!isEveryPointContained(*pt.getGeometryN(i), area)) {
+                return false;
+            }
+        }
+        return true;
+    }
+
+    const CoordinateXY* c = pt.getCoordinate();
+    if (c == nullptr) {
+        return true;
+    }
+    return isContained(*c, &area);
+}
+
+bool
+SimplePointInAreaLocator::isAreaContainingPoints(const geom::Geometry& pt, const geom::Geometry& area)
+{
+    // Contains requires II=T; empty point geometry cannot satisfy that.
+    if (pt.isEmpty()) {
+        return false;
+    }
+
+    bool anyInterior = false;
+    const std::size_t n = pt.getNumGeometries();
+
+    if (n > 1) {
+        for (std::size_t i = 0; i < n; i++) {
+            const geom::Geometry* gi = pt.getGeometryN(i);
+            if (gi->isEmpty()) {
+                continue;
+            }
+            if (gi->getNumGeometries() > 1) {
+                if (!isAreaContainingPoints(*gi, area)) {
+                    return false;
+                }
+                anyInterior = true;
+                continue;
+            }
+            const CoordinateXY* c = gi->getCoordinate();
+            if (c == nullptr) {
+                continue;
+            }
+            const Location loc = locate(*c, &area);
+            if (loc == Location::EXTERIOR) {
+                return false;
+            }
+            if (loc == Location::INTERIOR) {
+                anyInterior = true;
+            }
+        }
+        return anyInterior;
+    }
+
+    const CoordinateXY* c = pt.getCoordinate();
+    if (c == nullptr) {
+        return false;
+    }
+    return locate(*c, &area) == Location::INTERIOR;
+}
+
 geom::Location
 SimplePointInAreaLocator::locateInGeometry(const CoordinateXY& p, const Geometry* geom)
 {

src/geom/Geometry.cpp

@@ -261,6 +261,10 @@ Geometry::getEnvelope() const
 bool
 Geometry::disjoint(const Geometry* g) const
 {
+    if (hasCurvedComponents() || g->hasCurvedComponents()) {
+        return !intersects(g);
+    }
+
 #if USE_RELATENG
     return operation::relateng::RelateNG::disjoint(this, g);
 #else
@@ -323,13 +327,13 @@ Geometry::intersects(const Geometry* g) const
         return predicate::RectangleIntersects::intersects(*p, *this);
     }
 
-    auto typ = getGeometryTypeId();
-    if (typ == GEOS_CURVEPOLYGON && g->getGeometryTypeId() == GEOS_POINT) {
-        auto loc = locate::SimplePointInAreaLocator::locatePointInSurface(*g->getCoordinate(), *detail::down_cast<const Surface*>(this));
-        return loc != Location::EXTERIOR;
-    } else if (typ == GEOS_POINT && g->getGeometryTypeId() == GEOS_CURVEPOLYGON) {
-        auto loc = locate::SimplePointInAreaLocator::locatePointInSurface(*getCoordinate(), *detail::down_cast<const Surface*>(g));
-        return loc != Location::EXTERIOR;
+    if (hasCurvedComponents() || g->hasCurvedComponents()) {
+        if (getDimension() == Dimension::A && g->getDimension() == Dimension::P) {
+            return locate::SimplePointInAreaLocator::isAnyPointContained(*g, *this);
+        }
+        if (g->getDimension() == Dimension::A && getDimension() == Dimension::P) {
+            return locate::SimplePointInAreaLocator::isAnyPointContained(*this, *g);
+        }
     }
 
 #if USE_RELATENG
@@ -426,6 +430,15 @@ Geometry::within(const Geometry* g) const
 bool
 Geometry::contains(const Geometry* g) const
 {
+    if (hasCurvedComponents() || g->hasCurvedComponents()) {
+        if (getDimension() == Dimension::A && g->getDimension() == Dimension::P) {
+            return locate::SimplePointInAreaLocator::isAreaContainingPoints(*g, *this);
+        }
+        if (g->getDimension() == Dimension::A && getDimension() == Dimension::P) {
+            return locate::SimplePointInAreaLocator::isAreaContainingPoints(*this, *g);
+        }
+    }
+
 #if USE_RELATENG
     return operation::relateng::RelateNG::contains(this, g);
 #else

tests/unit/capi/GEOSContainsTest.cpp

@@ -207,5 +207,102 @@ void object::test<6>()
     ensure_equals(GEOSContains_r(ctxt_, geom2_, geom1_), 0);
 }
 
+template<>
+template<>
+void object::test<7>()
+{
+    set_test_name("Single-point multipoint contained by MultiSurface");
+
+    geom1_ = fromWKT("MULTISURFACE(POLYGON ((100 100, 200 100, 200 200, 100 100)), CURVEPOLYGON (COMPOUNDCURVE(CIRCULARSTRING(0 0, 1 1, 2 0), (2 0, 0 0))))");
+    geom2_ = fromWKT("MULTIPOINT ((0.1556955 0.5355459))");
+
+    ensure_equals(GEOSContains(geom1_, geom2_), 1);
+}
+
+template<>
+template<>
+void object::test<8>()
+{
+    set_test_name("Only 1 part of 2-point MultiPoint contained by MultiSurface");
+
+    geom1_ = fromWKT("MULTISURFACE(POLYGON ((100 100, 200 100, 200 200, 100 100)), CURVEPOLYGON (COMPOUNDCURVE(CIRCULARSTRING(0 0, 1 1, 2 0), (2 0, 0 0))))");
+    geom2_ = fromWKT("MULTIPOINT ((0.1556955 0.5355459), (500 500))");
+
+    ensure_equals(GEOSContains(geom1_, geom2_), 0);
+}
+
+template<>
+template<>
+void object::test<9>()
+{
+    set_test_name("MultiPoint contained by MultiSurface");
+
+    geom1_ = fromWKT("MULTISURFACE(POLYGON ((100 100, 200 100, 200 200, 100 100)), CURVEPOLYGON (COMPOUNDCURVE(CIRCULARSTRING(0 0, 1 1, 2 0), (2 0, 0 0))))");
+    geom2_ = fromWKT("MULTIPOINT ((0.1556955 0.5355459), (199 101))");
+
+    ensure_equals(GEOSContains(geom1_, geom2_), 1);
+}
+
+
+template<>
+template<>
+void object::test<10>()
+{
+    set_test_name("Empty MultiPoint is not contained by MultiSurface (II requires T)");
+
+    geom1_ = fromWKT("MULTISURFACE(CURVEPOLYGON (COMPOUNDCURVE(CIRCULARSTRING(0 0, 1 1, 2 0), (2 0, 0 0))))");
+    geom2_ = fromWKT("MULTIPOINT EMPTY");
+
+    ensure_equals(GEOSContains(geom1_, geom2_), 0);
+}
+
+template<>
+template<>
+void object::test<11>()
+{
+    set_test_name("Empty Point is not contained by CurvePolygon (II requires T)");
+
+    geom1_ = fromWKT("CURVEPOLYGON (COMPOUNDCURVE(CIRCULARSTRING(0 0, 1 1, 2 0), (2 0, 0 0)))");
+    geom2_ = fromWKT("POINT EMPTY");
+
+    ensure_equals(GEOSContains(geom1_, geom2_), 0);
+}
+
+template<>
+template<>
+void object::test<12>()
+{
+    set_test_name("Boundary point is not contained by CurvePolygon");
+
+    geom1_ = fromWKT("CURVEPOLYGON (COMPOUNDCURVE(CIRCULARSTRING(0 0, 1 1, 2 0), (2 0, 0 0)))");
+    geom2_ = fromWKT("POINT (0 0)");
+
+    ensure_equals(GEOSContains(geom1_, geom2_), 0);
+}
+
+template<>
+template<>
+void object::test<13>()
+{
+    set_test_name("Boundary + interior MultiPoint is contained");
+
+    geom1_ = fromWKT("CURVEPOLYGON (COMPOUNDCURVE(CIRCULARSTRING(0 0, 1 1, 2 0), (2 0, 0 0)))");
+    geom2_ = fromWKT("MULTIPOINT ((0 0), (0.1556955 0.5355459))");
+
+    ensure_equals(GEOSContains(geom1_, geom2_), 1);
+}
+
+template<>
+template<>
+void object::test<14>()
+{
+    set_test_name("Boundary-only MultiPoint is not contained");
+
+    geom1_ = fromWKT("CURVEPOLYGON (COMPOUNDCURVE(CIRCULARSTRING(0 0, 1 1, 2 0), (2 0, 0 0)))");
+    geom2_ = fromWKT("MULTIPOINT ((0 0), (2 0))");
+
+    ensure_equals(GEOSContains(geom1_, geom2_), 0);
+}
+
 } // namespace tut
 

tests/unit/capi/GEOSDisjointTest.cpp

@@ -58,5 +58,32 @@ void object::test<3>()
     ensure_equals(0, GEOSDisjoint_r(ctxt_, geom1_, geom2_));
 }
 
+template<>
+template<>
+void object::test<4>()
+{
+    set_test_name("MultiSurface / MultiPoint PIP");
+
+    geom1_ = fromWKT("MULTISURFACE(POLYGON ((100 100, 200 100, 200 200, 100 100)), CURVEPOLYGON (COMPOUNDCURVE(CIRCULARSTRING(0 0, 1 1, 2 0), (2 0, 0 0))))");
+    geom2_ = fromWKT("MULTIPOINT ((5000 5000), (0.1556955 0.5355459))");
+
+    ensure_equals(GEOSDisjoint(geom1_, geom2_), 0);
+    ensure_equals(GEOSDisjoint(geom2_, geom1_), 0);
+}
+
+
+template<>
+template<>
+void object::test<5>()
+{
+    set_test_name("Empty MultiPoint is disjoint from MultiSurface");
+
+    geom1_ = fromWKT("MULTISURFACE(CURVEPOLYGON (COMPOUNDCURVE(CIRCULARSTRING(0 0, 1 1, 2 0), (2 0, 0 0))))");
+    geom2_ = fromWKT("MULTIPOINT EMPTY");
+
+    ensure_equals(GEOSDisjoint(geom1_, geom2_), 1);
+    ensure_equals(GEOSDisjoint(geom2_, geom1_), 1);
+}
+
 } // namespace tut
 

tests/unit/capi/GEOSIntersectsTest.cpp

@@ -236,8 +236,8 @@ template<>
 template<>
 void object::test<12>()
 {
-    geom1_ = fromWKT("CURVEPOLYGON (COMPOUNDCURVE(CIRCULARSTRING(0 0, 1 1, 2 0), (2 0, 0 0)))");
-    geom2_ = fromWKT("POINT (0.1556955 0.5355459)");
+    geom1_ = fromWKT("MULTISURFACE(POLYGON ((100 100, 200 100, 200 200, 100 100)), CURVEPOLYGON (COMPOUNDCURVE(CIRCULARSTRING(0 0, 1 1, 2 0), (2 0, 0 0))))");
+    geom2_ = fromWKT("MULTIPOINT ((5000 5000), (0.1556955 0.5355459))");
 
     // PostGIS would return false here because geom2 is inside geom1
     // but outside the linearized form of geom1
@@ -272,5 +272,47 @@ void object::test<14>()
     ensure_equals(1, GEOSIntersects_r(ctxt_, geom1_, geom2_));
 }
 
+template<>
+template<>
+void object::test<15>()
+{
+    set_test_name("PostGIS ticket 5832");
+
+    geom1_ = fromWKT("MULTISURFACE(CURVEPOLYGON(COMPOUNDCURVE((25492739.7449 6677441.1816,25492771.1213 6677416.7832),CIRCULARSTRING(25492771.1213 6677416.7832,25492832.3384 6677400.8583,25492885.9851 6677434.3719),(25492885.9851 6677434.3719,25492900.7986 6677455.7498),CIRCULARSTRING(25492900.7986 6677455.7498,25492901.1068 6677457.3233,25492900.1601 6677458.6175,25492822.2626 6677477.2365,25492747.0232 6677449.7828),(25492747.0232 6677449.7828,25492739.7748 6677444.3615),CIRCULARSTRING(25492739.7748 6677444.3615,25492738.9731 6677442.7789,25492739.7449 6677441.1816))))");
+    geom2_ = fromWKT("POINT (25492818 6677399.98)");
+
+    ensure_equals(GEOSIntersects(geom1_, geom2_), 1);
+    ensure_equals(GEOSIntersects(geom2_, geom1_), 1);
+}
+
+
+
+template<>
+template<>
+void object::test<16>()
+{
+    set_test_name("Empty MultiPoint does not intersect MultiSurface");
+
+    geom1_ = fromWKT("MULTISURFACE(CURVEPOLYGON (COMPOUNDCURVE(CIRCULARSTRING(0 0, 1 1, 2 0), (2 0, 0 0))))");
+    geom2_ = fromWKT("MULTIPOINT EMPTY");
+
+    ensure_equals(GEOSIntersects(geom1_, geom2_), 0);
+    ensure_equals(GEOSIntersects(geom2_, geom1_), 0);
+}
+
+
+template<>
+template<>
+void object::test<17>()
+{
+    set_test_name("Boundary point intersects CurvePolygon");
+
+    geom1_ = fromWKT("CURVEPOLYGON (COMPOUNDCURVE(CIRCULARSTRING(0 0, 1 1, 2 0), (2 0, 0 0)))");
+    geom2_ = fromWKT("POINT (0 0)");
+
+    ensure_equals(GEOSIntersects(geom1_, geom2_), 1);
+    ensure_equals(GEOSIntersects(geom2_, geom1_), 1);
+}
+
 } // namespace tut