Convex hull

src/spatialjoin/GeometryCache.cpp

@@ -144,6 +144,9 @@ sj::Line sj::GeometryCache<sj::Line>::getFromDisk(size_t off,
     ret.cutouts[boxid] = cutout;
   }
 
+  // convex hull
+  readPoly(_geomsFReads[tid], ret.convexHull);
+
   // OBB
   readPoly(_geomsFReads[tid], ret.obb);
 
@@ -212,6 +215,9 @@ sj::Area sj::GeometryCache<sj::Area>::getFromDisk(size_t off,
     ret.cutouts[boxid] = cutout;
   }
 
+  // convex hull
+  readPoly(_geomsFReads[tid], ret.convexHull);
+
   // OBB
   readPoly(_geomsFReads[tid], ret.obb);
 
@@ -313,6 +319,9 @@ size_t sj::GeometryCache<sj::Line>::add(const sj::Line& val) {
     _geomsOffset += sizeof(int32_t) + sizeof(int32_t);
   }
 
+  // convex hull
+  writePoly(val.convexHull);
+
   // OBB
   writePoly(val.obb);
 
@@ -381,6 +390,9 @@ size_t sj::GeometryCache<sj::Area>::add(const sj::Area& val) {
     _geomsOffset += sizeof(int32_t) + sizeof(int32_t);
   }
 
+  // convex hull
+  writePoly(val.convexHull);
+
   // OBB
   writePoly(val.obb);
 

src/spatialjoin/GeometryCache.h

@@ -55,6 +55,9 @@ struct Area {
   // cutouts
   std::map<int32_t, size_t> cutouts;
 
+  // Convex hull
+  util::geo::I32XSortedPolygon convexHull;
+
   // OBB
   util::geo::I32XSortedPolygon obb;
 };
@@ -90,6 +93,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

