Intersection tests for tile bounding volumes

Cesium3DTilesSelection/include/Cesium3DTilesSelection/BoundingVolume.h

@@ -96,4 +96,9 @@ getOrientedBoundingBoxFromBoundingVolume(
     const BoundingVolume& boundingVolume,
     const CesiumGeospatial::Ellipsoid& ellipsoid CESIUM_DEFAULT_ELLIPSOID);
 
+/**
+ * @brief Returns true if two bounding volumes intersect.
+ */
+CESIUM3DTILESSELECTION_API bool
+testIntersection(const BoundingVolume& volume0, const BoundingVolume& volume1);
 } // namespace Cesium3DTilesSelection

Cesium3DTilesSelection/src/BoundingVolume.cpp

@@ -17,6 +17,7 @@
 
 #include <cstdint>
 #include <optional>
+#include <type_traits>
 #include <variant>
 
 using namespace CesiumGeometry;
@@ -277,4 +278,118 @@ OrientedBoundingBox getOrientedBoundingBoxFromBoundingVolume(
   return std::visit(Operation{ellipsoid}, boundingVolume);
 }
 
+namespace {
+
+template <typename BV>
+bool intersectsImpl(const BoundingSphere& bs, const BV& bv1) {
+  double d2 = bv1.computeDistanceSquaredToPosition(bs.getCenter());
+  return d2 < bs.getRadius() * bs.getRadius();
+}
+bool intersectsImpl(
+    const BoundingSphere& bs,
+    const BoundingRegionWithLooseFittingHeights& bv1) {
+  double d2 = bv1.computeConservativeDistanceSquaredToPosition(bs.getCenter());
+  return d2 < bs.getRadius() * bs.getRadius();
+}
+
+bool intersectsImpl(
+    const OrientedBoundingBox& obb0,
+    const OrientedBoundingBox& obb1) {
+  return intersects(obb0, obb1);
+}
+
+bool intersectsImpl(const OrientedBoundingBox& obb, const BoundingRegion& br) {
+  return intersects(obb, br.getBoundingBox());
+}
+
+bool intersectsImpl(
+    const OrientedBoundingBox& obb,
+    const BoundingCylinderRegion& bcr) {
+  return intersects(obb, bcr.toOrientedBoundingBox());
+}
+
+template <typename BV>
+bool intersectsImpl(const OrientedBoundingBox& obb, const BV& bv) {
+  return intersectsImpl(obb.toSphere(), bv);
+}
+
+template <typename BV>
+bool intersectsImpl(const BoundingRegion& br, const BV& bv) {
+  return intersectsImpl(br.getBoundingBox(), bv);
+}
+
+} // namespace
+
+// This function produces a lot of "unreachable code" on Windows, but it should
+// be ok.
+
+#ifdef _MSC_VER
+#pragma warning(push)
+#pragma warning(disable : 4702)
+#endif
+
+bool testIntersection(
+    const BoundingVolume& volume0,
+    const BoundingVolume& volume1) {
+  return visit(
+      [&](auto&& bv0, auto&& bv1) -> bool {
+        // If one of the volumes is a Bounding sphere, then don't bother to
+        // convert; just test against the sphere.
+        if constexpr (std::is_same_v<decltype(bv0), const BoundingSphere&>) {
+          return intersectsImpl(bv0, bv1);
+        } else if constexpr (std::is_same_v<
+                                 decltype(bv1),
+                                 const BoundingSphere&>) {
+          return intersectsImpl(bv1, bv0);
+          ;
+        }
+        if constexpr (std::is_same_v<
+                          decltype(bv0),
+                          const OrientedBoundingBox&>) {
+          return intersectsImpl(bv0, bv1);
+        } else if constexpr (std::is_same_v<
+                                 decltype(bv1),
+                                 const OrientedBoundingBox&>) {
+          return intersectsImpl(bv1, bv0);
+        }
+        if constexpr (std::is_same_v<decltype(bv0), const BoundingRegion&>) {
+          return intersectsImpl(bv0, bv1);
+        } else if constexpr (std::is_same_v<
+                                 decltype(bv1),
+                                 const BoundingRegion&>) {
+          return intersectsImpl(bv1, bv0);
+        }
+        if constexpr (std::is_same_v<
+                          decltype(bv0),
+                          const BoundingRegionWithLooseFittingHeights&>) {
+          return intersectsImpl(bv0.getBoundingRegion(), bv1);
+        } else if constexpr (
+            std::is_same_v<
+                decltype(bv1),
+                const BoundingRegionWithLooseFittingHeights&>) {
+          return intersectsImpl(bv1.getBoundingRegion(), bv0);
+        }
+        if constexpr (
+            std::is_same_v<decltype(bv0), const S2CellBoundingVolume&> ||
+            std::is_same_v<decltype(bv1), const S2CellBoundingVolume&>) {
+          // not dealing yet
+          return false;
+        }
+        if constexpr (std::is_same_v<
+                          decltype(bv0),
+                          const BoundingCylinderRegion&>) {
+          return intersectsImpl(bv0.toOrientedBoundingBox(), bv1);
+        } else if constexpr (std::is_same_v<
+                                 decltype(bv1),
+                                 const BoundingCylinderRegion&>) {
+          return intersectsImpl(bv1.toOrientedBoundingBox(), bv0);
+        }
+        return false;
+      },
+      volume0,
+      volume1);
+}
+#ifdef _MSC_VER
+#pragma warning(pop)
+#endif
 } // namespace Cesium3DTilesSelection

