unsigned char StorageIndex(unsigned char childIndex, unsigned char childMask)
{
  if (childIndex == 0) return 0;
  unsigned char storageIndex = (childMask & 1);
  if (childIndex == 1) return storageIndex; storageIndex += !!(childMask & 2);
  if (childIndex == 2) return storageIndex; storageIndex += !!(childMask & 4);
  if (childIndex == 3) return storageIndex; storageIndex += !!(childMask & 8);
  if (childIndex == 4) return storageIndex; storageIndex += !!(childMask & 16);
  if (childIndex == 5) return storageIndex; storageIndex += !!(childMask & 32);
  if (childIndex == 6) return storageIndex; storageIndex += !!(childMask & 64);
  // assume /*if (childIndex == 7) */
  return storageIndex;
}

float RayBoxIntersects(float3 rayPos, float3 rayDir, float3 boxMin, float3 boxMax)
{
  const float largestFloat = 9999999999.f;
  const float smallestFloat = 0.000000001f;
  float3 invDir = (float3)(largestFloat, largestFloat, largestFloat);
  if (fabs(rayDir.x) > smallestFloat) invDir.x = (float)1.0f / rayDir.x;
  if (fabs(rayDir.y) > smallestFloat) invDir.y = (float)1.0f / rayDir.y;
  if (fabs(rayDir.z) > smallestFloat) invDir.z = (float)1.0f / rayDir.z;
  float3 mins = ((float3)(boxMin) - rayPos) * invDir; float3 maxs = ((float3)(boxMax) - rayPos) * invDir;
  float tmax = min(min(max(mins[0], maxs[0]), max(mins[1], maxs[1])), max(mins[2], maxs[2]));
  if (tmax < 0) return largestFloat;
  float tmin = max(max(min(mins[0], maxs[0]), min(mins[1], maxs[1])), min(mins[2], maxs[2]));
  if (tmin <= tmax) return max(tmin, 0.f);
  return largestFloat;
}

kernel void Draw(__global unsigned int *pPixels, __global int *pTreeDepth, __global unsigned int *pChildData, __global unsigned int *pColorData, __global unsigned char *pMasksData, __global float *pRayDirs, __global float *pCamPos, __global float *pRotData)
{
  const bool colorBySteps = 0;

  const float largestFloat = 9999999999.f;
  const float smallestFloat = 0.000000001f;

  // Rotate Rays
  int id = get_global_id(0);
  float3 rayPos = (float3)(pCamPos[0], pCamPos[1], pCamPos[2]);
  float3 rayDir = (float3)(pRayDirs[id * 3 + 0], pRayDirs[id * 3 + 1], pRayDirs[id * 3 + 2]);
  rayDir = (float3)(rayDir.x, rayDir.y * pRotData[2] - rayDir.z * pRotData[3], rayDir.y * pRotData[3] + rayDir.z * pRotData[2]);
  rayDir = (float3)(rayDir.x * pRotData[0] - rayDir.y * pRotData[1], rayDir.x * pRotData[1] + rayDir.y * pRotData[0], rayDir.z);

  const int stackSize = 22;
  int stackPtr = 1;
  int descents = 0;

  float3 stackRegionPos[stackSize];
  unsigned int stackNodeID[stackSize];
  unsigned char stackDepth[stackSize];

  // Push Root Node
  stackNodeID[0] = 0;
  stackDepth[0] = *pTreeDepth;
  stackRegionPos[0] = (float3)(0, 0, 0);

  while (stackPtr > 0)
  {
    --stackPtr;
    float3 regionPos = stackRegionPos[stackPtr];
    unsigned int nodeID = stackNodeID[stackPtr];
    unsigned char depth = stackDepth[stackPtr];

    unsigned char childMask = pMasksData[nodeID];
    unsigned int parentID = pChildData[nodeID];

    for (int z = 0; z < 2; z++) for (int y = 0; y < 2; y++) for (int x = 0; x < 2; x++)
    {
      unsigned char c = 0;
      // Use HERO when on the last level of the trees
      if (depth == 0) c = (x ^ (rayDir.x < 0)) + (y ^ (rayDir.y < 0)) * 2 + (z ^ (rayDir.z < 0)) * 4; // HERO Algorithm
      if (depth > 0) c = (x ^ (rayDir.x > 0)) + (y ^ (rayDir.y > 0)) * 2 + (z ^ (rayDir.z > 0)) * 4; // Reverse HERO on higher layers (since were making a reversed stack of work todo)

      if (((childMask >> c) & 1) == 0) continue; // Skip empty children
      float3 ep = (float3)((c & 1) > 0, (c & 2) > 0, (c & 4) > 0); // Calculate Nodes Position
      float3 subRegionPos = (float3)(regionPos * 2) + ep;

      // Calculate Ray To Node Distance
      bool nearPixelSize = false;
      if (1)
      {
        float distToNode = RayBoxIntersects(rayPos, rayDir, (float3)(subRegionPos * (1 << depth)), (float3)((subRegionPos + (float3)(1,1,1)) * (1 << depth)));
        if (distToNode == largestFloat) continue; // Check we hit it
      }

      int childID = parentID + StorageIndex(c, childMask); // Calculate Child Index

      descents++;

      // Reached Bottom (color)
      if (depth == 0 || nearPixelSize)
      {
        if (colorBySteps) // Color By Steps
        {
          int dist = descents; dist = min(dist, 255); pPixels[id] = dist | (dist << 8) | (dist << 16);
        }
        else pPixels[id] = pColorData[childID]; // Color By Surface Hit
        return; // Finished!
      }

      stackRegionPos[stackPtr] = subRegionPos;
      stackNodeID[stackPtr] = childID;
      stackDepth[stackPtr] = depth - 1;
      ++stackPtr;
    }
  }
  if (colorBySteps) // Color By Steps
  {
    int dist = descents; dist = min(dist, 255); pPixels[id] = dist | (dist << 8) | (dist << 16);
  }
}