add gltf ray intersect function with tests

CHANGES.md

@@ -2,6 +2,11 @@
 
 ### ? - ?
 
+##### Additions :tada:
+
+- Added `rayTriangle` intersection function that returns the intersection point between a ray and a triangle.
+- Added `intersectRayGltfModel` intersection function that returns the first intersection point between a ray and a glTF model.
+
 ##### Fixes :wrench:
 
 - Fixed a crash in `SubtreeAvailability::loadSubtree`.

CesiumGeometry/CMakeLists.txt

@@ -12,6 +12,7 @@ set_target_properties(CesiumGeometry
     PROPERTIES
         TEST_SOURCES "${CESIUM_GEOMETRY_TEST_SOURCES}"
         TEST_HEADERS "${CESIUM_GEOMETRY_TEST_HEADERS}"
+        TEST_DATA_DIR ${CMAKE_CURRENT_LIST_DIR}/test/data
 )
 
 set_target_properties(CesiumGeometry

CesiumGeometry/include/CesiumGeometry/IntersectionTests.h

@@ -28,6 +28,44 @@ class CESIUMGEOMETRY_API IntersectionTests final {
   static std::optional<glm::dvec3>
   rayPlane(const Ray& ray, const Plane& plane) noexcept;
 
+  /**
+   * @brief Tests if a ray hits a triangle and returns the hit point, if any.
+   *
+   * @param ray The ray.
+   * @param p0 The first vertex of the triangle.
+   * @param p1 The second vertex of the triangle.
+   * @param p2 The third vertex of the triangle.
+   * @param cullBackFaces Whether to cull back faces or not.
+   * @return The hit point, if any.
+   */
+  static std::optional<glm::dvec3> rayTriangle(
+      const Ray& ray,
+      const glm::dvec3& p0,
+      const glm::dvec3& p1,
+      const glm::dvec3& p2,
+      bool cullBackFaces = false);
+
+  /**
+   * @brief Tests if a ray hits a triangle and outputs the distance to the hit
+   * point, if any.
+   *
+   * @param ray The ray.
+   * @param p0 The first vertex of the triangle.
+   * @param p1 The second vertex of the triangle.
+   * @param p2 The third vertex of the triangle.
+   * @param cullBackFaces Whether to cull back faces or not.
+   * @param[out] t The distance from the ray origin to the intersection point,
+   * if any.
+   * @return Whether the ray intersects the triangle.
+   */
+  static bool rayTriangleParametric(
+      const Ray& ray,
+      const glm::dvec3& p0,
+      const glm::dvec3& p1,
+      const glm::dvec3& p2,
+      double& t,
+      bool cullBackFaces = false);
+
   /**
    * @brief Determines whether the point is completely inside the triangle.
    *

CesiumGeometry/src/IntersectionTests.cpp

@@ -31,6 +31,81 @@ IntersectionTests::rayPlane(const Ray& ray, const Plane& plane) noexcept {
   return ray.getOrigin() + ray.getDirection() * t;
 }
 
+std::optional<glm::dvec3> IntersectionTests::rayTriangle(
+    const Ray& ray,
+    const glm::dvec3& V0,
+    const glm::dvec3& V1,
+    const glm::dvec3& V2,
+    bool cullBackFaces) {
+  double t;
+  if (rayTriangleParametric(ray, V0, V1, V2, t, cullBackFaces)) {
+    return std::make_optional<glm::dvec3>(
+        ray.getOrigin() + t * ray.getDirection());
+  } else {
+    return std::optional<glm::dvec3>();
+  }
+}
+
+bool IntersectionTests::rayTriangleParametric(
+    const Ray& ray,
+    const glm::dvec3& p0,
+    const glm::dvec3& p1,
+    const glm::dvec3& p2,
+    double& t,
+    bool cullBackFaces) {
+
+  const glm::dvec3& origin = ray.getOrigin();
+  const glm::dvec3& direction = ray.getDirection();
+
+  glm::dvec3 edge0 = p1 - p0;
+  glm::dvec3 edge1 = p2 - p0;
+
+  glm::dvec3 p = glm::cross(direction, edge1);
+  double det = glm::dot(edge0, p);
+  if (cullBackFaces) {
+    if (det < Math::Epsilon6)
+      return false;
+
+    glm::dvec3 tvec = origin - p0;
+    double u = glm::dot(tvec, p);
+    if (u < 0.0 || u > det)
+      return false;
+
+    glm::dvec3 q = glm::cross(tvec, edge0);
+    double v = glm::dot(direction, q);
+    if (v < 0.0 || u + v > det)
+      return false;
+
+    t = glm::dot(edge1, q);
+    if (t < 0) {
+      return false;
+    }
+    t /= det;
+
+  } else {
+
+    if (std::abs(det) < Math::Epsilon6)
+      return false;
+
+    double invDet = 1.0 / det;
+
+    glm::dvec3 tvec = origin - p0;
+    double u = glm::dot(tvec, p) * invDet;
+    if (u < 0.0 || u > 1.0)
+      return false;
+
+    glm::dvec3 q = glm::cross(tvec, edge0);
+    double v = glm::dot(direction, q) * invDet;
+    if (v < 0.0 || u + v > 1.0)
+      return false;
+
+    t = glm::dot(edge1, q) * invDet;
+    if (t < 0.0)
+      return false;
+  }
+  return true;
+}
+
 bool IntersectionTests::pointInTriangle2D(
     const glm::dvec2& point,
     const glm::dvec2& triangleVertA,

CesiumGeometry/test/TestIntersectionTests.cpp

@@ -35,3 +35,71 @@ TEST_CASE("IntersectionTests::rayPlane") {
       IntersectionTests::rayPlane(testCase.ray, testCase.plane);
   CHECK(intersectionPoint == testCase.expectedIntersectionPoint);
 }
+
+TEST_CASE("IntersectionTests::rayTriangle") {
+
+  glm::dvec3 V0 = glm::dvec3(-1.0, 0.0, 0.0);
+  glm::dvec3 V1 = glm::dvec3(1.0, 0.0, 0.0);
+  glm::dvec3 V2 = glm::dvec3(0.0, 1.0, 0.0);
+
+  struct TestCase {
+    Ray ray;
+    bool cullBackFaces;
+    std::optional<glm::dvec3> expectedIntersectionPoint;
+  };
+
+  auto testCase = GENERATE(
+      // rayTriangle intersects front face
+      TestCase{
+          Ray(glm::dvec3(0.0, 0.0, 1.0), glm::dvec3(0.0, 0.0, -1.0)),
+          false,
+          glm::dvec3(0.0, 0.0, 0.0)},
+      // rayTriangle intersects back face without culling
+      TestCase{
+          Ray(glm::dvec3(0.0, 0.0, -1.0), glm::dvec3(0.0, 0.0, 1.0)),
+          false,
+          glm::dvec3(0.0, 0.0, 0.0)},
+      // rayTriangle does not intersect back face with culling
+      TestCase{
+          Ray(glm::dvec3(0.0, 0.0, -1.0), glm::dvec3(0.0, 0.0, 1.0)),
+          true,
+          std::nullopt},
+      // rayTriangle does not intersect outside the 0-1 edge
+      TestCase{
+          Ray(glm::dvec3(0.0, -1.0, 1.0), glm::dvec3(0.0, 0.0, -1.0)),
+          false,
+          std::nullopt},
+      // rayTriangle does not intersect outside the 1-2 edge
+      TestCase{
+          Ray(glm::dvec3(1.0, 1.0, 10.0), glm::dvec3(0.0, 0.0, -1.0)),
+          false,
+          std::nullopt},
+      // rayTriangle does not intersect outside the 2-0 edge
+      TestCase{
+          Ray(glm::dvec3(2.0, 0.0, 0.0), glm::dvec3(0.0, 1.0, 0.0)),
+          false,
+          std::nullopt},
+      // rayTriangle does not intersect parallel ray and triangle
+      TestCase{
+          Ray(glm::dvec3(-1.0, 1.0, 1.0), glm::dvec3(0.0, 0.0, -1.0)),
+          false,
+          std::nullopt},
+      // rayTriangle does not intersect parallel ray and triangle
+      TestCase{
+          Ray(glm::dvec3(-1.0, 0.0, 1.0), glm::dvec3(1.0, 0.0, 0.0)),
+          false,
+          std::nullopt},
+      // rayTriangle does not intersect behind the ray origin
+      TestCase{
+          Ray(glm::dvec3(0.0, 0.0, 1.0), glm::dvec3(0.0, 0.0, 1.0)),
+          false,
+          std::nullopt});
+
+  std::optional<glm::dvec3> intersectionPoint = IntersectionTests::rayTriangle(
+      testCase.ray,
+      V0,
+      V1,
+      V2,
+      testCase.cullBackFaces);
+  CHECK(intersectionPoint == testCase.expectedIntersectionPoint);
+}

CesiumGltfContent/include/CesiumGltfContent/GltfUtilities.h

@@ -15,6 +15,10 @@ struct Buffer;
 struct Model;
 } // namespace CesiumGltf
 
+namespace CesiumGeometry {
+class Ray;
+} // namespace CesiumGeometry
+
 namespace CesiumGltfContent {
 /**
  * A collection of utility functions that are used to process and transform a
@@ -121,5 +125,23 @@ struct CESIUMGLTFCONTENT_API GltfUtilities {
       CesiumGltf::Model& gltf,
       CesiumGltf::Buffer& destination,
       CesiumGltf::Buffer& source);
+
+  /**
+   * @brief Intersects a ray with a glTF model and returns the first
+   * intersection point.
+   *
+   * This function only handles primitives with TRIANGLES or TRIANGLE_FAN mode.
+   * Other modes are ignored.
+   *
+   * @param ray The ray.
+   * @param gltf The glTF model.
+   * @param cullBackFaces An optional boolean flag to indicate whether to cull
+   * backfaces or not. Defaults to true.
+   * @param return The intersection point along the ray, if any.
+   */
+  static std::optional<glm::dvec3> intersectRayGltfModel(
+      const CesiumGeometry::Ray& ray,
+      CesiumGltf::Model& gltf,
+      bool cullBackFaces = true);
 };
 } // namespace CesiumGltfContent

