Added GPRTRayTracer::occluded

CMakeLists.txt

@@ -249,6 +249,8 @@ if (XDG_ENABLE_GPRT)
     ${device_code}
     HEADERS
     ${CMAKE_CURRENT_SOURCE_DIR}/include/xdg/gprt/shared_structs.h
+    ${CMAKE_CURRENT_SOURCE_DIR}/include/xdg/shared_enums.h
+    ${CMAKE_CURRENT_SOURCE_DIR}/include/xdg/geometry/dp_math.h
     SOURCES
     ${CMAKE_CURRENT_SOURCE_DIR}/src/gprt/${device_code}.slang
     )

include/xdg/geometry/dp_math.h

@@ -0,0 +1,43 @@
+#ifndef DP_MATH_H
+#define DP_MATH_H
+
+/*
+This header acts as a light wrapper to provide a common interface for vector math operations
+in both C++ and Slang compilation contexts. It defines a `dp::vec3` type and associated
+functions for dot product, cross product, and absolute value. In C++ compilation, it maps to `xdg::Vec3da`,
+while in Slang compilation, it maps to `double3`.
+*/
+
+#ifdef __SLANG__
+
+// Slang compilation, map dp::vec3 -> double3
+namespace dp { 
+  typedef double3 vec3;
+
+  inline double dot(vec3 a, vec3 b) { return ::dot(a, b); }
+  inline vec3 cross(vec3 a, vec3 b) { return ::cross(a, b); }
+  inline double abs(double a) { return ::abs(a); }
+
+  static const double DBL_ZERO_TOL = 20 * DBL_EPSILON;
+  static const double INFTY = 1.7976931348623157e+308; // std::numeric_limits<double>::max() is not available in slang
+}
+
+#else
+#include "xdg/vec3da.h"
+
+// C++ compilation map dp::vec3 -> xdg::Vec3da
+namespace dp {
+  typedef xdg::Vec3da vec3;
+
+  inline double dot(vec3 a, vec3 b) { return xdg::dot(a, b); }
+  inline vec3 cross(vec3 a, vec3 b) { return xdg::cross(a, b); }
+  inline double abs(double a) { return std::fabs(a); }
+
+  static constexpr double DBL_ZERO_TOL = 20.0 * std::numeric_limits<double>::epsilon();
+  constexpr double INFTY {std::numeric_limits<double>::max()};
+
+}
+
+#endif // end of ifdef __slang__
+
+#endif // DP_MATH_H

include/xdg/geometry/plucker.h

@@ -1,11 +1,10 @@
 #ifndef _XDG_PLUCKER_H
 #define _XDG_PLUCKER_H
 