Cesium3DTilesSelection/test/TestBoundingVolume.cpp

@@ -13,6 +13,8 @@
 
 #include <doctest/doctest.h>
 
+#include <numbers>
+
 using namespace Cesium3DTilesSelection;
 using namespace CesiumGeometry;
 using namespace CesiumGeospatial;
@@ -102,3 +104,46 @@ TEST_CASE("getOrientedBoundingBoxFromBoundingVolume") {
         newObb.getHalfAxes());
   }
 }
+
+TEST_CASE("testIntersection") {
+  SUBCASE("OrientedBoundingBox") {
+    // Simple tests with oriented bounding boxes
+    // Two boxes rotated 45 degrees around Z in opposite directions
+    using namespace std::numbers;
+    glm::dmat3 counterclockwise(
+        {sqrt2 / 2.0, sqrt2 / 2.0, 0.0},
+        {-sqrt2 / 2.0, sqrt2 / 2.0, 0.0},
+        {0.0, 0.0, 1.0});
+    glm::dmat3 clockwise(
+        {sqrt2 / 2.0, -sqrt2 / 2.0, 0.0},
+        {sqrt2 / 2.0, sqrt2 / 2.0, 0.0},
+        {0.0, 0.0, 1.0});
+    OrientedBoundingBox obb0({-2.0, 0.0, 0.0}, counterclockwise);
+    OrientedBoundingBox obb1({2.0, 0.0, 0.0}, clockwise);
+    CHECK(!testIntersection(obb0, obb1));
+    OrientedBoundingBox obb2({-1.0, 0.0, 0.0}, counterclockwise);
+    OrientedBoundingBox obb3({1.0, 0.0, 0.0}, clockwise);
+    CHECK(testIntersection(obb2, obb3));
+  }
+  SUBCASE("BoundingRegions") {
+    // A "nautical mile" square in Philadelphia
+    BoundingRegion phl(
+        GlobeRectangle(
+            -1.3120159199172432,
+            0.6969344194233807,
+            -1.311257041597562,
+            0.6975161958407122),
+        0.0,
+        300.0);
+    // A nautical mile square in NYC
+    BoundingRegion nyc(
+        GlobeRectangle(
+            -1.2921574402846652,
+            0.710289268896309,
+            -1.2913899085981675,
+            0.7108710453136404),
+        0.0,
+        300.0);
+    CHECK(!testIntersection(phl, nyc));
+  }
+}

CesiumGeometry/include/CesiumGeometry/AxisAlignedBox.h

@@ -7,6 +7,7 @@
 #include <vector>
 
 namespace CesiumGeometry {
+class OrientedBoundingBox;
 
 /**
  * @brief An Axis-Aligned Bounding Box (AABB), where the axes of the box are
@@ -133,4 +134,20 @@ struct CESIUMGEOMETRY_API AxisAlignedBox final {
   static AxisAlignedBox fromPositions(const std::vector<glm::dvec3>& positions);
 };
 
+/**
+ * @brief Test if two axis aligned boxes intersect.
+ */
+bool CESIUMGEOMETRY_API
+intersects(const AxisAlignedBox& b0, const AxisAlignedBox& b1);
+
+/**
+ * @brief Return the equivalent OrientedBoundingBox for an AxisAlignedBox. The
+ * axes will be the principal axes.
+ *
+ * @param box the source bounding box
+ *
+ * @returns An oriented bounding box
+ */
+OrientedBoundingBox toOrientedBoundingBox(const AxisAlignedBox& box);
+
 } // namespace CesiumGeometry

CesiumGeometry/include/CesiumGeometry/OrientedBoundingBox.h

@@ -147,4 +147,10 @@ class CESIUMGEOMETRY_API OrientedBoundingBox final {
   glm::dvec3 _lengths;
 };
 
