From 015fe83000ef77dae5eaf88d090e927891106ee5 Mon Sep 17 00:00:00 2001 From: waqar-ukaea Date: Fri, 3 Oct 2025 15:52:53 +0100 Subject: [PATCH 1/8] Attempting to add slang/c++ agnostic wrapper around plucker intersection tests --- include/xdg/geometry/plucker.h | 14 ++++++++---- src/geometry/plucker.cpp | 41 +++++++++++++++++----------------- src/mesh_manager_interface.cpp | 20 +++++++++-------- src/triangle_intersect.cpp | 12 +++++----- 4 files changed, 48 insertions(+), 39 deletions(-) diff --git a/include/xdg/geometry/plucker.h b/include/xdg/geometry/plucker.h index 6eb7ea7e..eef67eac 100644 --- a/include/xdg/geometry/plucker.h +++ b/include/xdg/geometry/plucker.h @@ -2,6 +2,7 @@ #define _XDG_PLUCKER_H #include "xdg/vec3da.h" +#include "xdg/geometry/dp_math.h" namespace xdg { @@ -20,10 +21,15 @@ 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 vertices, - const Position& origin, - const Direction& direction, - double& dist_out, + +struct PluckerIntersectionResult { + bool hit = false; // Whether an intersection occurred + double t = 0.0; // Distance along the ray to the intersection point +}; + +PluckerIntersectionResult plucker_ray_tri_intersect(const std::array vertices, + const dp::vec3& origin, + const dp::vec3& direction, const double nonneg_ray_len = INFTY, const double* neg_ray_len = nullptr, const int* orientation = nullptr); diff --git a/src/geometry/plucker.cpp b/src/geometry/plucker.cpp index 10a03f6a..d6cd5889 100644 --- a/src/geometry/plucker.cpp +++ b/src/geometry/plucker.cpp @@ -8,39 +8,38 @@ namespace xdg { -constexpr bool EXIT_EARLY = false; +constexpr PluckerIntersectionResult EXIT_EARLY = {false, 0.0}; -double plucker_edge_test(const Position& vertexa, const Position& vertexb, - const Position& ray, const Position& ray_normal) +double plucker_edge_test(const dp::vec3& vertexa, const dp::vec3& vertexb, + const dp::vec3& ray, const dp::vec3& ray_normal) { double pip; if (lower(vertexa, vertexb)) { - const Position edge = vertexb - vertexa; - const Position edge_normal = edge.cross(vertexa); - pip = ray.dot(edge_normal) + ray_normal.dot(edge); + 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 Position edge = vertexa - vertexb; - const Position edge_normal = edge.cross(vertexb); - pip = ray.dot(edge_normal) + ray_normal.dot(edge); + 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 (PLUCKER_ZERO_TOL > fabs(pip)) + if (PLUCKER_ZERO_TOL > dp::abs(pip)) // <-- absd pip = 0.0; return pip; } -bool plucker_ray_tri_intersect(const std::array vertices, - const Position& origin, - const Direction& direction, - double& dist_out, +PluckerIntersectionResult plucker_ray_tri_intersect(const std::array vertices, + const dp::vec3& origin, + const dp::vec3& direction, const double nonneg_ray_len, const double* neg_ray_len, const int* orientation) { - dist_out = INFTY; + double dist_out = INFTY; - const Position raya = direction; - const Position rayb = direction.cross(origin); + const dp::vec3 raya = direction; + const dp::vec3 rayb = direction.cross(origin); // Determine the value of the first Plucker coordinate from edge 0 double plucker_coord0 = @@ -98,7 +97,7 @@ bool plucker_ray_tri_intersect(const std::array vertices, 1.0 / (plucker_coord0 + plucker_coord1 + plucker_coord2); assert(0.0 != inverse_sum); - const Position intersection(plucker_coord0 * inverse_sum * vertices[2] + + const dp::vec3 intersection(plucker_coord0 * inverse_sum * vertices[2] + plucker_coord1 * inverse_sum * vertices[0] + plucker_coord2 * inverse_sum * vertices[1]); @@ -106,9 +105,9 @@ bool plucker_ray_tri_intersect(const std::array vertices, int idx = 0; double max_abs_dir = 0; for (unsigned int i = 0; i < 3; ++i) { - if (fabs(direction[i]) > max_abs_dir) { + if (dp::abs(direction[i]) > max_abs_dir) { idx = i; - max_abs_dir = fabs(direction[i]); + max_abs_dir = dp::abs(direction[i]); } } @@ -122,7 +121,7 @@ bool plucker_ray_tri_intersect(const std::array vertices, return EXIT_EARLY; } - return true; + return {true, dist_out}; } diff --git a/src/mesh_manager_interface.cpp b/src/mesh_manager_interface.cpp index 8a558354..a13cb765 100644 --- a/src/mesh_manager_interface.cpp +++ b/src/mesh_manager_interface.cpp @@ -145,7 +145,7 @@ MeshManager::next_element(MeshID current_element, const Position& r, const Position& u) const { - std::array dists = {INFTY, INFTY, INFTY, INFTY}; + std::array dists = {INFTY, INFTY, INFTY, INFTY}; std::array hit_types; auto element_face_accessor = ElementFaceAccessor::create(this, current_element); @@ -162,15 +162,17 @@ MeshManager::next_element(MeshID current_element, // with respect to the element int orientation = 1; // perform ray-triangle intersection - hit_types[i] = plucker_ray_tri_intersect(coords, - r, - u, - dists[i], - INFTY, - nullptr, - &orientation); + + auto result = plucker_ray_tri_intersect(coords, + r, + u, + INFTY, + nullptr, + &orientation); + + hit_types[i] = result.hit; // set distance and ensure it is non-negative - dists[i] = std::max(0.0, dists[i]); + dists[i] = result.hit ? std::max(0.0, result.t) : INFTY; } // determine the minimum distance to exit and the face number diff --git a/src/triangle_intersect.cpp b/src/triangle_intersect.cpp index 867592a3..e34e53ba 100644 --- a/src/triangle_intersect.cpp +++ b/src/triangle_intersect.cpp @@ -65,10 +65,10 @@ void TriangleIntersectionFunc(RTCIntersectFunctionNArguments* args) { Direction ray_direction = {ray.ddir[0], ray.ddir[1], ray.ddir[2]}; // local variable for distance to the triangle intersection - double plucker_dist; - bool hit_tri = plucker_ray_tri_intersect(vertices, ray_origin, ray_direction, plucker_dist); + auto result = plucker_ray_tri_intersect(vertices, ray_origin, ray_direction); + double plucker_dist = result.t; - if (!hit_tri) return; + if (!result.hit) return; if (plucker_dist > rayhit->ray.dtfar) return; @@ -141,8 +141,10 @@ void TriangleOcclusionFunc(RTCOccludedFunctionNArguments* args) { // get the double precision ray from the args RTCSurfaceDualRay* ray = (RTCSurfaceDualRay*) args->ray; - double plucker_dist; - if (plucker_ray_tri_intersect(vertices, ray->dorg, ray->ddir, plucker_dist)) { + auto result = plucker_ray_tri_intersect(vertices, ray->dorg, ray->ddir); + double plucker_dist = result.t; + + if (result.hit) { ray->set_tfar(-INFTY); } } From 8e4e88c930003761a20682958cbdb2c4ce3aa951 Mon Sep 17 00:00:00 2001 From: waqar-ukaea Date: Fri, 3 Oct 2025 17:42:30 +0100 Subject: [PATCH 2/8] Refactored EmbreeRayTracer to make use of language agnostic math wrapper --- include/xdg/geometry/plucker.h | 129 ++++++++++++++++++++++++++++++++- src/geometry/plucker.cpp | 115 ----------------------------- 2 files changed, 126 insertions(+), 118 deletions(-) diff --git a/include/xdg/geometry/plucker.h b/include/xdg/geometry/plucker.h index eef67eac..888085da 100644 --- a/include/xdg/geometry/plucker.h +++ b/include/xdg/geometry/plucker.h @@ -6,7 +6,6 @@ namespace xdg { - /* * Triangle vertex ordering convention: * @@ -27,12 +26,136 @@ struct PluckerIntersectionResult { double t = 0.0; // Distance along the ray to the intersection point }; -PluckerIntersectionResult plucker_ray_tri_intersect(const std::array vertices, +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(const dp::vec3& a, const 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(const dp::vec3& vertexa, const dp::vec3& vertexb, + const dp::vec3& ray, const dp::vec3& ray_normal) +{ + double pip; + if (lower(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 (PLUCKER_ZERO_TOL > dp::abs(pip)) // <-- absd + pip = 0.0; + return pip; +} + +inline PluckerIntersectionResult plucker_ray_tri_intersect(const std::array vertices, const dp::vec3& origin, const dp::vec3& direction, const double nonneg_ray_len = INFTY, const double* neg_ray_len = nullptr, - const int* orientation = nullptr); + const int* orientation = nullptr) +{ + double dist_out = INFTY; + + const dp::vec3 raya = direction; + const dp::vec3 rayb = direction.cross(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 (orientation && (*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 (orientation) { + 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 (orientation) { + 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); + assert(0.0 != inverse_sum); + + const dp::vec3 intersection(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 (unsigned int 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]; + + // is the intersection within distance limits? + if ((nonneg_ray_len && nonneg_ray_len < dist_out) || // intersection is beyond positive limit + (neg_ray_len && *neg_ray_len >= dist_out) || // intersection is behind negative limit + (!neg_ray_len && 0 > dist_out)) // unless neg_ray_len used, don't allow negative distances + { + return EXIT_EARLY; + } + + return {true, dist_out}; +} } // namespace xdg diff --git a/src/geometry/plucker.cpp b/src/geometry/plucker.cpp index d6cd5889..67ef1d05 100644 --- a/src/geometry/plucker.cpp +++ b/src/geometry/plucker.cpp @@ -8,121 +8,6 @@ namespace xdg { -constexpr PluckerIntersectionResult EXIT_EARLY = {false, 0.0}; - -double plucker_edge_test(const dp::vec3& vertexa, const dp::vec3& vertexb, - const dp::vec3& ray, const dp::vec3& ray_normal) -{ - double pip; - if (lower(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 (PLUCKER_ZERO_TOL > dp::abs(pip)) // <-- absd - pip = 0.0; - return pip; -} - -PluckerIntersectionResult plucker_ray_tri_intersect(const std::array vertices, - const dp::vec3& origin, - const dp::vec3& direction, - const double nonneg_ray_len, - const double* neg_ray_len, - const int* orientation) -{ - double dist_out = INFTY; - - const dp::vec3 raya = direction; - const dp::vec3 rayb = direction.cross(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 (orientation && (*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 (orientation) { - 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 (orientation) { - 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); - assert(0.0 != inverse_sum); - - const dp::vec3 intersection(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 (unsigned int 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]; - - // is the intersection within distance limits? - if ((nonneg_ray_len && nonneg_ray_len < dist_out) || // intersection is beyond positive limit - (neg_ray_len && *neg_ray_len >= dist_out) || // intersection is behind negative limit - (!neg_ray_len && 0 > dist_out)) // unless neg_ray_len used, don't allow negative distances - { - return EXIT_EARLY; - } - - return {true, dist_out}; -} } // namespace xdg \ No newline at end of file From 020b9f50e6a0bcd6965c2d516ad3ebbf01bb53dc Mon Sep 17 00:00:00 2001 From: waqar-ukaea Date: Mon, 6 Oct 2025 12:58:47 +0100 Subject: [PATCH 3/8] Refactored GPRTRayTracer to make use of language agnostic math wrapper --- CMakeLists.txt | 2 ++ include/xdg/gprt/shared_structs.h | 1 + src/gprt/dbl_deviceCode.slang | 52 +++++++++++++++---------------- 3 files changed, 29 insertions(+), 26 deletions(-) diff --git a/CMakeLists.txt b/CMakeLists.txt index 37b97734..9cbad767 100644 --- a/CMakeLists.txt +++ b/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 ) diff --git a/include/xdg/gprt/shared_structs.h b/include/xdg/gprt/shared_structs.h index 09192988..65af7b28 100644 --- a/include/xdg/gprt/shared_structs.h +++ b/include/xdg/gprt/shared_structs.h @@ -1,5 +1,6 @@ #include "gprt.h" #include "../shared_enums.h" +#include "xdg/geometry/dp_math.h" struct GPRTPrimitiveRef { diff --git a/src/gprt/dbl_deviceCode.slang b/src/gprt/dbl_deviceCode.slang index 2f069704..797f60ca 100644 --- a/src/gprt/dbl_deviceCode.slang +++ b/src/gprt/dbl_deviceCode.slang @@ -109,11 +109,11 @@ void populate_aabbs(uint3 DispatchThreadID: SV_DispatchThreadID, uniform DPTriangleGeomData record) { int primID = DispatchThreadID.x; int3 indices = record.index[primID]; - double3 A = record.vertex[indices[0]]; - double3 B = record.vertex[indices[1]]; - double3 C = record.vertex[indices[2]]; - double3 dpaabbmin = min(min(A, B), C); - double3 dpaabbmax = max(max(A, B), C); + dp::vec3 A = record.vertex[indices[0]]; + dp::vec3 B = record.vertex[indices[1]]; + dp::vec3 C = record.vertex[indices[2]]; + dp::vec3 dpaabbmin = min(min(A, B), C); + dp::vec3 dpaabbmax = max(max(A, B), C); float3 fpaabbmin = float3(dpaabbmin - float3(FLT_EPSILON, FLT_EPSILON, FLT_EPSILON)); float3 fpaabbmax = float3(dpaabbmax + float3(FLT_EPSILON, FLT_EPSILON, FLT_EPSILON)); @@ -151,19 +151,19 @@ void DPTrianglePluckerIntersection(uniform DPTriangleGeomData record) } int3 indices = record.index[primID]; - double3 v0 = record.vertex[indices[0]]; - double3 v1 = record.vertex[indices[1]]; - double3 v2 = record.vertex[indices[2]]; + dp::vec3 v0 = record.vertex[indices[0]]; + dp::vec3 v1 = record.vertex[indices[1]]; + dp::vec3 v2 = record.vertex[indices[2]]; - double3 origin = record.rayIn[rayID].origin; - double3 direction = record.rayIn[rayID].direction; + dp::vec3 origin = record.rayIn[rayID].origin; + dp::vec3 direction = record.rayIn[rayID].direction; // double tMin = record.rayIn[rayID].tMin; double tMin = record.rayIn[rayID].tMin; double tMax = record.rayIn[rayID].tMax; - const double3 raya = direction; - const double3 rayb = cross(direction, origin); + const dp::vec3 raya = direction; + const dp::vec3 rayb = dp::cross(direction, origin); double plucker_coord0 = plucker_edge_test(v0, v1, raya, rayb); if (useOrientation && orientation * plucker_coord0 > 0) { @@ -192,7 +192,7 @@ void DPTrianglePluckerIntersection(uniform DPTriangleGeomData record) } const double inverse_sum = 1.0 / (plucker_coord0 + plucker_coord1 + plucker_coord2); - const double3 intersection = double3(plucker_coord0 * inverse_sum * v2 + + const dp::vec3 intersection = dp::vec3(plucker_coord0 * inverse_sum * v2 + plucker_coord1 * inverse_sum * v0 + plucker_coord2 * inverse_sum * v1); @@ -227,7 +227,7 @@ void DPTrianglePluckerIntersection(uniform DPTriangleGeomData record) float f32t = float(t); if (double(f32t) < t) f32t = next_after(f32t); - double3 norm = record.normals[primID]; // recover double precision normal. TODO - Should we calculate from vertices instead? + dp::vec3 norm = record.normals[primID]; // recover double precision normal. TODO - Should we calculate from vertices instead? // sense adjustment of normal if (record.rayIn[rayID].volume_tree == record.reverse_tree) @@ -235,7 +235,7 @@ void DPTrianglePluckerIntersection(uniform DPTriangleGeomData record) norm = -norm; } - double norm_dot_dir = dot(norm, direction); + double norm_dot_dir = dp::dot(norm, direction); uint hit_kind = norm_dot_dir < 0 ? HIT_KIND_TRIANGLE_FRONT_FACE : HIT_KIND_TRIANGLE_BACK_FACE; @@ -259,30 +259,30 @@ void DPTrianglePluckerIntersection(uniform DPTriangleGeomData record) // ------------------------------------------------- Helper functions ------------------------------------------------- -bool orientation_cull(in double3 ray, in double3 normal, in xdg::HitOrientation orientation) { +bool orientation_cull(in dp::vec3 ray, in dp::vec3 normal, in HitOrientation orientation) { if (orientation == xdg::HitOrientation::ANY) return false; // No culling - if (orientation == xdg::HitOrientation::EXITING) return dot(ray, normal) < 0.0; // Cull exiting rays - if (orientation == xdg::HitOrientation::ENTERING) return dot(ray, normal) > 0.0; // Cull entering rays + if (orientation == xdg::HitOrientation::EXITING) return dp::dot(ray, normal) < 0.0; // Cull exiting rays + if (orientation == xdg::HitOrientation::ENTERING) return dp::dot(ray, normal) > 0.0; // Cull entering rays return false; // Default case, no culling } // Plucker coordinate -double plucker_edge_test(in double3 vertexa, in double3 vertexb, in double3 ray, in double3 ray_normal) +double plucker_edge_test(in dp::vec3 vertexa, in dp::vec3 vertexb, in dp::vec3 ray, in dp::vec3 ray_normal) { double pip; const double near_zero = 10 * DBL_EPSILON; if (first(vertexa, vertexb)) { - double3 edge = vertexb - vertexa; - double3 edge_normal = cross(edge, vertexa); - pip = dot(ray, edge_normal) + dot(ray_normal, edge); + dp::vec3 edge = vertexb - vertexa; + dp::vec3 edge_normal = dp::cross(edge, vertexa); + pip = dp::dot(ray, edge_normal) + dp::dot(ray_normal, edge); } else { - double3 edge = vertexa - vertexb; - double3 edge_normal = cross(edge, vertexb); - pip = dot(ray, edge_normal) + dot(ray_normal, edge); + dp::vec3 edge = vertexa - vertexb; + dp::vec3 edge_normal = dp::cross(edge, vertexb); + pip = dp::dot(ray, edge_normal) + dp::dot(ray_normal, edge); pip = -pip; } @@ -294,7 +294,7 @@ double plucker_edge_test(in double3 vertexa, in double3 vertexb, in double3 ray, 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(in double3 a, in double3 b) +inline bool first(in dp::vec3 a, in dp::vec3 b) { if (a[0] < b[0]) return true; From 0ad9d89116f209e54f603423b9afa4bb916411d8 Mon Sep 17 00:00:00 2001 From: waqar-ukaea Date: Mon, 6 Oct 2025 15:48:29 +0100 Subject: [PATCH 4/8] Refactored intersection shader to remove RT pipeline flags not used --- src/gprt/dbl_deviceCode.slang | 68 ++++++++++++----------------------- 1 file changed, 22 insertions(+), 46 deletions(-) diff --git a/src/gprt/dbl_deviceCode.slang b/src/gprt/dbl_deviceCode.slang index 797f60ca..26f56b78 100644 --- a/src/gprt/dbl_deviceCode.slang +++ b/src/gprt/dbl_deviceCode.slang @@ -133,63 +133,39 @@ void DPTrianglePluckerIntersection(uniform DPTriangleGeomData record) uint rayID = DispatchRaysIndex().x; uint nRays = DispatchRaysDimensions().x; - uint flags = RayFlags(); - if (rayID >= nRays) { - return; - } - - bool useOrientation = false; - int orientation = 0; - if ((flags & RAY_FLAG_CULL_BACK_FACING_TRIANGLES) != 0) { - orientation = -1; - useOrientation = true; - } - else if ((flags & RAY_FLAG_CULL_FRONT_FACING_TRIANGLES) != 0) { - orientation = 1; - useOrientation = true; - } + if (rayID >= nRays) return; + // Load vertices int3 indices = record.index[primID]; dp::vec3 v0 = record.vertex[indices[0]]; dp::vec3 v1 = record.vertex[indices[1]]; dp::vec3 v2 = record.vertex[indices[2]]; + dp::vec3[3] vertices = { v0, v1, v2 }; + // Load rays dp::vec3 origin = record.rayIn[rayID].origin; dp::vec3 direction = record.rayIn[rayID].direction; - - // double tMin = record.rayIn[rayID].tMin; double tMin = record.rayIn[rayID].tMin; double tMax = record.rayIn[rayID].tMax; + // ----- Perform plucker intersection test ----- + const dp::vec3 raya = direction; const dp::vec3 rayb = dp::cross(direction, origin); double plucker_coord0 = plucker_edge_test(v0, v1, raya, rayb); - if (useOrientation && orientation * plucker_coord0 > 0) { - return; - } - double plucker_coord1 = plucker_edge_test(v1, v2, raya, rayb); - if (useOrientation && orientation * plucker_coord1 > 0) { - return; - } - else if ((0.0 < plucker_coord0 && 0.0 > plucker_coord1) || (0.0 > plucker_coord0 && 0.0 < plucker_coord1)) { - return; - } + + // sign-consistency check + if ((0.0 < plucker_coord0 && 0.0 > plucker_coord1) || (0.0 > plucker_coord0 && 0.0 < plucker_coord1)) return; double plucker_coord2 = plucker_edge_test(v2, v0, raya, rayb); - if (useOrientation && orientation * plucker_coord2 > 0) { - return; - } - 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; - } + 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; - if (0.0 == plucker_coord0 && 0.0 == plucker_coord1 && 0.0 == plucker_coord2) { - return; - } + // coplanar check + if (0.0 == plucker_coord0 && 0.0 == plucker_coord1 && 0.0 == plucker_coord2) return; const double inverse_sum = 1.0 / (plucker_coord0 + plucker_coord1 + plucker_coord2); const dp::vec3 intersection = dp::vec3(plucker_coord0 * inverse_sum * v2 + @@ -200,7 +176,7 @@ void DPTrianglePluckerIntersection(uniform DPTriangleGeomData record) double max_abs_dir = 0; for(uint i = 0; i < 3; ++i) { if(abs(direction[i]) > max_abs_dir) { - idx = i; + idx = i; max_abs_dir = abs(direction[i]); } } @@ -210,16 +186,16 @@ void DPTrianglePluckerIntersection(uniform DPTriangleGeomData record) double u = plucker_coord2 * inverse_sum; double v = plucker_coord0 * inverse_sum; - - if(u < 0.0 || v < 0.0 || (u + v) > 1.0) { + // Barycentric coords check + if (u < 0.0 || v < 0.0 || (u + v) > 1.0) { t = -1.0; } - if (t > tMax) { - return; - } - if (t < tMin) { - return; - } + + // Check t range + if (t > tMax) return; + if (t < tMin) return; + + // ----- End of plucker intersection test ----- DPAttribute attr; attr.f64t = t; From 443f41debf9439f60fe5b35fe20065e464286d16 Mon Sep 17 00:00:00 2001 From: waqar-ukaea Date: Mon, 6 Oct 2025 17:41:32 +0100 Subject: [PATCH 5/8] Working plucker function on both CPU and GPU --- include/xdg/geometry/plucker.h | 64 ++++++++++++++++++---------------- src/gprt/dbl_deviceCode.slang | 61 ++++++++------------------------ src/mesh_manager_interface.cpp | 7 ++-- src/triangle_intersect.cpp | 16 +++++++-- 4 files changed, 66 insertions(+), 82 deletions(-) diff --git a/include/xdg/geometry/plucker.h b/include/xdg/geometry/plucker.h index 888085da..cae41596 100644 --- a/include/xdg/geometry/plucker.h +++ b/include/xdg/geometry/plucker.h @@ -1,7 +1,6 @@ #ifndef _XDG_PLUCKER_H #define _XDG_PLUCKER_H -#include "xdg/vec3da.h" #include "xdg/geometry/dp_math.h" namespace xdg { @@ -33,7 +32,7 @@ constexpr PluckerIntersectionResult EXIT_EARLY = {false, 0.0}; representation. This would be more simple with MOAB handles instead of coordinates... */ -inline bool first(const dp::vec3& a, const dp::vec3& b) { +inline bool first(dp::vec3 a, dp::vec3 b) { if (a[0] < b[0]) return true; if (a[0] > b[0]) return false; @@ -43,11 +42,11 @@ inline bool first(const dp::vec3& a, const dp::vec3& b) { return a[2] < b[2]; } -inline double plucker_edge_test(const dp::vec3& vertexa, const dp::vec3& vertexb, - const dp::vec3& ray, const dp::vec3& ray_normal) +inline double plucker_edge_test(dp::vec3 vertexa, dp::vec3 vertexb, + dp::vec3 ray, dp::vec3 ray_normal) { double pip; - if (lower(vertexa, vertexb)) { + 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); @@ -57,22 +56,23 @@ inline double plucker_edge_test(const dp::vec3& vertexa, const dp::vec3& vertexb pip = dp::dot(ray, edge_normal) + dp::dot(ray_normal, edge); pip = -pip; } - if (PLUCKER_ZERO_TOL > dp::abs(pip)) // <-- absd + if (dp::DBL_ZERO_TOL > dp::abs(pip)) // <-- absd pip = 0.0; return pip; } -inline PluckerIntersectionResult plucker_ray_tri_intersect(const std::array vertices, - const dp::vec3& origin, - const dp::vec3& direction, - const double nonneg_ray_len = INFTY, - const double* neg_ray_len = nullptr, - const int* orientation = nullptr) +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 = INFTY; + double dist_out = dp::INFTY; const dp::vec3 raya = direction; - const dp::vec3 rayb = direction.cross(origin); + const dp::vec3 rayb = dp::cross(direction, origin); // Determine the value of the first Plucker coordinate from edge 0 double plucker_coord0 = @@ -80,7 +80,7 @@ inline PluckerIntersectionResult plucker_ray_tri_intersect(const std::array 0) { + if (useOrientation && orientation * plucker_coord0 > 0) { return EXIT_EARLY; } @@ -90,8 +90,8 @@ inline PluckerIntersectionResult plucker_ray_tri_intersect(const std::array 0) { + if (useOrientation) { + if (orientation * plucker_coord1 > 0) { return EXIT_EARLY; } // If the orientation is not specified, all plucker_coords must be the same @@ -107,8 +107,8 @@ inline PluckerIntersectionResult plucker_ray_tri_intersect(const std::array 0) { + if (useOrientation) { + if (orientation * plucker_coord2 > 0) { return EXIT_EARLY; } // If the orientation is not specified, all plucker_coords must be the same @@ -128,16 +128,15 @@ inline PluckerIntersectionResult plucker_ray_tri_intersect(const std::array max_abs_dir) { idx = i; max_abs_dir = dp::abs(direction[i]); @@ -146,14 +145,19 @@ inline PluckerIntersectionResult plucker_ray_tri_intersect(const std::array= dist_out) || // intersection is behind negative limit - (!neg_ray_len && 0 > dist_out)) // unless neg_ray_len used, don't allow negative distances - { - return EXIT_EARLY; + // 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}; } diff --git a/src/gprt/dbl_deviceCode.slang b/src/gprt/dbl_deviceCode.slang index 26f56b78..12c56402 100644 --- a/src/gprt/dbl_deviceCode.slang +++ b/src/gprt/dbl_deviceCode.slang @@ -1,4 +1,5 @@ -#include "../../include/xdg/gprt/shared_structs.h" +#include "xdg/gprt/shared_structs.h" +#include "xdg/geometry/plucker.h" [[vk::push_constant]] dblRayFirePushConstants PC; @@ -141,7 +142,7 @@ void DPTrianglePluckerIntersection(uniform DPTriangleGeomData record) dp::vec3 v0 = record.vertex[indices[0]]; dp::vec3 v1 = record.vertex[indices[1]]; dp::vec3 v2 = record.vertex[indices[2]]; - dp::vec3[3] vertices = { v0, v1, v2 }; + dp::vec3 vertices[3] = { v0, v1, v2 }; // Load rays dp::vec3 origin = record.rayIn[rayID].origin; @@ -149,54 +150,20 @@ void DPTrianglePluckerIntersection(uniform DPTriangleGeomData record) double tMin = record.rayIn[rayID].tMin; double tMax = record.rayIn[rayID].tMax; - // ----- Perform plucker intersection test ----- + bool useOrientation = false; + int orientation = 0; - const dp::vec3 raya = direction; - const dp::vec3 rayb = dp::cross(direction, origin); + xdg::PluckerIntersectionResult result = xdg::plucker_ray_tri_intersect(vertices, + origin, + direction, + tMax, + tMin, + useOrientation, + orientation); + if (result.hit == false) return; // No intersection - double plucker_coord0 = plucker_edge_test(v0, v1, raya, rayb); - double plucker_coord1 = plucker_edge_test(v1, v2, raya, rayb); + double t = result.t; - // sign-consistency check - if ((0.0 < plucker_coord0 && 0.0 > plucker_coord1) || (0.0 > plucker_coord0 && 0.0 < plucker_coord1)) return; - - double plucker_coord2 = plucker_edge_test(v2, v0, raya, rayb); - 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; - - // coplanar check - if (0.0 == plucker_coord0 && 0.0 == plucker_coord1 && 0.0 == plucker_coord2) return; - - const double inverse_sum = 1.0 / (plucker_coord0 + plucker_coord1 + plucker_coord2); - const dp::vec3 intersection = dp::vec3(plucker_coord0 * inverse_sum * v2 + - plucker_coord1 * inverse_sum * v0 + - plucker_coord2 * inverse_sum * v1); - - int idx = 0; - double max_abs_dir = 0; - for(uint i = 0; i < 3; ++i) { - if(abs(direction[i]) > max_abs_dir) { - idx = i; - max_abs_dir = abs(direction[i]); - } - } - const double dist = (intersection[idx] - origin[idx]) / direction[idx]; - - double t = dist; - 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) { - t = -1.0; - } - - // Check t range - if (t > tMax) return; - if (t < tMin) return; - - // ----- End of plucker intersection test ----- - DPAttribute attr; attr.f64t = t; diff --git a/src/mesh_manager_interface.cpp b/src/mesh_manager_interface.cpp index a13cb765..7725c604 100644 --- a/src/mesh_manager_interface.cpp +++ b/src/mesh_manager_interface.cpp @@ -163,12 +163,13 @@ MeshManager::next_element(MeshID current_element, int orientation = 1; // perform ray-triangle intersection - auto result = plucker_ray_tri_intersect(coords, + auto result = plucker_ray_tri_intersect(coords.data(), r, u, INFTY, - nullptr, - &orientation); + 0.0, + true, + orientation); hit_types[i] = result.hit; // set distance and ensure it is non-negative diff --git a/src/triangle_intersect.cpp b/src/triangle_intersect.cpp index e34e53ba..79b5acfe 100644 --- a/src/triangle_intersect.cpp +++ b/src/triangle_intersect.cpp @@ -65,7 +65,13 @@ void TriangleIntersectionFunc(RTCIntersectFunctionNArguments* args) { Direction ray_direction = {ray.ddir[0], ray.ddir[1], ray.ddir[2]}; // local variable for distance to the triangle intersection - auto result = plucker_ray_tri_intersect(vertices, ray_origin, ray_direction); + auto result = plucker_ray_tri_intersect(vertices.data(), + ray_origin, + ray_direction, + rayhit->ray.dtfar, + 0.0, + false, + 0); double plucker_dist = result.t; if (!result.hit) return; @@ -141,7 +147,13 @@ void TriangleOcclusionFunc(RTCOccludedFunctionNArguments* args) { // get the double precision ray from the args RTCSurfaceDualRay* ray = (RTCSurfaceDualRay*) args->ray; - auto result = plucker_ray_tri_intersect(vertices, ray->dorg, ray->ddir); + auto result = plucker_ray_tri_intersect(vertices.data(), + ray->dorg, + ray->ddir, + ray->dtfar, + 0.0, + false, + 0); double plucker_dist = result.t; if (result.hit) { From 7edc25b91e58dea1ef2121de84a94cc9ce504606 Mon Sep 17 00:00:00 2001 From: waqar-ukaea Date: Tue, 7 Oct 2025 13:44:18 +0100 Subject: [PATCH 6/8] Fixed empty lines --- include/xdg/geometry/plucker.h | 1 - src/gprt/dbl_deviceCode.slang | 1 - 2 files changed, 2 deletions(-) diff --git a/include/xdg/geometry/plucker.h b/include/xdg/geometry/plucker.h index cae41596..725d4f43 100644 --- a/include/xdg/geometry/plucker.h +++ b/include/xdg/geometry/plucker.h @@ -157,7 +157,6 @@ inline PluckerIntersectionResult plucker_ray_tri_intersect(dp::vec3 vertices[3], // is the intersection within distance limits? if (dist_out < tMin || dist_out > tMax) return EXIT_EARLY; - return {true, dist_out}; } diff --git a/src/gprt/dbl_deviceCode.slang b/src/gprt/dbl_deviceCode.slang index 12c56402..2a5266cb 100644 --- a/src/gprt/dbl_deviceCode.slang +++ b/src/gprt/dbl_deviceCode.slang @@ -161,7 +161,6 @@ void DPTrianglePluckerIntersection(uniform DPTriangleGeomData record) useOrientation, orientation); if (result.hit == false) return; // No intersection - double t = result.t; DPAttribute attr; From 8e954e0cc99352282d0974131660eefea00e6361 Mon Sep 17 00:00:00 2001 From: waqar-ukaea Date: Mon, 13 Oct 2025 16:06:55 +0100 Subject: [PATCH 7/8] Reverted back to relative include paths for slang code --- include/xdg/geometry/dp_math.h | 43 +++++++++++++++++++++++++++++++ include/xdg/geometry/plucker.h | 2 +- include/xdg/gprt/shared_structs.h | 2 +- src/gprt/dbl_deviceCode.slang | 10 ++++++- src/triangle_intersect.cpp | 3 +-- 5 files changed, 55 insertions(+), 5 deletions(-) create mode 100644 include/xdg/geometry/dp_math.h diff --git a/include/xdg/geometry/dp_math.h b/include/xdg/geometry/dp_math.h new file mode 100644 index 00000000..5b162d3f --- /dev/null +++ b/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::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::epsilon(); + constexpr double INFTY {std::numeric_limits::max()}; + +} + +#endif // end of ifdef __slang__ + +#endif // DP_MATH_H \ No newline at end of file diff --git a/include/xdg/geometry/plucker.h b/include/xdg/geometry/plucker.h index 725d4f43..b35140a4 100644 --- a/include/xdg/geometry/plucker.h +++ b/include/xdg/geometry/plucker.h @@ -1,7 +1,7 @@ #ifndef _XDG_PLUCKER_H #define _XDG_PLUCKER_H -#include "xdg/geometry/dp_math.h" +#include "dp_math.h" namespace xdg { diff --git a/include/xdg/gprt/shared_structs.h b/include/xdg/gprt/shared_structs.h index 65af7b28..cf87964e 100644 --- a/include/xdg/gprt/shared_structs.h +++ b/include/xdg/gprt/shared_structs.h @@ -1,6 +1,6 @@ #include "gprt.h" #include "../shared_enums.h" -#include "xdg/geometry/dp_math.h" +#include "../geometry/dp_math.h" struct GPRTPrimitiveRef { diff --git a/src/gprt/dbl_deviceCode.slang b/src/gprt/dbl_deviceCode.slang index 2a5266cb..23d3943d 100644 --- a/src/gprt/dbl_deviceCode.slang +++ b/src/gprt/dbl_deviceCode.slang @@ -1,5 +1,13 @@ +#include "../../include/xdg/gprt/shared_structs.h" +#include "../../include/xdg/geometry/plucker.h" + +/* +For now we have to use relative paths for includes, which is not ideal. If https://github.com/gprt-org/GPRT/pull/82 gets +merged into GPRT we will be able to more robustly include these headers in the manner below: + #include "xdg/gprt/shared_structs.h" #include "xdg/geometry/plucker.h" +*/ [[vk::push_constant]] dblRayFirePushConstants PC; @@ -201,7 +209,7 @@ void DPTrianglePluckerIntersection(uniform DPTriangleGeomData record) // ------------------------------------------------- Helper functions ------------------------------------------------- -bool orientation_cull(in dp::vec3 ray, in dp::vec3 normal, in HitOrientation orientation) { +bool orientation_cull(in dp::vec3 ray, in dp::vec3 normal, in xdg::HitOrientation orientation) { if (orientation == xdg::HitOrientation::ANY) return false; // No culling if (orientation == xdg::HitOrientation::EXITING) return dp::dot(ray, normal) < 0.0; // Cull exiting rays if (orientation == xdg::HitOrientation::ENTERING) return dp::dot(ray, normal) > 0.0; // Cull entering rays diff --git a/src/triangle_intersect.cpp b/src/triangle_intersect.cpp index 79b5acfe..c3bc6905 100644 --- a/src/triangle_intersect.cpp +++ b/src/triangle_intersect.cpp @@ -72,9 +72,8 @@ void TriangleIntersectionFunc(RTCIntersectFunctionNArguments* args) { 0.0, false, 0); - double plucker_dist = result.t; - if (!result.hit) return; + double plucker_dist = result.t; if (plucker_dist > rayhit->ray.dtfar) return; From 08455090f00305c3ff0c025cce63fe91d038eb07 Mon Sep 17 00:00:00 2001 From: waqar-ukaea Date: Wed, 29 Oct 2025 17:59:15 +0000 Subject: [PATCH 8/8] Implemented GPRTRayTracer::occluded() --- include/xdg/gprt/ray_tracer.h | 10 ++- include/xdg/gprt/shared_structs.h | 1 + include/xdg/shared_enums.h | 5 ++ src/gprt/dbl_deviceCode.slang | 109 ++++++++++++++++++++++++++---- src/gprt/ray_tracer.cpp | 57 ++++++++++++++-- tests/test_occluded.cpp | 51 +++++++++----- 6 files changed, 190 insertions(+), 43 deletions(-) diff --git a/include/xdg/gprt/ray_tracer.h b/include/xdg/gprt/ray_tracer.h index 0805af27..592e766c 100644 --- a/include/xdg/gprt/ray_tracer.h +++ b/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 rayGenProgram_; GPRTRayGenOf rayGenPointInVolProgram_; - GPRTMissOf missProgram_; + GPRTRayGenOf rayGenOccludedProgram_; + GPRTMissOf missProgram_; GPRTComputeOf aabbPopulationProgram_; //> surface_to_geometry_map_; // surface_volume_tree_to_accel_map; // Map from XDG::TreeID to GPRTAccel for volume TLAS + std::unordered_map surface_volume_tree_to_accel_map_; // Map from XDG::TreeID to GPRTAccel for volume TLAS std::vector blas_handles_; // Store BLAS handles so that they can be explicitly referenced in destructor // Global Tree IDs diff --git a/include/xdg/gprt/shared_structs.h b/include/xdg/gprt/shared_structs.h index cf87964e..465e2a82 100644 --- a/include/xdg/gprt/shared_structs.h +++ b/include/xdg/gprt/shared_structs.h @@ -27,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 */ diff --git a/include/xdg/shared_enums.h b/include/xdg/shared_enums.h index f6198f60..7aa91445 100644 --- a/include/xdg/shared_enums.h +++ b/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, diff --git a/src/gprt/dbl_deviceCode.slang b/src/gprt/dbl_deviceCode.slang index 23d3943d..7e99ad65 100644 --- a/src/gprt/dbl_deviceCode.slang +++ b/src/gprt/dbl_deviceCode.slang @@ -12,24 +12,25 @@ merged into GPRT we will be able to more robustly include these headers in the m [[vk::push_constant]] dblRayFirePushConstants PC; -struct RayFirePayload { +struct Payload { double distance; // Distance to intersection int surf_id; // ID of the surface hit SurfaceAccelerationStructure tlas; int primitive_id; // ID of the primitive hit - xdg::PointInVolume piv; // Point in volume check (0 for outside, 1 for inside) -}; - + xdg::PointInVolume piv; // Point in volume check (0 for outside, 1 for inside) + xdg::Visibility visibility; // 1 if visible, 1 if occluded +}; struct DPAttribute { double f64t; // double precision hit distance int global_prim_id; }; +// ------------------------------------------------- Hit Shaders ------------------------------------------------- [shader("closesthit")] -void ray_fire_hit(uniform DPTriangleGeomData record, inout RayFirePayload payload, in DPAttribute attr) { +void ray_fire_hit(uniform DPTriangleGeomData record, inout Payload payload, in DPAttribute attr) { // Distance from the ray origin to the hit point uint hit_kind = HitKind(); uint rayID = DispatchRaysIndex().x; @@ -48,19 +49,31 @@ void ray_fire_hit(uniform DPTriangleGeomData record, inout RayFirePayload payloa payload.primitive_id = attr.global_prim_id; } +[shader("anyhit")] +void occlusion_anyhit(uniform DPTriangleGeomData record, inout Payload payload) { + payload.visibility = xdg::Visibility::OCCLUDED; + AcceptHitAndEndSearch(); +} + +// ------------------------------------------------- Miss Shaders ------------------------------------------------- + +// GPRT doesn't currently support setting multiple miss shaders so only a single miss shader is possible in a single GPRT context + [shader("miss")] -void ray_fire_miss(inout RayFirePayload payload) { +void miss(inout Payload payload) { // Set the miss payload to default values payload.distance = -1.0f; payload.surf_id = -1; payload.primitive_id = -1; + + payload.visibility = xdg::Visibility::VISIBLE; } -// This ray generation program will kick off the ray tracing process, -// generating rays and tracing them into the world. +// ------------------------------------------------- RayGen Shaders ------------------------------------------------- + [shader("raygeneration")] void ray_fire(uniform dblRayGenData record, uniform DPTriangleGeomData mesh) { - RayFirePayload payload; + Payload payload; uint rayID = DispatchRaysIndex().x; // Trace the ray into the scene @@ -77,7 +90,9 @@ void ray_fire(uniform dblRayGenData record, uniform DPTriangleGeomData mesh) { payload.surf_id = -1; payload.tlas = world; - TraceRay(world, RAY_FLAG_NONE, 0xff, 0, 1, rayDesc, payload); + int missIndex = 0; // must be 0 as we only have a single miss shader + int hitGroupIndex = 0; // 0 index used for ray hitgroup [ray_fire] + [point_in_volume] + TraceRay(world, RAY_FLAG_NONE, 0xff, hitGroupIndex, missIndex, rayDesc, payload); // Store the distance to the hit point and the surface ID in buffers for CPU record.out[rayID].distance = payload.distance; @@ -87,7 +102,7 @@ void ray_fire(uniform dblRayGenData record, uniform DPTriangleGeomData mesh) { [shader("raygeneration")] void point_in_volume(uniform dblRayGenData record, uniform DPTriangleGeomData mesh) { - RayFirePayload payload; + Payload payload; uint rayID = DispatchRaysIndex().x; // Trace the ray into the scene @@ -104,12 +119,38 @@ void point_in_volume(uniform dblRayGenData record, uniform DPTriangleGeomData me payload.tlas = world; payload.piv = xdg::PointInVolume::OUTSIDE; // Initialize point in volume check result to outside (0) - TraceRay(world, RAY_FLAG_NONE, 0xff, 0, 1, rayDesc, payload); + int missIndex = 0; // must be 0 as we only have a single miss shader + int hitGroupIndex = 0; // 0 index used for ray hitgroup [ray_fire] + [point_in_volume] + + TraceRay(world, RAY_FLAG_NONE, 0xff, hitGroupIndex, missIndex, rayDesc, payload); record.out.surf_id = payload.surf_id; record.out[rayID].piv = payload.piv; // Point in volume check result } +[shader("raygeneration")] +void occluded(uniform dblRayGenData record, uniform DPTriangleGeomData mesh) { + Payload payload; + uint rayID = DispatchRaysIndex().x; + + // Trace the ray into the scene + RayDesc rayDesc; + rayDesc.Origin = float3(record.ray[rayID].origin); + rayDesc.Direction = normalize(float3(record.ray[rayID].direction)); + rayDesc.TMin = float(record.ray[rayID].tMin); + rayDesc.TMax = float(record.ray[rayID].tMax); + + SurfaceAccelerationStructure world = record.ray[rayID].volume_accel; + payload.visibility = xdg::Visibility::VISIBLE; // Initialize visibility to visible (0) + + int missIndex = 0; // must be 0 as we only have a single miss shader + int hitGroupIndex = 1; // 1 index used for ray hitgroup [occluded] + + TraceRay(world, RAY_FLAG_NONE, 0xff, hitGroupIndex, missIndex, rayDesc, payload); + + record.out[rayID].visibility = payload.visibility; +} + // ------------------------------------------------- Compute Shaders ------------------------------------------------- /* A shader to compute and store AABB min/maxes in single precision using double precision coords*/ [shader("compute")] @@ -206,6 +247,50 @@ void DPTrianglePluckerIntersection(uniform DPTriangleGeomData record) ReportHit(f32t, hit_kind, attr); } +// DP plucker intersection shader for occlusion queries +[shader("intersection")] +void DPTrianglePluckerIntersection_occluded(uniform DPTriangleGeomData record) +{ + int primID = PrimitiveIndex(); + int global_prim_id = record.primitive_refs[primID].id; + + uint rayID = DispatchRaysIndex().x; + uint nRays = DispatchRaysDimensions().x; + + if (rayID >= nRays) return; + + // Load vertices + int3 indices = record.index[primID]; + dp::vec3 v0 = record.vertex[indices[0]]; + dp::vec3 v1 = record.vertex[indices[1]]; + dp::vec3 v2 = record.vertex[indices[2]]; + dp::vec3 vertices[3] = { v0, v1, v2 }; + + // Load rays + dp::vec3 origin = record.rayIn[rayID].origin; + dp::vec3 direction = record.rayIn[rayID].direction; + double tMin = record.rayIn[rayID].tMin; + double tMax = record.rayIn[rayID].tMax; + + xdg::PluckerIntersectionResult result = xdg::plucker_ray_tri_intersect(vertices, + origin, + direction, + tMax, + tMin, + false, + 0); + if (result.hit == false) return; // No intersection + double t = result.t; + + DPAttribute attr; + attr.f64t = t; + + float f32t = float(t); + if (double(f32t) < t) f32t = next_after(f32t); + + // hit kind, distance, attribute not actually needed. But we still need to report a hit to end the search + ReportHit(f32t, HIT_KIND_TRIANGLE_FRONT_FACE, attr); +} // ------------------------------------------------- Helper functions ------------------------------------------------- diff --git a/src/gprt/ray_tracer.cpp b/src/gprt/ray_tracer.cpp index c5778f97..c5022bfc 100644 --- a/src/gprt/ray_tracer.cpp +++ b/src/gprt/ray_tracer.cpp @@ -5,6 +5,7 @@ namespace xdg { GPRTRayTracer::GPRTRayTracer() { + numRayTypes_ = 2; // ray_fire/point_in_volume + occluded gprtRequestRayTypeCount(numRayTypes_); // Set the number of shaders which can be set to the same geometry context_ = gprtContextCreate(); module_ = gprtModuleCreate(context_, dbl_deviceCode); @@ -26,6 +27,11 @@ GPRTRayTracer::GPRTRayTracer() rayGenPIVData->ray = gprtBufferGetDevicePointer(rayInputBuffer_); rayGenPIVData->out = gprtBufferGetDevicePointer(rayOutputBuffer_); + + dblRayGenData* rayGenOccludedData = gprtRayGenGetParameters(rayGenOccludedProgram_); + rayGenOccludedData->ray = gprtBufferGetDevicePointer(rayInputBuffer_); + rayGenOccludedData->out = gprtBufferGetDevicePointer(rayOutputBuffer_); + // Set up build parameters for acceleration structures buildParams_.buildMode = GPRT_BUILD_MODE_FAST_BUILD_NO_UPDATE; } @@ -38,7 +44,7 @@ GPRTRayTracer::~GPRTRayTracer() // Destroy TLAS structures - for (const auto& [tree, accel] : surface_volume_tree_to_accel_map) { + for (const auto& [tree, accel] : surface_volume_tree_to_accel_map_) { gprtAccelDestroy(accel); } @@ -68,13 +74,18 @@ void GPRTRayTracer::setup_shaders() // Set up ray generation and miss programs rayGenProgram_ = gprtRayGenCreate(context_, module_, "ray_fire"); rayGenPointInVolProgram_ = gprtRayGenCreate(context_, module_, "point_in_volume"); - missProgram_ = gprtMissCreate(context_, module_, "ray_fire_miss"); + rayGenOccludedProgram_ = gprtRayGenCreate(context_, module_, "occluded"); + missProgram_ = gprtMissCreate(context_, module_, "miss"); aabbPopulationProgram_ = gprtComputeCreate(context_, module_, "populate_aabbs"); - // Create a "triangle" geometry type and set its closest-hit program trianglesGeomType_ = gprtGeomTypeCreate(context_, GPRT_AABBS); + // Set its closest-hit program and intersection program for ray-fire and point-in-volume queries gprtGeomTypeSetClosestHitProg(trianglesGeomType_, 0, module_, "ray_fire_hit"); // closesthit for ray queries gprtGeomTypeSetIntersectionProg(trianglesGeomType_, 0, module_, "DPTrianglePluckerIntersection"); // set intersection program for double precision rays + + // Set any-hit program and intersection program for occlusion queries + gprtGeomTypeSetAnyHitProg(trianglesGeomType_, 1, module_, "occlusion_anyhit"); // anyhit for occlusion queries + gprtGeomTypeSetIntersectionProg(trianglesGeomType_, 1, module_, "DPTrianglePluckerIntersection_occluded"); // set intersection program for double precision rays for occlusion queries } void GPRTRayTracer::init() @@ -196,7 +207,7 @@ GPRTRayTracer::create_surface_tree(const std::shared_ptr& mesh_mana auto instanceBuffer = gprtDeviceBufferCreate(context_, surfaceBlasInstances.size(), surfaceBlasInstances.data()); GPRTAccel volume_tlas = gprtInstanceAccelCreate(context_, surfaceBlasInstances.size(), instanceBuffer); gprtAccelBuild(context_, volume_tlas, buildParams_); - surface_volume_tree_to_accel_map[tree] = volume_tlas; + surface_volume_tree_to_accel_map_[tree] = volume_tlas; return tree; } @@ -213,7 +224,7 @@ bool GPRTRayTracer::point_in_volume(SurfaceTreeID tree, const Direction* direction, const std::vector* exclude_primitives) const { - GPRTAccel volume = surface_volume_tree_to_accel_map.at(tree); + GPRTAccel volume = surface_volume_tree_to_accel_map_.at(tree); dblRayGenData* rayGenPIVData = gprtRayGenGetParameters(rayGenPointInVolProgram_); // Use provided direction or if Direction == nulptr use default direction @@ -273,7 +284,7 @@ std::pair GPRTRayTracer::ray_fire(SurfaceTreeID tree, HitOrientation orientation, std::vector* const exclude_primitives) { - GPRTAccel volume = surface_volume_tree_to_accel_map.at(tree); + GPRTAccel volume = surface_volume_tree_to_accel_map_.at(tree); dblRayGenData* rayGenData = gprtRayGenGetParameters(rayGenProgram_); gprtBufferMap(rayInputBuffer_); // Update the ray input buffer @@ -320,6 +331,38 @@ std::pair GPRTRayTracer::ray_fire(SurfaceTreeID tree, return {distance, surface}; } + +bool GPRTRayTracer::occluded(SurfaceTreeID tree, + const Position& origin, + const Direction& direction, + double& distance) const +{ + GPRTAccel volume = surface_volume_tree_to_accel_map_.at(tree); + dblRayGenData* rayGenOccludedData = gprtRayGenGetParameters(rayGenOccludedProgram_); + + gprtBufferMap(rayInputBuffer_); // Update the ray input buffer + dblRayInput* rayInput = gprtBufferGetHostPointer(rayInputBuffer_); + rayInput[0].volume_accel = gprtAccelGetDeviceAddress(volume); + rayInput[0].origin = {origin.x, origin.y, origin.z}; + rayInput[0].direction = {direction.x, direction.y, direction.z}; + rayInput[0].tMax = INFTY; // Set a large distance limit + rayInput[0].tMin = 0.0; + rayInput[0].volume_tree = tree; // Set the TreeID of the volume being queried + rayInput[0].hitOrientation = HitOrientation::ANY; // No orientation culling for occlusion check + gprtBufferUnmap(rayInputBuffer_); // required to sync buffer back on GPU? + + gprtRayGenLaunch1D(context_, rayGenOccludedProgram_, 1); // Launch raygen shader (entry point to RT pipeline) + gprtGraphicsSynchronize(context_); // Ensure all GPU operations are complete before returning control flow to CPU + + // Retrieve the output from the ray output buffer + gprtBufferMap(rayOutputBuffer_); + dblRayOutput* rayOutput = gprtBufferGetHostPointer(rayOutputBuffer_); + auto visibility = rayOutput->visibility; + gprtBufferUnmap(rayOutputBuffer_); // required to sync buffer back on GPU? Maybe this second unmap isn't actually needed since we dont need to resyncrhonize after retrieving the data from device + + return visibility; +} + void GPRTRayTracer::create_global_surface_tree() { // Create a TLAS (Top-Level Acceleration Structure) for all the volumes @@ -329,7 +372,7 @@ void GPRTRayTracer::create_global_surface_tree() SurfaceTreeID tree = next_surface_tree_id(); surface_trees_.push_back(tree); - surface_volume_tree_to_accel_map[tree] = global_accel; + surface_volume_tree_to_accel_map_[tree] = global_accel; global_surface_tree_ = tree; global_surface_accel_ = global_accel; } diff --git a/tests/test_occluded.cpp b/tests/test_occluded.cpp index 4201ced9..9d328dcb 100644 --- a/tests/test_occluded.cpp +++ b/tests/test_occluded.cpp @@ -1,11 +1,12 @@ - // for testing #include +#include // xdg includes #include "xdg/mesh_manager_interface.h" #include "xdg/xdg.h" -#include "xdg/embree/ray_tracer.h" + +#include "util.h" #include "mesh_mock.h" using namespace xdg; @@ -14,20 +15,34 @@ TEST_CASE("Test Occluded") { std::shared_ptr mm = std::make_shared(); mm->init(); // this should do nothing, just good practice to call it - std::shared_ptr xdg = std::make_shared(mm); - xdg->prepare_raytracer(); - auto rti = xdg->ray_tracing_interface(); - auto [volume_tree, element_tree] = rti->register_volume(mm, mm->volumes()[0]); - - // setup ray to fire that won't hit the mock model - Position r {-100.0, 0.0, 0.0}; - Direction u {1.0, 0.0, 0.0}; - double dist {0.0}; - - bool result = rti->occluded(volume_tree, r, u, dist); - REQUIRE(result == true); - - u = {-1.0, 0.0, 0.0}; - result = rti->occluded(volume_tree, r, u, dist); - REQUIRE(result == false); + // Generate one test run per enabled backend + auto rt_backend = GENERATE(RTLibrary::EMBREE, RTLibrary::GPRT); + + DYNAMIC_SECTION(fmt::format("Backend = {}", rt_backend)) { + check_ray_tracer_supported(rt_backend); // skip if backend not enabled at configuration time + auto rti = create_raytracer(rt_backend); + REQUIRE(rti); + + auto mm = std::make_shared(false); + mm->init(); + REQUIRE(mm->mesh_library() == MeshLibrary::MOCK); + + auto [volume_tree, element_tree] = rti->register_volume(mm, mm->volumes()[0]); + REQUIRE(volume_tree != ID_NONE); + REQUIRE(element_tree == ID_NONE); + + rti->init(); // Ensure ray tracer is initialized (e.g. build SBT for GPRT) + + // setup ray to fire that won't hit the mock model + Position r {-100.0, 0.0, 0.0}; + Direction u {1.0, 0.0, 0.0}; + double dist {0.0}; + + bool result = rti->occluded(volume_tree, r, u, dist); + REQUIRE(result == true); + + u = {-1.0, 0.0, 0.0}; + result = rti->occluded(volume_tree, r, u, dist); + REQUIRE(result == false); + } } \ No newline at end of file