Add pointInTriangle for 3D coordinates

CesiumGeometry/include/CesiumGeometry/IntersectionTests.h

@@ -42,6 +42,38 @@ public:
       const glm::dvec2& triangleVertA,
       const glm::dvec2& triangleVertB,
       const glm::dvec2& triangleVertC) noexcept;
+
+  /**
+   * @brief Determines whether the point is completely inside the triangle.
+   *
+   * @param point The point to check.
+   * @param triangleVertA The first vertex of the triangle.
+   * @param triangleVertB The second vertex of the triangle.
+   * @param triangleVertC The third vertex of the triangle.
+   * @return Whether the point is within the triangle.
+   */
+  static bool pointInTriangle(
+      const glm::dvec3& point,
+      const glm::dvec3& triangleVertA,
+      const glm::dvec3& triangleVertB,
+      const glm::dvec3& triangleVertC) noexcept;
+
+  /**
+   * @brief Determines whether the point is completely inside the triangle and
+   * outputs the barycentric coordinates for the point.
+   *
+   * @param point The point to check.
+   * @param triangleVertA The first vertex of the triangle.
+   * @param triangleVertB The second vertex of the triangle.
+   * @param triangleVertC The third vertex of the triangle.
+   * @return Whether the point is within the triangle.
+   */
+  static bool pointInTriangle(
+      const glm::dvec3& point,
+      const glm::dvec3& triangleVertA,
+      const glm::dvec3& triangleVertB,
+      const glm::dvec3& triangleVertC,
+      glm::dvec3& barycentricCoordinates) noexcept;
 };
 
 } // namespace CesiumGeometry

CesiumGeometry/src/IntersectionTests.cpp

@@ -6,6 +6,7 @@
 #include <CesiumUtility/Math.h>
 
 #include <glm/geometric.hpp>
+#include <glm/gtx/norm.hpp>
 
 using namespace CesiumUtility;
 
@@ -36,17 +37,21 @@ bool IntersectionTests::pointInTriangle2D(
     const glm::dvec2& triangleVertA,
     const glm::dvec2& triangleVertB,
     const glm::dvec2& triangleVertC) noexcept {
-
+  // Define the triangle edges.
   const glm::dvec2 ab = triangleVertB - triangleVertA;
-  const glm::dvec2 ab_perp(-ab.y, ab.x);
   const glm::dvec2 bc = triangleVertC - triangleVertB;
-  const glm::dvec2 bc_perp(-bc.y, bc.x);
   const glm::dvec2 ca = triangleVertA - triangleVertC;
+
+  // Get the vector perpendicular to each triangle edge in the same 2D plane).
+  const glm::dvec2 ab_perp(-ab.y, ab.x);
+  const glm::dvec2 bc_perp(-bc.y, bc.x);
   const glm::dvec2 ca_perp(-ca.y, ca.x);
 
+  // Get vectors from triangle corners to point.
   const glm::dvec2 av = point - triangleVertA;
   const glm::dvec2 cv = point - triangleVertC;
 
+  // Compute the area of the sub-triangle.
   const double v_proj_ab_perp = glm::dot(av, ab_perp);
   const double v_proj_bc_perp = glm::dot(cv, bc_perp);
   const double v_proj_ca_perp = glm::dot(cv, ca_perp);
@@ -58,4 +63,61 @@ bool IntersectionTests::pointInTriangle2D(
           v_proj_bc_perp <= 0.0);
 }
 
+bool IntersectionTests::pointInTriangle(
+    const glm::dvec3& point,
+    const glm::dvec3& triangleVertA,
+    const glm::dvec3& triangleVertB,
+    const glm::dvec3& triangleVertC) noexcept {
+  glm::dvec3 unused;
+  return IntersectionTests::pointInTriangle(
+      point,
+      triangleVertA,
+      triangleVertB,
+      triangleVertC,
+      unused);
+}
+
+bool IntersectionTests::pointInTriangle(
+    const glm::dvec3& point,
+    const glm::dvec3& triangleVertA,
+    const glm::dvec3& triangleVertB,
+    const glm::dvec3& triangleVertC,
+    glm::dvec3& barycentricCoordinates) noexcept {
+  // Define the triangle edges.
+  const glm::dvec3 ab = triangleVertB - triangleVertA;
+  const glm::dvec3 bc = triangleVertC - triangleVertB;
+
+  const glm::dvec3 triangleNormal = glm::cross(ab, bc);
+  const double lengthSquared = glm::length2(triangleNormal);
+  if (lengthSquared < CesiumUtility::Math::Epsilon8) {
+    // Don't handle invalid triangles.
+    return false;
+  }
+
+  // Technically this is the area of the parallelogram between the vectors,
+  // but triangles are 1/2 of the area, and the division will cancel out later.
+  const double triangleAreaInv = 1.0 / glm::sqrt(lengthSquared);
+
+  const glm::dvec3 ap = point - triangleVertA;
+  const double triangleABPRatio =
+      glm::length(glm::cross(ab, ap)) * triangleAreaInv;
+  if (triangleABPRatio < 0 || triangleABPRatio > 1) {
+    return false;
+  }
+
+  const glm::dvec3 bp = point - triangleVertB;
+  const double triangleBCPRatio =
+      glm::length(glm::cross(bc, bp)) * triangleAreaInv;
+  if (triangleBCPRatio < 0 || triangleBCPRatio > 1) {
+    return false;
+  }
+
+  const double triangleCAPRatio = 1.0 - triangleABPRatio - triangleBCPRatio;
+  barycentricCoordinates.x = triangleBCPRatio;
+  barycentricCoordinates.y = triangleCAPRatio;
+  barycentricCoordinates.z = triangleABPRatio;
+
+  return true;
+}
+
 } // namespace CesiumGeometry

