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/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 6eb7ea7e..b35140a4 100644 --- a/include/xdg/geometry/plucker.h +++ b/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 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 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 09192988..465e2a82 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 "../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 */ 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/geometry/plucker.cpp b/src/geometry/plucker.cpp index 10a03f6a..67ef1d05 100644 --- a/src/geometry/plucker.cpp +++ b/src/geometry/plucker.cpp @@ -8,122 +8,6 @@ namespace xdg { -constexpr bool EXIT_EARLY = false; - -double plucker_edge_test(const Position& vertexa, const Position& vertexb, - const Position& ray, const Position& 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); - } else { - const Position edge = vertexa - vertexb; - const Position edge_normal = edge.cross(vertexb); - pip = ray.dot(edge_normal) + ray_normal.dot(edge); - pip = -pip; - } - if (PLUCKER_ZERO_TOL > fabs(pip)) - pip = 0.0; - return pip; -} - -bool plucker_ray_tri_intersect(const std::array vertices, - const Position& origin, - const Direction& direction, - double& dist_out, - const double nonneg_ray_len, - const double* neg_ray_len, - const int* orientation) -{ - dist_out = INFTY; - - const Position raya = direction; - const Position 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 Position 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 (fabs(direction[i]) > max_abs_dir) { - idx = i; - max_abs_dir = fabs(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; -} } // namespace xdg \ No newline at end of file diff --git a/src/gprt/dbl_deviceCode.slang b/src/gprt/dbl_deviceCode.slang index 2f069704..7e99ad65 100644 --- a/src/gprt/dbl_deviceCode.slang +++ b/src/gprt/dbl_deviceCode.slang @@ -1,26 +1,36 @@ #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; -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; @@ -39,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 @@ -68,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; @@ -78,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 @@ -95,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")] @@ -109,11 +159,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)); @@ -133,101 +183,42 @@ 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]; - double3 v0 = record.vertex[indices[0]]; - double3 v1 = record.vertex[indices[1]]; - double3 v2 = record.vertex[indices[2]]; - - double3 origin = record.rayIn[rayID].origin; - double3 direction = record.rayIn[rayID].direction; - - // double tMin = record.rayIn[rayID].tMin; + 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; - const double3 raya = direction; - const double3 rayb = 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; - } - - 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_coord0 && 0.0 == plucker_coord1 && 0.0 == plucker_coord2) { - return; - } - - const double inverse_sum = 1.0 / (plucker_coord0 + plucker_coord1 + plucker_coord2); - const double3 intersection = double3(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; - - - if(u < 0.0 || v < 0.0 || (u + v) > 1.0) { - t = -1.0; - } - if (t > tMax) { - return; - } - if (t < tMin) { - return; - } + bool useOrientation = false; + int orientation = 0; + xdg::PluckerIntersectionResult result = xdg::plucker_ray_tri_intersect(vertices, + origin, + direction, + tMax, + tMin, + useOrientation, + orientation); + 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); - 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 +226,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; @@ -256,33 +247,77 @@ 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 ------------------------------------------------- -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 xdg::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 +329,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; 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/src/mesh_manager_interface.cpp b/src/mesh_manager_interface.cpp index 8a558354..7725c604 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,18 @@ 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.data(), + r, + u, + INFTY, + 0.0, + true, + 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..c3bc6905 100644 --- a/src/triangle_intersect.cpp +++ b/src/triangle_intersect.cpp @@ -65,10 +65,15 @@ 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); - - if (!hit_tri) return; + auto result = plucker_ray_tri_intersect(vertices.data(), + ray_origin, + ray_direction, + rayhit->ray.dtfar, + 0.0, + false, + 0); + if (!result.hit) return; + double plucker_dist = result.t; if (plucker_dist > rayhit->ray.dtfar) return; @@ -141,8 +146,16 @@ 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.data(), + ray->dorg, + ray->ddir, + ray->dtfar, + 0.0, + false, + 0); + double plucker_dist = result.t; + + if (result.hit) { ray->set_tfar(-INFTY); } } 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