-#include "xdg/vec3da.h"
+#include "dp_math.h"
 
 namespace xdg {
 
-
 /*
  * Triangle vertex ordering convention:
  *
@@ -20,13 +19,146 @@ namespace xdg {
  * (normal pointing out of the plane). This ordering is based on the reference:
  * https://doi.org/10.1002/cnm.1237
  */
-bool plucker_ray_tri_intersect(const std::array<Position, 3> vertices,
-                               const Position& origin,
-                               const Direction& direction,
-                               double& dist_out,
-                               const double nonneg_ray_len = INFTY,
-                               const double* neg_ray_len = nullptr,
-                               const int* orientation = nullptr);
+
+struct PluckerIntersectionResult {
+  bool hit = false;    // Whether an intersection occurred
+  double t = 0.0;      // Distance along the ray to the intersection point
+};
+
+constexpr PluckerIntersectionResult EXIT_EARLY = {false, 0.0};
+
+/* Function to return the vertex with the lowest coordinates. To force the same
+  ray-edge computation, the Plücker test needs to use consistent edge
+  representation. This would be more simple with MOAB handles instead of
+  coordinates... 
+*/
+inline bool first(dp::vec3 a, dp::vec3 b) {
+  if (a[0] < b[0]) return true;  
+  if (a[0] > b[0]) return false;
+
+  if (a[1] < b[1]) return true;  
+  if (a[1] > b[1]) return false;
+
+  return a[2] < b[2];
+}
+
+inline double plucker_edge_test(dp::vec3 vertexa, dp::vec3 vertexb,
+                                dp::vec3 ray, dp::vec3 ray_normal)
+{
+  double pip;
+  if (first(vertexa, vertexb)) {
+    const dp::vec3 edge = vertexb - vertexa;
+    const dp::vec3 edge_normal = dp::cross(edge, vertexa);
+    pip = dp::dot(ray, edge_normal) + dp::dot(ray_normal, edge);
+  } else {
+    const dp::vec3 edge = vertexa - vertexb;
+    const dp::vec3 edge_normal = dp::cross(edge, vertexb);
+    pip = dp::dot(ray, edge_normal) + dp::dot(ray_normal, edge);
+    pip = -pip;
+  }
+  if (dp::DBL_ZERO_TOL > dp::abs(pip))  // <-- absd
+    pip = 0.0;
+  return pip;
+}
+
+inline PluckerIntersectionResult plucker_ray_tri_intersect(dp::vec3 vertices[3],
+                                dp::vec3 origin,
+                                dp::vec3 direction,
+                                double tMax,
+                                double tMin,
+                                bool useOrientation,
+                                int orientation)
+{
+  double dist_out = dp::INFTY;
+
+  const dp::vec3 raya = direction;
+  const dp::vec3 rayb = dp::cross(direction, origin);
+
+  // Determine the value of the first Plucker coordinate from edge 0
+  double plucker_coord0 =
+    plucker_edge_test(vertices[0], vertices[1], raya, rayb);
+
+  // If orientation is set, confirm that sign of plucker_coordinate indicate
+  // correct orientation of intersection
+  if (useOrientation && orientation * plucker_coord0 > 0) {
+    return EXIT_EARLY;
+  }
+
+  // Determine the value of the second Plucker coordinate from edge 1
+  double plucker_coord1 =
+    plucker_edge_test(vertices[1], vertices[2], raya, rayb);
+
+  // If orientation is set, confirm that sign of plucker_coordinate indicate
+  // correct orientation of intersection
+  if (useOrientation) {
+    if (orientation * plucker_coord1 > 0) {
+      return EXIT_EARLY;
+    }
+    // If the orientation is not specified, all plucker_coords must be the same
+    // sign or zero.
+  } else if ((0.0 < plucker_coord0 && 0.0 > plucker_coord1) ||
+             (0.0 > plucker_coord0 && 0.0 < plucker_coord1)) {
+    return EXIT_EARLY;
+  }
+
+  // Determine the value of the third Plucker coordinate from edge 2
+  double plucker_coord2 =
+    plucker_edge_test(vertices[2], vertices[0], raya, rayb);
+
+  // If orientation is set, confirm that sign of plucker_coordinate indicate
+  // correct orientation of intersection
+  if (useOrientation) {
+    if (orientation * plucker_coord2 > 0) {
+      return EXIT_EARLY;
+    }
+    // If the orientation is not specified, all plucker_coords must be the same
+    // sign or zero.
+  } else if ((0.0 < plucker_coord1 && 0.0 > plucker_coord2) ||
+             (0.0 > plucker_coord1 && 0.0 < plucker_coord2) ||
+             (0.0 < plucker_coord0 && 0.0 > plucker_coord2) ||
+             (0.0 > plucker_coord0 && 0.0 < plucker_coord2)) {
+    return EXIT_EARLY;
+  }
+
+  // check for coplanar case to avoid dividing by zero
+  if (0.0 == plucker_coord0 && 0.0 == plucker_coord1 && 0.0 == plucker_coord2) {
+    return EXIT_EARLY;
+  }
+
+  // get the distance to intersection
+  const double inverse_sum =
+    1.0 / (plucker_coord0 + plucker_coord1 + plucker_coord2);
+
+  const dp::vec3 intersection = dp::vec3(plucker_coord0 * inverse_sum * vertices[2] +
+                                         plucker_coord1 * inverse_sum * vertices[0] +
+                                         plucker_coord2 * inverse_sum * vertices[1]);
+
+  // To minimize numerical error, get index of largest magnitude direction.
+  int idx = 0;
+  double max_abs_dir = 0;
+  for (uint i = 0; i < 3; ++i) {
+    if (dp::abs(direction[i]) > max_abs_dir) {
+      idx = i;
+      max_abs_dir = dp::abs(direction[i]);
+    }
+  }
+
+  dist_out = (intersection[idx] - origin[idx]) / direction[idx];
+
+  // Barycentric coords check
+  double u = plucker_coord2 * inverse_sum;
+  double v = plucker_coord0 * inverse_sum;
+
+  // Barycentric coords check
+  if (u < 0.0 || v < 0.0 || (u + v) > 1.0) {
+      dist_out = -1.0;
+  }
+
+  // is the intersection within distance limits?
+  if (dist_out < tMin || dist_out > tMax) return EXIT_EARLY;
+
+  return {true, dist_out};
+}
 
 } // namespace xdg
 

