Repair GEO_INTERSECTS

CHANGELOG

@@ -1,6 +1,11 @@
 devel
 -----
 
+* Fixed various problems in GEO_INTERSECTS: wrong results, not implemented
+  cases and numerically unstable behaviour. In particular, the case of
+  the intersection of two polygons in which one is an S2LngLatRect
+  is fixed (BTS-475).
+
 * Fixed ES-867 and ES-922: removed eslint from NPM packages descriptions and
   updated netmask package to non-vulnerable version.
 

arangod/Aql/Functions.cpp

@@ -1200,8 +1200,15 @@ AqlValue geoContainsIntersect(ExpressionContext* expressionContext,
     return AqlValue(AqlValueHintNull());
   }
 
-  bool result = contains ? outer.contains(&inner) : outer.intersects(&inner);
-  return AqlValue(AqlValueHintBool(result));
+  bool result;
+  try {
+    result = contains ? outer.contains(&inner) : outer.intersects(&inner);
+    return AqlValue(AqlValueHintBool(result));
+  } catch (arangodb::basics::Exception const& ex) {
+    res.reset(ex.code(), ex.what());
+    registerWarning(expressionContext, func, res);
+    return AqlValue(AqlValueHintBool(false));
+  }
 }
 
 static Result parseGeoPolygon(VPackSlice polygon, VPackBuilder& b) {

lib/Geo/ShapeContainer.cpp

@@ -55,6 +55,24 @@
 using namespace arangodb;
 using namespace arangodb::geo;
 
+namespace {
+
+S2Polygon latLngRectToPolygon(S2LatLngRect const* rect) {
+  // Construct polygon from rect:
+  std::vector<S2Point> v;
+  v.reserve(5);
+  v.emplace_back(rect->GetVertex(0).ToPoint());
+  v.emplace_back(rect->GetVertex(1).ToPoint());
+  v.emplace_back(rect->GetVertex(2).ToPoint());
+  v.emplace_back(rect->GetVertex(3).ToPoint());
+  v.emplace_back(rect->GetVertex(0).ToPoint());
+  std::unique_ptr<S2Loop> loop;
+  loop = std::make_unique<S2Loop>(std::move(v), S2Debug::DISABLE);
+  return S2Polygon{std::move(loop), S2Debug::DISABLE};
+}
+
+}
+
 Result ShapeContainer::parseCoordinates(VPackSlice const& json, bool geoJson) {
   if (!json.isArray() || json.length() < 2) {
     return Result(TRI_ERROR_BAD_PARAMETER, "Invalid coordinate pair");
@@ -557,41 +575,57 @@ bool ShapeContainer::equals(ShapeContainer const* cc) const {
 bool ShapeContainer::intersects(S2Polyline const* other) const {
   switch (_type) {
     case ShapeContainer::Type::S2_POINT: {
-      S2Point const& p = static_cast<S2PointRegion*>(_data)->point();
-      // containment is only numerically defined on the endpoints
-      for (int k = 0; k < other->num_vertices(); k++) {
-        if (other->vertex(k) == p) {
-          return true;
-        }
-      }
-      return false;
+      THROW_ARANGO_EXCEPTION_MESSAGE(TRI_ERROR_NOT_IMPLEMENTED,
+          "The case GEO_INTERSECTS(<point>, <polyline>) is numerically instable and thus not supported.");
     }
     case ShapeContainer::Type::S2_POLYLINE: {
       S2Polyline const* ll = static_cast<S2Polyline const*>(_data);
       return ll->Intersects(other);
     }
     case ShapeContainer::Type::S2_LATLNGRECT: {
-      S2LatLngRect const* rect = static_cast<S2LatLngRect const*>(_data);
-      for (int k = 0; k < other->num_vertices(); k++) {
-        if (rect->Contains(other->vertex(k))) {
-          return true;
-        }
-      }
-      return false;
+      auto rectPoly = ::latLngRectToPolygon(static_cast<S2LatLngRect const*>(_data));
+      auto cuts = rectPoly.IntersectWithPolyline(*other);
+      return !cuts.empty();
     }
     case ShapeContainer::Type::S2_POLYGON: {
       S2Polygon const* poly = static_cast<S2Polygon const*>(_data);
       auto cuts = poly->IntersectWithPolyline(*other);
       return !cuts.empty();
     }
     case ShapeContainer::Type::S2_MULTIPOINT:
-    case ShapeContainer::Type::S2_MULTIPOLYLINE:
+    case ShapeContainer::Type::S2_MULTIPOLYLINE: {
+      THROW_ARANGO_EXCEPTION_MESSAGE(TRI_ERROR_NOT_IMPLEMENTED,
+          "The case GEO_INTERSECTS(<multipoint or multipolyline>, <polyline>) is not yet implemented.");
+    }
     case ShapeContainer::Type::EMPTY:
       TRI_ASSERT(false);
   }
   return false;
 }
 
+namespace {
+bool intersectRectPolygon(S2LatLngRect const* rect, S2Polygon const* poly) {
+  if (rect->is_full()) {
+    return true;  // rectangle spans entire sphere
+  } else if (rect->is_point()) {
+    return poly->Contains(rect->lo().ToPoint());  // easy case
+  } else if (!rect->Intersects(poly->GetRectBound())) {
+    return false;  // cheap rejection
+  }
+  auto rectPoly = ::latLngRectToPolygon(rect);
+  return poly->Intersects(&rectPoly);
+}
+
+bool insersectMultiPointsRegion(S2MultiPointRegion const* points, S2Region const* region) {
+  for (int i = 0; i < points->num_points(); ++i) {
+    if (region->Contains(points->point(i))) {
+      return true;
+    }
+  }
+  return false;
+}
+}
+
 bool ShapeContainer::intersects(S2LatLngRect const* other) const {
   switch (_type) {
     case ShapeContainer::Type::S2_POINT: {
@@ -600,7 +634,10 @@ bool ShapeContainer::intersects(S2LatLngRect const* other) const {
     }
 
     case ShapeContainer::Type::S2_POLYLINE: {
-      return contains(other);
+      auto rectPoly = ::latLngRectToPolygon(static_cast<S2LatLngRect const*>(other));
+      S2Polyline const* self = static_cast<S2Polyline const*>(_data);
+      auto cuts = rectPoly.IntersectWithPolyline(*self);
+      return !cuts.empty();
     }
 
     case ShapeContainer::Type::S2_LATLNGRECT: {
@@ -610,22 +647,17 @@ bool ShapeContainer::intersects(S2LatLngRect const* other) const {
 
     case ShapeContainer::Type::S2_POLYGON: {
       S2Polygon const* self = static_cast<S2Polygon const*>(_data);
-      if (other->is_full()) {
-        return true;  // rectangle spans entire sphere
-      } else if (other->is_point()) {
-        return self->Contains(other->lo().ToPoint());  // easy case
-      } else if (!other->Intersects(self->GetRectBound())) {
-        return false;  // cheap rejection
-      }
-      // construct bounding polyline of rect
-      S2Polyline rectBound({other->GetVertex(0), other->GetVertex(1),
-                            other->GetVertex(2), other->GetVertex(3)});
-      return self->Intersects(rectBound);
+      return intersectRectPolygon(other, self);
+    }
+
+    case ShapeContainer::Type::S2_MULTIPOINT: {
+      S2MultiPointRegion* self = static_cast<S2MultiPointRegion*>(_data);
+      return insersectMultiPointsRegion(self, other);
     }
 
-    case ShapeContainer::Type::S2_MULTIPOINT:
     case ShapeContainer::Type::S2_MULTIPOLYLINE: {
-      return contains(other);  // same
+      THROW_ARANGO_EXCEPTION_MESSAGE(TRI_ERROR_NOT_IMPLEMENTED,
+          "The case GEO_INTERSECTS(<multiline>, <latlngrect>) is not yet implemented.");
     }
 
     case ShapeContainer::Type::EMPTY:
@@ -641,31 +673,24 @@ bool ShapeContainer::intersects(S2Polygon const* other) const {
       return other->Contains(p);
     }
     case ShapeContainer::Type::S2_POLYLINE: {
-      LOG_TOPIC("2cb3c", ERR, Logger::FIXME)
-          << "intersection with polyline is not well defined";
-      return false;  // numerically not well defined
+      S2Polyline* line = static_cast<S2Polyline*>(_data);
+      auto cuts = other->IntersectWithPolyline(*line);
+      return !cuts.empty();
     }
     case ShapeContainer::Type::S2_LATLNGRECT: {
       S2LatLngRect const* self = static_cast<S2LatLngRect const*>(_data);
-      if (self->is_full()) {
-        return true;  // rectangle spans entire sphere
-      } else if (self->is_point()) {
-        return other->Contains(self->lo().ToPoint());  // easy case
-      } else if (!self->Intersects(other->GetRectBound())) {
-        return false;  // cheap rejection
-      }
-      // construct bounding polyline of rect
-      S2Polyline rectBound({self->GetVertex(0), self->GetVertex(1),
-                            self->GetVertex(2), self->GetVertex(3)});
-      return other->Intersects(rectBound);
+      return intersectRectPolygon(self, other);
     }
     case ShapeContainer::Type::S2_POLYGON: {
       S2Polygon const* self = static_cast<S2Polygon const*>(_data);
       return self->Intersects(other);
     }
-    case ShapeContainer::Type::EMPTY:
     case ShapeContainer::Type::S2_MULTIPOINT:
-    case ShapeContainer::Type::S2_MULTIPOLYLINE:
+    case ShapeContainer::Type::S2_MULTIPOLYLINE: {
+      THROW_ARANGO_EXCEPTION_MESSAGE(TRI_ERROR_NOT_IMPLEMENTED,
+          "The case GEO_INTERSECTS(<multipoint or multipolyline>, <polygon>) is not yet implemented.");
+    }
+    case ShapeContainer::Type::EMPTY:
       TRI_ASSERT(false);
   }
   return false;
@@ -674,6 +699,11 @@ bool ShapeContainer::intersects(S2Polygon const* other) const {
 bool ShapeContainer::intersects(ShapeContainer const* cc) const {
   switch (cc->_type) {
     case ShapeContainer::Type::S2_POINT: {
+      if (_type == ShapeContainer::Type::S2_POLYLINE ||
+          _type == ShapeContainer::Type::S2_MULTIPOLYLINE) {
+        THROW_ARANGO_EXCEPTION_MESSAGE(TRI_ERROR_NOT_IMPLEMENTED,
+            "The case GEO_INTERSECTS(<polyline>, <point>) is numerically instable and thus not supported.");
+      }
       S2Point const& p = static_cast<S2PointRegion*>(cc->_data)->point();
       return _data->Contains(p);  // same
     }
@@ -687,13 +717,13 @@ bool ShapeContainer::intersects(ShapeContainer const* cc) const {
       return intersects(static_cast<S2LatLngRect const*>(cc->_data));
     }
     case ShapeContainer::Type::S2_MULTIPOINT: {
-      auto pts = static_cast<S2MultiPointRegion const*>(cc->_data);
-      for (int k = 0; k < pts->num_points(); k++) {
-        if (_data->Contains(pts->point(k))) {
-          return true;
-        }
+      if (_type == ShapeContainer::Type::S2_POLYLINE ||
+          _type == ShapeContainer::Type::S2_MULTIPOLYLINE) {
+        THROW_ARANGO_EXCEPTION_MESSAGE(TRI_ERROR_NOT_IMPLEMENTED,
+            "The case GEO_INTERSECTS(<polyline>, <multipoint>) is numerically instable and thus not supported.");
       }
-      return false;
+      auto pts = static_cast<S2MultiPointRegion const*>(cc->_data);
+      return insersectMultiPointsRegion(pts, _data);
     }
     case ShapeContainer::Type::S2_MULTIPOLYLINE: {
       auto lines = static_cast<S2MultiPolyline const*>(cc->_data);