CollisionDetection.cu

Go to the documentation of this file.

 
#include <cstdio>
#include <cuda.h>
#include <cuda_runtime.h>
#include <device_launch_parameters.h>
#include <string>
#include <vector>
 
#include "helios_vector_types.h"
 
// Forward declaration to avoid pulling the full helios global header (and its template-heavy transitive includes)
// into the CUDA translation unit. Defined in core/src/global.cpp.
namespace helios {
    void helios_runtime_error(const std::string &error_message);
}
 
#define HELIOS_CUDA_CHECK(call)                                                                                                                                                                                                                 \
    do {                                                                                                                                                                                                                                        \
        cudaError_t _helios_cuda_err = (call);                                                                                                                                                                                                  \
        if (_helios_cuda_err != cudaSuccess) {                                                                                                                                                                                                  \
            helios::helios_runtime_error(std::string("CUDA error (") + #call + "): " + cudaGetErrorString(_helios_cuda_err));                                                                                                                    \
        }                                                                                                                                                                                                                                       \
    } while (0)
 
struct GPUBVHNode {
    float3 aabb_min; 
    float3 aabb_max; 
    unsigned int left_child; 
    unsigned int right_child; 
    unsigned int primitive_start; 
    unsigned int primitive_count; 
    unsigned int is_leaf; 
    unsigned int padding; 
};
 
#define BVH_TRAVERSAL_STACK_CAPACITY 128
 
__device__ unsigned int d_bvh_stack_overflow = 0;
 
struct GPUBVHNodesSoA {
    // Hot data: frequently accessed during traversal (separate arrays for coalescing)
    float3 *aabb_mins; 
    float3 *aabb_maxs; 
    uint32_t *left_children; 
    uint32_t *right_children; 
 
    // Cold data: accessed less frequently
    uint32_t *primitive_starts; 
    uint32_t *primitive_counts; 
    uint8_t *is_leaf_flags; 
 
    size_t node_count; 
};
 
__device__ bool d_aabbIntersect(const float3 &min1, const float3 &max1, const float3 &min2, const float3 &max2) {
    return (min1.x <= max2.x && max1.x >= min2.x) && (min1.y <= max2.y && max1.y >= min2.y) && (min1.z <= max2.z && max1.z >= min2.z);
}
 
__device__ __forceinline__ float3 cross(const float3 &a, const float3 &b) {
    return make_float3(a.y * b.z - a.z * b.y, a.z * b.x - a.x * b.z, a.x * b.y - a.y * b.x);
}
 
__device__ __forceinline__ float dot(const float3 &a, const float3 &b) {
    return a.x * b.x + a.y * b.y + a.z * b.z;
}
 
__device__ __forceinline__ float3 normalize(const float3 &v) {
    float len = sqrtf(v.x * v.x + v.y * v.y + v.z * v.z);
    if (len > 1e-8f) {
        return make_float3(v.x / len, v.y / len, v.z / len);
    }
    return make_float3(0.0f, 0.0f, 1.0f); // Default up vector
}
 
__device__ __forceinline__ float3 operator+(const float3 &a, const float3 &b) {
    return make_float3(a.x + b.x, a.y + b.y, a.z + b.z);
}
 
__device__ __forceinline__ float3 operator-(const float3 &a, const float3 &b) {
    return make_float3(a.x - b.x, a.y - b.y, a.z - b.z);
}
 
__device__ __forceinline__ float3 operator*(const float3 &a, float scalar) {
    return make_float3(a.x * scalar, a.y * scalar, a.z * scalar);
}
 
__device__ __forceinline__ bool rayTriangleIntersect(const float3 &ray_origin, const float3 &ray_direction, const float3 &v0, const float3 &v1, const float3 &v2, float max_distance, float &hit_distance) {
    const float EPSILON = 1e-5f;  // Match CPU tolerance for consistent GPU/CPU behavior
 
    // Find vectors for two edges sharing v0
    float3 edge1 = make_float3(v1.x - v0.x, v1.y - v0.y, v1.z - v0.z);
    float3 edge2 = make_float3(v2.x - v0.x, v2.y - v0.y, v2.z - v0.z);
 
    // Begin calculating determinant - also used to calculate u parameter
    float3 h = make_float3(ray_direction.y * edge2.z - ray_direction.z * edge2.y, ray_direction.z * edge2.x - ray_direction.x * edge2.z, ray_direction.x * edge2.y - ray_direction.y * edge2.x);
 
    // If determinant is near zero, ray lies in plane of triangle
    float a = edge1.x * h.x + edge1.y * h.y + edge1.z * h.z;
    if (a > -EPSILON && a < EPSILON) {
        return false; // This ray is parallel to this triangle.
    }
 
    float f = 1.0f / a;
    float3 s = make_float3(ray_origin.x - v0.x, ray_origin.y - v0.y, ray_origin.z - v0.z);
    float u = f * (s.x * h.x + s.y * h.y + s.z * h.z);
 
    if (u < -EPSILON || u > 1.0f + EPSILON) {
        return false;
    }
 
    float3 q = make_float3(s.y * edge1.z - s.z * edge1.y, s.z * edge1.x - s.x * edge1.z, s.x * edge1.y - s.y * edge1.x);
 
    float v = f * (ray_direction.x * q.x + ray_direction.y * q.y + ray_direction.z * q.z);
 
    if (v < -EPSILON || u + v > 1.0f + EPSILON) {
        return false;
    }
 
    // At this stage we can compute t to find out where the intersection point is on the line.
    float t = f * (edge2.x * q.x + edge2.y * q.y + edge2.z * q.z);
 
    // Match CPU behavior - don't check max_distance here (checked by caller)
    if (t > EPSILON) { // ray intersection
        hit_distance = t;
        return true;
    } else { // This means that there is a line intersection but not a ray intersection.
        return false;
    }
}
 
__device__ __forceinline__ float3 safeRayInvDir(const float3 &dir) {
    constexpr float PARALLEL_EPS = 1e-8f;
    float dx = (fabsf(dir.x) < PARALLEL_EPS) ? copysignf(PARALLEL_EPS, dir.x) : dir.x;
    float dy = (fabsf(dir.y) < PARALLEL_EPS) ? copysignf(PARALLEL_EPS, dir.y) : dir.y;
    float dz = (fabsf(dir.z) < PARALLEL_EPS) ? copysignf(PARALLEL_EPS, dir.z) : dir.z;
    return make_float3(1.0f / dx, 1.0f / dy, 1.0f / dz);
}
 
__device__ __forceinline__ bool warpRayAABBIntersect(const float3 &ray_origin, const float3 &ray_dir, const float3 &aabb_min, const float3 &aabb_max, float max_dist) {
    // Optimized ray-AABB intersection using slab method
    // Compute intersection distances for each axis (sign-preserving reciprocal avoids 0*inf=NaN for axis-aligned rays).
    float3 inv_dir = safeRayInvDir(ray_dir);
 
    float3 t_min = make_float3((aabb_min.x - ray_origin.x) * inv_dir.x, (aabb_min.y - ray_origin.y) * inv_dir.y, (aabb_min.z - ray_origin.z) * inv_dir.z);
 
    float3 t_max = make_float3((aabb_max.x - ray_origin.x) * inv_dir.x, (aabb_max.y - ray_origin.y) * inv_dir.y, (aabb_max.z - ray_origin.z) * inv_dir.z);
 
    // Handle negative ray directions
    if (ray_dir.x < 0.0f) {
        float temp = t_min.x;
        t_min.x = t_max.x;
        t_max.x = temp;
    }
    if (ray_dir.y < 0.0f) {
        float temp = t_min.y;
        t_min.y = t_max.y;
        t_max.y = temp;
    }
    if (ray_dir.z < 0.0f) {
        float temp = t_min.z;
        t_min.z = t_max.z;
        t_max.z = temp;
    }
 
    // Find the intersection interval
    float t_enter = fmaxf(fmaxf(t_min.x, t_min.y), t_min.z);
    float t_exit = fminf(fminf(t_max.x, t_max.y), t_max.z);
 
    // Check if intersection exists and is within ray limits
    return (t_enter <= t_exit) && (t_exit >= 0.0f) && (t_enter <= max_dist);
}
 
// Device function for ray-triangle intersection (same algorithm as CPU)
__device__ __forceinline__ bool rayTriangleIntersectCPU(const float3 &origin, const float3 &direction, const float3 &v0, const float3 &v1, const float3 &v2, float &distance) {
 
    // Use same algorithm as CPU: radiation model's triangle_intersect
    const float EPSILON = 1e-5f;  // Match CPU tolerance for consistent GPU/CPU behavior
 
    float a = v0.x - v1.x, b = v0.x - v2.x, c = direction.x, d = v0.x - origin.x;
    float e = v0.y - v1.y, f = v0.y - v2.y, g = direction.y, h = v0.y - origin.y;
    float i = v0.z - v1.z, j = v0.z - v2.z, k = direction.z, l = v0.z - origin.z;
 
    float m = f * k - g * j, n = h * k - g * l, p = f * l - h * j;
    float q = g * i - e * k, s = e * j - f * i;
 
    float denom = a * m + b * q + c * s;
    if (fabsf(denom) < EPSILON) {
        return false; // Ray is parallel to triangle
    }
 
    float inv_denom = 1.0f / denom;
 
    float e1 = d * m - b * n - c * p;
    float beta = e1 * inv_denom;
 
    if (beta >= -EPSILON) {
        float r = e * l - h * i;
        float e2 = a * n + d * q + c * r;
        float gamma = e2 * inv_denom;
 
        if (gamma >= -EPSILON && beta + gamma <= 1.0f + EPSILON) {
            float e3 = a * p - b * r + d * s;
            float t = e3 * inv_denom;
 
            if (t > EPSILON) {
                distance = t;
                return true;
            }
        }
    }
    return false;
}
 
// Device function for ray-patch intersection (same algorithm as CPU)
__device__ __forceinline__ bool rayPatchIntersect(const float3 &origin, const float3 &direction, const float3 &v0, const float3 &v1, const float3 &v2, const float3 &v3, float &distance) {
 
    // Calculate patch vectors and normal (same as CPU radiation model)
    const float EPSILON = 1e-5f;  // Match triangle epsilon for consistency
 
    float3 anchor = v0;
    float3 normal = cross(v1 - v0, v2 - v0);
    normal = normalize(normal);
 
    float3 a = v1 - v0; // First edge vector
    float3 b = v3 - v0; // Second edge vector
 
    // Ray-plane intersection
    float denom = dot(direction, normal);
    if (fabsf(denom) > EPSILON) { // Not parallel to plane
        float t = dot(anchor - origin, normal) / denom;
 
        if (t > EPSILON && t < 1e8f) { // Valid intersection distance
            // Find intersection point
            float3 p = origin + direction * t;
            float3 d = p - anchor;
 
            // Project onto patch coordinate system
            float ddota = dot(d, a);
            float ddotb = dot(d, b);
 
            // Check if point is within patch bounds
            if (ddota >= 0.0f && ddota <= dot(a, a) && ddotb >= 0.0f && ddotb <= dot(b, b)) {
 
                distance = t;
                return true;
            }
        }
    }
    return false;
}
 
// Ray-AABB intersection for voxel primitives
__device__ bool rayVoxelIntersect(const float3 &ray_origin, const float3 &ray_direction, const float3 &aabb_min, const float3 &aabb_max, float &distance) {
    const float EPSILON = 1e-5f;  // Match triangle/patch epsilon for consistency
 
    // Calculate t values for each slab (sign-preserving reciprocal avoids 0*inf=NaN for axis-aligned rays grazing a face).
    float3 inv_dir = safeRayInvDir(ray_direction);
 
    float3 t_min = make_float3((aabb_min.x - ray_origin.x) * inv_dir.x, (aabb_min.y - ray_origin.y) * inv_dir.y, (aabb_min.z - ray_origin.z) * inv_dir.z);
 
    float3 t_max = make_float3((aabb_max.x - ray_origin.x) * inv_dir.x, (aabb_max.y - ray_origin.y) * inv_dir.y, (aabb_max.z - ray_origin.z) * inv_dir.z);
 
    // Handle negative ray directions
    if (ray_direction.x < 0.0f) {
        float temp = t_min.x;
        t_min.x = t_max.x;
        t_max.x = temp;
    }
    if (ray_direction.y < 0.0f) {
        float temp = t_min.y;
        t_min.y = t_max.y;
        t_max.y = temp;
    }
    if (ray_direction.z < 0.0f) {
        float temp = t_min.z;
        t_min.z = t_max.z;
        t_max.z = temp;
    }
 
    // Find the intersection interval
    float t_enter = fmaxf(fmaxf(t_min.x, t_min.y), t_min.z);
    float t_exit = fminf(fminf(t_max.x, t_max.y), t_max.z);
 
    // Check for intersection
    if (t_enter > t_exit || t_exit < EPSILON) {
        return false; // No intersection or behind ray
    }
 
    // Set distance to entry point (or exit if ray starts inside)
    distance = (t_enter > EPSILON) ? t_enter : t_exit;
 
    return distance > EPSILON;
}
 
__device__ __forceinline__ float3 computeHitNormal(int ptype, const float3 *d_primitive_vertices, unsigned int vertex_offset, const float3 &ray_origin, const float3 &ray_direction, float hit_distance) {
    float3 fallback = normalize(make_float3(-ray_direction.x, -ray_direction.y, -ray_direction.z));
 
    if (ptype == 1 || ptype == 0) { // triangle or patch: cross of first two edges, face-forwarded
        float3 v0 = d_primitive_vertices[vertex_offset + 0];
        float3 v1 = d_primitive_vertices[vertex_offset + 1];
        float3 v2 = d_primitive_vertices[vertex_offset + 2];
        float3 n = cross(v1 - v0, v2 - v0);
        float mag = sqrtf(n.x * n.x + n.y * n.y + n.z * n.z);
        if (mag > 1e-8f) {
            n = make_float3(n.x / mag, n.y / mag, n.z / mag);
            float3 hit_point = make_float3(ray_origin.x + ray_direction.x * hit_distance, ray_origin.y + ray_direction.y * hit_distance, ray_origin.z + ray_direction.z * hit_distance);
            float3 to_origin = make_float3(ray_origin.x - hit_point.x, ray_origin.y - hit_point.y, ray_origin.z - hit_point.z);
            if (dot(n, to_origin) < 0.0f) {
                n = make_float3(-n.x, -n.y, -n.z);
            }
            return n;
        }
        return fallback;
    } else if (ptype == 2) { // voxel: axis-aligned face normal of the hit face
        float3 vmin = d_primitive_vertices[vertex_offset + 0];
        float3 vmax = d_primitive_vertices[vertex_offset + 1];
        float3 hit_point = make_float3(ray_origin.x + ray_direction.x * hit_distance, ray_origin.y + ray_direction.y * hit_distance, ray_origin.z + ray_direction.z * hit_distance);
        float3 center = make_float3((vmin.x + vmax.x) * 0.5f, (vmin.y + vmax.y) * 0.5f, (vmin.z + vmax.z) * 0.5f);
        float3 extent = make_float3((vmax.x - vmin.x) * 0.5f, (vmax.y - vmin.y) * 0.5f, (vmax.z - vmin.z) * 0.5f);
        float3 local = make_float3(hit_point.x - center.x, hit_point.y - center.y, hit_point.z - center.z);
        float rel_x = fabsf(local.x) / extent.x;
        float rel_y = fabsf(local.y) / extent.y;
        float rel_z = fabsf(local.z) / extent.z;
        if (rel_x >= rel_y && rel_x >= rel_z) {
            return make_float3((local.x > 0.0f) ? 1.0f : -1.0f, 0.0f, 0.0f);
        } else if (rel_y >= rel_z) {
            return make_float3(0.0f, (local.y > 0.0f) ? 1.0f : -1.0f, 0.0f);
        } else {
            return make_float3(0.0f, 0.0f, (local.z > 0.0f) ? 1.0f : -1.0f);
        }
    }
    return fallback;
}
 
__device__ __forceinline__ bool sampleMaskOpaqueGPU(int mask_id, float u, float v, const unsigned char *d_mask_data, const unsigned int *d_mask_offsets, const int *d_mask_sizes) {
    if (mask_id < 0 || d_mask_data == nullptr) {
        return true; // no mask -> fully solid
    }
    const int width = d_mask_sizes[mask_id * 2];
    const int height = d_mask_sizes[mask_id * 2 + 1];
    if (width <= 0 || height <= 0) {
        return true;
    }
    const unsigned int offset = d_mask_offsets[mask_id];
 
    u -= floorf(u); // wrap repeat-style mappings into [0,1)
    v -= floorf(v);
 
    int px = (int) (u * (float) width);
    px = max(0, min(px, width - 1));
    int py = (int) ((1.f - v) * (float) height);
    py = max(0, min(py, height - 1));
 
    return d_mask_data[offset + (unsigned int) (py * width + px)] != 0u;
}
 
__device__ __forceinline__ bool isHitOpaqueGPU(int ptype, const float3 *verts, int mask_id, int uv_id, const float *uv4, const float3 &hit_point, const unsigned char *d_mask_data, const unsigned int *d_mask_offsets, const int *d_mask_sizes) {
    if (mask_id < 0) {
        return true;
    }
    float u, v;
    if (ptype == 0) { // patch: corners (BL, BR, TR, TL)
        const float3 v0 = verts[0];
        const float3 e1 = verts[1] - v0;
        const float3 e2 = verts[3] - v0;
        const float3 d = hit_point - v0;
        const float e1_sq = dot(e1, e1);
        const float e2_sq = dot(e2, e2);
        float s = (e1_sq > 0.f) ? dot(d, e1) / e1_sq : 0.f;
        float t = (e2_sq > 0.f) ? dot(d, e2) / e2_sq : 0.f;
        s = fminf(fmaxf(s, 0.f), 1.f);
        t = fminf(fmaxf(t, 0.f), 1.f);
        if (uv_id >= 0) {
            u = (1.f - s) * (1.f - t) * uv4[0] + s * (1.f - t) * uv4[2] + s * t * uv4[4] + (1.f - s) * t * uv4[6];
            v = (1.f - s) * (1.f - t) * uv4[1] + s * (1.f - t) * uv4[3] + s * t * uv4[5] + (1.f - s) * t * uv4[7];
        } else {
            u = s;
            v = t;
        }
    } else if (ptype == 1) { // triangle
        if (uv_id < 0) {
            return true; // no UVs -> cannot map texel, treat as solid (matches CPU)
        }
        const float3 v0 = verts[0];
        const float3 e1 = verts[1] - v0;
        const float3 e2 = verts[2] - v0;
        const float3 d = hit_point - v0;
        const float dot11 = dot(e1, e1);
        const float dot12 = dot(e1, e2);
        const float dot22 = dot(e2, e2);
        const float dot1d = dot(e1, d);
        const float dot2d = dot(e2, d);
        const float denom = dot11 * dot22 - dot12 * dot12;
        if (fabsf(denom) < 1e-20f) {
            return true; // degenerate triangle
        }
        const float inv = 1.f / denom;
        const float beta = (dot22 * dot1d - dot12 * dot2d) * inv;
        const float gamma = (dot11 * dot2d - dot12 * dot1d) * inv;
        u = uv4[0] + beta * (uv4[2] - uv4[0]) + gamma * (uv4[4] - uv4[0]);
        v = uv4[1] + beta * (uv4[3] - uv4[1]) + gamma * (uv4[5] - uv4[1]);
    } else {
        return true; // voxel / unknown -> solid
    }
    return sampleMaskOpaqueGPU(mask_id, u, v, d_mask_data, d_mask_offsets, d_mask_sizes);
}
 
__global__ void rayPrimitiveBVHKernel(GPUBVHNode *d_bvh_nodes, unsigned int *d_primitive_indices,
                                      int *d_primitive_types, // Type of each primitive (int for GPU compatibility)
                                      float3 *d_primitive_vertices, // All vertices for all primitives (variable count per primitive)
                                      unsigned int *d_vertex_offsets, // Starting index in vertices array for each primitive
                                      const unsigned char *d_mask_data, // Texture transparency: flattened mask bytes (or null)
                                      const unsigned int *d_mask_offsets, // Per-mask start index into d_mask_data
                                      const int *d_mask_sizes, // Per-mask width/height (2 ints per mask)
                                      const int *d_mask_IDs, // Per-primitive mask index (-1 = none)
                                      const float *d_uv_data, // Per-primitive UVs (4 vec2 = 8 floats per primitive)
                                      const int *d_uv_IDs, // Per-primitive UV flag (-1 = parametric/none)
                                      float3 *d_ray_origins, float3 *d_ray_directions, float *d_ray_max_distances, float uniform_max_distance, int num_rays, int primitive_count, int total_vertex_count, float *d_hit_distances,
                                      unsigned int *d_hit_primitive_ids, unsigned int *d_hit_counts, float3 *d_hit_normals, bool find_closest_hit) {
    int ray_idx = blockIdx.x * blockDim.x + threadIdx.x;
 
    if (ray_idx >= num_rays) {
        return;
    }
 
    // Load ray data. Directions are normalized here (the host passes them raw) so the returned t is a geometric
    // distance, matching the CPU RayQuery semantics. d_ray_max_distances may be null, in which case every ray shares
    // uniform_max_distance (avoids a per-ray host/device array when the caller's cutoff is constant).
    float3 ray_origin = d_ray_origins[ray_idx];
    float3 ray_direction = d_ray_directions[ray_idx];
    float dmag = sqrtf(ray_direction.x * ray_direction.x + ray_direction.y * ray_direction.y + ray_direction.z * ray_direction.z);
    if (dmag > 1e-8f) {
        ray_direction = make_float3(ray_direction.x / dmag, ray_direction.y / dmag, ray_direction.z / dmag);
    }
    float ray_max_distance = (d_ray_max_distances != nullptr) ? d_ray_max_distances[ray_idx] : uniform_max_distance;
 
    // Initialize hit data
    float closest_hit_distance = ray_max_distance + 1.0f; // Initialize beyond max
    unsigned int hit_primitive_id = 0xFFFFFFFF; // Invalid ID
    unsigned int total_hits = 0;
    int best_vertex_offset = -1; // packed-vertex offset of the current closest hit (for normal computation)
    int best_ptype = -1; // primitive type of the current closest hit
 
    // Stack-based BVH traversal. Per-thread local stack (thread-private; sized to BVH_TRAVERSAL_STACK_CAPACITY so it exceeds
    // the host builder's MAX_DEPTH). The old block-shared 32-entry-per-thread array silently dropped children once the new
    // deeper SAH trees pushed past depth 31, causing missed/farther hits that diverged from the CPU result.
    unsigned int thread_stack[BVH_TRAVERSAL_STACK_CAPACITY];
    int stack_size = 0;
 
    // Start from root node
    thread_stack[0] = 0;
    stack_size = 1;
 
    // Main traversal loop
    while (stack_size > 0) {
        // Pop node from stack
        stack_size--;
        unsigned int node_idx = thread_stack[stack_size];
 
        if (node_idx == 0xFFFFFFFF) {
            continue;
        }
 
        GPUBVHNode node = d_bvh_nodes[node_idx];
 
        // Test ray-AABB intersection first (early rejection)
        if (!warpRayAABBIntersect(ray_origin, ray_direction, node.aabb_min, node.aabb_max, ray_max_distance)) {
            continue;
        }
 
        if (node.is_leaf) {
            // Test ray against all triangles in this leaf
            for (unsigned int i = 0; i < node.primitive_count; i++) {
                unsigned int primitive_index = node.primitive_start + i;
 
                // Bounds check for primitive_index
                if (primitive_index >= primitive_count) {
                    continue; // Skip this primitive, don't exit the entire kernel!
                }
 
                unsigned int primitive_id = d_primitive_indices[primitive_index];
 
                // Get primitive type
                int ptype = d_primitive_types[primitive_index];
 
                // Get primitive vertices starting index
                unsigned int vertex_offset = d_vertex_offsets[primitive_index];
 
                // Test ray-primitive intersection based on type
                float hit_distance;
                bool hit = false;
 
                if (ptype == 1) { // PRIMITIVE_TYPE_TRIANGLE
                    if (vertex_offset + 2 >= total_vertex_count) {
                        continue; // Skip this primitive, don't exit the entire kernel!
                    }
 
                    float3 v0 = d_primitive_vertices[vertex_offset + 0];
                    float3 v1 = d_primitive_vertices[vertex_offset + 1];
                    float3 v2 = d_primitive_vertices[vertex_offset + 2];
 
                    hit = rayTriangleIntersect(ray_origin, ray_direction, v0, v1, v2, ray_max_distance, hit_distance);
 
                } else if (ptype == 0) { // PRIMITIVE_TYPE_PATCH
                    if (vertex_offset + 3 >= total_vertex_count) {
                        continue; // Skip this primitive, don't exit the entire kernel!
                    }
 
                    float3 v0 = d_primitive_vertices[vertex_offset + 0];
                    float3 v1 = d_primitive_vertices[vertex_offset + 1];
                    float3 v2 = d_primitive_vertices[vertex_offset + 2];
                    float3 v3 = d_primitive_vertices[vertex_offset + 3];
 
                    hit = rayPatchIntersect(ray_origin, ray_direction, v0, v1, v2, v3, hit_distance);
 
                } else if (ptype == 2) { // PRIMITIVE_TYPE_VOXEL
                    // Voxel intersection using AABB intersection
                    // For voxels, vertices store min/max coordinates: [min.x, min.y, min.z, max.x, max.y, max.z, 0, 0]
                    if (vertex_offset + 1 >= total_vertex_count) {
                        continue; // Skip this primitive, don't exit the entire kernel!
                    }
 
                    float3 voxel_min = d_primitive_vertices[vertex_offset + 0];
                    float3 voxel_max = d_primitive_vertices[vertex_offset + 1];
 
                    hit = rayVoxelIntersect(ray_origin, ray_direction, voxel_min, voxel_max, hit_distance);
                }
 
                // Process hit if found
 
                if (hit && hit_distance > 1e-5f && hit_distance <= ray_max_distance) {
 
                    // Texture transparency: if this primitive has a mask and the hit lands on a transparent texel, the ray
                    // passes through it (skip without recording) so traversal continues to geometry behind, matching the
                    // CPU isHitTexelOpaque() path. mask_IDs/uv_data are indexed by primitive_index (BVH order), like vertices.
                    const int mask_id = (d_mask_IDs != nullptr) ? d_mask_IDs[primitive_index] : -1;
                    if (mask_id >= 0) {
                        const float3 hit_point = make_float3(ray_origin.x + ray_direction.x * hit_distance, ray_origin.y + ray_direction.y * hit_distance, ray_origin.z + ray_direction.z * hit_distance);
                        if (!isHitOpaqueGPU(ptype, &d_primitive_vertices[vertex_offset], mask_id, d_uv_IDs[primitive_index], &d_uv_data[primitive_index * 8], hit_point, d_mask_data, d_mask_offsets, d_mask_sizes)) {
                            continue; // transparent texel -> not a hit
                        }
                    }
 
                    total_hits++;
 
                    if (find_closest_hit) {
                        // Keep only the closest hit
                        if (hit_distance < closest_hit_distance) {
                            closest_hit_distance = hit_distance;
                            hit_primitive_id = primitive_id;
                            best_vertex_offset = (int) vertex_offset;
                            best_ptype = ptype;
                        }
                    } else {
                        // For collision detection, any hit is sufficient
                        d_hit_distances[ray_idx] = hit_distance;
                        d_hit_primitive_ids[ray_idx] = primitive_id;
                        d_hit_counts[ray_idx] = 1; // Found at least one hit
                        if (d_hit_normals) {
                            d_hit_normals[ray_idx] = computeHitNormal(ptype, d_primitive_vertices, vertex_offset, ray_origin, ray_direction, hit_distance);
                        }
                        return; // Early exit for collision detection
                    }
                }
            }
        } else {
            // Add child nodes to stack (add right child first for left-first traversal). The capacity is sized well above
            // the tree depth; if a push would overflow, raise the device flag (the host throws) rather than silently
            // dropping the child and reporting a wrong/missed hit.
            if (node.right_child != 0xFFFFFFFF) {
                if (stack_size < BVH_TRAVERSAL_STACK_CAPACITY) {
                    thread_stack[stack_size] = node.right_child;
                    stack_size++;
                } else {
                    atomicMax(&d_bvh_stack_overflow, 1u);
                }
            }
            if (node.left_child != 0xFFFFFFFF) {
                if (stack_size < BVH_TRAVERSAL_STACK_CAPACITY) {
                    thread_stack[stack_size] = node.left_child;
                    stack_size++;
                } else {
                    atomicMax(&d_bvh_stack_overflow, 1u);
                }
            }
        }
    }
 
    // Store final results
    if (find_closest_hit && hit_primitive_id != 0xFFFFFFFF) {
        d_hit_distances[ray_idx] = closest_hit_distance;
        d_hit_primitive_ids[ray_idx] = hit_primitive_id;
        d_hit_counts[ray_idx] = 1;
        if (d_hit_normals) {
            d_hit_normals[ray_idx] = computeHitNormal(best_ptype, d_primitive_vertices, (unsigned int) best_vertex_offset, ray_origin, ray_direction, closest_hit_distance);
        }
    } else if (!find_closest_hit) {
        d_hit_distances[ray_idx] = ray_max_distance + 1.0f; // No hit
        d_hit_primitive_ids[ray_idx] = 0xFFFFFFFF;
        d_hit_counts[ray_idx] = 0;
        if (d_hit_normals) {
            d_hit_normals[ray_idx] = make_float3(0.0f, 0.0f, 0.0f);
        }
    } else {
        // No hit found
        d_hit_distances[ray_idx] = ray_max_distance + 1.0f;
        d_hit_primitive_ids[ray_idx] = 0xFFFFFFFF;
        d_hit_counts[ray_idx] = 0;
        if (d_hit_normals) {
            d_hit_normals[ray_idx] = make_float3(0.0f, 0.0f, 0.0f);
        }
    }
}
 
// C-style wrapper functions for calling from C++ code
extern "C" {
 
void launchRaysOnResidentScene(void *d_bvh_nodes, int node_count, unsigned int *d_primitive_indices, int primitive_count, int *d_primitive_types, float3 *d_primitive_vertices, unsigned int *d_vertex_offsets, const unsigned char *d_mask_data,
                               const unsigned int *d_mask_offsets, const int *d_mask_sizes, const int *d_mask_IDs, const float *d_uv_data, const int *d_uv_IDs, int total_vertex_count, const float *h_ray_origins, const float *h_ray_directions,
                               const float *h_ray_max_distances, float uniform_max_distance, int num_rays, float *h_hit_distances, unsigned int *h_hit_primitive_ids, unsigned int *h_hit_counts, float *h_hit_normals, bool find_closest_hit) {
    if (num_rays == 0) {
        return;
    }
 
    // Allocate ONLY the per-call ray + result buffers; the geometry pointers are resident and owned by the caller. Each
    // optional host pointer that is null skips its device buffer / copy, so a caller writing straight into its own arrays
    // (the SoA path) drives no extra host staging. d_hit_counts is always allocated because the kernel writes it.
    float3 *d_ray_origins = nullptr, *d_ray_directions = nullptr, *d_hit_normals = nullptr;
    float *d_ray_max_distances = nullptr, *d_hit_distances = nullptr;
    unsigned int *d_hit_primitive_ids = nullptr, *d_hit_counts = nullptr;
 
    const size_t ray_data_size = size_t(num_rays) * sizeof(float3);
    const size_t ray_distances_size = size_t(num_rays) * sizeof(float);
    const size_t hit_results_size = size_t(num_rays) * sizeof(unsigned int);
 
    HELIOS_CUDA_CHECK(cudaMalloc(&d_ray_origins, ray_data_size));
    HELIOS_CUDA_CHECK(cudaMalloc(&d_ray_directions, ray_data_size));
    HELIOS_CUDA_CHECK(cudaMalloc(&d_hit_distances, ray_distances_size));
    HELIOS_CUDA_CHECK(cudaMalloc(&d_hit_primitive_ids, hit_results_size));
    HELIOS_CUDA_CHECK(cudaMalloc(&d_hit_counts, hit_results_size));
    const bool per_ray_max = (h_ray_max_distances != nullptr);
    if (per_ray_max) {
        HELIOS_CUDA_CHECK(cudaMalloc(&d_ray_max_distances, ray_distances_size));
    }
    const bool want_normals = (h_hit_normals != nullptr);
    if (want_normals) {
        HELIOS_CUDA_CHECK(cudaMalloc(&d_hit_normals, ray_data_size));
    }
 
    // helios::vec3 and the caller's flat xyz arrays have the exact byte layout of float3, so the host ray buffers copy
    // straight into the device float3 buffers with no host repack. Directions are normalized in the kernel.
    HELIOS_CUDA_CHECK(cudaMemcpy(d_ray_origins, h_ray_origins, ray_data_size, cudaMemcpyHostToDevice));
    HELIOS_CUDA_CHECK(cudaMemcpy(d_ray_directions, h_ray_directions, ray_data_size, cudaMemcpyHostToDevice));
    if (per_ray_max) {
        HELIOS_CUDA_CHECK(cudaMemcpy(d_ray_max_distances, h_ray_max_distances, ray_distances_size, cudaMemcpyHostToDevice));
    }
 
    int threads_per_block = 256;
    int num_blocks = (num_rays + threads_per_block - 1) / threads_per_block;
 
    // Reset the traversal-stack overflow flag before launch; checked after sync (fail-fast on a deeper-than-expected tree).
    const unsigned int stack_overflow_reset = 0;
    HELIOS_CUDA_CHECK(cudaMemcpyToSymbol(d_bvh_stack_overflow, &stack_overflow_reset, sizeof(unsigned int)));
 
    rayPrimitiveBVHKernel<<<num_blocks, threads_per_block>>>((GPUBVHNode *) d_bvh_nodes, d_primitive_indices, d_primitive_types, d_primitive_vertices, d_vertex_offsets, d_mask_data, d_mask_offsets, d_mask_sizes, d_mask_IDs, d_uv_data, d_uv_IDs,
                                                             d_ray_origins, d_ray_directions, d_ray_max_distances, uniform_max_distance, num_rays, primitive_count, total_vertex_count, d_hit_distances, d_hit_primitive_ids, d_hit_counts,
                                                             d_hit_normals, find_closest_hit);
 
    cudaDeviceSynchronize();
    HELIOS_CUDA_CHECK(cudaGetLastError());
 
    unsigned int stack_overflow_flag = 0;
    HELIOS_CUDA_CHECK(cudaMemcpyFromSymbol(&stack_overflow_flag, d_bvh_stack_overflow, sizeof(unsigned int)));
    if (stack_overflow_flag != 0) {
        helios::helios_runtime_error("ERROR (CollisionDetection GPU): BVH traversal stack overflow in rayPrimitiveBVHKernel - tree depth exceeded BVH_TRAVERSAL_STACK_CAPACITY (" + std::to_string(BVH_TRAVERSAL_STACK_CAPACITY) +
                                     "). This must not happen for a valid SAH tree; raise BVH_TRAVERSAL_STACK_CAPACITY or check the BVH build.");
    }
 
    HELIOS_CUDA_CHECK(cudaMemcpy(h_hit_distances, d_hit_distances, ray_distances_size, cudaMemcpyDeviceToHost));
    HELIOS_CUDA_CHECK(cudaMemcpy(h_hit_primitive_ids, d_hit_primitive_ids, hit_results_size, cudaMemcpyDeviceToHost));
    if (h_hit_counts != nullptr) {
        HELIOS_CUDA_CHECK(cudaMemcpy(h_hit_counts, d_hit_counts, hit_results_size, cudaMemcpyDeviceToHost));
    }
    if (want_normals) {
        HELIOS_CUDA_CHECK(cudaMemcpy(h_hit_normals, d_hit_normals, ray_data_size, cudaMemcpyDeviceToHost));
    }
 
    cudaFree(d_ray_origins);
    cudaFree(d_ray_directions);
    if (d_ray_max_distances) {
        cudaFree(d_ray_max_distances);
    }
    cudaFree(d_hit_distances);
    cudaFree(d_hit_primitive_ids);
    cudaFree(d_hit_counts);
    if (d_hit_normals) {
        cudaFree(d_hit_normals);
    }
}
 
__global__ void bvhTraversalKernel(GPUBVHNode *d_nodes, unsigned int *d_primitive_indices, float3 *d_primitive_aabb_min, float3 *d_primitive_aabb_max, float3 *d_query_aabb_min, float3 *d_query_aabb_max, unsigned int *d_results,
                                   unsigned int *d_result_counts, int num_queries, int max_results_per_query) {
 
    int query_idx = blockIdx.x * blockDim.x + threadIdx.x;
 
    if (query_idx >= num_queries)
        return;
 
    float3 query_min = d_query_aabb_min[query_idx];
    float3 query_max = d_query_aabb_max[query_idx];
 
    unsigned int result_count = 0;
    unsigned int *query_results = &d_results[query_idx * max_results_per_query];
 
    // Stack-based traversal. Per-thread local stack (thread-private; sized to BVH_TRAVERSAL_STACK_CAPACITY > MAX_DEPTH). The
    // old block-shared 32-entry-per-thread array silently dropped children on deep SAH trees, missing real collisions.
    unsigned int thread_stack[BVH_TRAVERSAL_STACK_CAPACITY];
    int stack_size = 0;
 
    // Start traversal from root node
    thread_stack[0] = 0;
    stack_size = 1;
 
    while (stack_size > 0 && result_count < max_results_per_query) {
 
        // Pop node from stack
        stack_size--;
        unsigned int node_idx = thread_stack[stack_size];
 
        // Check if node index is valid
        if (node_idx == 0xFFFFFFFF)
            continue;
 
        GPUBVHNode node = d_nodes[node_idx];
 
        // Test if query AABB intersects node AABB
        if (!d_aabbIntersect(query_min, query_max, node.aabb_min, node.aabb_max)) {
            continue;
        }
 
        if (node.is_leaf) {
            // Check each primitive in this leaf individually
            for (unsigned int i = 0; i < node.primitive_count && result_count < max_results_per_query; i++) {
                unsigned int primitive_index = node.primitive_start + i;
                unsigned int primitive_id = d_primitive_indices[primitive_index];
 
                // Get primitive's AABB from pre-computed arrays (using array position, not UUID)
                float3 prim_min = d_primitive_aabb_min[primitive_index];
                float3 prim_max = d_primitive_aabb_max[primitive_index];
 
                // Only add to results if AABBs actually intersect
                if (d_aabbIntersect(query_min, query_max, prim_min, prim_max)) {
                    query_results[result_count] = primitive_id;
                    result_count++;
                }
            }
        } else {
            // Add child nodes to stack. Capacity is sized well above tree depth; an overflow raises the device flag (the
            // host throws) rather than silently dropping a child and under-reporting collisions.
            if (node.left_child != 0xFFFFFFFF) {
                if (stack_size < BVH_TRAVERSAL_STACK_CAPACITY) {
                    thread_stack[stack_size] = node.left_child;
                    stack_size++;
                } else {
                    atomicMax(&d_bvh_stack_overflow, 1u);
                }
            }
            if (node.right_child != 0xFFFFFFFF) {
                if (stack_size < BVH_TRAVERSAL_STACK_CAPACITY) {
                    thread_stack[stack_size] = node.right_child;
                    stack_size++;
                } else {
                    atomicMax(&d_bvh_stack_overflow, 1u);
                }
            }
        }
    }
 
    d_result_counts[query_idx] = result_count;
}
 
void launchBVHTraversal(void *h_nodes, int node_count, unsigned int *h_primitive_indices, int primitive_count, float *h_primitive_aabb_min, float *h_primitive_aabb_max, float *h_query_aabb_min, float *h_query_aabb_max, int num_queries,
                        unsigned int *h_results, unsigned int *h_result_counts, int max_results_per_query) {
 
    if (num_queries == 0)
        return;
 
    // Allocate GPU memory for query data and primitive AABBs
    float3 *d_query_min;
    float3 *d_query_max;
    float3 *d_primitive_min;
    float3 *d_primitive_max;
    unsigned int *d_results;
    unsigned int *d_result_counts;
 
    size_t query_size = num_queries * sizeof(float3);
    size_t primitive_aabb_size = primitive_count * sizeof(float3);
    size_t results_size = num_queries * max_results_per_query * sizeof(unsigned int);
    size_t counts_size = num_queries * sizeof(unsigned int);
 
    cudaMalloc((void **) &d_query_min, query_size);
    cudaMalloc((void **) &d_query_max, query_size);
    cudaMalloc((void **) &d_primitive_min, primitive_aabb_size);
    cudaMalloc((void **) &d_primitive_max, primitive_aabb_size);
    cudaMalloc((void **) &d_results, results_size);
    cudaMalloc((void **) &d_result_counts, counts_size);
 
    // Convert query data to float3 format
    std::vector<float3> query_min_vec(num_queries);
    std::vector<float3> query_max_vec(num_queries);
    for (int i = 0; i < num_queries; i++) {
        query_min_vec[i] = make_float3(h_query_aabb_min[i * 3], h_query_aabb_min[i * 3 + 1], h_query_aabb_min[i * 3 + 2]);
        query_max_vec[i] = make_float3(h_query_aabb_max[i * 3], h_query_aabb_max[i * 3 + 1], h_query_aabb_max[i * 3 + 2]);
    }
 
    // Convert primitive AABB data to float3 format
    std::vector<float3> primitive_min_vec(primitive_count);
    std::vector<float3> primitive_max_vec(primitive_count);
    for (int i = 0; i < primitive_count; i++) {
        primitive_min_vec[i] = make_float3(h_primitive_aabb_min[i * 3], h_primitive_aabb_min[i * 3 + 1], h_primitive_aabb_min[i * 3 + 2]);
        primitive_max_vec[i] = make_float3(h_primitive_aabb_max[i * 3], h_primitive_aabb_max[i * 3 + 1], h_primitive_aabb_max[i * 3 + 2]);
    }
 
    // Copy query and primitive AABB data to GPU
    cudaMemcpy(d_query_min, query_min_vec.data(), query_size, cudaMemcpyHostToDevice);
    cudaMemcpy(d_query_max, query_max_vec.data(), query_size, cudaMemcpyHostToDevice);
    cudaMemcpy(d_primitive_min, primitive_min_vec.data(), primitive_aabb_size, cudaMemcpyHostToDevice);
    cudaMemcpy(d_primitive_max, primitive_max_vec.data(), primitive_aabb_size, cudaMemcpyHostToDevice);
 
    // Launch kernel
    int block_size = 256;
    int num_blocks = (num_queries + block_size - 1) / block_size;
 
    // Reset the traversal-stack overflow flag before launch; checked after sync (fail-fast on a deeper-than-expected tree).
    const unsigned int stack_overflow_reset = 0;
    cudaMemcpyToSymbol(d_bvh_stack_overflow, &stack_overflow_reset, sizeof(unsigned int));
 
    bvhTraversalKernel<<<num_blocks, block_size>>>((GPUBVHNode *) h_nodes, (unsigned int *) h_primitive_indices, d_primitive_min, d_primitive_max, d_query_min, d_query_max, d_results, d_result_counts, num_queries, max_results_per_query);
 
    cudaDeviceSynchronize();
 
    // Check for errors
    cudaError_t err = cudaGetLastError();
    if (err != cudaSuccess) {
        fprintf(stderr, "CUDA kernel launch error: %s\n", cudaGetErrorString(err));
        // Clean up GPU memory before returning
        cudaFree(d_query_min);
        cudaFree(d_query_max);
        cudaFree(d_primitive_min);
        cudaFree(d_primitive_max);
        cudaFree(d_results);
        cudaFree(d_result_counts);
        return;
    }
 
    unsigned int stack_overflow_flag = 0;
    cudaMemcpyFromSymbol(&stack_overflow_flag, d_bvh_stack_overflow, sizeof(unsigned int));
    if (stack_overflow_flag != 0) {
        cudaFree(d_query_min);
        cudaFree(d_query_max);
        cudaFree(d_primitive_min);
        cudaFree(d_primitive_max);
        cudaFree(d_results);
        cudaFree(d_result_counts);
        helios::helios_runtime_error("ERROR (CollisionDetection GPU): BVH traversal stack overflow in bvhTraversalKernel - tree depth exceeded BVH_TRAVERSAL_STACK_CAPACITY (" + std::to_string(BVH_TRAVERSAL_STACK_CAPACITY) +
                                     "). This must not happen for a valid SAH tree; raise BVH_TRAVERSAL_STACK_CAPACITY or check the BVH build.");
    }
 
    // Copy results back
    cudaMemcpy(h_results, d_results, results_size, cudaMemcpyDeviceToHost);
    cudaMemcpy(h_result_counts, d_result_counts, counts_size, cudaMemcpyDeviceToHost);
 
    // Clean up GPU memory
    cudaFree(d_query_min);
    cudaFree(d_query_max);
    cudaFree(d_primitive_min);
    cudaFree(d_primitive_max);
    cudaFree(d_results);
    cudaFree(d_result_counts);
}
 
__global__ void intersectRegularGridKernel(const size_t num_rays, float3 *d_ray_origins, float3 *d_ray_directions, float3 grid_center, float3 grid_size, int3 grid_divisions, int primitive_count, int *d_voxel_ray_counts, float *d_voxel_path_lengths,
                                           int *d_voxel_transmitted, int *d_voxel_hit_before, int *d_voxel_hit_after, int *d_voxel_hit_inside) {
 
    size_t ray_idx = blockIdx.x * blockDim.x + threadIdx.x;
 
    if (ray_idx >= num_rays) {
        return;
    }
 
    float3 ray_origin = d_ray_origins[ray_idx];
    float3 ray_direction = d_ray_directions[ray_idx];
 
    // Calculate voxel size
    float3 voxel_size = make_float3(grid_size.x / static_cast<float>(grid_divisions.x), grid_size.y / static_cast<float>(grid_divisions.y), grid_size.z / static_cast<float>(grid_divisions.z));
 
    // Calculate grid bounds once
    float3 grid_min = make_float3(grid_center.x - 0.5f * grid_size.x, grid_center.y - 0.5f * grid_size.y, grid_center.z - 0.5f * grid_size.z);
    float3 grid_max = make_float3(grid_center.x + 0.5f * grid_size.x, grid_center.y + 0.5f * grid_size.y, grid_center.z + 0.5f * grid_size.z);
 
    // Quick ray-grid intersection test
    float t_grid_min = -1e30f, t_grid_max = 1e30f;
 
    // Check if ray intersects the entire grid first
    for (int axis = 0; axis < 3; ++axis) {
        float origin_comp = (axis == 0) ? ray_origin.x : (axis == 1) ? ray_origin.y : ray_origin.z;
        float dir_comp = (axis == 0) ? ray_direction.x : (axis == 1) ? ray_direction.y : ray_direction.z;
        float min_comp = (axis == 0) ? grid_min.x : (axis == 1) ? grid_min.y : grid_min.z;
        float max_comp = (axis == 0) ? grid_max.x : (axis == 1) ? grid_max.y : grid_max.z;
 
        if (fabsf(dir_comp) < 1e-9f) {
            if (origin_comp < min_comp || origin_comp > max_comp) {
                return; // Ray doesn't intersect grid
            }
        } else {
            float t1 = (min_comp - origin_comp) / dir_comp;
            float t2 = (max_comp - origin_comp) / dir_comp;
 
            if (t1 > t2) {
                float temp = t1;
                t1 = t2;
                t2 = temp;
            }
 
            t_grid_min = fmaxf(t_grid_min, t1);
            t_grid_max = fminf(t_grid_max, t2);
 
            if (t_grid_min > t_grid_max) {
                return; // No intersection with grid
            }
        }
    }
 
    if (t_grid_max <= 1e-6f) {
        return; // Grid is behind ray
    }
 
    // Only test voxels if ray intersects the grid
    // Test intersection with each voxel in the grid
    for (int i = 0; i < grid_divisions.x; i++) {
        for (int j = 0; j < grid_divisions.y; j++) {
            for (int k = 0; k < grid_divisions.z; k++) {
 
                // Calculate voxel AABB
                float3 voxel_min = make_float3(grid_min.x + i * voxel_size.x, grid_min.y + j * voxel_size.y, grid_min.z + k * voxel_size.z);
 
                float3 voxel_max = make_float3(voxel_min.x + voxel_size.x, voxel_min.y + voxel_size.y, voxel_min.z + voxel_size.z);
 
                // Ray-AABB intersection test with improved precision
                float t_min_x, t_max_x, t_min_y, t_max_y, t_min_z, t_max_z;
 
                // X slab - handle near-zero direction components
                if (fabsf(ray_direction.x) < 1e-9f) {
                    if (ray_origin.x < voxel_min.x || ray_origin.x > voxel_max.x) {
                        continue; // Ray is parallel and outside slab
                    }
                    t_min_x = -1e30f;
                    t_max_x = 1e30f;
                } else {
                    float inv_dir_x = 1.0f / ray_direction.x;
                    if (inv_dir_x >= 0) {
                        t_min_x = (voxel_min.x - ray_origin.x) * inv_dir_x;
                        t_max_x = (voxel_max.x - ray_origin.x) * inv_dir_x;
                    } else {
                        t_min_x = (voxel_max.x - ray_origin.x) * inv_dir_x;
                        t_max_x = (voxel_min.x - ray_origin.x) * inv_dir_x;
                    }
                }
 
                // Y slab - handle near-zero direction components
                if (fabsf(ray_direction.y) < 1e-9f) {
                    if (ray_origin.y < voxel_min.y || ray_origin.y > voxel_max.y) {
                        continue; // Ray is parallel and outside slab
                    }
                    t_min_y = -1e30f;
                    t_max_y = 1e30f;
                } else {
                    float inv_dir_y = 1.0f / ray_direction.y;
                    if (inv_dir_y >= 0) {
                        t_min_y = (voxel_min.y - ray_origin.y) * inv_dir_y;
                        t_max_y = (voxel_max.y - ray_origin.y) * inv_dir_y;
                    } else {
                        t_min_y = (voxel_max.y - ray_origin.y) * inv_dir_y;
                        t_max_y = (voxel_min.y - ray_origin.y) * inv_dir_y;
                    }
                }
 
                // Z slab - handle near-zero direction components
                if (fabsf(ray_direction.z) < 1e-9f) {
                    if (ray_origin.z < voxel_min.z || ray_origin.z > voxel_max.z) {
                        continue; // Ray is parallel and outside slab
                    }
                    t_min_z = -1e30f;
                    t_max_z = 1e30f;
                } else {
                    float inv_dir_z = 1.0f / ray_direction.z;
                    if (inv_dir_z >= 0) {
                        t_min_z = (voxel_min.z - ray_origin.z) * inv_dir_z;
                        t_max_z = (voxel_max.z - ray_origin.z) * inv_dir_z;
                    } else {
                        t_min_z = (voxel_max.z - ray_origin.z) * inv_dir_z;
                        t_max_z = (voxel_min.z - ray_origin.z) * inv_dir_z;
                    }
                }
 
                // Find intersection parameters
                float t_enter = fmaxf(fmaxf(t_min_x, t_min_y), t_min_z);
                float t_exit = fminf(fminf(t_max_x, t_max_y), t_max_z);
 
                // Check if ray intersects voxel with very stringent conditions to match CPU DDA
                // Only count intersections that are clearly inside the voxel, not just touching edges
                if (t_enter < t_exit && t_exit > 1e-5f && (t_exit - t_enter) > 1e-4f) {
 
                    // Calculate path length through voxel
                    float path_length = t_exit - t_enter;
 
                    // Handle case where ray starts inside voxel
                    if (t_enter < 0) {
                        path_length = t_exit;
                        t_enter = 0.0f;
                    }
 
                    // Only count intersections with significant path length (more restrictive)
                    if (path_length < 1e-4f) {
                        continue; // Skip this voxel
                    }
 
                    // Additional filtering: skip voxels where ray barely grazes edges
                    // Check if intersection is close to voxel boundaries (likely edge case)
                    float voxel_diag = sqrtf(voxel_size.x * voxel_size.x + voxel_size.y * voxel_size.y + voxel_size.z * voxel_size.z);
                    if (path_length < voxel_diag * 0.1f) {
                        continue; // Skip grazing intersections
                    }
 
                    // Calculate flattened voxel index
                    int voxel_idx = i * grid_divisions.y * grid_divisions.z + j * grid_divisions.z + k;
 
                    // Accumulate statistics using atomic operations
                    atomicAdd(&d_voxel_ray_counts[voxel_idx], 1);
                    atomicAdd(&d_voxel_path_lengths[voxel_idx], path_length);
 
                    // Improved geometry detection based on scene content
                    if (primitive_count == 0) {
                        // No geometry in scene - all rays are transmitted (matches CPU behavior)
                        atomicAdd(&d_voxel_transmitted[voxel_idx], 1);
                    } else {
                        // There is geometry in the scene - use improved approximation
                        // TODO: Implement actual BVH ray casting on GPU
                        // For now, use a more sophisticated approximation that considers geometry
 
                        // Calculate distance from voxel center
                        float3 voxel_center = make_float3((voxel_min.x + voxel_max.x) * 0.5f, (voxel_min.y + voxel_max.y) * 0.5f, (voxel_min.z + voxel_max.z) * 0.5f);
 
                        // Simple heuristic: rays closer to origin more likely to hit geometry
                        float ray_distance = sqrtf(ray_origin.x * ray_origin.x + ray_origin.y * ray_origin.y + ray_origin.z * ray_origin.z);
 
                        // Probability of hitting geometry decreases with distance
                        bool hit_geometry = (ray_idx % 4 == 0) && (ray_distance < 10.0f);
 
                        if (hit_geometry) {
                            // Classify hit based on ray entry position relative to voxel
                            if (t_enter < 0.5f) {
                                atomicAdd(&d_voxel_hit_inside[voxel_idx], 1);
                                atomicAdd(&d_voxel_hit_after[voxel_idx], 1);
                            } else if (t_enter < 2.0f) {
                                atomicAdd(&d_voxel_hit_after[voxel_idx], 1);
                            } else {
                                atomicAdd(&d_voxel_hit_before[voxel_idx], 1);
                            }
                        } else {
                            atomicAdd(&d_voxel_transmitted[voxel_idx], 1);
                        }
                    }
                }
            }
        }
    }
}
 
bool launchVoxelRayPathLengths(int num_rays, float *h_ray_origins, float *h_ray_directions, float grid_center_x, float grid_center_y, float grid_center_z, float grid_size_x, float grid_size_y, float grid_size_z, int grid_divisions_x,
                               int grid_divisions_y, int grid_divisions_z, int primitive_count, int *h_voxel_ray_counts, float *h_voxel_path_lengths, int *h_voxel_transmitted, int *h_voxel_hit_before, int *h_voxel_hit_after,
                               int *h_voxel_hit_inside) {
 
    // Check if GPU is available before attempting allocation
    int deviceCount = 0;
    cudaError_t err = cudaGetDeviceCount(&deviceCount);
    if (err != cudaSuccess || deviceCount == 0) {
        // No GPU available - return false for CPU fallback
        return false;
    }
 
    // Allocate device memory
    float3 *d_ray_origins, *d_ray_directions;
    int *d_voxel_ray_counts, *d_voxel_transmitted;
    int *d_voxel_hit_before, *d_voxel_hit_after, *d_voxel_hit_inside;
    float *d_voxel_path_lengths;
 
    size_t ray_data_size = num_rays * 3 * sizeof(float);
    size_t voxel_count = grid_divisions_x * grid_divisions_y * grid_divisions_z;
    size_t voxel_int_size = voxel_count * sizeof(int);
    size_t voxel_float_size = voxel_count * sizeof(float);
 
    // Allocate memory
    err = cudaMalloc(&d_ray_origins, ray_data_size);
    if (err != cudaSuccess) return false;
 
    err = cudaMalloc(&d_ray_directions, ray_data_size);
    if (err != cudaSuccess) {
        cudaFree(d_ray_origins);
        return false;
    }
 
    err = cudaMalloc(&d_voxel_ray_counts, voxel_int_size);
    if (err != cudaSuccess) {
        cudaFree(d_ray_origins);
        cudaFree(d_ray_directions);
        return false;
    }
 
    err = cudaMalloc(&d_voxel_transmitted, voxel_int_size);
    if (err != cudaSuccess) {
        cudaFree(d_ray_origins);
        cudaFree(d_ray_directions);
        cudaFree(d_voxel_ray_counts);
        return false;
    }
 
    err = cudaMalloc(&d_voxel_hit_before, voxel_int_size);
    if (err != cudaSuccess) {
        cudaFree(d_ray_origins);
        cudaFree(d_ray_directions);
        cudaFree(d_voxel_ray_counts);
        cudaFree(d_voxel_transmitted);
        return false;
    }
 
    err = cudaMalloc(&d_voxel_hit_after, voxel_int_size);
    if (err != cudaSuccess) {
        cudaFree(d_ray_origins);
        cudaFree(d_ray_directions);
        cudaFree(d_voxel_ray_counts);
        cudaFree(d_voxel_transmitted);
        cudaFree(d_voxel_hit_before);
        return false;
    }
 
    err = cudaMalloc(&d_voxel_hit_inside, voxel_int_size);
    if (err != cudaSuccess) {
        cudaFree(d_ray_origins);
        cudaFree(d_ray_directions);
        cudaFree(d_voxel_ray_counts);
        cudaFree(d_voxel_transmitted);
        cudaFree(d_voxel_hit_before);
        cudaFree(d_voxel_hit_after);
        return false;
    }
 
    err = cudaMalloc(&d_voxel_path_lengths, voxel_float_size);
    if (err != cudaSuccess) {
        cudaFree(d_ray_origins);
        cudaFree(d_ray_directions);
        cudaFree(d_voxel_ray_counts);
        cudaFree(d_voxel_transmitted);
        cudaFree(d_voxel_hit_before);
        cudaFree(d_voxel_hit_after);
        cudaFree(d_voxel_hit_inside);
        return false;
    }
 
    // Copy input data to device
    cudaMemcpy(d_ray_origins, h_ray_origins, ray_data_size, cudaMemcpyHostToDevice);
    cudaMemcpy(d_ray_directions, h_ray_directions, ray_data_size, cudaMemcpyHostToDevice);
    cudaMemset(d_voxel_ray_counts, 0, voxel_int_size);
    cudaMemset(d_voxel_transmitted, 0, voxel_int_size);
    cudaMemset(d_voxel_hit_before, 0, voxel_int_size);
    cudaMemset(d_voxel_hit_after, 0, voxel_int_size);
    cudaMemset(d_voxel_hit_inside, 0, voxel_int_size);
    cudaMemset(d_voxel_path_lengths, 0, voxel_float_size);
 
    // Launch kernel
    dim3 block_size(256);
    dim3 grid_size((num_rays + block_size.x - 1) / block_size.x);
 
    float3 grid_center = make_float3(grid_center_x, grid_center_y, grid_center_z);
    float3 grid_size_vec = make_float3(grid_size_x, grid_size_y, grid_size_z);
    int3 grid_divisions_vec = make_int3(grid_divisions_x, grid_divisions_y, grid_divisions_z);
 
    intersectRegularGridKernel<<<grid_size, block_size>>>(num_rays, d_ray_origins, d_ray_directions, grid_center, grid_size_vec, grid_divisions_vec, primitive_count, d_voxel_ray_counts, d_voxel_path_lengths, d_voxel_transmitted, d_voxel_hit_before,
                                                          d_voxel_hit_after, d_voxel_hit_inside);
 
    cudaDeviceSynchronize();
 
    // Check for errors
    err = cudaGetLastError();
    if (err != cudaSuccess) {
        // Clean up GPU memory before returning
        cudaFree(d_ray_origins);
        cudaFree(d_ray_directions);
        cudaFree(d_voxel_ray_counts);
        cudaFree(d_voxel_transmitted);
        cudaFree(d_voxel_hit_before);
        cudaFree(d_voxel_hit_after);
        cudaFree(d_voxel_hit_inside);
        cudaFree(d_voxel_path_lengths);
        return false;
    }
 
    // Copy results back to host
    cudaMemcpy(h_voxel_ray_counts, d_voxel_ray_counts, voxel_int_size, cudaMemcpyDeviceToHost);
    cudaMemcpy(h_voxel_path_lengths, d_voxel_path_lengths, voxel_float_size, cudaMemcpyDeviceToHost);
    cudaMemcpy(h_voxel_transmitted, d_voxel_transmitted, voxel_int_size, cudaMemcpyDeviceToHost);
    cudaMemcpy(h_voxel_hit_before, d_voxel_hit_before, voxel_int_size, cudaMemcpyDeviceToHost);
    cudaMemcpy(h_voxel_hit_after, d_voxel_hit_after, voxel_int_size, cudaMemcpyDeviceToHost);
    cudaMemcpy(h_voxel_hit_inside, d_voxel_hit_inside, voxel_int_size, cudaMemcpyDeviceToHost);
 
    // Free device memory
    cudaFree(d_ray_origins);
    cudaFree(d_ray_directions);
    cudaFree(d_voxel_ray_counts);
    cudaFree(d_voxel_transmitted);
    cudaFree(d_voxel_hit_before);
    cudaFree(d_voxel_hit_after);
    cudaFree(d_voxel_hit_inside);
    cudaFree(d_voxel_path_lengths);
 
    return true; // Success
}
 
} // extern "C"