CesiumGeometry/test/TestIntersectionTests.cpp

@@ -35,3 +35,214 @@ TEST_CASE("IntersectionTests::rayPlane") {
       IntersectionTests::rayPlane(testCase.ray, testCase.plane);
   CHECK(intersectionPoint == testCase.expectedIntersectionPoint);
 }
+
+TEST_CASE("IntersectionTests::pointInTriangle2D") {
+  struct TestCase {
+    glm::dvec2 point;
+    glm::dvec2 triangleVert1;
+    glm::dvec2 triangleVert2;
+    glm::dvec2 triangleVert3;
+    bool expected;
+  };
+
+  auto testCase = GENERATE(
+      // is in triangle (corner)
+      TestCase{
+          glm::dvec2(1.0, 0.0),
+          glm::dvec2(-1.0, 0.0),
+          glm::dvec2(0.0, 1.0),
+          glm::dvec2(1.0, 0.0),
+          true},
+      // is in triangle (edge)
+      TestCase{
+          glm::dvec2(0.0, 0.0),
+          glm::dvec2(-1.0, 0.0),
+          glm::dvec2(0.0, 1.0),
+          glm::dvec2(1.0, 0.0),
+          true},
+      // is in triangle (middle)
+      TestCase{
+          glm::dvec2(0.2, 0.5),
+          glm::dvec2(-1.0, 0.0),
+          glm::dvec2(0.0, 1.0),
+          glm::dvec2(1.0, 0.0),
+          true},
+      // is in triangle (regardless of winding order)
+      TestCase{
+          glm::dvec2(0.2, 0.5),
+          glm::dvec2(-1.0, 0.0),
+          glm::dvec2(1.0, 0.0),
+          glm::dvec2(0.0, 1.0),
+          true},
+      // is not in triangle
+      TestCase{
+          glm::dvec2(-2.0, 0.5),
+          glm::dvec2(-1.0, 0.0),
+          glm::dvec2(1.0, 0.0),
+          glm::dvec2(0.0, 1.0),
+          false},
+      // is not in triangle (regardless of winding order)
+      TestCase{
+          glm::dvec2(-2.0, 0.5),
+          glm::dvec2(-1.0, 0.0),
+          glm::dvec2(1.0, 0.0),
+          glm::dvec2(0.0, 1.0),
+          false});
+
+  bool result = IntersectionTests::pointInTriangle2D(
+      testCase.point,
+      testCase.triangleVert1,
+      testCase.triangleVert2,
+      testCase.triangleVert3);
+  CHECK(result == testCase.expected);
+}
+
+TEST_CASE("IntersectionTests::pointInTriangle") {
+  struct TestCase {
+    glm::dvec3 point;
+    glm::dvec3 triangleVert1;
+    glm::dvec3 triangleVert2;
+    glm::dvec3 triangleVert3;
+    bool expected;
+  };
+
+  auto testCase = GENERATE(
+      // is in triangle (corner)
+      TestCase{
+          glm::dvec3(1.0, 0.0, 0.0),
+          glm::dvec3(-1.0, 0.0, 0.0),
+          glm::dvec3(0.0, 1.0, 0.0),
+          glm::dvec3(1.0, 0.0, 0.0),
+          true},
+      // is in triangle (edge)
+      TestCase{
+          glm::dvec3(0.0, 0.0, 0.0),
+          glm::dvec3(-1.0, 0.0, 0.0),
+          glm::dvec3(0.0, 1.0, 0.0),
+          glm::dvec3(1.0, 0.0, 0.0),
+          true},
+      // is in triangle (middle)
+      TestCase{
+          glm::dvec3(0.2, 0.5, 0.0),
+          glm::dvec3(-1.0, 0.0, 0.0),
+          glm::dvec3(0.0, 1.0, 0.0),
+          glm::dvec3(1.0, 0.0, 0.0),
+          true},
+      // is in triangle (regardless of winding order)
+      TestCase{
+          glm::dvec3(0.2, 0.5, 0.0),
+          glm::dvec3(-1.0, 0.0, 0.0),
+          glm::dvec3(1.0, 0.0, 0.0),
+          glm::dvec3(0.0, 1.0, 0.0),
+          true},
+      // is not in triangle (same plane)
+      TestCase{
+          glm::dvec3(-2.0, 0.5, 0.0),
+          glm::dvec3(-1.0, 0.0, 0.0),
+          glm::dvec3(1.0, 0.0, 0.0),
+          glm::dvec3(0.0, 1.0, 0.0),
+          false},
+      // is not in triangle (different plane)
+      TestCase{
+          glm::dvec3(0.2, 0.5, 1.0),
+          glm::dvec3(-1.0, 0.0, 0.0),
+          glm::dvec3(1.0, 0.0, 0.0),
+          glm::dvec3(0.0, 1.0, 0.0),
+          false},
+      // is not in triangle (regardless of winding order)
+      TestCase{
+          glm::dvec3(-2.0, 0.5, 0.0),
+          glm::dvec3(-1.0, 0.0, 0.0),
+          glm::dvec3(1.0, 0.0, 0.0),
+          glm::dvec3(0.0, 1.0, 0.0),
+          false});
+
+  bool result = IntersectionTests::pointInTriangle(
+      testCase.point,
+      testCase.triangleVert1,
+      testCase.triangleVert2,
+      testCase.triangleVert3);
+  CHECK(result == testCase.expected);
+}
+
+TEST_CASE("IntersectionTests::pointInTriangle with barycentric coordinates") {
+  struct TestCase {
+    glm::dvec3 point;
+    glm::dvec3 triangleVert1;
+    glm::dvec3 triangleVert2;
+    glm::dvec3 triangleVert3;
+    bool expected;
+    glm::dvec3 expectedCoordinates;
+  };
+
+  auto testCase = GENERATE(
+      // is in triangle (corner)
+      TestCase{
+          glm::dvec3(1.0, 0.0, 0.0),
+          glm::dvec3(-1.0, 0.0, 0.0),
+          glm::dvec3(0.0, 1.0, 0.0),
+          glm::dvec3(1.0, 0.0, 0.0),
+          true,
+          glm::dvec3(0.0, 0.0, 1.0)},
+      // is in triangle (edge)
+      TestCase{
+          glm::dvec3(0.0, 0.0, 0.0),
+          glm::dvec3(-1.0, 0.0, 0.0),
+          glm::dvec3(0.0, 1.0, 0.0),
+          glm::dvec3(1.0, 0.0, 0.0),
+          true,
+          glm::dvec3(0.5, 0.0, 0.5)},
+      // is in triangle (middle)
+      TestCase{
+          glm::dvec3(0.0, 0.5, 0.0),
+          glm::dvec3(-1.0, 0.0, 0.0),
+          glm::dvec3(0.0, 1.0, 0.0),
+          glm::dvec3(1.0, 0.0, 0.0),
+          true,
+          glm::dvec3(0.25, 0.5, 0.25)},
+      // is in triangle (regardless of winding order)
+      TestCase{
+          glm::dvec3(0.0, 0.5, 0.0),
+          glm::dvec3(-1.0, 0.0, 0.0),
+          glm::dvec3(1.0, 0.0, 0.0),
+          glm::dvec3(0.0, 1.0, 0.0),
+          true,
+          glm::dvec3(0.25, 0.25, 0.5)},
+      // is not in triangle (same plane)
+      TestCase{
+          glm::dvec3(-2.0, 0.5, 0.0),
+          glm::dvec3(-1.0, 0.0, 0.0),
+          glm::dvec3(1.0, 0.0, 0.0),
+          glm::dvec3(0.0, 1.0, 0.0),
+          false,
+          glm::dvec3()},
+      // is not in triangle (different plane)
+      TestCase{
+          glm::dvec3(0.2, 0.5, 1.0),
+          glm::dvec3(-1.0, 0.0, 0.0),
+          glm::dvec3(1.0, 0.0, 0.0),
+          glm::dvec3(0.0, 1.0, 0.0),
+          false,
+          glm::dvec3()},
+      // is not in triangle (regardless of winding order)
+      TestCase{
+          glm::dvec3(-2.0, 0.5, 0.0),
+          glm::dvec3(-1.0, 0.0, 0.0),
+          glm::dvec3(1.0, 0.0, 0.0),
+          glm::dvec3(0.0, 1.0, 0.0),
+          false,
+          glm::dvec3()});
+
+  glm::dvec3 barycentricCoordinates;
+  bool result = IntersectionTests::pointInTriangle(
+      testCase.point,
+      testCase.triangleVert1,
+      testCase.triangleVert2,
+      testCase.triangleVert3,
+      barycentricCoordinates);
+
+  REQUIRE(result == testCase.expected);
+  if (result) {
+    CHECK(barycentricCoordinates == testCase.expectedCoordinates);
+  }
+}