Shadow caster mapのときに、やっぱガチのGPUレイトレじゃないとだめかも。って書いたのでとりあえず試しにCUDAでレイとシーンのポリゴンとの交差判定をするコードを書いてみた。
（メインの部分だけだけどカーネルのコードを最後の方に載せる。）

その動作確認として、視点からレイを飛ばして交差したところまでの距離（いわゆる深度値）を直接可視化してみた。

1024 x 1024で計算したものを縮小してある。
左右逆になっとりますね。多分レイトレ部分ではなく最終的な画像出力処理部分の問題。

シーンのデータからBVHを構築（CPU側）→BVHをトラバース（GPU側）→レイとポリゴンの交差判定（GPU側）という流れ。
上の画像の場合は、シーンのポリゴン数32、レイ数1024 x 1024 ≒ 100万で、1.7秒ほど掛かったからあまり速くない。
Geforce 9400Mで贅沢言い過ぎかな？
やはり可能ならばGLSLとかで普通のグラフィックス処理に近いやりかたで処理するようにしたほうがGPUの性能は引き出せるみたい。
まだ全然最適化はしてないけど、今の構造のままで大きく高速化する余地があるのかは分からない。
分岐しまくりなのでそれを減らせれば多少は速くなるんじゃなかろうか。
もっと新しいGPUなら特性も違うだろうしもっと色々やりようがある感じはする。

ただ、CUDAってうちの環境ではなんか不安定で、ちゃんとメモリを解放してるつもりなんだけど何回か実行してるとメモリが足りなくなったり、たまにcuInit()がCUDA_ERROR_NOT_INITIALIZEDを返してきたり（おまえがそれ返してどうすんねんw）、微妙感がある。
OpenGLでやってるときはそういう事は全くない。
やはりWindowsで開発したほうが安定性は高いのかもしれない。

ちなみにGPUのメモリが足りなくなったら、Macを再起動すれば当然解放されるんだけど、スリープしてすぐ解除するだけでも解放されるみたい。

以下、今回のカーネルのコード。
恥を忍んで全体を公開したかったんだけど、Xcodeを使ってるもんで、個人情報的にそのまま公開して大丈夫かいまいち確信がもてないので、肝となる部分だけ。
BVHのデータ構造としては、ノードはAABBそのもので、全体を二分ヒープ化して1次元配列化したものなので、インデックス計算だけで親や子を辿れるようになっている。
CPU側のコードではtraverseとintersectのキックで1セットで、これを全てのレイの処理が終わるまで複数回繰り返すようにしている。

#define IN_CUDA_KERNEL
#include "CUDARay.hpp"
#include "CUDARayHitInfo.hpp"
#include "CUDAAabb.hpp"

texture<float4, cudaTextureType1D, cudaReadModeElementType> bvhTexRef;
texture<float4, cudaTextureType1D, cudaReadModeElementType> triTexRef;

// 親ノードのIDを得る
__device__ int parent(int nodeId) {
    return (nodeId - (2 - nodeId % 2)) / 2;
}

// 左の子ノードのIDを得る
__device__ int leftChild(int nodeId) {
    return nodeId * 2 + 1;
}

// 右の子ノードのIDを得る
__device__ int rightChild(int nodeId) {
    return nodeId * 2 + 2;
}

// 近い側の子と遠い側の子を得る
__device__ void getNearAndFarChild(CUDARay &ray, int nodeId, int &near, int &far) {
    int left = leftChild(nodeId);
    int right = rightChild(nodeId);
    
    CUDAAabb left_aabb(bvhTexRef, left);
    CUDAAabb right_aabb(bvhTexRef, right);
    
    // AABBの中心とレイの始点で比較する
    float ld = (left_aabb.getCenter() - ray.origin).length2();
    float rd = (right_aabb.getCenter() - ray.origin).length2();
    
    near = (ld < rd)?left:right;
    far = (ld < rd)?right:left;
}

// AABBとレイの交差判定
__device__ bool boxtest(CUDARay &ray, float maxT, int nodeId)
{
    CUDAAabb aabb(bvhTexRef, nodeId);
    return aabb.isIntersect(ray, maxT);
}

// 葉ノードかどうか
__device__ bool isLeaf(int nodeId, int aabb_count)
{    
    return (((aabb_count + 1) / 2 - 1) <= nodeId && nodeId <= (aabb_count - 1));
}

extern "C"
{
    // トラバース処理
    __global__ void traverse(CUDARay *rays, CUDARayHitInfo *hitInfos, int *stagedCount, int n, int aabb_count)
    {
        int threadIndex = blockDim.x * blockIdx.x + threadIdx.x;
        if (!(threadIndex < n)) return;
        
        CUDARay ray = rays[threadIndex];
        float maxT = hitInfos[threadIndex].t;
        
        while (ray.current_node > -1) {
            int near, far;
            getNearAndFarChild(ray, ray.current_node, near, far);
            
            // 遠いほうの子ノードから戻ってきた場合はさらに親に戻る
            if (ray.last_node == far) {
                ray.last_node = ray.current_node;
                ray.current_node = parent(ray.current_node);
                continue;
            }
            
            // 子ノードとレイ-AABBの交差テストをする。
            // 親から降りてきた場合はnear, そうじゃない場合はfarを処理する
            int tryChild = (ray.last_node == parent(ray.current_node))?near:far;
            if (boxtest(ray, maxT, tryChild)) {
                ray.last_node = ray.current_node;
                ray.current_node = tryChild;
                if (isLeaf(ray.current_node, aabb_count)) {
                    // 葉ノードなら終了
                    atomicAdd(stagedCount, 1);
                    break;
                }
            } else {
                if (tryChild == near) {
                    // 交差判定の失敗がnearなら次はfarを処理する
                    ray.last_node = near;
                } else {
                    // 2つの子との交差判定が終わったら親に戻る
                    ray.last_node = ray.current_node;
                    ray.current_node = parent(ray.current_node);
                }
                //ray.last_node = (tryChild == near)?near:ray.current_node;
                //ray.current_node = (tryChild == near)?ray.current_node:parent(ray.current_node);
            }
        }
        
        //__syncthreads();
        rays[threadIndex] = ray;
        // このフェーズではhitInfoに書き込みはない
    }
    
    // レイ-三角形交差判定
    __global__ void intersect(CUDARay *rays, CUDARayHitInfo *hitInfos, int n, int leaf_base_id)
    {
        int threadIndex = blockDim.x * blockIdx.x + threadIdx.x;
        if (!(threadIndex < n)) return;
        
        CUDARay ray = rays[threadIndex];
        
        if (ray.current_node >= leaf_base_id) {
            CUDARayHitInfo hitInfo = hitInfos[threadIndex];
            // ray.current_nodeから三角形を取得
            CUDATriangle tri(triTexRef, ray.current_node - leaf_base_id);
            if (tri.isIntersect(ray, hitInfo)) {
                // hitInfo.tより近くで交差
                hitInfo.triangle_id = ray.current_node - leaf_base_id;
                hitInfos[threadIndex] = hitInfo;
            }
        }
    }
}