include/xdg/gprt/ray_tracer.h

@@ -82,10 +82,7 @@ class GPRTRayTracer : public RayTracer {
     bool occluded(TreeID scene,
                   const Position& origin,
                   const Direction& direction,
-                  double& dist) const override {
-      fatal_error("Occlusion queries are not currently supported with GPRT ray tracer");
-      return false;
-    }
+                  double& dist) const override;
 
   private:
     // GPRT objects 
@@ -98,7 +95,8 @@ class GPRTRayTracer : public RayTracer {
     // Shader programs
     GPRTRayGenOf<dblRayGenData> rayGenProgram_; 
     GPRTRayGenOf<dblRayGenData> rayGenPointInVolProgram_;
-    GPRTMissOf<void> missProgram_; 
+    GPRTRayGenOf<dblRayGenData> rayGenOccludedProgram_;
+    GPRTMissOf<void> missProgram_;
     GPRTComputeOf<DPTriangleGeomData> aabbPopulationProgram_; //<! AABB population program for double precision rays
 
     // Buffers 
@@ -118,7 +116,7 @@ class GPRTRayTracer : public RayTracer {
     std::map<MeshID, GPRTGeomOf<DPTriangleGeomData>> surface_to_geometry_map_; //<! Map from mesh surface to embree geometry
 
     // Internal GPRT Mappings
-    std::unordered_map<SurfaceTreeID, GPRTAccel> surface_volume_tree_to_accel_map; // Map from XDG::TreeID to GPRTAccel for volume TLAS
+    std::unordered_map<SurfaceTreeID, GPRTAccel> surface_volume_tree_to_accel_map_; // Map from XDG::TreeID to GPRTAccel for volume TLAS
     std::vector<GPRTAccel> blas_handles_; // Store BLAS handles so that they can be explicitly referenced in destructor
 
     // Global Tree IDs

include/xdg/gprt/shared_structs.h

@@ -1,5 +1,6 @@
 #include "gprt.h"
 #include "../shared_enums.h"
+#include "../geometry/dp_math.h"
 
 struct GPRTPrimitiveRef
 {
@@ -26,6 +27,7 @@ struct dblRayOutput
   int surf_id;
   int primitive_id;
   xdg::PointInVolume piv; // Point in volume check result (0 for outside, 1 for inside)
+  xdg::Visibility visibility; // 0 if visible, 1 if occluded
 };
 
 /* variables for double precision triangle mesh geometry */

include/xdg/shared_enums.h

@@ -8,6 +8,11 @@ namespace xdg {
     INSIDE = 1 
   };
 
+  enum Visibility : int { 
+    VISIBLE = 0, 
+    OCCLUDED = 1 
+  };
+
   enum HitOrientation : int {
     ANY = -1,
     EXITING = 0,