__constant sampler_t sampler = CLK_FILTER_NEAREST;

// Inline definition of horizontal max
inline float max4(float a, float b, float c, float d)
{
    return max(max(max(a, b), c), d);
}

// Inline definition of horizontal min
inline float min4(float a, float b, float c, float d)
{
    return min(min(min(a, b), c), d);
}

// Traversal kernel
__kernel void traverse( __read_only image2d_t nodes,
            __global const float4* triangles,
            __global const float4* rays,
            __global float4* result,
            const int num,
            const int w,
            const int h)
{
    // Ray index
    int idx = get_global_id(0);

    if(idx < num)
    {
        // Stack
            int todo[32];
            int todoOffset = 0;

        // Current node
        int nodeNum = 0;

            float tmin = 0.0f;
            float depth = 2e30f;

            // Fetch ray origin, direction and compute invdirection
            float4 origin = rays[2 * idx + 0];
            float4 direction = rays[2 * idx + 1];
            float4 invdir = native_recip(direction);

            float4 temp = (float4)(0.0f, 0.0f, 0.0f, 1.0f);

            // Traversal loop
            while(true)
            {
            // Fetch node information
                    int2 nodeCoord = (int2)((nodeNum << 2) % w, (nodeNum << 2) / w);
                    int4 specs = read_imagei(nodes, sampler, nodeCoord + (int2)(3, 0));

            // While node isn't leaf
                    while(specs.z == 0)
                    {
                // Fetch child bounding boxes
                        float4 n0xy = read_imagef(nodes, sampler, nodeCoord);
                        float4 n1xy = read_imagef(nodes, sampler, nodeCoord + (int2)(1, 0));
                        float4 nz = read_imagef(nodes, sampler, nodeCoord + (int2)(2, 0));

                // Test ray against child bounding boxes
                        float oodx = origin.x * invdir.x;
                        float oody = origin.y * invdir.y;
                        float oodz = origin.z * invdir.z;
                        float c0lox = n0xy.x * invdir.x - oodx;
                        float c0hix = n0xy.y * invdir.x - oodx;
                        float c0loy = n0xy.z * invdir.y - oody;
                        float c0hiy = n0xy.w * invdir.y - oody;
                        float c0loz = nz.x * invdir.z - oodz;
                        float c0hiz = nz.y * invdir.z - oodz;
                        float c1loz = nz.z * invdir.z - oodz;
                        float c1hiz = nz.w * invdir.z - oodz;
                        float c0min = max4(min(c0lox, c0hix), min(c0loy, c0hiy), min(c0loz, c0hiz), tmin);
                        float c0max = min4(max(c0lox, c0hix), max(c0loy, c0hiy), max(c0loz, c0hiz), depth);
                        float c1lox = n1xy.x * invdir.x - oodx;
                        float c1hix = n1xy.y * invdir.x - oodx;
                        float c1loy = n1xy.z * invdir.y - oody;
                        float c1hiy = n1xy.w * invdir.y - oody;
                        float c1min = max4(min(c1lox, c1hix), min(c1loy, c1hiy), min(c1loz, c1hiz), tmin);
                        float c1max = min4(max(c1lox, c1hix), max(c1loy, c1hiy), max(c1loz, c1hiz), depth);

                        bool traverseChild0 = (c0max >= c0min);
                        bool traverseChild1 = (c1max >= c1min);

                nodeNum = specs.x;
                        int nodeAbove = specs.y;

                // We hit just one out of 2 childs
                        if(traverseChild0 != traverseChild1)
                        {
                                if(traverseChild1)
                                {
                                    nodeNum = nodeAbove;
                                }
                        }
                // We hit either both or none
                        else
                        {
                    // If we hit none, pop node from stack (or exit traversal, if stack is empty)
                                if (!traverseChild0)
                                {
                                    if(todoOffset == 0)
                                    {
                            break;
                                    }
                        nodeNum = todo[--todoOffset];
                                }
                    // If we hit both
                                else
                                {
                        // Sort them (so nearest goes 1st, further 2nd)
                                    if(c1min < c0min)
                                    {
                                            unsigned int tmp = nodeNum;
                                            nodeNum = nodeAbove;
                                            nodeAbove = tmp;
                                    }

                        // Push further on stack
                                    todo[todoOffset++] = nodeAbove;
                                }
                        }

                // Fetch next node information
                        nodeCoord = (int2)((nodeNum << 2) % w, (nodeNum << 2) / w);
                        specs = read_imagei(nodes, sampler, nodeCoord + (int2)(3, 0));
                    }

            // If node is leaf & has some primitives
                    if(specs.z > 0)
                    {
                // Loop through primitives & perform intersection with them (Woop triangles)
                        #pragma nounroll
                for(int i = specs.x; i < specs.y; i++)
                        {
                    // Fetch first point from global memory
                                float4 v0 = triangles[i * 4 + 0];

                                float o_z = v0.w - origin.x * v0.x - origin.y * v0.y - origin.z * v0.z;
                                float i_z = 1.0f / (direction.x * v0.x + direction.y * v0.y + direction.z * v0.z);
                                float t = o_z * i_z;

                                if(t > 0.0f && t < depth)
                                {
                        // Fetch second point from global memory
                                    float4 v1 = triangles[i * 4 + 1];

                                    float o_x = v1.w + origin.x * v1.x + origin.y * v1.y + origin.z * v1.z;
                                    float d_x = direction.x * v1.x + direction.y * v1.y + direction.z * v1.z;
                                    float u = o_x + t * d_x;

                                    if(u >= 0.0f && u <= 1.0f)
                                    {
                            // Fetch third point from global memory
                                            float4 v2 = triangles[i * 4 + 2];

                                            float o_y = v2.w + origin.x * v2.x + origin.y * v2.y + origin.z * v2.z;
                                            float d_y = direction.x * v2.x + direction.y * v2.y + direction.z * v2.z;
                                            float v = o_y + t * d_y;

                                            if(v >= 0.0f && u + v <= 1.0f)
                                            {
                                // We got successful hit, store the information
                                                depth = t;
                                                temp.x = u;
                                                temp.y = v;
                                                temp.z = t;
                                                temp.w = as_float(i);
                                            }
                                    }
                                }
                        }
                    }

            // Pop node from stack (if empty, finish traversal)
                    if(todoOffset == 0)
                    {
                        break;
                    }

                    nodeNum = todo[--todoOffset];
            }

        // Store the ray traversal result in global memory
            result[idx] = temp;
        }
}