

layout(binding=0, rgba16ui)
#ifdef LOGIC_MODIFIER
  uniform uimage3D outputData;
#else
  uniform writeonly uimage3D outputData;
#endif

  
#define NTHREAD 32

#ifdef SPACELEAPING_TILE_SIZE

#define STORAGE_ELEMENTS (SPACELEAPING_TILE_SIZE * 64)

#else
#error SPACELEAPING_TILE_SIZE not defined
#endif

layout(local_size_x = NTHREAD, local_size_y = 4, local_size_z = 4) in;


// Specifics to this shader
uniform float thresholdValue;


//shared float arMin [STORAGE_ELEMENTS];
//shared float arMax [STORAGE_ELEMENTS];
shared bool  IsVisibleBit [STORAGE_ELEMENTS];


bool
IsVisible(uint uStorageIndex)
{
/*
  float
    rMin = arMin[uStorageIndex],
    rMax = arMax[uStorageIndex];

  return (rMin < thresholdValue);
*/
  return IsVisibleBit[uStorageIndex];

}

void
ProcessVoxelValue(uint uStorageIndex, float rValue)
{
  bool
    isAlreadyVisible = IsVisibleBit[uStorageIndex];

  bool
    isThisVisible = (rValue < thresholdValue);

  IsVisibleBit[uStorageIndex] = isAlreadyVisible || isThisVisible;

/*
  float
    rMin = arMin [uStorageIndex],
    rMax = arMax [uStorageIndex];
      
  arMin [uStorageIndex] = min (rMin, rValue);
  arMax [uStorageIndex] = max (rMax, rValue);
*/
}

void
InitStorage (uint uStorageIndex)
{
  //arMin[uStorageIndex] = FLT_MAX;
  //arMax[uStorageIndex] = -FLT_MAX;

  IsVisibleBit[uStorageIndex] = false;
}


void
ProcessCPUAlignedXKernel(uint uXIn, uint uYIn, uint uZIn, uint j, uint k, uvec3 cMin, uvec3 cMax)
{

  uint
    uBitCubeSize = uSubTileSize / (4u * uSpaceLeapingTileSize);
  
  uint
    uXReadPosition = uXIn - cMin.x;

  uint
    uCubeOutXCenter = uXReadPosition / (4u * uBitCubeSize);

  uint
    iCenter = (uXReadPosition / uBitCubeSize) % 4u;

  int
    border = 1;

  uint
    uCubeBeforeX = 0u,
    iBefore = 0u,
    uCubeAfterX = 0u,
    iAfter = 0u;

  bool
    isBefore = false,
    isAfter = false,
    isCenter = true;    // Always center

  uint
    positionInCubeX = uXReadPosition % uBitCubeSize;

  if (positionInCubeX == 0u)
  {
    if (iCenter == 0u)
    {
      if (uCubeOutXCenter > 0u)
      {
        uCubeBeforeX = uCubeOutXCenter - 1u;
        iBefore = 4u -1u;
        isBefore = true;
      }
    }
    else
    {
      uCubeBeforeX = uCubeOutXCenter;
      iBefore = iCenter - 1;
      isBefore = true;
    }
  }
          
  else if (positionInCubeX == uBitCubeSize -1u)
  {
    if (uXIn < cMax.x -1u)
    {
      if (iCenter == 4u -1u)
      {
        if (uCubeOutXCenter < uSpaceLeapingTileSize -1u)
        {
          uCubeAfterX = uCubeOutXCenter + 1u;
          iAfter = 0u;
          isAfter = true;
        }
      }
      else
      {
        iAfter = iCenter +1;
        uCubeAfterX = uCubeOutXCenter;
        isAfter = true;
      }
    }
    else
    {
      ;      
    }
  }

  bool
    isBeforeProcessing = false,
    isCenterProcessing = false,
    isAfterProcessing = false;

  if (isBefore)
  {
    uint
      uMinMaxBefore = k * 4 * 4 + j * 4 + iBefore;

    uint
      uSharedDataIndex = uCubeBeforeX * 64u + uMinMaxBefore;

    isBeforeProcessing = ! IsVisible(uSharedDataIndex);
  }

  if (isCenter)
  {
    uint
      uMinMaxCenter = k * 4 * 4 + j * 4 + iCenter;

    uint
      uSharedDataIndex = uCubeOutXCenter * 64u + uMinMaxCenter;

    isCenterProcessing =  ! IsVisible(uSharedDataIndex);

  }

  if (isAfter)
  {
    uint
      uMinMaxAfter = k * 4 * 4 + j * 4 + iAfter;

    uint
      uSharedDataIndex = uCubeAfterX * 64u + uMinMaxAfter;

    isAfterProcessing = ! IsVisible(uSharedDataIndex);

  }

  bool
    isAnyProcessing = isBeforeProcessing || isCenterProcessing || isAfterProcessing;

  if (!isAnyProcessing)
  {
    return;
  }

  ivec3
    samplePosition = ivec3 (uXIn, uYIn, uZIn);

  vec4
    sampleValue = texelFetch(samplerData, samplePosition, 0);

  float
  	rValue = sampleValue.x;

  if (isBeforeProcessing)
  {
    uint
      uMinMaxBefore = k * 4 * 4 + j * 4 + iBefore;

    uint
      uSharedDataIndex = uCubeBeforeX * 64u + uMinMaxBefore;

    ProcessVoxelValue(uSharedDataIndex, rValue);
  }

  if (isCenterProcessing)
  {
    uint
      uMinMaxCenter = k * 4 * 4 + j * 4 + iCenter;

    uint
      uSharedDataIndex = uCubeOutXCenter * 64u + uMinMaxCenter;

    ProcessVoxelValue(uSharedDataIndex, rValue);
  }

  if (isAfterProcessing)
  {
    uint
      uMinMaxAfter = k * 4 * 4 + j * 4 + iAfter;

    uint
      uSharedDataIndex = uCubeAfterX * 64u + uMinMaxAfter;

    ProcessVoxelValue(uSharedDataIndex, rValue);
  }
}

