From 781e09db84e11ceef107a55ce60a455d2bde8cb3 Mon Sep 17 00:00:00 2001 From: waqar-ukaea Date: Wed, 1 Jul 2026 10:09:08 +0100 Subject: [PATCH 1/5] Avoid full intersection point construction in Plucker distance calculation --- include/xdg/geometry/plucker.h | 19 ++++++++++--------- 1 file changed, 10 insertions(+), 9 deletions(-) diff --git a/include/xdg/geometry/plucker.h b/include/xdg/geometry/plucker.h index 8b7581c9..1916cee3 100644 --- a/include/xdg/geometry/plucker.h +++ b/include/xdg/geometry/plucker.h @@ -123,14 +123,6 @@ inline PluckerIntersectionResult plucker_ray_tri_intersect(dp::vec3 vertices[3], return {false, 0.0}; } - // 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; @@ -141,7 +133,16 @@ inline PluckerIntersectionResult plucker_ray_tri_intersect(dp::vec3 vertices[3], } } - dist_out = (intersection[idx] - origin[idx]) / direction[idx]; + // get the distance to intersection + const double inverse_sum = + 1.0 / (plucker_coord0 + plucker_coord1 + plucker_coord2); + + const double intersection_component = + (plucker_coord0 * vertices[2][idx] + + plucker_coord1 * vertices[0][idx] + + plucker_coord2 * vertices[1][idx]) * inverse_sum; + + dist_out = (intersection_component - origin[idx]) / direction[idx]; // Barycentric coords check double u = plucker_coord2 * inverse_sum; From cec22288f24efd81cbab37a6d484ba8673e292f7 Mon Sep 17 00:00:00 2001 From: waqar-ukaea Date: Wed, 1 Jul 2026 15:13:22 +0100 Subject: [PATCH 2/5] Scalarize Plucker edge test arithmetic - Replace temporary vec3 edge/cross/dot operations with explicit component wise checks - Not sure if this is robust enough. Changes suggested by codex --- include/xdg/geometry/plucker.h | 28 +++++++++++++++++++--------- 1 file changed, 19 insertions(+), 9 deletions(-) diff --git a/include/xdg/geometry/plucker.h b/include/xdg/geometry/plucker.h index 1916cee3..e94b86e0 100644 --- a/include/xdg/geometry/plucker.h +++ b/include/xdg/geometry/plucker.h @@ -43,15 +43,25 @@ inline bool first(dp::vec3 a, dp::vec3 b) { 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); + const bool vertexa_is_first = first(vertexa, vertexb); + + const double edge_x = vertexa_is_first ? vertexb[0] - vertexa[0] : vertexa[0] - vertexb[0]; + const double edge_y = vertexa_is_first ? vertexb[1] - vertexa[1] : vertexa[1] - vertexb[1]; + const double edge_z = vertexa_is_first ? vertexb[2] - vertexa[2] : vertexa[2] - vertexb[2]; + + const double origin_x = vertexa_is_first ? vertexa[0] : vertexb[0]; + const double origin_y = vertexa_is_first ? vertexa[1] : vertexb[1]; + const double origin_z = vertexa_is_first ? vertexa[2] : vertexb[2]; + + const double edge_normal_x = edge_y * origin_z - edge_z * origin_y; + const double edge_normal_y = edge_z * origin_x - edge_x * origin_z; + const double edge_normal_z = edge_x * origin_y - edge_y * origin_x; + + double pip = + ray[0] * edge_normal_x + ray[1] * edge_normal_y + ray[2] * edge_normal_z + + ray_normal[0] * edge_x + ray_normal[1] * edge_y + ray_normal[2] * edge_z; + + if (!vertexa_is_first) { pip = -pip; } if (dp::DBL_ZERO_TOL > dp::abs(pip)) // <-- absd From f1aa845085acfcd1d986d6310a59a301661247fa Mon Sep 17 00:00:00 2001 From: waqar-ukaea Date: Mon, 20 Jul 2026 11:16:11 +0100 Subject: [PATCH 3/5] Reduce live hit information during traversal by storing only best hit - Attempt to reduce register pressure - Recover full hit information once at the end after best hit found --- include/xdg/cuBQL/triangles.h | 4 +++ src/cuBQL/intersection.cpp | 59 ++++++++++++++++++++++++++--------- 2 files changed, 49 insertions(+), 14 deletions(-) diff --git a/include/xdg/cuBQL/triangles.h b/include/xdg/cuBQL/triangles.h index 07016816..6bb1f051 100644 --- a/include/xdg/cuBQL/triangles.h +++ b/include/xdg/cuBQL/triangles.h @@ -2,6 +2,7 @@ #define _XDG_CUBQL_TRIANGLES_H #include +#include // Guards to prevent CUDA headers from being included in host code, which causes // failed compilation with LLVM-clang. @@ -67,6 +68,9 @@ struct CuBQLSurfaceMesh { type is the compact, non-owning device-data view used during traversal. */ struct CuBQLVolumeGroup { + // Constant for invalid primitive index, used to indicate no hit in traversal kernel + static constexpr std::uint32_t INVALID_BVH_PRIMITIVE = std::numeric_limits::max(); + struct SurfaceDD { CuBQLSurfaceMesh::DD mesh; diff --git a/src/cuBQL/intersection.cpp b/src/cuBQL/intersection.cpp index e469803f..9e7fd8f2 100644 --- a/src/cuBQL/intersection.cpp +++ b/src/cuBQL/intersection.cpp @@ -17,6 +17,38 @@ static inline float reject_candidate(const cuBQL::ray3f& traversal_ray) return traversal_ray.tMax; } +static inline void store_surface_hit(CuBQLVolumeGroup::DD volume_group, + std::uint32_t prim_ref_index, + double distance, + const cuBQL::vec3d& direction, + CuBQLSurfaceHit* hit) +{ + const auto ref = volume_group.prim_refs[prim_ref_index]; + const auto surface = volume_group.surfaces[ref.surface_index]; + const auto mesh = surface.mesh; + const auto local_index = ref.primitive_index; + const cuBQL::vec3i vertex_indices = mesh.indices[local_index]; + + const cuBQL::vec3d vertex_a = mesh.vertices[vertex_indices.x]; + const cuBQL::vec3d vertex_b = mesh.vertices[vertex_indices.y]; + const cuBQL::vec3d vertex_c = mesh.vertices[vertex_indices.z]; + const cuBQL::vec3d normal = cuBQL::cross(vertex_b - vertex_a, + vertex_c - vertex_a); + + double normal_dot_direction = dot(normal, direction); + if (surface.reverse_sense) { + normal_dot_direction = -normal_dot_direction; + } + + hit->distance = distance; + hit->surface = mesh.surface_id; + hit->primitive = mesh.primitive_ids[local_index]; + hit->piv = normal_dot_direction > 0.0 ? INSIDE : OUTSIDE; + hit->next_volume = surface.next_volume; + hit->boundary_condition = surface.boundary_condition; + hit->normal = normal; +} + static inline void intersect_surface_tree(CuBQLVolumeGroup::DD volume_group, CuBQLRay intersection_ray, CuBQLSurfaceHit* hit, @@ -34,7 +66,11 @@ static inline void intersect_surface_tree(CuBQLVolumeGroup::DD volume_group, traversal_ray.tMin = static_cast(intersection_ray.tMin); traversal_ray.tMax = static_cast(hit->distance); - auto intersect_prim = [=, &traversal_ray] + // Nearest hit state. + double best_distance = hit->distance; + std::uint32_t best_prim_ref_index = CuBQLVolumeGroup::INVALID_BVH_PRIMITIVE; + + auto intersect_prim = [=, &traversal_ray, &best_distance, &best_prim_ref_index] (std::uint32_t bvh_primitive_index) -> float { const auto ref = volume_group.prim_refs[bvh_primitive_index]; @@ -43,13 +79,13 @@ static inline void intersect_surface_tree(CuBQLVolumeGroup::DD volume_group, const auto local_index = ref.primitive_index; const MeshID primitive_id = mesh.primitive_ids[local_index]; - // Reject the previously hit primitive to avoid immediate self-intersection. + // Reject the previously hit primitive to avoid immediate self-intersection if (primitive_id == last_hit_primitive) { return reject_candidate(traversal_ray); } - // Scalar queries may provide an arbitrary primitive exclusion history. - // TODO - Think about how to provide arbitrary history checks for batch queries. + // Scalar queries may provide an arbitrary primitive exclusion history + // TODO - Think about how to provide arbitrary history checks for batch queries for (int i = 0; i < exclude_count; ++i) { if (exclude_primitives[i] == primitive_id) { return reject_candidate(traversal_ray); @@ -80,25 +116,20 @@ static inline void intersect_surface_tree(CuBQLVolumeGroup::DD volume_group, auto intersection = plucker_ray_tri_intersect(vertices, intersection_ray.origin, intersection_ray.direction, - hit->distance, + best_distance, intersection_ray.tMin, false, 0); - // store ray payload if hit found + // Store only the best hit state needed to materialize the final hit after traversal if (intersection.hit) { - hit->distance = intersection.t; - hit->surface = mesh.surface_id; - hit->primitive = primitive_id; - hit->piv = normal_dot_direction > 0.0 ? INSIDE : OUTSIDE; - hit->next_volume = surface.next_volume; - hit->boundary_condition = surface.boundary_condition; - hit->normal = normal; + best_distance = intersection.t; + best_prim_ref_index = bvh_primitive_index; traversal_ray.tMax = static_cast(intersection.t); } // Return value is only the FP32 traversal shrink distance. The accepted hit - // distance stored above remains the FP64 Plucker result. + // distance stored above remains the FP64 Plucker result return reject_candidate(traversal_ray); }; From 51c8650d5074871f1dc7af957fc669297d9c6662 Mon Sep 17 00:00:00 2001 From: waqar-ukaea Date: Thu, 23 Jul 2026 17:10:16 +0100 Subject: [PATCH 4/5] Ensure the best hit is stored in cuBQL traversal/intersection --- src/cuBQL/intersection.cpp | 8 ++++++++ 1 file changed, 8 insertions(+) diff --git a/src/cuBQL/intersection.cpp b/src/cuBQL/intersection.cpp index 9e7fd8f2..8c317c37 100644 --- a/src/cuBQL/intersection.cpp +++ b/src/cuBQL/intersection.cpp @@ -135,6 +135,14 @@ static inline void intersect_surface_tree(CuBQLVolumeGroup::DD volume_group, // Single level traversal call for a shrinking ray query against the flattened BVH of the volume group. cuBQL::shrinkingRayQuery::forEachPrim(intersect_prim, volume_group.bvh, traversal_ray); + + if (best_prim_ref_index != CuBQLVolumeGroup::INVALID_BVH_PRIMITIVE) { + store_surface_hit(volume_group, + best_prim_ref_index, + best_distance, + intersection_ray.direction, + hit); + } } #pragma omp end declare target From 1c177259b4e67cdb3b12eba946d9f2052ec7156e Mon Sep 17 00:00:00 2001 From: waqar-ukaea Date: Wed, 5 Aug 2026 17:53:31 +0100 Subject: [PATCH 5/5] Ensure hit normal returned in rayhit buffer is normalised --- include/xdg/cuBQL/triangles.h | 4 ---- src/cuBQL/intersection.cpp | 18 +++++++++++++----- 2 files changed, 13 insertions(+), 9 deletions(-) diff --git a/include/xdg/cuBQL/triangles.h b/include/xdg/cuBQL/triangles.h index 6bb1f051..07016816 100644 --- a/include/xdg/cuBQL/triangles.h +++ b/include/xdg/cuBQL/triangles.h @@ -2,7 +2,6 @@ #define _XDG_CUBQL_TRIANGLES_H #include -#include // Guards to prevent CUDA headers from being included in host code, which causes // failed compilation with LLVM-clang. @@ -68,9 +67,6 @@ struct CuBQLSurfaceMesh { type is the compact, non-owning device-data view used during traversal. */ struct CuBQLVolumeGroup { - // Constant for invalid primitive index, used to indicate no hit in traversal kernel - static constexpr std::uint32_t INVALID_BVH_PRIMITIVE = std::numeric_limits::max(); - struct SurfaceDD { CuBQLSurfaceMesh::DD mesh; diff --git a/src/cuBQL/intersection.cpp b/src/cuBQL/intersection.cpp index 8c317c37..89301e71 100644 --- a/src/cuBQL/intersection.cpp +++ b/src/cuBQL/intersection.cpp @@ -67,8 +67,9 @@ static inline void intersect_surface_tree(CuBQLVolumeGroup::DD volume_group, traversal_ray.tMax = static_cast(hit->distance); // Nearest hit state. + constexpr std::uint32_t invalid_bvh_primitive = static_cast(-1); double best_distance = hit->distance; - std::uint32_t best_prim_ref_index = CuBQLVolumeGroup::INVALID_BVH_PRIMITIVE; + std::uint32_t best_prim_ref_index = invalid_bvh_primitive; auto intersect_prim = [=, &traversal_ray, &best_distance, &best_prim_ref_index] (std::uint32_t bvh_primitive_index) -> float @@ -136,7 +137,7 @@ static inline void intersect_surface_tree(CuBQLVolumeGroup::DD volume_group, // Single level traversal call for a shrinking ray query against the flattened BVH of the volume group. cuBQL::shrinkingRayQuery::forEachPrim(intersect_prim, volume_group.bvh, traversal_ray); - if (best_prim_ref_index != CuBQLVolumeGroup::INVALID_BVH_PRIMITIVE) { + if (best_prim_ref_index != invalid_bvh_primitive) { store_surface_hit(volume_group, best_prim_ref_index, best_distance, @@ -268,9 +269,16 @@ intersect_surface_tree_batch(const cubql::Context& context, d_ray_hits[ray_id].point_in_volume = static_cast(hit.piv); d_ray_hits[ray_id].next_volume = hit.next_volume; d_ray_hits[ray_id].boundary_condition = static_cast(hit.boundary_condition); - d_ray_hits[ray_id].normal[0] = hit.normal.x; - d_ray_hits[ray_id].normal[1] = hit.normal.y; - d_ray_hits[ray_id].normal[2] = hit.normal.z; + if (hit.primitive != ID_NONE) { + const cuBQL::vec3d normal = cuBQL::normalize(hit.normal); + d_ray_hits[ray_id].normal[0] = normal.x; + d_ray_hits[ray_id].normal[1] = normal.y; + d_ray_hits[ray_id].normal[2] = normal.z; + } else { + d_ray_hits[ray_id].normal[0] = 0.0; + d_ray_hits[ray_id].normal[1] = 0.0; + d_ray_hits[ray_id].normal[2] = 0.0; + } } }