@@ -103,6 +103,7 @@ int main(int argc, char** argv) {
   bool useOBB = true;
   bool useCutouts = true;
   bool useDiagBox = true;
+  bool useConvexHulls = true;
 
   for (int i = 1; i < argc; i++) {
     std::string cur = argv[i];
@@ -147,6 +148,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 {
           std::cerr << "Unknown option '" << cur << "', see -h" << std::endl;
           exit(1);
@@ -208,7 +211,7 @@ int main(int argc, char** argv) {
 
   Sweeper sweeper({NUM_THREADS, prefix, intersects, contains, covers, touches,
                    equals, overlaps, crosses, suffix, useBoxIds, useArea,
-                   useOBB, useCutouts, useDiagBox},
+                   useOBB, useCutouts, useDiagBox, useConvexHulls},
                   useCache, cache, output);
 
   if (!useCache) {

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;
@@ -56,7 +61,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 +
@@ -115,6 +122,22 @@ 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(timeFullGeoCheckAreaArea) / 1000000000.0;
   ss << "time for " << fullGeoChecksAreaArea
      << " full geom checks AREA/AREA: " << t << " s (" << ((t / sum) * 100.0)

src/spatialjoin/Sweeper.cpp

@@ -145,6 +145,11 @@ void Sweeper::add(const I32Polygon& poly, const std::string& gid,
     }
   }
 
+  util::geo::I32Polygon convexHull;
+  if (_cfg.useConvexHull) {
+    convexHull = util::geo::convexHull(poly);
+  }
+
   util::geo::I32Polygon obb;
   obb = util::geo::convexHull(util::geo::getOrientedEnvelope(poly));
 
@@ -161,7 +166,8 @@ void Sweeper::add(const I32Polygon& poly, const std::string& gid,
   if (!_cfg.useArea) areaSize = 0;
 
   size_t id = _areaCache.add(
-      {spoly, box, gid, subid, areaSize, outerAreaSize, boxIds, cutouts, obb});
+      {spoly, box, gid, subid, areaSize, outerAreaSize, boxIds, cutouts,
+       convexHull, obb});
 
   diskAdd({id, box.getLowerLeft().getY(), box.getUpperRight().getY(),
            box.getLowerLeft().getX(), false, POLYGON, areaSize, box45});
@@ -193,6 +199,11 @@ void Sweeper::add(const I32Line& line, const std::string& gid, size_t subid) {
 
   double len = util::geo::len(line);
 
+  util::geo::I32Polygon convexHull;
+  if (_cfg.useConvexHull) {
+    convexHull = util::geo::convexHull(line);
+  }
+
   util::geo::I32Polygon obb;
   obb = util::geo::convexHull(util::geo::getOrientedEnvelope(line));
 
@@ -222,7 +233,8 @@ void Sweeper::add(const I32Line& line, const std::string& gid, size_t subid) {
     const I32XSortedLine sline(line);
 
     size_t id =
-        _lineCache.add({sline, box, gid, subid, len, boxIds, cutouts, obb});
+        _lineCache.add({sline, box, gid, subid, len, boxIds, cutouts,
+                        convexHull, obb});
 
     diskAdd({id, box.getLowerLeft().getY(), box.getUpperRight().getY(),
              box.getLowerLeft().getX(), false, LINE, len, box45});
@@ -641,6 +653,15 @@ std::tuple<bool, bool, bool, bool, bool> Sweeper::check(const Area* a,
     }
   }
 
+  if (_cfg.useConvexHull) {
+    auto ts = TIME();
+    auto r = util::geo::intersectsContainsCovers(
+        a->convexHull, a->box, a->outerArea,
+        b->convexHull, b->box, b->outerArea);
+    _stats[t].timeConvexHullIsectAreaArea += TOOK(ts);
+    if (!std::get<0>(r)) return {0, 0, 0, 0, 0};
+  }
+
   auto ts = TIME();
   auto res = intersectsContainsCovers(a->geom, a->box, a->outerArea, b->geom,
                                       b->box, b->outerArea);
@@ -699,6 +720,14 @@ std::tuple<bool, bool, bool, bool, bool> Sweeper::check(const Line* a,
     }
   }
 
+  if (_cfg.useConvexHull) {
+    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};
+  }
+
   auto ts = TIME();
   auto res = intersectsContainsCovers(a->geom, a->box, b->geom, b->box);
   _stats[t].timeFullGeoCheckAreaLine += TOOK(ts);
@@ -761,6 +790,14 @@ std::tuple<bool, bool, bool, bool, bool> Sweeper::check(const Line* a,
     }
   }
 
+  if (_cfg.useConvexHull) {
+    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();
   auto res = intersectsCovers(a->geom, b->geom, a->box, b->box);
   _stats[t].timeFullGeoCheckLineLine += TOOK(ts);
@@ -983,6 +1020,13 @@ std::pair<bool, bool> Sweeper::check(const I32Point& a, const Area* b,
     if (!std::get<1>(r)) return {0, 0};
   }
 
+  if (_cfg.useConvexHull) {
+    auto ts = TIME();
+    auto r = containsCovers(a, b->convexHull);
+    _stats[t].timeConvexHullIsectAreaPoint += TOOK(ts);
+    if (!std::get<0>(r)) return {0, 0};
+  }
+
   auto ts = TIME();
   auto res = containsCovers(a, b->geom);
   _stats[t].timeFullGeoCheckAreaPoint += TOOK(ts);

src/spatialjoin/Sweeper.h

@@ -84,6 +84,7 @@ struct SweeperCfg {
   bool useOBB;
   bool useCutouts;
   bool useDiagBox;
+  bool useConvexHull;
 };
 
 // buffer sizes _must_ be multiples of sizeof(BoxVal)

src/spatialjoin/tests/TestMain.cpp

@@ -17,7 +17,7 @@ using sj::Sweeper;
 std::string fullRun(const std::string& file) {
   Sweeper sweeper({1, "$", " intersects ", " contains ", " covers ",
                    " touches ", " equals ", " overlaps ", " crosses ", "$\n",
-                   true, true, true, true, true},
+                   true, true, true, true, true, true},
                   false, ".", ".resTmp");
 
   size_t gid = 0;