Mess of a kernel

extern "C" __global__ void otrace_kernel(void)
{
    // Get information from constant memory
    const OtraceInput& in = c_OtraceInput;
    int numRays = in.numRays;
    bool anyHit = in.anyHit;
    float4* rays = (float4*)in.rays;
    int4* results = (int4*)in.results;
    float4* nodesA = (float4*)in.nodesA;
    float4* nodesB = (float4*)in.nodesB;
    float4* nodesC = (float4*)in.nodesC;
    float4* nodesD = (float4*)in.nodesD;
    float4* trisA = (float4*)in.trisA;
    float4* trisB = (float4*)in.trisB;
    float4* trisC = (float4*)in.trisC;
    int* triIndices = (int*)in.triIndices;
    float2* texCoords = (float2*)in.texCoords;
    float3* normals = (float3*)in.normals;
    int3* vertIdx = (int3*)in.triVertIndex;
    float4* atlasInfo = (float4*)in.atlasInfo;
    int* matId = (int*)in.matId;
    float4* matInfo = (float4*)in.matInfo;
    int emissiveCount = in.emissiveNum;
    int3* emissive = (int3*)in.emissive;
    int trisCount = in.trisCount;
    int3* tris = (int3*)in.tris;
    int vertsCount = in.vertsCount;
    float3* verts = (float3*)in.verts;
    int randomSeed = in.randomSeed;
    unsigned int* triMaterialColor    = (unsigned int*)in.matColor;

    int rayidx = threadIdx.x + blockDim.x * (threadIdx.y + blockDim.y * (blockIdx.x + gridDim.x * blockIdx.y));
    if (rayidx >= numRays)
        return;

    // Cast primary ray
    int hitIndex;
    float hitU;
    float hitV;
    float hitT;
    float4 result = make_float4(0.0f, 0.0f, 0.0f, 1.0f);
    float4 primO = rays[rayidx * 2 + 0];
    float4 primD = rays[rayidx * 2 + 1];
    float4 tempDir;
    float mult_val = 1.0f;

    /*result = make_float4(-random(in.randomSeed, rayidx),
        -random(in.randomSeed, rayidx),
        -random(in.randomSeed, rayidx),
        1.0f);
    STORE_RESULT(rayidx, hitIndex, result.x, result.y, result.z);
    return;*/

    // For each sample per pixel
    for (int spp = 0; spp < SAMPLES_PER_PIXEL; spp++)
    {
        // Grab primitive origin and direction
        rays[rayidx * 2 + 0] = primO;
        rays[rayidx * 2 + 1] = primD;

        int depth = 0;
        float4 color = make_float4(0.0f, 0.0f, 0.0f, 1.0f);             // Accumulated color
        float4 reflectance = make_float4(1.0f, 1.0f, 1.0f, 1.0f);       // Accumulated reflectance
        while(true)
        {
            // If we miss, return black
            bool hit = traversal(rayidx, numRays, rays, nodesA, anyHit, results, &hitIndex, &hitT, &hitU, &hitV);
            if (hit == false || hitIndex == -1)
            {
                break;
            }

            // Get hit triangle index, compute barycentric coordinates, texture coordinates, hitpoint and normal; break on hitting light
            int tri = hitIndex;
            if (matInfo[matId[tri]].x > 0.0f)
            {
                color = color + reflectance * make_float4(matInfo[matId[tri]].x, matInfo[matId[tri]].x, matInfo[matId[tri]].x, 1.0f);
                break;
            }
            float3 bary = make_float3(hitU, hitV, 1.0f - hitU - hitV);
            float4 hitPoint = rays[rayidx * 2 + 0] + rays[rayidx * 2 + 1] * hitT * 0.99f;
            float4 hitPointBehind = rays[rayidx * 2 + 0] + rays[rayidx * 2 + 1] * hitT * 1.01f;
            float3 normal = interpolateAttribute3f(bary, normals[vertIdx[tri].x], normals[vertIdx[tri].y], normals[vertIdx[tri].z]);
            float2 tc = interpolateAttribute2f(bary, texCoords[vertIdx[tri].x], texCoords[vertIdx[tri].y], texCoords[vertIdx[tri].z]);
            // Calculate normal (always directing towards ray origin)
            float3 normalDir = dot(normal, make_float3(rays[rayidx * 2 + 1])) < 0.0f ? normal : normal * -1;
            // Sample color
            float4 f = matInfo[matId[tri]].x > 0.0f ? make_float4(matInfo[matId[tri]].x, matInfo[matId[tri]].x, matInfo[matId[tri]].x, 1.0f) :
                matInfo[matId[tri]].w == 0.0f ? fromABGR(triMaterialColor[tri]) : sample2D(bary, tc, atlasInfo[tri], t_textureAtlas);
            f.x = powf(f.x, 2.2f);
            f.y = powf(f.y, 2.2f);
            f.z = powf(f.z, 2.2f);
            // Calculate maximum reflectance
            float p = f.x > f.y && f.x > f.z ? f.x : f.y > f.z ? f.y : f.z;

            // Store previous ray dir
            float4 raydir = rays[rayidx * 2 + 1];

            // Explicit step
            if (matInfo[matId[tri]].y == 0.0f && matInfo[matId[tri]].z == 0.0f)
            {
                int emiId = (uint)randomInt(in.randomSeed, rayidx + (spp + 127) * (rayidx + (depth + 334) * (rayidx + 72))) % emissiveCount;
                float3 baryRand = make_float3(random(in.randomSeed, rayidx + (spp + 17) * (rayidx + (depth + 68) * (rayidx + 62))),
                    random(in.randomSeed, rayidx + (spp + 63) * (rayidx + (depth + 68) * (rayidx + 628))),
                    0.0f);
                if(baryRand.x + baryRand.y > 1.0f)
                {
                    baryRand.x = 1.0f - baryRand.x;
                    baryRand.y = 1.0f - baryRand.y;
                }
                baryRand.z = 1.0f - baryRand.x - baryRand.y;
                float3 sp = interpolateAttribute3f(baryRand, verts[emissive[emiId].x], verts[emissive[emiId].y], verts[emissive[emiId].z]);
                float3 l = sp - make_float3(hitPoint);
                float atten = sqrtf(1.0f / (length(l) + 1.0f));
                l = normalize(l);
                rays[rayidx * 2 + 0] = hitPoint;
                rays[rayidx * 2 + 0].w = 0.01f;
                rays[rayidx * 2 + 1] = make_float4(l, 10000.0f);
                hit = traversal(rayidx, numRays, rays, nodesA, anyHit, results, &hitIndex, &hitT, &hitU, &hitV);
                float4 expl = make_float4(0.0f, 0.0f, 0.0f, 0.0f);
                if (hit && matInfo[matId[hitIndex]].x > 0.0f)
                {
                    expl = matInfo[matId[hitIndex]].x * f * dot(normalDir, l) * atten;
                }

                color = color + reflectance * (matInfo[matId[tri]].x > 0.0f ? matInfo[matId[tri]].x : 0.0f) + reflectance * expl;
            }
            else
            {
                color = color + reflectance * (matInfo[matId[tri]].x > 0.0f ? matInfo[matId[tri]].x : 0.0f);
            }

            // Russian roulette
            if (++depth > 3)
            {
                if (random(in.randomSeed, rayidx + spp * (rayidx + depth * (rayidx + 256))) < p)
                {
                    f = f * (1.0f / p);

                    if (depth > MAX_PATH_LENGTH)
                    {
                        break;
                    }
                }
                else
                {
                    break;
                }
            }

            reflectance = reflectance * f;
            // Diffuse
            if (matInfo[matId[tri]].y == 0.0f && matInfo[matId[tri]].z == 0.0f)
            {
                float r1 = 2.0f * 3.141592654f * random(in.randomSeed, rayidx + (spp + 37) * (rayidx + (depth + 17) * (rayidx + 171)));
                float r2 = random(in.randomSeed, rayidx + (spp + 63) * (rayidx + (depth + 94) * (rayidx + 317)));
                float r2s = sqrtf(r2);
                float3 w = normalDir;
                float3 u = normalize(cross((fabsf(w.x) > 0.1f ? make_float3(0.0f, 1.0f, 0.0f) : make_float3(1.0f, 0.0f, 0.0f)), w));
                float3 v = cross(w, u);
                float3 d = normalize(u * cosf(r1) * r2s + v * sinf(r1) * r2s + w * sqrtf(1.0f - r2));
                rays[rayidx * 2 + 0] = hitPoint;
                rays[rayidx * 2 + 0].w = 0.01f;
                rays[rayidx * 2 + 1] = make_float4(d, 10000.0f);
            }
            // Ideal dielectric reflection
            else if (matInfo[matId[tri]].z == 0.0f)
            {
                rays[rayidx * 2 + 0] = hitPoint;
                rays[rayidx * 2 + 0].w = 0.01f;
                rays[rayidx * 2 + 1] = raydir - 2.0f * make_float4(normal, 0.0f) * dot(make_float4(normal, 0.0f), raydir);
            }
            // Ideal dielectric refraction
            else if (matInfo[matId[tri]].y == 0.0f)
            {
                rays[rayidx * 2 + 0] = hitPoint;
                rays[rayidx * 2 + 0].w = 0.01f;
                rays[rayidx * 2 + 1] = raydir - 2.0f * make_float4(normal, 0.0f) * dot(make_float4(normal, 0.0f), raydir);

                bool into = dot(normal, normalDir) > 0.0f;
                float nc = 1.0f;
                float nt = 1.5f;
                float nnt = into ? nc / nt : nt / nc;
                float ddn = dot(make_float3(raydir), normalDir);
                float cos2t = 1.0f - nnt * nnt * (1.0f - ddn * ddn);
                if (cos2t < 0.0f)
                {
                    continue;
                }
                float3 tdir = normalize(make_float3(raydir) * nnt - normal * ((into ? 1.0f : -1.f) * (ddn * nnt + sqrtf(cos2t))));
                //float3 tdir = normalize(make_float3(raydir));
                float a = nt - nc;
                float b = nt + nc;
                float c = 1.0f - (into ? -ddn : dot(tdir, normal));
                float Re = matInfo[matId[tri]].y + (1.0f - matInfo[matId[tri]].y) * c * c * c * c * c;
                float Tr = 1.0f - Re;
                float P = 0.25f + 0.5f * Re;
                float RP = Re / P;
                float TP = Tr / (1.0f - P);

                if (random(in.randomSeed, rayidx + (spp + 63) * (rayidx + (depth + 94) * (rayidx + 317))) < P)
                {
                    reflectance = reflectance * RP;
                }
                else
                {
                    reflectance = reflectance * TP;
                    rays[rayidx * 2 + 0] = hitPointBehind;
                    rays[rayidx * 2 + 0].w = 0.01f;
                    rays[rayidx * 2 + 1] = make_float4(tdir, 10000.0f);
                }
            }
        }

        result = result + color;

        // Until path is terminated or path maximum length is reached
        /*for (int path = 0; path < MAX_PATH_LENGTH; path++)
        {
            float4 color = make_float4(0.0f, 0.0f, 0.0f, 1.0f);
            float3 normal = make_float3(0.0f, 0.0f, 0.0f);
            float len;
            bool hit;

            tempDir = rays[rayidx * 2 + 1];

            // Check hit
            hit = traversal(rayidx, numRays, rays, nodesA, anyHit, results, &hitIndex, &hitT, &hitU, &hitV);
            if(hit)
            {
                int tri = hitIndex;
                if (tri == -1)
                {
                    // Invalid triangle hit - return black
                    color = make_float4(0.0f, 0.0f, 0.0f, 1.0f);
                }
                else
                {
                    if(matInfo[matId[tri]].w == 0.0f)
                    {
                        // We hit a light
                        color = make_float4(1.0f, 1.0f, 1.0f, 1.0f);
                        normal = make_float3(rays[rayidx * 2 + 1]);
                    }
                    else
                    {
                        // We hit a triangle
                        float3 bary = make_float3(hitU, hitV, 1.0f - hitU - hitV);

                        // Sample diffuse texture for color
                        float2 tc = interpolateAttribute2f(bary, texCoords[vertIdx[tri].x], texCoords[vertIdx[tri].y], texCoords[vertIdx[tri].z]);
                        color = sample2D(bary, tc, atlasInfo[tri], t_textureAtlas);

                        // Interpolate normal
                        normal = interpolateAttribute3f(bary, normals[vertIdx[tri].x], normals[vertIdx[tri].y], normals[vertIdx[tri].z]);

                        // Generate random value for emissive triangle ID
                        int hashA = in.randomSeed + path * spp + spp + path + rayidx * (path + threadIdx.x + 17);
                        int hashB = 0x9e3779b9u;
                        int hashC = 0x9e3779b9u;
                        jenkinsMix(hashA, hashB, hashC);
                        int emiId = (uint)hashC % emissiveCount;

                        // Generate random barycentric coordinate
                        float3 baryRand = make_float3((float)hashA * exp2f(-32), (float)hashB * exp2f(-32), 0.0f);
                        if(baryRand.x + baryRand.y > 1.0f)
                        {
                            baryRand.x = 1.0f - baryRand.x;
                            baryRand.y = 1.0f - baryRand.y;
                        }
                        baryRand.z = 1.0f - baryRand.x - baryRand.y;

                        // Calculate new origin and direction to light
                        rays[rayidx * 2 + 0] = rays[rayidx * 2 + 0] + rays[rayidx * 2 + 1] * hitT * 0.99f;
                        rays[rayidx * 2 + 0].w = 0.01f;
                        float3 pos = interpolateAttribute3f(baryRand, verts[emissive[emiId].x], verts[emissive[emiId].y], verts[emissive[emiId].z]);
                        float3 dir = pos - make_float3(rays[rayidx * 2 + 0]);
                        len = length(dir);
                        dir = normalize(dir);
                        rays[rayidx * 2 + 1] = make_float4(dir, 10000.0f);

                        // Calculate BRDF

                        color = color * max(dot(dir, normal), 0.0f);
                    }
                }
            }

            // Cast a ray against random light source position
            int hitIdx;
            hit = traversal(rayidx, numRays, rays, nodesA, anyHit, results, &hitIdx, &hitT, &hitU, &hitV);
            if (hit)
            {
                // Hit means shadow, unless emissive material
                if (hitIdx != -1)
                {
                    if(matInfo[matId[hitIdx]].w != 0.0f)
                    {
                        color *= 0.0f;
                    }
                }
            }

            // Generate random value for next step
            int hashA = in.randomSeed + path * spp + spp + path * 4 + 27 + rayidx * (path + threadIdx.x + 34);
            int hashB = 0x9e3779b9u;
            int hashC = 0x9e3779b9u;
            jenkinsMix(hashA, hashB, hashC);
            rays[rayidx * 2 + 0].w = 0.01f;
            rays[rayidx * 2 + 1] = normalize(make_float4((float)hashA * exp2f(-32), (float)hashB * exp2f(-32), (float)hashC * exp2f(-32), 0.0f));
            rays[rayidx * 2 + 1].w = 10000.0f;

            mult_val = dot(tempDir, make_float4(-normal, 0.0f));
            color.x = max(color.x, 0.0f);
            color.y = max(color.y, 0.0f);
            color.z = max(color.z, 0.0f);
            result += color * mult_val;
        }*/
    }
    result /= SAMPLES_PER_PIXEL;

    STORE_RESULT(rayidx, hitIndex, result.x, result.y, result.z);
}