
layout(binding=0)
   uniform sampler3D samplerData;

layout(binding=0, rgba16ui)
  uniform writeonly uimage3D outputData;

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

uniform uvec3 cubeCoordinateIn[128];
uniform uvec3 cubeCoordinateOut[128];

uniform int thresholdValue;

shared bool bit_local [4*4*4];

/////////////////////////////////////////////////////////////////////////////
// IsVoxelInRange
//

bool
IsVoxelInRange(int x, int y, int z, uvec3 cMin)
{
  if (x < 0 || y < 0 || z < 0)
  {
    return false;
  }

  ivec3
    iMin = ivec3 (cMin.x, cMin.y, cMin.z);

  if (x < iMin.x || y < iMin.y || z < iMin.z)
  {
    return false;
  }

  uvec3
    iMax = ivec3(iMin.x + 64,  iMin.y + 64,  iMin.z + 64);

  if (x >= iMax.x || y >= iMax.y || z >= iMax.z)
  {
    return false;
  }

  return true;
}

/////////////////////////////////////////////////////////////////////////////
// calcMultiKernel
//

bool
calcMultiKernel(uvec3 cCurrentCoordIn, uint x, uint y, uint z, uint i, uint j, uint k)
{

  uvec3
    cGlobalCoordIn = uvec3 (cCurrentCoordIn.x, cCurrentCoordIn.y, cCurrentCoordIn.z);

  uvec3
    cCornerTexelMin = uvec3 (cGlobalCoordIn.x + 16 * x + 4 *i, cGlobalCoordIn.y + 16 * y + 4 * j, cGlobalCoordIn.z + 16 * z + 4 * k);

  int
    border = 1;

  float
    accumulator = 0.f;

  ivec3
    iMin = ivec3 ( int (cCornerTexelMin.x) - border, int (cCornerTexelMin.y) - border, int (cCornerTexelMin.z) - border);

  ivec3   
	iMax = ivec3 ( iMin.x + 4 +  2 * border,		 iMin.y + 4 +  2 * border,		   iMin.z + 4 +  2 * border);

  // Can calculate loop limits once here, not needed in inner loop for every iteration.

  for (int z_in = iMin.z; z_in < iMax.z; z_in++)
  {
    for (int y_in = iMin.y; y_in < iMax.y; y_in++)
    {
      for (int x_in = iMin.x; x_in < iMax.x ; x_in++)
      {
        if (IsVoxelInRange(x_in, y_in, z_in, cGlobalCoordIn))
        {
		  ivec3
			samplePosition = ivec3 (x_in, y_in, z_in);

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

		  accumulator += samplerValue.x;
		  accumulator += samplerValue.y;
		  accumulator += samplerValue.z;
		  accumulator += samplerValue.w;
        }
        else
        {
          // Outside texture
        }
      }
    }
  }

  int
     a = int(accumulator * 3.4f);
 
  return (a < thresholdValue);
}

/////////////////////////////////////////////////////////////////////////////
// highPass
//

bool
highPass(uvec3 cCurrentCoordIn, uint x, uint y, uint z, uint i, uint j, uint k)
{

  uvec3
    cGlobalCoordIn = uvec3 (cCurrentCoordIn.x, cCurrentCoordIn.y, cCurrentCoordIn.z);

  uvec3
    cCornerTexelMin = uvec3 (cGlobalCoordIn.x + 16 * x + 4 *i, cGlobalCoordIn.y + 16 * y + 4 * j, cGlobalCoordIn.z + 16 * z + 4 * k);

  int
    border = 1;

  float
    highValue = thresholdValue/2000.f;

  ivec3
    iMin = ivec3 ( int (cCornerTexelMin.x) - border, int (cCornerTexelMin.y) - border, int (cCornerTexelMin.z) - border);

  ivec3   
	iMax = ivec3 (iMin.x + 4 +  2 * border, iMin.y + 4 +  2 * border, iMin.z + 4 +  2 * border);

  iMin.x = max(iMin.x, int(cGlobalCoordIn.x));
  iMax.x = min(iMax.x, int(cGlobalCoordIn.x + 64));

  iMin.y = max(iMin.y, int(cGlobalCoordIn.y));
  iMax.y = min(iMax.y, int(cGlobalCoordIn.y + 64));

  iMin.z = max(iMin.z, int(cGlobalCoordIn.z));
  iMax.z = min(iMax.z, int(cGlobalCoordIn.z + 64));


  for (int z_in = iMin.z; z_in < iMax.z; z_in++)
  {
    for (int y_in = iMin.y; y_in < iMax.y; y_in++)
    {
      for (int x_in = iMin.x; x_in < iMax.x ; x_in++)
      {
		ivec3
		  samplePosition = ivec3 (x_in, y_in, z_in);

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

		float
		  red = samplerValue.x * samplerValue.w;	

		if (red > highValue )
		{
		  return false;
		}
      }
    }
  }

  return true;
}

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

void main()
{

  // 64, 1, 1
  uint
    id = gl_WorkGroupID.x;

  uint
    iCube = gl_WorkGroupID.y;

  uvec3
     cubeIn = cubeCoordinateIn[iCube];

  uvec3
     cubeOut = cubeCoordinateOut[iCube];


  // id (0..63) => xyz (0..3, 0..3, 0..3)
  uint
    x =  id % 4u,
	y = (id /4u) % 4u,
	z = (id /16u) %4u;

  uint
    i = gl_LocalInvocationID.x,
  	j = gl_LocalInvocationID.y,
  	k = gl_LocalInvocationID.z;

  uint
    bit_offset = k * 4 * 4 + j * 4 + i;

  bit_local[bit_offset] = highPass( cubeIn, x, y, z, i, j, k); //                  calcMultiKernel(cubeIn, x, y, z, i, j, k);

  barrier();

  // Gather set bit_local into the singel 4x4x4 texel using a single thread per workgroup

  if (i == 0 && j == 0 && k == 0)
  {
	uint
	  vec[4];

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

	for (k = 0; k < 4; k++)
	{
	  for (j = 0; j < 4; j++)
	  {
		for (i = 0; i < 4; i++)
		{
		  uint
		    ibitindex = k * 4 * 4 + j * 4 + i;

          if (bit_local[ibitindex])
		  {
			uint
			  thisbit = 4 * j + i;

		    uint
			  mask = 1 << thisbit;

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

	ivec3
	  cGlobalCoordOut = ivec3 (cubeOut.x + x, cubeOut.y + y, cubeOut.z + z);

	uvec4
	  outValue = uvec4 (vec[0], vec[1], vec[2], vec[3]);
	
	imageStore(outputData, cGlobalCoordOut, outValue);
  }
}

