diff --git a/CMakeLists.txt b/CMakeLists.txt index 8d129141..a9f1afc1 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -149,7 +149,14 @@ src/gprt/ray_tracer.cpp list(APPEND xdg_device_codes dbl_deviceCode ) - +list(APPEND xdg_slang_headers +src/gprt/triangle_rt_shaders.slang +src/gprt/tetrahedron_rt_shaders.slang +include/xdg/gprt/shared_structs.h +include/xdg/gprt/rt_common.slangh +include/xdg/shared_enums.h +include/xdg/geometry/dp_math.h +) endif() if (XDG_ENABLE_LIBMESH) @@ -267,9 +274,7 @@ if (XDG_ENABLE_GPRT) OUTPUT_TARGET ${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 + ${CMAKE_CURRENT_SOURCE_DIR}/${xdg_slang_headers} 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 index 5b162d3f..a2fece8f 100644 --- a/include/xdg/geometry/dp_math.h +++ b/include/xdg/geometry/dp_math.h @@ -8,18 +8,19 @@ functions for dot product, cross product, and absolute value. In C++ compilation while in Slang compilation, it maps to `double3`. */ -#ifdef __SLANG__ +#if defined(__SLANG__) || defined(__SLANG_COMPILER__) // Slang compilation, map dp::vec3 -> double3 namespace dp { typedef double3 vec3; + static const double DBL_EPS = 2.2204460492503131e-016; + static const double PLUCKER_ZERO_TOL = 20.0 * DBL_EPS; + static const double INFTY = 1.7976931348623157e+308; + 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 @@ -29,15 +30,14 @@ namespace dp { namespace dp { typedef xdg::Vec3da vec3; + static constexpr double PLUCKER_ZERO_TOL = 20.0 * std::numeric_limits::epsilon(); + constexpr double INFTY {std::numeric_limits::max()}; + 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 +#endif // DP_MATH_H diff --git a/include/xdg/geometry/plucker.h b/include/xdg/geometry/plucker.h index b35140a4..200bb4db 100644 --- a/include/xdg/geometry/plucker.h +++ b/include/xdg/geometry/plucker.h @@ -25,7 +25,7 @@ struct PluckerIntersectionResult { double t = 0.0; // Distance along the ray to the intersection point }; -constexpr PluckerIntersectionResult EXIT_EARLY = {false, 0.0}; +static const 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 @@ -56,7 +56,7 @@ inline double plucker_edge_test(dp::vec3 vertexa, dp::vec3 vertexb, pip = dp::dot(ray, edge_normal) + dp::dot(ray_normal, edge); pip = -pip; } - if (dp::DBL_ZERO_TOL > dp::abs(pip)) // <-- absd + if (dp::PLUCKER_ZERO_TOL > dp::abs(pip)) // <-- absd pip = 0.0; return pip; } @@ -160,6 +160,37 @@ inline PluckerIntersectionResult plucker_ray_tri_intersect(dp::vec3 vertices[3], return {true, dist_out}; } +// Plücker containment test for tetrahedra. Returns true if the point is inside the tetrahedron defined by vertices v0, v1, v2, v3. +inline bool plucker_tet_containment_test(const dp::vec3 point, + const dp::vec3 v0, + const dp::vec3 v1, + const dp::vec3 v2, + const dp::vec3 v3) { + + // TODO - I've decided to use Cramer's rule here instead of matrix inversion to make it easier to implement a cross-compilable version of this function + // Not sure if that is necessarily the best choice however. + const dp::vec3 e0 = v1 - v0; + const dp::vec3 e1 = v2 - v0; + const dp::vec3 e2 = v3 - v0; + const dp::vec3 rhs = point - v0; + const double det = dp::dot(e0, dp::cross(e1, e2)); // scalar triple product of matrix [e0 e1 e2] + + const double inv_det = 1.0 / det; + const double lambda1 = dp::dot(rhs, dp::cross(e1, e2)) * inv_det; + const double lambda2 = dp::dot(e0, dp::cross(rhs, e2)) * inv_det; + const double lambda3 = dp::dot(e0, dp::cross(e1, rhs)) * inv_det; + const double lambda0 = 1.0 - (lambda1 + lambda2 + lambda3); + + const double barycentric_min = -dp::PLUCKER_ZERO_TOL; + const double barycentric_max = 1.0 + dp::PLUCKER_ZERO_TOL; + + // Check all λ_i in [0, 1] + return (lambda0 >= barycentric_min && lambda0 <= barycentric_max) && + (lambda1 >= barycentric_min && lambda1 <= barycentric_max) && + (lambda2 >= barycentric_min && lambda2 <= barycentric_max) && + (lambda3 >= barycentric_min && lambda3 <= barycentric_max); +} + } // namespace xdg -#endif // include guard \ No newline at end of file +#endif // include guard diff --git a/include/xdg/gprt/ray_tracer.h b/include/xdg/gprt/ray_tracer.h index 8d24d107..13ad8956 100644 --- a/include/xdg/gprt/ray_tracer.h +++ b/include/xdg/gprt/ray_tracer.h @@ -22,7 +22,8 @@ enum class RayGenType { RAY_FIRE, POINT_IN_VOLUME, OCCLUDED, - CLOSEST + CLOSEST, + FIND_ELEMENT }; struct gprtRayHit { @@ -49,16 +50,9 @@ class GPRTRayTracer : public RayTracer { // Setup the different shader programs for use with this ray tracer void setup_shaders(); - MeshID find_element(const Position& point) const override - { - fatal_error("Element trees not currently supported with GPRT ray tracer"); - return ID_NONE; - }; + MeshID find_element(const Position& point) const override; - MeshID find_element(TreeID tree, const Position& point) const override { - fatal_error("Element trees not currently supported with GPRT ray tracer"); - return ID_NONE; - }; + MeshID find_element(TreeID tree, const Position& point) const override; std::pair register_volume(const std::shared_ptr& mesh_manager, @@ -114,9 +108,11 @@ class GPRTRayTracer : public RayTracer { // Shader programs std::map> rayGenPrograms_; - GPRTMissOf missProgram_; - GPRTComputeOf aabbPopulationProgram_; // triangleMissProgram_; + GPRTMissOf tetMissProgram_; + GPRTComputeOf aabbTriPopulationProgram_; // aabbTetPopulationProgram_; // excludePrimitivesBuffer_; // globalBlasInstances_; // trianglesGeomType_; // tetrahedraGeomType_; //> surface_to_geometry_map_; // surface_volume_tree_to_accel_map; // Map from XDG::TreeID to GPRTAccel for volume TLAS + std::unordered_map element_volume_tree_to_accel_map; // Map from XDG::TreeID to GPRTAccel for element TLAS + std::vector blas_handles_; // Store BLAS handles so that they can be explicitly referenced in destructor // Global Tree IDs @@ -143,4 +142,4 @@ class GPRTRayTracer : public RayTracer { } // namespace xdg -#endif // include guard \ No newline at end of file +#endif // include guard diff --git a/include/xdg/gprt/rt_common.slangh b/include/xdg/gprt/rt_common.slangh new file mode 100644 index 00000000..cddec57d --- /dev/null +++ b/include/xdg/gprt/rt_common.slangh @@ -0,0 +1,79 @@ +#ifndef XDG_GPRT_RT_COMMON_SLANGH +#define XDG_GPRT_RT_COMMON_SLANGH + +#include "shared_structs.h" +#include "../geometry/dp_math.h" // dp math shared between C++ and slang +#include "../geometry/plucker.h" // Plücker ray-edge intersection test + +// Shared types for both Triangle and Tetrahedron ray tracing pipelines. +[[vk::push_constant]] +dblRayFirePushConstants PC; + +static const int32_t ID_NONE = -1; // slang copy of the C++ constant for invalid IDs +struct DPAttribute +{ + double f64t; // double precision hit distance + int global_prim_id; +}; + +// ------------------------------------------------- Helper functions ------------------------------------------------- + +// Templated/generic function to compute AABBs in double precision and fill the float AABB buffer. Generic to handle triangles, tets and future primitives with different numbers of vertices +__generic +inline void populate_aabb(uint primID, double3 *vertex, float3 *aabb, uint indices[N]) +{ + double3 dpAabbMin = vertex[indices[0]]; + double3 dpAabbMax = dpAabbMin; + + // unroll the loop since N is known to be small(ish) at compile time to keep performance + [unroll] + for (uint i = 1; i < N; ++i) + { + double3 vert = vertex[indices[i]]; + dpAabbMin = min(dpAabbMin, vert); + dpAabbMax = max(dpAabbMax, vert); + } + + float3 fpAabbMin = float3(dpAabbMin - double3(FLT_EPSILON, FLT_EPSILON, FLT_EPSILON)); + float3 fpAabbMax = float3(dpAabbMax + double3(FLT_EPSILON, FLT_EPSILON, FLT_EPSILON)); + + aabb[2 * primID] = fpAabbMin; + aabb[2 * primID + 1] = fpAabbMax; +} + + +bool orientation_cull(in double3 ray, in double3 normal, in xdg::HitOrientation orientation) { + if (orientation == xdg::HitOrientation::ANY) return false; // No culling + double dot_product = dot(ray, normal); + if (orientation == xdg::HitOrientation::EXITING) return dot_product < 0.0; // Cull exiting rays + if (orientation == xdg::HitOrientation::ENTERING) return dot_product > 0.0; // Cull entering rays + return false; // Default case, no culling +} + + +/* 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(in double3 a, in double3 b) +{ + if (a[0] < b[0]) return true; + + if (a[0] == b[0] && a[1] < b[1]) return true; + + if (a[1] == b[1] && a[2] < b[2]) return true; + + return false; +} + +float next_after(float a) { + uint a_ = asuint(a); + if (a < 0) { + a_--; + } else { + a_++; + } + return asfloat(a_); +} + +#endif // XDG_GPRT_RT_COMMON_SLANGH diff --git a/include/xdg/gprt/shared_structs.h b/include/xdg/gprt/shared_structs.h index caa9d553..3ec4cba4 100644 --- a/include/xdg/gprt/shared_structs.h +++ b/include/xdg/gprt/shared_structs.h @@ -2,6 +2,14 @@ #include "../shared_enums.h" #include "../geometry/dp_math.h" + +namespace xdg { + static const uint RT_SURFACE_RAY_INDEX = 0; // SBT ray index for surface ray queries + static const uint RT_VOLUME_RAY_INDEX = 1; // SBT ray index for volumetric ray queries + static const uint RT_NUM_RAY_TYPES = 2; // total number of ray types (surface and volume) + static const uint RT_SURFACE_MISS_INDEX = 0; // SBT ray index for surface ray miss shader + static const uint RT_VOLUME_MISS_INDEX = 1; // SBT ray index for +} struct GPRTPrimitiveRef { int id; // ID of the primitive @@ -18,7 +26,8 @@ struct dblRay int32_t exclude_count; // Number of excluded primitives xdg::HitOrientation hitOrientation; int volume_tree; // TreeID of the volume being queried - SurfaceAccelerationStructure volume_accel; // The volume accel + SurfaceAccelerationStructure volume_accel_surf; // The volume accel for surface acceleration structure + SolidAccelerationStructure volume_accel_solid; // The volume accel for solid acceleration structure }; struct dblHit @@ -47,6 +56,17 @@ struct DPTriangleGeomData { int num_faces; // Number of faces in the geometry }; +struct DPTetrahedronGeomData { + double3 *vertex; // vertex buffer + float3 *aabbs; // AABB buffer + uint4 *index; // index buffer + int32_t vol_id; + dblRay *ray; // double precision rays + xdg::HitOrientation hitOrientation; + GPRTPrimitiveRef* primitive_refs; + int num_tets; // Number of tetrahedra in the geometry +}; + struct dblRayGenData { dblRay *ray; dblHit *hit; diff --git a/include/xdg/shared_enums.h b/include/xdg/shared_enums.h index f6198f60..a20c27fa 100644 --- a/include/xdg/shared_enums.h +++ b/include/xdg/shared_enums.h @@ -4,6 +4,7 @@ namespace xdg { enum PointInVolume : int { + UNSET = -1, OUTSIDE = 0, INSIDE = 1 }; diff --git a/include/xdg/tetrahedron_contain.h b/include/xdg/tetrahedron_contain.h index 51f1a338..03250447 100644 --- a/include/xdg/tetrahedron_contain.h +++ b/include/xdg/tetrahedron_contain.h @@ -3,6 +3,7 @@ #include "xdg/vec3da.h" +#include "xdg/geometry/plucker.h" namespace xdg { @@ -30,11 +31,6 @@ namespace xdg * @return `true` if the point is inside or on the boundary of the tetrahedron, * `false` otherwise. */ -bool plucker_tet_containment_test(const Position& point, - const Vertex& v0, - const Vertex& v1, - const Vertex& v2, - const Vertex& v3); // Embree call back functions for element search void VolumeElementBoundsFunc(RTCBoundsFunctionArguments* args); diff --git a/src/gprt/dbl_deviceCode.slang b/src/gprt/dbl_deviceCode.slang index 4e643fd9..49323e98 100644 --- a/src/gprt/dbl_deviceCode.slang +++ b/src/gprt/dbl_deviceCode.slang @@ -1,266 +1,7 @@ -#include "../../include/xdg/gprt/shared_structs.h" -#include "../../include/xdg/geometry/plucker.h" +#ifndef DBL_DEVICE_CODE_SLANGH +#define DBL_DEVICE_CODE_SLANGH -/* -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 "triangle_rt_shaders.slang" +#include "tetrahedron_rt_shaders.slang" -#include "xdg/gprt/shared_structs.h" -#include "xdg/geometry/plucker.h" -*/ - -[[vk::push_constant]] -dblRayFirePushConstants PC; - -struct RayFirePayload { - 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) -}; - - -struct DPAttribute -{ - double f64t; // double precision hit distance - int global_prim_id; -}; - - -[shader("closesthit")] -void ray_fire_hit(uniform DPTriangleGeomData record, inout RayFirePayload payload, in DPAttribute attr) { - // Distance from the ray origin to the hit point - uint hit_kind = HitKind(); - uint rayID = DispatchRaysIndex().x; - - // There is some logic for handling next volumes inside the h5m-reader which I could make use of too - // TODO : Should the dblHit struct return the next volume ID for the ray back to the host - - payload.piv = (hit_kind == HIT_KIND_TRIANGLE_FRONT_FACE) - ? xdg::PointInVolume::OUTSIDE - : xdg::PointInVolume::INSIDE; - - int instanceID = InstanceID(); - - payload.distance = attr.f64t; - payload.surf_id = record.surf_id; - payload.primitive_id = attr.global_prim_id; -} - -[shader("miss")] -void ray_fire_miss(inout RayFirePayload payload) { - // Set the miss payload to default values - payload.distance = -1.0f; - payload.surf_id = -1; - payload.primitive_id = -1; -} - -// This ray generation program will kick off the ray tracing process, -// generating rays and tracing them into the world. -[shader("raygeneration")] -void ray_fire(uniform dblRayGenData record, uniform DPTriangleGeomData mesh) { - RayFirePayload 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; - - // Pass the ray's origin and direction to the payload - payload.distance = -1.0f; - payload.surf_id = -1; - payload.tlas = world; - - TraceRay(world, RAY_FLAG_NONE, 0xff, 0, 1, rayDesc, payload); - - // Store the distance to the hit point and the surface ID in buffers for CPU - record.hit[rayID].distance = payload.distance; - record.hit[rayID].surf_id = payload.surf_id; - record.hit[rayID].primitive_id = payload.primitive_id; -} - -[shader("raygeneration")] -void point_in_volume(uniform dblRayGenData record, uniform DPTriangleGeomData mesh) { - RayFirePayload payload; - uint rayID = DispatchRaysIndex().x; - - // Trace the ray into the scene - RayDesc rayDesc; - rayDesc.Origin = float3(record.ray[rayID].origin); - rayDesc.Direction = float3(normalize(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; - - // Pass the ray's origin and direction to the payload - payload.surf_id = -1; - 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); - - record.hit[rayID].surf_id = payload.surf_id; - record.hit[rayID].piv = payload.piv; // Point in volume check result -} - -// ------------------------------------------------- Compute Shaders ------------------------------------------------- -/* A shader to compute and store AABB min/maxes in single precision using double precision coords*/ -[shader("compute")] -[numthreads(1, 1, 1)] -void -populate_aabbs(uint3 DispatchThreadID: SV_DispatchThreadID, uniform DPTriangleGeomData record) { - int primID = DispatchThreadID.x; - int3 indices = record.index[primID]; - 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)); - - record.aabbs[2 * primID] = fpaabbmin; - record.aabbs[2 * primID + 1] = fpaabbmax; -} - -// ------------------------------------------------ CUSTOM INTERSECTION SHADERS ------------------------------------------------ - - -/* 1D ray generation intersection with a double precision triangle using the Plucker intersection algorithm*/ -[shader("intersection")] -void DPTrianglePluckerIntersection(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.ray[rayID].origin; - dp::vec3 direction = record.ray[rayID].direction; - double tMin = record.ray[rayID].tMin; - double tMax = record.ray[rayID].tMax; - - 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); - - dp::vec3 norm = record.normals[primID]; // recover double precision normal. TODO - Should we calculate from vertices instead? - - // sense adjustment of normal - if (record.ray[rayID].volume_tree == record.reverse_tree) - { - norm = -norm; - } - - 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; - - xdg::HitOrientation hitOrientation = record.ray[rayID].hitOrientation; - - if (orientation_cull(direction, norm, hitOrientation)) - { - return; - } - - for (int i = 0; i < record.ray[rayID].exclude_count; ++i) - { - if (record.ray[rayID].exclude_primitives[i] == global_prim_id) { - return; - } - } - attr.global_prim_id = global_prim_id; - ReportHit(f32t, hit_kind, attr); -} - - -// ------------------------------------------------- Helper functions ------------------------------------------------- - -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 - return false; // Default case, no culling -} - -// Plucker coordinate -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)) - { - 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 - { - 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; - } - - if (near_zero > abs(pip)) pip = 0.0; - return pip; -} - -/* 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(in dp::vec3 a, in dp::vec3 b) -{ - if (a[0] < b[0]) return true; - - if (a[0] == b[0] && a[1] < b[1]) return true; - - if (a[1] == b[1] && a[2] < b[2]) return true; - - return false; -} - -float next_after(float a) { - uint a_ = asuint(a); - if (a < 0) { - a_--; - } else { - a_++; - } - return asfloat(a_); -} \ No newline at end of file +#endif // DBL_DEVICE_CODE_SLANGH \ No newline at end of file diff --git a/src/gprt/ray_tracer.cpp b/src/gprt/ray_tracer.cpp index 5ad08e72..6b6db475 100644 --- a/src/gprt/ray_tracer.cpp +++ b/src/gprt/ray_tracer.cpp @@ -5,7 +5,7 @@ namespace xdg { GPRTRayTracer::GPRTRayTracer() { - gprtRequestRayTypeCount(numRayTypes_); // Set the number of shaders which can be set to the same geometry + gprtRequestRayTypeCount(RT_NUM_RAY_TYPES); context_ = gprtContextCreate(); module_ = gprtModuleCreate(context_, dbl_deviceCode); @@ -30,6 +30,10 @@ GPRTRayTracer::GPRTRayTracer() rayGenPIVData->ray = gprtBufferGetDevicePointer(rayHitBuffers_.ray); rayGenPIVData->hit = gprtBufferGetDevicePointer(rayHitBuffers_.hit); + dblRayGenData* rayGenFindElementData = gprtRayGenGetParameters(rayGenPrograms_.at(RayGenType::FIND_ELEMENT)); + rayGenFindElementData->ray = gprtBufferGetDevicePointer(rayHitBuffers_.ray); + rayGenFindElementData->hit = gprtBufferGetDevicePointer(rayHitBuffers_.hit); + // Set up build parameters for acceleration structures buildParams_.buildMode = GPRT_BUILD_MODE_FAST_BUILD_NO_UPDATE; } @@ -41,11 +45,16 @@ GPRTRayTracer::~GPRTRayTracer() gprtComputeSynchronize(context_); - // Destroy TLAS structures + // Destroy TLAS structures for surface trees for (const auto& [tree, accel] : surface_volume_tree_to_accel_map) { gprtAccelDestroy(accel); } + // Destroy TLAS structures for element trees + for (const auto& [tree, accel] : element_volume_tree_to_accel_map) { + gprtAccelDestroy(accel); + } + // Destroy BLAS structures for (const auto& blas : blas_handles_) { gprtAccelDestroy(blas); @@ -56,6 +65,7 @@ GPRTRayTracer::~GPRTRayTracer() gprtGeomDestroy(geom); } gprtGeomTypeDestroy(trianglesGeomType_); + gprtGeomTypeDestroy(tetrahedraGeomType_); // Destroy Buffers gprtBufferDestroy(rayHitBuffers_.ray); @@ -72,15 +82,24 @@ void GPRTRayTracer::setup_shaders() // Set up ray generation and miss programs rayGenPrograms_[RayGenType::RAY_FIRE] = gprtRayGenCreate(context_, module_, "ray_fire"); rayGenPrograms_[RayGenType::POINT_IN_VOLUME] = gprtRayGenCreate(context_, module_, "point_in_volume"); + rayGenPrograms_[RayGenType::FIND_ELEMENT] = gprtRayGenCreate(context_, module_, "find_element"); // TODO: Add Occluded and closest raygen entry points - missProgram_ = gprtMissCreate(context_, module_, "ray_fire_miss"); - aabbPopulationProgram_ = gprtComputeCreate(context_, module_, "populate_aabbs"); + triangleMissProgram_ = gprtMissCreate(context_, module_, "ray_fire_miss"); + tetMissProgram_ = gprtMissCreate(context_, module_, "tet_miss"); + + aabbTriPopulationProgram_ = gprtComputeCreate(context_, module_, "populate_tri_aabbs"); + aabbTetPopulationProgram_ = gprtComputeCreate(context_, module_, "populate_tet_aabbs"); // Create a "triangle" geometry type and set its closest-hit program trianglesGeomType_ = gprtGeomTypeCreate(context_, GPRT_AABBS); - gprtGeomTypeSetClosestHitProg(trianglesGeomType_, 0, module_, "ray_fire_hit"); // closesthit for ray queries - gprtGeomTypeSetIntersectionProg(trianglesGeomType_, 0, module_, "DPTrianglePluckerIntersection"); // set intersection program for double precision rays + gprtGeomTypeSetClosestHitProg(trianglesGeomType_, RT_SURFACE_RAY_INDEX, module_, "ray_fire_hit"); // closesthit for ray queries + gprtGeomTypeSetIntersectionProg(trianglesGeomType_, RT_SURFACE_RAY_INDEX, module_, "DPTrianglePluckerIntersection"); // set intersection program for double precision rays against triangles + + // Create a "tetrahedron" geometry type and set its closest-hit program + tetrahedraGeomType_ = gprtGeomTypeCreate(context_, GPRT_AABBS); + gprtGeomTypeSetClosestHitProg(tetrahedraGeomType_, RT_VOLUME_RAY_INDEX, module_, "tet_contain_hit"); // closesthit for point-in-volume queries + gprtGeomTypeSetIntersectionProg(tetrahedraGeomType_, RT_VOLUME_RAY_INDEX, module_, "DPTetrahedronPluckerIntersection"); // set intersection program for double precision rays against tetrahedra } void GPRTRayTracer::init() @@ -109,45 +128,39 @@ GPRTRayTracer::create_surface_tree(const std::shared_ptr& mesh_mana std::vector surfaceBlasInstances; // BLAS for each (surface) geometry in this volume for (const auto &surf : volume_surfaces) { - auto num_faces = mesh_manager->num_surface_faces(surf); - - // get the sense of this surface with respect to the volume - Sense triangle_sense {Sense::UNSET}; - auto surf_to_vol_senses = mesh_manager->get_parent_volumes(surf); - if (volume_id == surf_to_vol_senses.first) triangle_sense = Sense::FORWARD; - else if (volume_id == surf_to_vol_senses.second) triangle_sense = Sense::REVERSE; - DPTriangleGeomData* geom_data = nullptr; auto triangleGeom = gprtGeomCreate(context_, trianglesGeomType_); geom_data = gprtGeomGetParameters(triangleGeom); // pointer to assign data to - // Get storage for vertices + // Get storage for vertices and indices auto vertices = mesh_manager->get_surface_vertices(surf); auto indices = mesh_manager->get_surface_connectivity(surf); - std::vector dbl3Vertices; - dbl3Vertices.reserve(vertices.size()); - for (const auto &vertex : vertices) { - dbl3Vertices.push_back({vertex.x, vertex.y, vertex.z}); + if (indices.size() % 3 != 0) { + fatal_error("Surface {} connectivity size ({}) is not divisible by 3 for triangles", surf, indices.size()); + } + std::vector dbl3Vertices(vertices.size()); + for (size_t i = 0; i < vertices.size(); ++i) { + const auto& vertex = vertices[i]; + dbl3Vertices[i] = {vertex.x, vertex.y, vertex.z}; } // Get storage for indices - std::vector ui3Indices; - ui3Indices.reserve(indices.size() / 3); - for (size_t i = 0; i < indices.size(); i += 3) { - ui3Indices.emplace_back(indices[i], indices[i + 1], indices[i + 2]); + std::vector ui3Indices(indices.size() / 3); + for (size_t i = 0, triIdx = 0; i < indices.size(); i += 3, ++triIdx) { + ui3Indices[triIdx] = uint3(indices[i], indices[i + 1], indices[i + 2]); } // Get storage for normals - std::vector normals; - std::vector primitive_refs; - primitive_refs.reserve(num_faces); - normals.reserve(num_faces); - for (const auto &face : mesh_manager->get_surface_faces(surf)) { + auto surface_faces = mesh_manager->get_surface_faces(surf); + const size_t num_faces = surface_faces.size(); + + std::vector normals(surface_faces.size()); + std::vector primitive_refs(surface_faces.size()); + for (size_t i = 0; i < surface_faces.size(); ++i) { + const auto face = surface_faces[i]; auto norm = mesh_manager->face_normal(face); - normals.push_back({norm.x, norm.y, norm.z}); - GPRTPrimitiveRef prim_ref; - prim_ref.id = face; - primitive_refs.push_back(prim_ref); + normals[i] = {norm.x, norm.y, norm.z}; + primitive_refs[i].id = face; } auto vertex_buffer = gprtDeviceBufferCreate(context_, dbl3Vertices.size(), dbl3Vertices.data()); @@ -166,22 +179,18 @@ GPRTRayTracer::create_surface_tree(const std::shared_ptr& mesh_mana geom_data->primitive_refs = gprtBufferGetDevicePointer(primitive_refs_buffer); geom_data->num_faces = num_faces; - gprtComputeLaunch(aabbPopulationProgram_, {num_faces, 1, 1}, {1, 1, 1}, *geom_data); + gprtComputeLaunch(aabbTriPopulationProgram_, {num_faces, 1, 1}, {1, 1, 1}, *geom_data); GPRTAccel blas = gprtAABBAccelCreate(context_, triangleGeom, buildParams_.buildMode); gprtAccelBuild(context_, blas, buildParams_); - gprt::Instance instance; - instance = gprtAccelGetInstance(blas); // create instance of BLAS to be added to TLAS + gprt::Instance instance = gprtAccelGetInstance(blas); // create instance of BLAS to be added to TLAS instance.mask = 0xff; // mask can be used to filter instances during ray traversal. 0xff ensures no filtering // Store in maps surface_to_geometry_map_[surf] = triangleGeom; - geom_data = gprtGeomGetParameters(triangleGeom); - instance = gprtAccelGetInstance(blas); - instance.mask = 0xff; surfaceBlasInstances.push_back(instance); globalBlasInstances_.push_back(instance); @@ -210,8 +219,65 @@ GPRTRayTracer::create_surface_tree(const std::shared_ptr& mesh_mana ElementTreeID GPRTRayTracer::create_element_tree(const std::shared_ptr& mesh_manager, MeshID volume_id) { - warning("Element trees not currently supported with GPRT ray tracer"); - return TREE_NONE; + auto volume_elements = mesh_manager->get_volume_elements(volume_id); + if (volume_elements.empty()) return TREE_NONE; // No elements in this volume, so no tree to create + + ElementTreeID tree = next_element_tree_id(); + element_trees_.push_back(tree); + + DPTetrahedronGeomData* geom_data = nullptr; + auto tetrahedraGeom = gprtGeomCreate(context_, tetrahedraGeomType_); + geom_data = gprtGeomGetParameters(tetrahedraGeom); // pointer to assign data to + + auto vertices = mesh_manager->get_volume_vertices(volume_id); + auto indices = mesh_manager->get_volume_connectivity(volume_id); + + std::vector dbl3Vertices(vertices.size()); + for (size_t i = 0; i < vertices.size(); ++i) { + const auto& vertex = vertices[i]; + dbl3Vertices[i] = {vertex.x, vertex.y, vertex.z}; + } + + // Get storage for indices + std::vector ui4Indices(indices.size() / 4); + for (size_t i = 0, tetIdx = 0; i < indices.size(); i += 4, ++tetIdx) { + ui4Indices[tetIdx] = uint4(indices[i], indices[i + 1], indices[i + 2], indices[i + 3]); + } + + // Get storage for prim IDs + std::vector primitive_refs(volume_elements.size()); + for (size_t i = 0; i < volume_elements.size(); ++i) { + primitive_refs[i].id = volume_elements[i]; + } + + auto vertex_buffer = gprtDeviceBufferCreate(context_, dbl3Vertices.size(), dbl3Vertices.data()); + auto connectivity_buffer = gprtDeviceBufferCreate(context_, ui4Indices.size(), ui4Indices.data()); + auto primitive_refs_buffer = gprtDeviceBufferCreate(context_, primitive_refs.size(), primitive_refs.data()); + auto aabb_buffer = gprtDeviceBufferCreate(context_, 2*volume_elements.size(), 0); // AABBs for each tetrahedron + gprtAABBsSetPositions(tetrahedraGeom, aabb_buffer, volume_elements.size(), 2*sizeof(float3), 0); + + geom_data->aabbs = gprtBufferGetDevicePointer(aabb_buffer); + geom_data->vertex = gprtBufferGetDevicePointer(vertex_buffer); + geom_data->index = gprtBufferGetDevicePointer(connectivity_buffer); + geom_data->num_tets = volume_elements.size(); + geom_data->vol_id = volume_id; + geom_data->ray = gprtBufferGetDevicePointer(rayHitBuffers_.ray); + geom_data->primitive_refs = gprtBufferGetDevicePointer(primitive_refs_buffer); + + gprtComputeLaunch(aabbTetPopulationProgram_, {volume_elements.size(), 1, 1}, {1, 1, 1}, *geom_data); + + GPRTAccel blas = gprtAABBAccelCreate(context_, tetrahedraGeom, buildParams_.buildMode); + gprtAccelBuild(context_, blas, buildParams_); + gprt::Instance instance = gprtAccelGetInstance(blas); // create instance of BLAS to be added to TLAS + instance.mask = 0xff; // mask can be used to filter instances during ray traversal. 0xff ensures no filtering + + auto instanceBuffer = gprtDeviceBufferCreate(context_, 1, &instance); + GPRTAccel volume_tlas = gprtInstanceAccelCreate(context_, 1, instanceBuffer); + gprtAccelBuild(context_, volume_tlas, buildParams_); + + element_volume_tree_to_accel_map[tree] = volume_tlas; + + return tree; }; bool GPRTRayTracer::point_in_volume(SurfaceTreeID tree, @@ -229,7 +295,7 @@ bool GPRTRayTracer::point_in_volume(SurfaceTreeID tree, gprtBufferMap(rayHitBuffers_.ray); // Update the ray input buffer dblRay* ray = gprtBufferGetHostPointer(rayHitBuffers_.ray); - ray[0].volume_accel = gprtAccelGetDeviceAddress(volume); + ray[0].volume_accel_surf = gprtAccelGetDeviceAddress(volume); ray[0].origin = {point.x, point.y, point.z}; ray[0].direction = {directionUsed.x, directionUsed.y, directionUsed.z}; ray[0].tMax = INFTY; // Set a large distance limit @@ -286,7 +352,7 @@ std::pair GPRTRayTracer::ray_fire(SurfaceTreeID tree, gprtBufferMap(rayHitBuffers_.ray); // Update the ray input buffer dblRay* ray = gprtBufferGetHostPointer(rayHitBuffers_.ray); - ray[0].volume_accel = gprtAccelGetDeviceAddress(volume); + ray[0].volume_accel_surf = gprtAccelGetDeviceAddress(volume); ray[0].origin = {origin.x, origin.y, origin.z}; ray[0].direction = {direction.x, direction.y, direction.z}; ray[0].tMax = dist_limit; @@ -369,4 +435,36 @@ void GPRTRayTracer::check_ray_buffer_capacity(size_t N) gprtBuildShaderBindingTable(context_, static_cast(GPRT_SBT_GEOM | GPRT_SBT_RAYGEN)); } +MeshID GPRTRayTracer::find_element(const Position& point) const +{ + return find_element(global_element_tree_, point); +} + +MeshID GPRTRayTracer::find_element(TreeID tree, const Position& point) const +{ + + GPRTAccel volume = element_volume_tree_to_accel_map.at(tree); + auto rayGen = rayGenPrograms_.at(RayGenType::FIND_ELEMENT); + dblRayGenData* rayGenData = gprtRayGenGetParameters(rayGen); + + gprtBufferMap(rayHitBuffers_.ray); // Update the ray input buffer + dblRay* ray = gprtBufferGetHostPointer(rayHitBuffers_.ray); + ray[0].volume_accel_solid = gprtAccelGetDeviceAddress(volume); + ray[0].origin = {point.x, point.y, point.z}; + ray[0].volume_tree = tree; // Set the TreeID of the volume being queried + + gprtBufferUnmap(rayHitBuffers_.ray); // required to sync buffer back on GPU? + + gprtRayGenLaunch1D(context_, rayGen, 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 hit from the dblHit buffer + gprtBufferMap(rayHitBuffers_.hit); + dblHit* hit = gprtBufferGetHostPointer(rayHitBuffers_.hit); + auto primitive_id = hit[0].primitive_id; + gprtBufferUnmap(rayHitBuffers_.hit); // 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 primitive_id; +} + } // namespace xdg diff --git a/src/gprt/tetrahedron_rt_shaders.slang b/src/gprt/tetrahedron_rt_shaders.slang new file mode 100644 index 00000000..47df4479 --- /dev/null +++ b/src/gprt/tetrahedron_rt_shaders.slang @@ -0,0 +1,87 @@ +#ifndef XDG_GPRT_TETRAHEDRON_RT_SHADERS_SLANG +#define XDG_GPRT_TETRAHEDRON_RT_SHADERS_SLANG + +#include "../../include/xdg/gprt/rt_common.slangh" + +struct SolidRayFirePayload { + int primitive_id; // ID of the primitive the point is contained within +}; + + +// ------------------------------------------------- Compute Shaders ------------------------------------------------- +/* A shader to compute and store AABB min/maxes in single precision using double precision coords for Tetrahedra*/ +[shader("compute")] +[numthreads(1, 1, 1)] +void +populate_tet_aabbs(uint3 DispatchThreadID: SV_DispatchThreadID, uniform DPTetrahedronGeomData record) { + uint primID = DispatchThreadID.x; + uint4 indices = record.index[primID]; + uint vertexIndices[4] = { indices.x, indices.y, indices.z, indices.w }; + populate_aabb<4>(primID, record.vertex, record.aabbs, vertexIndices); +} + +[shader("raygeneration")] +void find_element(uniform dblRayGenData record) { + // This shader is only used to find the element containing a point, so we can ignore the raygen data and just use the ray's origin as the point to test for containment + SolidRayFirePayload payload; + uint rayID = DispatchRaysIndex().x; + PointDesc pointDesc; + pointDesc.Origin = float3(record.ray[rayID].origin); + + SolidAccelerationStructure world = record.ray[rayID].volume_accel_solid; + + // Pass the ray's origin and direction to the payload + payload.primitive_id = ID_NONE; // Initialize to -1 to indicate the point is + + TracePoint(world, RAY_FLAG_NONE, 0xff, xdg::RT_VOLUME_RAY_INDEX, xdg::RT_VOLUME_MISS_INDEX, pointDesc, payload); + + // Store the ID of the containing element in the hit buffer for CPU readback + record.hit[rayID].primitive_id = payload.primitive_id; // we're reusing dblHit's primitive_id field to store the ID of the containing element since we don't need hit distance or surf_id for this shader +} + +[shader("miss")] +void tet_miss(inout SolidRayFirePayload payload) { + payload.primitive_id = ID_NONE; // Mark the ray as outside all tetrahedra +} + +[shader("closesthit")] +void tet_contain_hit(uniform DPTetrahedronGeomData record, inout SolidRayFirePayload payload, in DPAttribute attr) { + payload.primitive_id = attr.global_prim_id; // Store the ID of the tet that contains the point +} + + +/* 1D ray generation intersection with a double precision Tetrahedron using the Plucker intersection algorithm*/ +[shader("intersection")] +void DPTetrahedronPluckerIntersection(uniform DPTetrahedronGeomData record) +{ + int primID = PrimitiveIndex(); + int global_prim_id = record.primitive_refs[primID].id; + + uint rayID = DispatchRaysIndex().x; + uint nRays = DispatchRaysDimensions().x; + uint flags = RayFlags(); + + if (rayID >= nRays) { + return; + } + + int4 indices = record.index[primID]; + double3 v0 = record.vertex[indices[0]]; + double3 v1 = record.vertex[indices[1]]; + double3 v2 = record.vertex[indices[2]]; + double3 v3 = record.vertex[indices[3]]; + double3 origin = record.ray[rayID].origin; + + bool inside = xdg::plucker_tet_containment_test(origin, v0, v1, v2, v3); + if (!inside) return; + + // If we are inside the tet we need to report a hit to move through the rest of the RT pipeline + float f32t = 0.0f; // We dont care about intersection distance + uint hit_kind = 0; // No hit kind since it doesn't matter + + DPAttribute attr; + attr.global_prim_id = global_prim_id; + ReportHit(f32t, hit_kind, attr); +} + +#endif // XDG_GPRT_TETRAHEDRON_RT_SHADERS_SLANG diff --git a/src/gprt/triangle_rt_shaders.slang b/src/gprt/triangle_rt_shaders.slang new file mode 100644 index 00000000..f8e63ab3 --- /dev/null +++ b/src/gprt/triangle_rt_shaders.slang @@ -0,0 +1,181 @@ +#ifndef XDG_GPRT_TRIANGLE_RT_SHADERS_SLANG +#define XDG_GPRT_TRIANGLE_RT_SHADERS_SLANG + +#include "../../include/xdg/gprt/rt_common.slangh" + +struct SurfaceRayFirePayload { + 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) +}; + +[shader("closesthit")] +void ray_fire_hit(uniform DPTriangleGeomData record, inout SurfaceRayFirePayload payload, in DPAttribute attr) { + // Distance from the ray origin to the hit point + uint hit_kind = HitKind(); + uint rayID = DispatchRaysIndex().x; + + payload.piv = (hit_kind == HIT_KIND_TRIANGLE_FRONT_FACE) + ? xdg::PointInVolume::OUTSIDE + : xdg::PointInVolume::INSIDE; + + int instanceID = InstanceID(); + + payload.distance = attr.f64t; + payload.surf_id = record.surf_id; + payload.primitive_id = attr.global_prim_id; +} + +[shader("miss")] +void ray_fire_miss(inout SurfaceRayFirePayload payload) { + // Set the miss payload to default values + payload.distance = -1.0f; + payload.surf_id = ID_NONE; + payload.primitive_id = ID_NONE; + payload.piv = xdg::PointInVolume::UNSET; +} + +// This ray generation program will kick off the ray tracing process, +// generating rays and tracing them into the world. +[shader("raygeneration")] +void ray_fire(uniform dblRayGenData record) { + SurfaceRayFirePayload 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_surf; + + // Pass the ray's origin and direction to the payload + payload.distance = -1.0f; + payload.surf_id = ID_NONE; + payload.tlas = world; + + TraceRay(world, RAY_FLAG_NONE, 0xff, xdg::RT_SURFACE_RAY_INDEX, xdg::RT_SURFACE_MISS_INDEX, rayDesc, payload); + + // Store the distance to the hit point and the surface ID in buffers for CPU + record.hit[rayID].distance = payload.distance; + record.hit[rayID].surf_id = payload.surf_id; + record.hit[rayID].primitive_id = payload.primitive_id; +} + +[shader("raygeneration")] +void point_in_volume(uniform dblRayGenData record) { + SurfaceRayFirePayload payload; + uint rayID = DispatchRaysIndex().x; + + // Trace the ray into the scene + RayDesc rayDesc; + rayDesc.Origin = float3(record.ray[rayID].origin); + rayDesc.Direction = float3(normalize(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_surf; + + // Pass the ray's origin and direction to the payload + payload.surf_id = ID_NONE; + payload.tlas = world; + payload.piv = xdg::PointInVolume::OUTSIDE; // Initialize point in volume check result to outside (0) + + TraceRay(world, RAY_FLAG_NONE, 0xff, xdg::RT_SURFACE_RAY_INDEX, xdg::RT_SURFACE_MISS_INDEX, rayDesc, payload); + + record.hit[rayID].surf_id = payload.surf_id; + record.hit[rayID].piv = payload.piv; // Point in volume check result +} + +// ------------------------------------------------- Compute Shaders ------------------------------------------------- +/* A shader to compute and store AABB min/maxes in single precision using double precision coords for Triangles*/ +[shader("compute")] +[numthreads(1, 1, 1)] +void +populate_tri_aabbs(uint3 DispatchThreadID: SV_DispatchThreadID, uniform DPTriangleGeomData record) { + uint primID = DispatchThreadID.x; + uint3 indices = record.index[primID]; + uint vertexIndices[3] = {indices.x, indices.y, indices.z}; + populate_aabb<3>(primID, record.vertex, record.aabbs, vertexIndices); +} + +// ------------------------------------------------ CUSTOM INTERSECTION SHADERS ------------------------------------------------ + +// Custom FP64 Plucker ray-triangle intersection algorithm for each ray-triangle pair +[shader("intersection")] +void DPTrianglePluckerIntersection(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.ray[rayID].origin; + dp::vec3 direction = record.ray[rayID].direction; + double tMin = record.ray[rayID].tMin; + double tMax = record.ray[rayID].tMax; + + 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); + + dp::vec3 norm = record.normals[primID]; // recover double precision normal. TODO - Should we calculate from vertices instead? + + // sense adjustment of normal + if (record.ray[rayID].volume_tree == record.reverse_tree) + { + norm = -norm; + } + + 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; + + xdg::HitOrientation hitOrientation = record.ray[rayID].hitOrientation; + + if (orientation_cull(direction, norm, hitOrientation)) + { + return; + } + + for (int i = 0; i < record.ray[rayID].exclude_count; ++i) + { + if (record.ray[rayID].exclude_primitives[i] == global_prim_id) { + return; + } + } + attr.global_prim_id = global_prim_id; + ReportHit(f32t, hit_kind, attr); +} + +#endif // XDG_GPRT_TRIANGLE_RT_SHADERS_SLANG diff --git a/src/tetrahedron_contain.cpp b/src/tetrahedron_contain.cpp index 743ea1d0..95dfc0b3 100644 --- a/src/tetrahedron_contain.cpp +++ b/src/tetrahedron_contain.cpp @@ -2,46 +2,12 @@ #include "xdg/ray_tracing_interface.h" #include "xdg/ray.h" #include "xdg/vec3da.h" - +#include "xdg/geometry/plucker.h" #include "xdg/util/linalg.h" namespace xdg { -bool plucker_tet_containment_test(const Position& point, - const Position& v0, - const Position& v1, - const Position& v2, - const Position& v3) { - using namespace linalg::aliases; - // Create matrix T = [v1 - v0, v2 - v0, v3 - v0] - Vec3da e0 = v1 - v0; - Vec3da e1 = v2 - v0; - Vec3da e2 = v3 - v0; - double3x3 T = { {e0.x, e0.y, e0.z}, - {e1.x, e1.y, e1.z}, - {e2.x, e2.y, e2.z}}; - - // Vector from v0 to point - Vec3da rhs = point - v0; - - // Solve T * [λ1, λ2, λ3] = rhs - double3 lambda123 = mul(inverse(T),{rhs.x, rhs.y, rhs.z}); - - // Compute λ0 - double lambda0 = 1.0f - (lambda123.x + lambda123.y + lambda123.z); - - // Final barycentric coordinate vector - double4 bary = { lambda0, lambda123.x, lambda123.y, lambda123.z }; - - // Check all λ_i in [0, 1] - for (int i = 0; i < 4; ++i) { - if (bary[i] < -PLUCKER_ZERO_TOL || bary[i] > 1.0f + PLUCKER_ZERO_TOL) - return false; - } - return true; -} - // Embree callbacks void VolumeElementBoundsFunc(RTCBoundsFunctionArguments* args) diff --git a/tests/test_find_element.cpp b/tests/test_find_element.cpp index 17023f2b..72ab64b3 100644 --- a/tests/test_find_element.cpp +++ b/tests/test_find_element.cpp @@ -1,45 +1,60 @@ // for testing #include +#include + // xdg includes #include "xdg/constants.h" #include "xdg/mesh_manager_interface.h" #include "xdg/embree/ray_tracer.h" +#include "util.h" #include "mesh_mock.h" using namespace xdg; +using namespace xdg::test; -TEST_CASE("Test Find Volumetric Element") +TEMPLATE_TEST_CASE("Test Find Volumetric Element", "[find_element][mock][volumetric]", + Embree_Raytracer, + GPRT_Raytracer) { - // create a mock mesh manager without volumetric elements - std::shared_ptr mm = std::make_shared(); - mm->init(); // this should do nothing - - REQUIRE(mm->num_volumes() == 1); - REQUIRE(mm->num_surfaces() == 6); - REQUIRE(mm->num_volume_elements(1) == 12); // should return 12 volumetric elements - - std::shared_ptr rti = std::make_shared(); - std::unordered_map volume_to_scene_map; - std::unordered_map element_to_scene_map; - for (auto volume: mm->volumes()) { - auto [volume_tree, element_tree] = rti->register_volume(mm, volume); - volume_to_scene_map[volume] = volume_tree; - element_to_scene_map[volume_tree] = element_tree; + // Generate one test run per enabled backend + constexpr auto rt_backend = TestType::value; + check_ray_tracer_supported(rt_backend); // skip if backend not enabled at configuration time + + DYNAMIC_SECTION(fmt::format("Backend = {}", rt_backend)) + { + auto rti = create_raytracer(rt_backend); + REQUIRE(rti); + // create a mock mesh manager without volumetric elements + std::shared_ptr mm = std::make_shared(); + mm->init(); // this should do nothing + + REQUIRE(mm->num_volumes() == 1); + REQUIRE(mm->num_surfaces() == 6); + REQUIRE(mm->num_volume_elements(1) == 12); // should return 12 volumetric elements + + std::unordered_map volume_to_scene_map; + std::unordered_map element_to_scene_map; + for (auto volume: mm->volumes()) { + auto [volume_tree, element_tree] = rti->register_volume(mm, volume); + volume_to_scene_map[volume] = volume_tree; + element_to_scene_map[volume_tree] = element_tree; + } + REQUIRE(rti->num_registered_trees() == 2); + rti->init(); // Ensure ray tracer is initialized (e.g. build SBT for GPRT after volumes registered) + + MeshID volume = 1; + + // test finding a volumetric element within the volume + Position point_inside {0.0, 0.0, 0.0}; // point inside the volume + MeshID element_id = rti->find_element(element_to_scene_map[volume], point_inside); + REQUIRE(element_id != ID_NONE); // should find an element + REQUIRE(element_id == 7); // Added this hardcode check to ensure GPRT and Embree return the same element ID for the same point + + Position point_outside {10.0, 10.0, 10.0}; // point outside the volume + element_id = rti->find_element(element_to_scene_map[volume], point_outside); + REQUIRE(element_id == ID_NONE); // should not find an element since the point is outside the volume } - REQUIRE(rti->num_registered_trees() == 2); - - MeshID volume = 1; - - // test finding a volumetric element within the volume - Position point_inside {0.0, 0.0, 0.0}; // point inside the volume - MeshID element_id = rti->find_element(element_to_scene_map[volume], point_inside); - REQUIRE(element_id != ID_NONE); // should find an element - - Position point_outside {10.0, 10.0, 10.0}; // point outside the volume - element_id = rti->find_element(element_to_scene_map[volume], point_outside); - REQUIRE(element_id == ID_NONE); // should not find an element since the point is outside the volume } -