+/**
+ * @brief Test if two oriented bounding boxes intersect.
+ */
+bool CESIUMGEOMETRY_API
+intersects(const OrientedBoundingBox& b1, const OrientedBoundingBox& b2);
+
 } // namespace CesiumGeometry

CesiumGeometry/src/AxisAlignedBox.cpp

@@ -1,4 +1,5 @@
 #include <CesiumGeometry/AxisAlignedBox.h>
+#include <CesiumGeometry/OrientedBoundingBox.h>
 
 #include <glm/common.hpp>
 #include <glm/ext/vector_double3.hpp>
@@ -23,4 +24,20 @@ AxisAlignedBox::fromPositions(const std::vector<glm::dvec3>& positions) {
 
   return AxisAlignedBox(min.x, min.y, min.z, max.x, max.y, max.z);
 }
+
+bool intersects(const AxisAlignedBox& b0, const AxisAlignedBox& b1) {
+  // Do all the axes overlap?
+  return b0.minimumX <= b1.maximumX && b0.maximumX >= b1.minimumX &&
+         b0.minimumY <= b1.maximumY && b0.maximumY >= b1.minimumY &&
+         b0.minimumZ <= b1.maximumZ && b0.maximumZ >= b1.minimumZ;
+}
+
+OrientedBoundingBox toOrientedBoundingBox(const AxisAlignedBox& box) {
+  return OrientedBoundingBox(
+      box.center,
+      glm::dmat3(
+          {box.lengthX * .5, 0.0, 0.0},
+          {0.0, box.lengthY * .5, 0.0},
+          {0.0, 0.0, box.lengthZ * .5}));
+}
 } // namespace CesiumGeometry

CesiumGeometry/src/OrientedBoundingBox.cpp

@@ -11,6 +11,8 @@
 #include <glm/ext/vector_double3.hpp>
 #include <glm/geometric.hpp>
 
+#include <cmath>
+
 namespace CesiumGeometry {
 CullingResult
 OrientedBoundingBox::intersectPlane(const Plane& plane) const noexcept {
@@ -225,4 +227,53 @@ OrientedBoundingBox::fromSphere(const BoundingSphere& sphere) noexcept {
   return OrientedBoundingBox(center, halfAxes);
 }
 
+bool intersects(const OrientedBoundingBox& b0, const OrientedBoundingBox& b1) {
+  const OrientedBoundingBox* boxes[2] = {&b0, &b1};
+  const glm::dvec3 T = b1.getCenter() - b0.getCenter();
+  const glm::dvec3 lengths[2] = {b0.getLengths() / 2.0, b1.getLengths() / 2.0};
+  // Normalized axes
+  glm::dmat3 axes[2] = {b0.getHalfAxes(), b1.getHalfAxes()};
+  for (int j = 0; j < 2; ++j) {
+    for (int i = 0; i < 3; ++i) {
+      axes[j][i] = axes[j][i] / lengths[j][i];
+    }
+  }
+  auto faceTest = [&](int box, int face, int other) {
+    const auto& otherHalfAxes = boxes[other]->getHalfAxes();
+    double projectedBoxLengths = lengths[box][face];
+    for (int i = 0; i < 3; ++i) {
+      projectedBoxLengths +=
+          std::abs(glm::dot(otherHalfAxes[i], axes[box][face]));
+    }
+    return std::abs(glm::dot(T, axes[box][face])) > projectedBoxLengths;
+  };
+  // Test the three face planes of the boxes against each box
+  for (int box = 0; box < 2; ++box) {
+    for (int face = 0; face < 3; ++face) {
+      if (faceTest(box, face, (box + 1) % 2)) {
+        return false;
+      }
+    }
+  }
+  // Test planes made by edges of box 0 and box 1
+  auto edgeTest = [&](int edge0, int edge1) {
+    glm::dvec3 planeNormal = glm::cross(axes[0][edge0], axes[1][edge1]);
+    double projectedBoxLengths = 0.0;
+    for (int box = 0; box < 2; ++box) {
+      for (int i = 0; i < 3; ++i) {
+        projectedBoxLengths +=
+            std::abs(glm::dot(boxes[box]->getHalfAxes()[i], planeNormal));
+      }
+    }
+    return std::abs(glm::dot(T, planeNormal)) > projectedBoxLengths;
+  };
+  for (int a = 0; a < 3; ++a) {
+    for (int b = a; b < 3; ++b) {
+      if (edgeTest(a, b)) {
+        return false;
+      }
+    }
+  }
+  return true;
+}
 } // namespace CesiumGeometry