Convex hull

src/spatialjoin/GeometryCache.cpp

@@ -183,6 +183,9 @@ sj::Line sj::GeometryCache<sj::Line>::getFromDisk(size_t off,
   // OBB
   readPoly(_geomsFReads[tid], ret.obb);
 
+  // convex hull
+  readPoly(_geomsFReads[tid], ret.convexHull);
+
   return ret;
 }
 
@@ -244,6 +247,9 @@ sj::Area sj::GeometryCache<sj::Area>::getFromDisk(size_t off,
   // OBB
   readPoly(_geomsFReads[tid], ret.obb);
 
+  // convex hull
+  readPoly(_geomsFReads[tid], ret.convexHull);
+
   // simplified inner
   readPoly(_geomsFReads[tid], ret.inner);
 
@@ -414,6 +420,9 @@ size_t sj::GeometryCache<sj::Line>::writeTo(const sj::Line& val,
   // OBB
   ret += writePoly(val.obb, str);
 
+  // convex hull
+  ret += writePoly(val.convexHull, str);
+
   return ret;
 }
 
@@ -476,6 +485,9 @@ size_t sj::GeometryCache<sj::Area>::writeTo(const sj::Area& val,
   // OBB
   ret += writePoly(val.obb, str);
 
+  // convex hull
+  ret += writePoly(val.convexHull, str);
+
   // innerGeom
   ret += writePoly(val.inner, str);
 

src/spatialjoin/GeometryCache.h

@@ -50,6 +50,9 @@ struct Area {
   // cutouts
   std::map<int32_t, size_t> cutouts;
 
+  // Convex hull
+  util::geo::I32XSortedPolygon convexHull;
+
   // OBB
   util::geo::I32XSortedPolygon obb;
 
@@ -103,6 +106,9 @@ struct Line {
   // cutouts
   std::map<int32_t, size_t> cutouts;
 
+  // Convex hull
+  util::geo::I32XSortedPolygon convexHull;
+
   // OBB
   util::geo::I32XSortedPolygon obb;
 };

src/spatialjoin/SpatialJoinMain.cpp

@@ -74,6 +74,8 @@ void printHelp(int argc, char** argv) {
       << "disable cutouts\n"
       << std::setw(41) << "  --no-diag-box"
       << "disable diagonal bounding-box based pre-filter\n"
+      << std::setw(41) << "  --no-convex-hulls"
+      << "disable convex hulls\n"
       << std::setw(41) << "  --no-fast-sweep-skip"
       << "disable fast sweep skip using binary search\n"
       << std::setw(41) << "  --use-inner-outer"
@@ -118,6 +120,7 @@ int main(int argc, char** argv) {
   bool useOBB = true;
   bool useCutouts = true;
   bool useDiagBox = true;
+  bool useConvexHulls = true;
   bool useFastSweepSkip = true;
   bool useInnerOuter = false;
 
@@ -170,6 +173,8 @@ int main(int argc, char** argv) {
           useCutouts = false;
         } else if (cur == "--no-diag-box") {
           useDiagBox = false;
+        } else if (cur == "--no-convex-hulls") {
+          useConvexHulls = false;
         } else if (cur == "--no-fast-sweep-skip") {
           useFastSweepSkip = false;
         } else if (cur == "--use-inner-outer") {
@@ -245,7 +250,7 @@ int main(int argc, char** argv) {
   Sweeper sweeper(
       {numThreads, numCaches, prefix, intersects, contains, covers, touches,
        equals, overlaps, crosses, suffix, useBoxIds, useArea, useOBB,
-       useCutouts, useDiagBox, useFastSweepSkip, useInnerOuter},
+       useCutouts, useDiagBox, useConvexHulls, useFastSweepSkip, useInnerOuter},
       cache, output);
 
   LOGTO(INFO, std::cerr) << "Parsing input geometries...";

src/spatialjoin/Stats.h

@@ -23,6 +23,11 @@ struct Stats {
   uint64_t timeOBBIsectAreaPoint = 0;
   uint64_t timeOBBIsectLineLine = 0;
 
+  uint64_t timeConvexHullIsectAreaArea = 0;
+  uint64_t timeConvexHullIsectAreaLine = 0;
+  uint64_t timeConvexHullIsectAreaPoint = 0;
+  uint64_t timeConvexHullIsectLineLine = 0;
+
   uint64_t timeFullGeoCheckAreaArea = 0;
   uint64_t timeFullGeoCheckAreaLine = 0;
   uint64_t timeFullGeoCheckAreaPoint = 0;
@@ -69,7 +74,9 @@ inline std::string Stats::toString() {
              timeGeoCacheRetrievalPoint + timeWrite + timeBoxIdIsectAreaArea +
              timeBoxIdIsectAreaLine + timeOBBIsectAreaArea +
              timeOBBIsectAreaLine + timeOBBIsectAreaPoint +
-             timeOBBIsectLineLine + timeBoxIdIsectAreaPoint +
+             timeOBBIsectLineLine + timeConvexHullIsectAreaArea +
+             timeConvexHullIsectAreaLine + timeConvexHullIsectAreaPoint +
+             timeConvexHullIsectLineLine + timeBoxIdIsectAreaPoint +
              timeBoxIdIsectLineLine + timeBoxIdIsectLinePoint +
              timeFullGeoCheckAreaArea + timeFullGeoCheckAreaLine +
              timeFullGeoCheckAreaPoint + timeFullGeoCheckLineLine +
@@ -129,6 +136,38 @@ inline std::string Stats::toString() {
   ss << "time for obb intersections LINE/LINE: " << t << " s ("
      << ((t / sum) * 100.0) << "%)\n";
 
+  t = double(timeConvexHullIsectAreaArea) / 1000000000.0;
+  ss << "time for convex hull intersections AREA/AREA: " << t << " s ("
+     << ((t / sum) * 100.0) << "%)\n";
+
+  t = double(timeConvexHullIsectAreaLine) / 1000000000.0;
+  ss << "time for convex hull intersections AREA/LINE: " << t << " s ("
+     << ((t / sum) * 100.0) << "%)\n";
+
+  t = double(timeConvexHullIsectAreaPoint) / 1000000000.0;
+  ss << "time for convex hull intersections AREA/POINT: " << t << " s ("
+     << ((t / sum) * 100.0) << "%)\n";
+
+  t = double(timeConvexHullIsectLineLine) / 1000000000.0;
+  ss << "time for convex hull intersections LINE/LINE: " << t << " s ("
+     << ((t / sum) * 100.0) << "%)\n";
+
+  t = double(timeConvexHullIsectAreaArea) / 1000000000.0;
+  ss << "time for convex hull intersections AREA/AREA: " << t << " s ("
+     << ((t / sum) * 100.0) << "%)\n";
+
+  t = double(timeConvexHullIsectAreaLine) / 1000000000.0;
+  ss << "time for convex hull intersections AREA/LINE: " << t << " s ("
+     << ((t / sum) * 100.0) << "%)\n";
+
+  t = double(timeConvexHullIsectAreaPoint) / 1000000000.0;
+  ss << "time for convex hull intersections AREA/POINT: " << t << " s ("
+     << ((t / sum) * 100.0) << "%)\n";
+
+  t = double(timeConvexHullIsectLineLine) / 1000000000.0;
+  ss << "time for convex hull intersections LINE/LINE: " << t << " s ("
+     << ((t / sum) * 100.0) << "%)\n";
+
   t = double(timeFullGeoCheckAreaArea) / 1000000000.0;
   ss << "time for " << fullGeoChecksAreaArea
      << " full geom checks AREA/AREA: " << t << " s (" << ((t / sum) * 100.0)
@@ -280,6 +319,10 @@ inline Stats operator+(const Stats& a, const Stats& b) {
                a.timeOBBIsectAreaLine + b.timeOBBIsectAreaLine,
                a.timeOBBIsectAreaPoint + b.timeOBBIsectAreaPoint,
                a.timeOBBIsectLineLine + b.timeOBBIsectLineLine,
+               a.timeConvexHullIsectAreaArea + b.timeConvexHullIsectAreaArea,
+               a.timeConvexHullIsectAreaLine + b.timeConvexHullIsectAreaLine,
+               a.timeConvexHullIsectAreaPoint + b.timeConvexHullIsectAreaPoint,
+               a.timeConvexHullIsectLineLine + b.timeConvexHullIsectLineLine,
                a.timeFullGeoCheckAreaArea + b.timeFullGeoCheckAreaArea,
                a.timeFullGeoCheckAreaLine + b.timeFullGeoCheckAreaLine,
                a.timeFullGeoCheckAreaPoint + b.timeFullGeoCheckAreaPoint,

src/spatialjoin/Sweeper.cpp

@@ -46,6 +46,7 @@ using util::geo::webMercToLatLng;
 
 const static size_t CUTOUTS_MIN_SIZE = 100;
 const static size_t OBB_MIN_SIZE = 100;
+const static size_t CONVEXHULL_MIN_SIZE = 100;
 
 const static double sin45 = 1.0 / sqrt(2);
 const static double cos45 = 1.0 / sqrt(2);
@@ -136,6 +137,7 @@ void Sweeper::add(const I32Polygon& poly, const std::string& gid, size_t subid,
                   WriteBatch& batch) {
   WriteCand cur;
   const auto& box = getBoundingBox(poly);
+  const auto& hull = util::geo::convexHull(poly);
   I32XSortedPolygon spoly(poly);
 
   double areaSize = area(poly);
@@ -214,7 +216,6 @@ void Sweeper::add(const I32Polygon& poly, const std::string& gid, size_t subid,
     }
 
     util::geo::I32Polygon obb;
-
     if (_cfg.useOBB && poly.getOuter().size() >= OBB_MIN_SIZE) {
       obb = util::geo::convexHull(
           util::geo::pad(util::geo::getOrientedEnvelope(poly), 10));
@@ -223,11 +224,29 @@ void Sweeper::add(const I32Polygon& poly, const std::string& gid, size_t subid,
       if (obb.getOuter().size() >= poly.getOuter().size()) obb = {};
     }
 
+    util::geo::I32Polygon convexHull;
+    if (_cfg.useConvexHull && poly.getOuter().size() >= CONVEXHULL_MIN_SIZE) {
+      convexHull = util::geo::convexHull(poly);
+
+      // drop redundant convex hull
+      if (convexHull.getOuter().size() >= poly.getOuter().size())
+        convexHull = {};
+      if (!obb.getOuter().empty() && convexHull.getOuter().size() == obb.getOuter().size() && outerArea(convexHull) == outerArea(obb))
+        convexHull = {};
+      if (!outer.empty() && convexHull.getOuter().size() == outer.getOuter().rawRing().size())
+        convexHull = {};
+    }
+
+    // if (subid > 0) {
+    // std::unique_lock<std::mutex> lock(_multiAddMtx);
+    // multiAdd(gid, box.getLowerLeft().getX(), box.getUpperRight().getX());
+    // }
+
     std::stringstream str;
     GeometryCache<Area>::writeTo(
         {spoly, box, gid, subid, areaSize, _cfg.useArea ? outerAreaSize : 0,
-         boxIds, cutouts, obb, inner, innerBox, innerOuterAreaSize, outer,
-         outerBox, outerOuterAreaSize},
+         boxIds, cutouts, obb, convexHull, inner, innerBox, innerOuterAreaSize,
+         outer, outerBox, outerOuterAreaSize},
         str);
     ;
 
@@ -327,13 +346,23 @@ void Sweeper::add(const I32Line& line, const std::string& gid, size_t subid,
       if (obb.getOuter().size() >= line.size()) obb = {};
     }
 
+    util::geo::I32Polygon convexHull;
+    if (_cfg.useConvexHull && line.size() >= CONVEXHULL_MIN_SIZE) {
+      convexHull = util::geo::convexHull(line);
+
+      // drop redundant hull
+      if (convexHull.getOuter().size() >= line.size()) convexHull = {};
+      if (!obb.getOuter().empty() && convexHull.getOuter().size() == obb.getOuter().size() && outerArea(convexHull) == outerArea(obb))
+        convexHull = {};
+    }
+
     if (!_cfg.useFastSweepSkip) {
       sline.setMaxSegLen(std::numeric_limits<int32_t>::max());
     }
 
     std::stringstream str;
     GeometryCache<Line>::writeTo(
-        {sline, box, gid, subid, len, boxIds, cutouts, obb}, str);
+        {sline, box, gid, subid, len, boxIds, cutouts, obb, convexHull}, str);
     cur.raw = str.str();
 
     cur.boxvalIn = {0,  // placeholder, will be overwritten later on
@@ -1003,6 +1032,7 @@ sj::Area Sweeper::areaFromSimpleArea(const SimpleArea* sa) const {
           {},
           {},
           {},
+          {},
           0,
           {},
           {},
@@ -1054,13 +1084,19 @@ std::tuple<bool, bool, bool, bool, bool> Sweeper::check(const Area* a,
     }
   }
 
-  if (_cfg.useOBB) {
-    if (!a->obb.empty() && !b->obb.empty()) {
-      auto ts = TIME();
-      auto r = util::geo::intersectsContainsCovers(a->obb, b->obb);
-      _stats[t].timeOBBIsectAreaArea += TOOK(ts);
-      if (!std::get<0>(r)) return {0, 0, 0, 0, 0};
-    }
+  if (_cfg.useOBB && !a->obb.empty() && !b->obb.empty()) {
+    auto ts = TIME();
+    auto r = util::geo::intersectsContainsCovers(a->obb, b->obb);
+    _stats[t].timeOBBIsectAreaArea += TOOK(ts);
+    if (!std::get<0>(r)) return {0, 0, 0, 0, 0};
+  }
+
+  if (_cfg.useConvexHull && !a->convexHull.empty() && !b->convexHull.empty()) {
+    auto ts = TIME();
+    auto r = util::geo::intersectsContainsCovers(a->convexHull, a->box, 0,
+                                                 b->convexHull, b->box, 0);
+    _stats[t].timeConvexHullIsectAreaArea += TOOK(ts);
+    if (!std::get<0>(r)) return {0, 0, 0, 0, 0};
   }
 
   if (_cfg.useInnerOuter && !a->outer.empty() && !b->outer.empty()) {
@@ -1160,6 +1196,14 @@ std::tuple<bool, bool, bool, bool, bool> Sweeper::check(const Line* a,
     }
   }
 
+  if (_cfg.useConvexHull && !a->convexHull.empty() && !b->convexHull.empty()) {
+    auto ts = TIME();
+    auto r = intersectsContainsCovers(a->convexHull, a->box, 0, b->convexHull,
+                                      b->box, 0);
+    _stats[t].timeConvexHullIsectAreaLine += TOOK(ts);
+    if (!std::get<0>(r)) return {0, 0, 0, 0, 0};
+  }
+
   if (_cfg.useInnerOuter && !b->outer.empty()) {
     auto ts = TIME();
     auto r = util::geo::intersectsContainsCovers(a->geom, a->box, b->outer,
@@ -1229,13 +1273,19 @@ std::tuple<bool, bool, bool, bool, bool> Sweeper::check(const Line* a,
     }
   }
 
-  if (_cfg.useOBB) {
-    if (!a->obb.empty() && !b->obb.empty()) {
-      auto ts = TIME();
-      auto r = util::geo::intersectsContainsCovers(a->obb, b->obb);
-      _stats[t].timeOBBIsectLineLine += TOOK(ts);
-      if (!std::get<0>(r)) return {0, 0, 0, 0, 0};
-    }
+  if (_cfg.useOBB && !a->obb.empty() && !b->obb.empty()) {
+    auto ts = TIME();
+    auto r = util::geo::intersectsContainsCovers(a->obb, b->obb);
+    _stats[t].timeOBBIsectLineLine += TOOK(ts);
+    if (!std::get<0>(r)) return {0, 0, 0, 0, 0};
+  }
+
+  if (_cfg.useConvexHull && !a->convexHull.empty() && !b->convexHull.empty()) {
+    auto ts = TIME();
+    auto r = intersectsContainsCovers(a->convexHull, a->box, 0, b->convexHull,
+                                      b->box, 0);
+    _stats[t].timeConvexHullIsectLineLine += TOOK(ts);
+    if (!std::get<0>(r)) return {0, 0, 0, 0, 0};
   }
 
   auto ts = TIME();
@@ -1288,6 +1338,14 @@ std::tuple<bool, bool, bool, bool, bool> Sweeper::check(const SimpleLine* a,
     if (!std::get<0>(r)) return {0, 0, 0, 0, 0};
   }
 
+  if (_cfg.useConvexHull && !b->convexHull.empty()) {
+    auto ts = TIME();
+    auto r = intersectsContainsCovers(
+        I32XSortedLine(LineSegment<int32_t>(a->a, a->b)), b->convexHull);
+    _stats[t].timeConvexHullIsectLineLine += TOOK(ts);
+    if (!std::get<0>(r)) return {0, 0, 0, 0, 0};
+  }
+
   if (_cfg.useInnerOuter && !b->outer.empty()) {
     auto ts = TIME();
     auto r = util::geo::intersectsContainsCovers(
@@ -1476,6 +1534,13 @@ std::pair<bool, bool> Sweeper::check(const I32Point& a, const Area* b,
     if (!std::get<1>(r)) return {0, 0};
   }
 
+  if (_cfg.useConvexHull && !b->convexHull.empty()) {
+    auto ts = TIME();
+    auto r = containsCovers(a, b->convexHull);
+    _stats[t].timeConvexHullIsectAreaPoint += TOOK(ts);
+    if (!std::get<1>(r)) return {0, 0};
+  }
+
   if (_cfg.useInnerOuter && !b->outer.empty()) {
     auto ts = TIME();
     auto r = containsCovers(a, b->outer);

src/spatialjoin/Sweeper.h

@@ -109,6 +109,7 @@ struct SweeperCfg {
   bool useDiagBox;
   bool useFastSweepSkip;
   bool useInnerOuter;
+  bool useConvexHull;
 };
 
 // buffer size _must_ be multiples of sizeof(BoxVal)