/////////////////////////////////////////////////////////////////////////////
// ProcessOutputTexel
//

void
ProcessOutputTexel(uvec3 cubeOutMin, uint uTexel, uint uY, uint uZ, uvec3 cubeInMin, uvec3 cubeInMax)
{
  uvec3
    cCoordOutTexel = uvec3 ( cubeOutMin.x + uTexel, cubeOutMin.y + uY, cubeOutMin.z + uZ);

  uint
    uBitCubeSize = uSubTileSize / (4u * uSpaceLeapingTileSize);

  uint
    vec[4];

  vec[0] = 0u;
  vec[1] = 0u;
  vec[2] = 0u;
  vec[3] = 0u;

  for (uint k = 0u; k < 4u; k++)
  {
	for (uint j = 0u; j < 4u; j++)
	{
	  for (uint i = 0u; i < 4u; i++)
	  {
		uint
		  ibitindex = k * 4u * 4u + j * 4u + i;
          
        // Check if this bit correspond to non null input volume.
        uint
          uX = uTexel;

        uvec3
           cCoordInMin = cubeInMin + uvec3 (4u * uBitCubeSize * uX + uBitCubeSize * i,
                                            4u * uBitCubeSize * uY + uBitCubeSize * j,
                                            4u * uBitCubeSize * uZ + uBitCubeSize * k);

        bool
          isDataInRange = (cCoordInMin.x < cubeInMax.x) && (cCoordInMin.y < cubeInMax.y) && (cCoordInMin.z < cubeInMax.z);

        uint
          uStorageIndex = uTexel * 64u + ibitindex;

        bool
          isVisible = isDataInRange && IsVisible(uStorageIndex);     

        bool
          isSpaceLeap = !isVisible;

        if (isSpaceLeap)
        {
		  uint
			thisbit = 4u * j + i;

		  uint
			mask = 1u << thisbit;

          vec[k] = vec[k] | mask;
        }
      }
    }
  }

  uvec4
    outValue = uvec4 (vec[0], vec[1], vec[2], vec[3]);

  ivec3
     cubeOut = ivec3 ( int(cCoordOutTexel.x), int (cCoordOutTexel.y), int (cCoordOutTexel.z));


#ifdef LOGIC_MODIFIER

    uvec4
      readValue = imageLoad(outputData, cubeOut);

    outValue.x = readValue.x LOGIC_MODIFIER outValue.x;
    outValue.y = readValue.y LOGIC_MODIFIER outValue.y;
    outValue.z = readValue.z LOGIC_MODIFIER outValue.z;
    outValue.w = readValue.w LOGIC_MODIFIER outValue.w;

#endif
	
  imageStore(outputData, cubeOut, outValue);
}


/////////////////////////////////////////////////////////////////////////////
// main
//