CesiumGltfContent/src/GltfUtilities.cpp

@@ -1,4 +1,6 @@
 #include <CesiumGeometry/Axis.h>
+#include <CesiumGeometry/IntersectionTests.h>
+#include <CesiumGeometry/Ray.h>
 #include <CesiumGeometry/Transforms.h>
 #include <CesiumGeospatial/BoundingRegionBuilder.h>
 #include <CesiumGltf/AccessorView.h>
@@ -213,4 +215,126 @@ GltfUtilities::parseGltfCopyright(const CesiumGltf::Model& gltf) {
   }
 }
 
+namespace {
+template <class T>
+static void intersectRayPrimitive(
+    const CesiumGeometry::Ray& ray,
+    const CesiumGltf::Model& model,
+    const CesiumGltf::MeshPrimitive& primitive,
+    bool cullBackFaces,
+    double& tMin) {
+
+  if (primitive.mode != MeshPrimitive::Mode::TRIANGLES &&
+      primitive.mode != MeshPrimitive::Mode::TRIANGLE_STRIP) {
+    return;
+  }
+
+  auto positionAccessorIt = primitive.attributes.find("POSITION");
+  if (positionAccessorIt == primitive.attributes.end()) {
+    return;
+  }
+
+  int positionAccessorID = positionAccessorIt->second;
+  const Accessor* pPositionAccessor =
+      Model::getSafe(&model.accessors, positionAccessorID);
+  if (!pPositionAccessor) {
+    return;
+  }
+
+  AccessorView<T> indicesView(model, primitive.indices);
+  AccessorView<glm::vec3> positionView(model, *pPositionAccessor);
+  double tCurr;
+
+  if (primitive.mode == CesiumGltf::MeshPrimitive::Mode::TRIANGLES) {
+    for (int32_t i = 0; i < indicesView.size(); i += 3) {
+      if (CesiumGeometry::IntersectionTests::rayTriangleParametric(
+              ray,
+              glm::dvec3(positionView[static_cast<int32_t>(indicesView[i])]),
+              glm::dvec3(
+                  positionView[static_cast<int32_t>(indicesView[i + 1])]),
+              glm::dvec3(
+                  positionView[static_cast<int32_t>(indicesView[i + 2])]),
+              tCurr,
+              cullBackFaces)) {
+        if (tCurr < tMin) {
+          tMin = tCurr;
+        }
+      }
+    }
+  } else {
+    for (int32_t i = 0; i < indicesView.size() - 2; ++i) {
+      if (i % 2) {
+        CesiumGeometry::IntersectionTests::rayTriangleParametric(
+            ray,
+            glm::dvec3(positionView[static_cast<int32_t>(indicesView[i])]),
+            glm::dvec3(positionView[static_cast<int32_t>(indicesView[i + 2])]),
+            glm::dvec3(positionView[static_cast<int32_t>(indicesView[i + 1])]),
+            tCurr,
+            true);
+      } else {
+        CesiumGeometry::IntersectionTests::rayTriangleParametric(
+            ray,
+            glm::dvec3(positionView[static_cast<int32_t>(indicesView[i])]),
+            glm::dvec3(positionView[static_cast<int32_t>(indicesView[i + 1])]),
+            glm::dvec3(positionView[static_cast<int32_t>(indicesView[i + 2])]),
+            tCurr,
+            true);
+      }
+    }
+  }
+}
+} // namespace
+
+std::optional<glm::dvec3> GltfUtilities::intersectRayGltfModel(
+    const CesiumGeometry::Ray& ray,
+    CesiumGltf::Model& gltf,
+    bool cullBackFaces) {
+  double t = DBL_MAX;
+  gltf.forEachPrimitiveInScene(
+      -1,
+      [ray, cullBackFaces, &t](
+          const CesiumGltf::Model& model,
+          const CesiumGltf::Node& /*node*/,
+          const CesiumGltf::Mesh& /*mesh*/,
+          const CesiumGltf::MeshPrimitive& primitive,
+          const glm::dmat4& /*nodeTransform*/) {
+        if (primitive.mode < CesiumGltf::MeshPrimitive::Mode::TRIANGLES)
+          return;
+
+        switch (model.accessors[primitive.indices].componentType) {
+        case Accessor::ComponentType::UNSIGNED_BYTE:
+          intersectRayPrimitive<uint8_t>(
+              ray,
+              model,
+              primitive,
+              cullBackFaces,
+              t);
+          break;
+        case Accessor::ComponentType::UNSIGNED_SHORT:
+          intersectRayPrimitive<uint16_t>(
+              ray,
+              model,
+              primitive,
+              cullBackFaces,
+              t);
+          break;
+        case Accessor::ComponentType::UNSIGNED_INT:
+          intersectRayPrimitive<uint32_t>(
+              ray,
+              model,
+              primitive,
+              cullBackFaces,
+              t);
+          break;
+        }
+      });
+
+  if (t != DBL_MAX) {
+    return std::make_optional<glm::dvec3>(
+        ray.getOrigin() + t * ray.getDirection());
+  } else {
+    return std::optional<glm::dvec3>();
+  }
+}
+
 } // namespace CesiumGltfContent