void main()
{
  
  // Initialize accumulators

  int
    nChunk = (STORAGE_ELEMENTS / NTHREAD) + ((STORAGE_ELEMENTS % NTHREAD) == 0? 0: 1);

  int
    iThreadID = int(gl_LocalInvocationID.x);

  for (int iChunk = 0; iChunk < nChunk; iChunk++)
  {
    int
      nOffset = NTHREAD * iChunk;
    
    if (nOffset + iThreadID < STORAGE_ELEMENTS)
    {
      uint
        uStorage = uint (nOffset + iThreadID);

      InitStorage(uStorage);      
    }
    else
    {
      ;
    }
  }

  barrier();

  uint
    uCube = gl_WorkGroupID.y;

  uint
    u1;

  uint
    u2;

  if (uCube % 2 == 0)
  {
    u1 = coordinateSetValue[uCube/2].x;
    u2  = coordinateSetValue[uCube/2].y;
  }
  else
  {
    u1 = coordinateSetValue[uCube/2].z;
    u2  = coordinateSetValue[uCube/2].w;
  }

  uvec3
    uvecInMin = uvec3 ( (u1 >> 0u) & ((1u << 9u) - 1u), (u1 >> 9u) & ((1u << 9u) - 1u),  (u1 >> (9u + 9u)) & ((1u << 9u) - 1u) );

  uint
    uLOD = (u1 >> (9u + 9u + 9u)) &  ((1u << 4u) - 1u);

  uint
    uEdgeX = (u1 >> (9u + 9u + 9u  + 4u)) & ((1u << 1u) - 1u);
  
  uint
    uEdgeYZ = (u2 >> 30u)  & (( 1u << 2u) -1u);

  uint
    uEdge = uEdgeX | (uEdgeYZ << 1u);

  uvec3
    uvecOutMin = uvec3 ( (u2 >> 0u) & (( 1 << 10u) -1u), (u2 >> 10u)  & (( 1u << 10u) -1u), (u2 >> 20u)  & (( 1u << 10u) - 1u) );

  uvec3                           
    uvecCornerSize = uvec3 (cornerTileSize[uLOD].x, cornerTileSize[uLOD].y, cornerTileSize[uLOD].z); 

  bool
    isEdgeX = (uEdge & 0x01) != 0,
    isEdgeY = (uEdge & 0x02) != 0,
    isEdgeZ = (uEdge & 0x04) != 0;

  uvec3                           
    uvecInSize = uvec3 (isEdgeX ? uvecCornerSize.x : uSubTileSize, isEdgeY ? uvecCornerSize.y : uSubTileSize, isEdgeZ ? uvecCornerSize.z : uSubTileSize);

  uvec3
    cubeInMin = uSubTileSize * uvecInMin,
    cubeInMax = cubeInMin + uvecInSize;

  uvec3
    cubeOutMin = uSpaceLeapingTileSize * uvecOutMin;

  uint
    uBitCubeSize = uSubTileSize / (4u * uSpaceLeapingTileSize);

  uint
    uID = gl_WorkGroupID.x;

  uint
    uYOut =  uID % uSpaceLeapingTileSize,
	uZOut = (uID /uSpaceLeapingTileSize) % uSpaceLeapingTileSize;

  int
    border = 1;

  // Calculate

  //for (uint k = 0; k < 4u; k++)
  uint j = gl_LocalInvocationID.y;
  //uint j = (gl_WorkGroupID.z %4u);
  {
    
    //for (uint j = 0; j < 4u; j++)
    uint k = gl_LocalInvocationID.z;
    //uint k = (gl_WorkGroupID.z/4u) % 4u;
    {
      ivec3
        cThisMin = ivec3 (0, cubeInMin.y + 4 * uBitCubeSize * uYOut + uBitCubeSize * j, cubeInMin.z + 4 * uBitCubeSize * uZOut + uBitCubeSize * k);

      uint
        nXMin = cubeInMin.x,
        nXMax = cubeInMax.x;

      uint
        nYMin = uint ( max ( cThisMin.y - border,                      int (cubeInMin.y))),
        nYMax = uint ( min ( cThisMin.y + int(uBitCubeSize) + border,  int (cubeInMax.y)));

      uint
        nZMin = uint ( max ( cThisMin.z - border,                      int (cubeInMin.z))),
        nZMax = uint ( min ( cThisMin.z + int(uBitCubeSize) + border,  int (cubeInMax.z)));

      uvec3
        cMin = uvec3 (nXMin, nYMin, nZMin),
        cMax = uvec3 (nXMax, nYMax, nZMax);
      
      uint
        uXIn = cMin.x;
      
      for (uint uZIn = cMin.z; uZIn < cMax.z; uZIn++)
      {
        for (uint uYIn = cMin.y; uYIn < cMax.y; uYIn++)
        {
          int
            nX = int(cMax.x - cMin.x);

          int
            nChunk = (nX / NTHREAD) + ((nX % NTHREAD) == 0 ? 0: 1);

          for (int iChunk = 0; iChunk < nChunk; iChunk++)
          {
            uint
              uXIn = cMin.x + uint (iChunk * NTHREAD) + uint(iThreadID);

            if (uXIn < cMax.x)
            {
              ProcessCPUAlignedXKernel(uXIn, uYIn, uZIn, j, k, cMin, cMax);
            }
          }
        }
      }
    }
  }

  barrier();

  // Output

  if ( uint(iThreadID) < uSpaceLeapingTileSize)
  {
    uint
      uMinMaxTexel = uint (iThreadID);

    ProcessOutputTexel(cubeOutMin, uMinMaxTexel, uYOut, uZOut, cubeInMin, cubeInMax);
  }
}



