#version 430
#include <shaders/materials/commons.glsl>
#include <shaders/materials/raytrace_buffers.glsl>
#include <shaders/materials/raytrace_commons.glsl>

// config

#define GRID_RES_PARTITION_Z 8
#define GRID_RES_PARTITION_Y 8
#define GRID_RES_PARTITION_X 8
layout (local_size_x = GRID_RES_PARTITION_X, local_size_y = GRID_RES_PARTITION_Y, local_size_z = GRID_RES_PARTITION_Z) in;

// pass1: count triangles. we don't have that much so this will be cheap
//        for each bucket store number of triangles that qualify
//        this is control structure with just num_triangles + buffer_offset into main buffers
// pass2: do another sweep through this. just use the offsets from previous pass
//


// Start with the simple grid for now. Even dont share calcs between lanes...
// USe 16x16x16 grid for this crap. Max triangle count 5k. Each bucket 
//

shared int grid_counters[GRID_RES_PARTITION_Z][GRID_RES_PARTITION_Y][GRID_RES_PARTITION_X];
shared int total_intersections;
shared uint shared_voxelization_group_sizes[16 * 16 * 16];

uniform int numFaces;

#define TRIBOX_VEC

#ifndef TRIBOX_VEC
#define X 0
#define Y 1
#define Z 2

#define CROSS(dest,v1,v2) \
          dest[0]=v1[1]*v2[2]-v1[2]*v2[1]; \
          dest[1]=v1[2]*v2[0]-v1[0]*v2[2]; \
          dest[2]=v1[0]*v2[1]-v1[1]*v2[0]; 

#define DOT(v1,v2) (v1[0]*v2[0]+v1[1]*v2[1]+v1[2]*v2[2])

#define SUB(dest,v1,v2) \
          dest[0]=v1[0]-v2[0]; \
          dest[1]=v1[1]-v2[1]; \
          dest[2]=v1[2]-v2[2]; 

#define FINDMINMAX(x0,x1,x2,min,max) \
  min = max = x0;   \
  if(x1<min) min=x1;\
  if(x1>max) max=x1;\
  if(x2<min) min=x2;\
  if(x2>max) max=x2;

bool planeBoxOverlap(float normal[3], float vert[3], float maxbox[3])	// -NJMP-
{
  int q;
  float vmin[3],vmax[3],v;
  for(q=X;q<=Z;q++)
  {
    v=vert[q];					// -NJMP-
    if(normal[q]>0.0f)
    {
      vmin[q]=-maxbox[q] - v;	// -NJMP-
      vmax[q]= maxbox[q] - v;	// -NJMP-
    }
    else
    {
      vmin[q]= maxbox[q] - v;	// -NJMP-
      vmax[q]=-maxbox[q] - v;	// -NJMP-
    }
  }
  if(DOT(normal,vmin)>0.0f) return false;	// -NJMP-
  if(DOT(normal,vmax)>=0.0f) return true;	// -NJMP-
  
  return false;
}


/*======================== X-tests ========================*/
#define AXISTEST_X01(a, b, fa, fb)			   \
	p0 = a*v0[Y] - b*v0[Z];			       	   \
	p2 = a*v2[Y] - b*v2[Z];			       	   \
        if(p0<p2) {min=p0; max=p2;} else {min=p2; max=p0;} \
	rad = fa * boxhalfsize[Y] + fb * boxhalfsize[Z];   \
	if(min>rad || max<-rad) return 0;

#define AXISTEST_X2(a, b, fa, fb)			   \
	p0 = a*v0[Y] - b*v0[Z];			           \
	p1 = a*v1[Y] - b*v1[Z];			       	   \
        if(p0<p1) {min=p0; max=p1;} else {min=p1; max=p0;} \
	rad = fa * boxhalfsize[Y] + fb * boxhalfsize[Z];   \
	if(min>rad || max<-rad) return 0;

/*======================== Y-tests ========================*/
#define AXISTEST_Y02(a, b, fa, fb)			   \
	p0 = -a*v0[X] + b*v0[Z];		      	   \
	p2 = -a*v2[X] + b*v2[Z];	       	       	   \
        if(p0<p2) {min=p0; max=p2;} else {min=p2; max=p0;} \
	rad = fa * boxhalfsize[X] + fb * boxhalfsize[Z];   \
	if(min>rad || max<-rad) return 0;

#define AXISTEST_Y1(a, b, fa, fb)			   \
	p0 = -a*v0[X] + b*v0[Z];		      	   \
	p1 = -a*v1[X] + b*v1[Z];	     	       	   \
        if(p0<p1) {min=p0; max=p1;} else {min=p1; max=p0;} \
	rad = fa * boxhalfsize[X] + fb * boxhalfsize[Z];   \
	if(min>rad || max<-rad) return 0;

/*======================== Z-tests ========================*/

#define AXISTEST_Z12(a, b, fa, fb)			   \
	p1 = a*v1[X] - b*v1[Y];			           \
	p2 = a*v2[X] - b*v2[Y];			       	   \
        if(p2<p1) {min=p2; max=p1;} else {min=p1; max=p2;} \
	rad = fa * boxhalfsize[X] + fb * boxhalfsize[Y];   \
	if(min>rad || max<-rad) return 0;

#define AXISTEST_Z0(a, b, fa, fb)			   \
	p0 = a*v0[X] - b*v0[Y];				   \
	p1 = a*v1[X] - b*v1[Y];			           \
        if(p0<p1) {min=p0; max=p1;} else {min=p1; max=p0;} \
	rad = fa * boxhalfsize[X] + fb * boxhalfsize[Y];   \
	if(min>rad || max<-rad) return 0;

int triBoxOverlap(float boxcenter[3],float boxhalfsize[3],float triverts[3][3])
{

  /*    use separating axis theorem to test overlap between triangle and box */
  /*    need to test for overlap in these directions: */
  /*    1) the {x,y,z}-directions (actually, since we use the AABB of the triangle */
  /*       we do not even need to test these) */
  /*    2) normal of the triangle */
  /*    3) crossproduct(edge from tri, {x,y,z}-directin) */
  /*       this gives 3x3=9 more tests */
   float v0[3],v1[3],v2[3];
//   float axis[3];
   float min,max,p0,p1,p2,rad,fex,fey,fez;		// -NJMP- "d" local variable removed
   float normal[3],e0[3],e1[3],e2[3];

   /* This is the fastest branch on Sun */
   /* move everything so that the boxcenter is in (0,0,0) */
   SUB(v0,triverts[0],boxcenter);
   SUB(v1,triverts[1],boxcenter);
   SUB(v2,triverts[2],boxcenter);

   /* compute triangle edges */
   SUB(e0,v1,v0);      /* tri edge 0 */
   SUB(e1,v2,v1);      /* tri edge 1 */
   SUB(e2,v0,v2);      /* tri edge 2 */

   /* Bullet 3:  */
   /*  test the 9 tests first (this was faster) */
   fex = abs(e0[X]);
   fey = abs(e0[Y]);
   fez = abs(e0[Z]);
   AXISTEST_X01(e0[Z], e0[Y], fez, fey);
   AXISTEST_Y02(e0[Z], e0[X], fez, fex);
   AXISTEST_Z12(e0[Y], e0[X], fey, fex);

   fex = abs(e1[X]);
   fey = abs(e1[Y]);
   fez = abs(e1[Z]);
   AXISTEST_X01(e1[Z], e1[Y], fez, fey);
   AXISTEST_Y02(e1[Z], e1[X], fez, fex);
   AXISTEST_Z0(e1[Y], e1[X], fey, fex);

   fex = abs(e2[X]);
   fey = abs(e2[Y]);
   fez = abs(e2[Z]);
   AXISTEST_X2(e2[Z], e2[Y], fez, fey);
   AXISTEST_Y1(e2[Z], e2[X], fez, fex);
   AXISTEST_Z12(e2[Y], e2[X], fey, fex);

   /* Bullet 1: */
   /*  first test overlap in the {x,y,z}-directions */
   /*  find min, max of the triangle each direction, and test for overlap in */
   /*  that direction -- this is equivalent to testing a minimal AABB around */
   /*  the triangle against the AABB */

   /* test in X-direction */
   FINDMINMAX(v0[X],v1[X],v2[X],min,max);
   if(min>boxhalfsize[X] || max<-boxhalfsize[X]) return 0;

   /* test in Y-direction */
   FINDMINMAX(v0[Y],v1[Y],v2[Y],min,max);
   if(min>boxhalfsize[Y] || max<-boxhalfsize[Y]) return 0;

   /* test in Z-direction */
   FINDMINMAX(v0[Z],v1[Z],v2[Z],min,max);
   if(min>boxhalfsize[Z] || max<-boxhalfsize[Z]) return 0;

   /* Bullet 2: */
   /*  test if the box intersects the plane of the triangle */
   /*  compute plane equation of triangle: normal*x+d=0 */
   CROSS(normal,e0,e1);
   // -NJMP- (line removed here)
   if(!planeBoxOverlap(normal,v0,boxhalfsize)) return 0;	// -NJMP-

   return 1;   /* box and triangle overlaps */
}

#else

int planeBoxOverlap(vec3 normal, vec3 vert, vec3 maxbox)	// -NJMP-
{
  vec3 vmin, vmax;

  #if 0
  {
    v=vert.x;					// -NJMP-
    if(normal.x>0.0f)
    {
      vmin.x=-maxbox.x - v;	// -NJMP-
      vmax.x= maxbox.x - v;	// -NJMP-
    }
    else
    {
      vmin.x= maxbox.x - v;	// -NJMP-
      vmax.x=-maxbox.x - v;	// -NJMP-
    }
  }
  #endif
  
  vec3 s = step(vec3(0.0), normal);          // 0.0 or 1.0
  vec3 sgn = s * 2.0 - 1.0;                  // -1.0 or 1.0 
  vmin = -sgn * maxbox - vert;
  vmax =  sgn * maxbox - vert;
  
  if(dot(normal, vmin) > 0.0f) return 0;	// -NJMP-
  if(dot(normal, vmax) >=0.0f) return 1;		// -NJMP-
  
  return 0;
}

/*======================== X-tests ========================*/
#define AXISTEST_X01(a, b, fa, fb)			   \
	p0 = a*v0.y - b*v0.z;			       	   \
	p2 = a*v2.y - b*v2.z;			       	   \
    if(p0<p2) {mi=p0; ma=p2;} else {mi=p2; ma=p0;} \
	rad = fa * boxhalfsize.y + fb * boxhalfsize.z;   \
	if(mi>rad || ma<-rad) return 0;

#define AXISTEST_X2(a, b, fa, fb)			   \
	p0 = a*v0.y - b*v0.z;			           \
	p1 = a*v1.y - b*v1.z;			       	   \
    if(p0<p1) {mi=p0; ma=p1;} else {mi=p1; ma=p0;} \
	rad = fa * boxhalfsize.y + fb * boxhalfsize.z;   \
	if(mi>rad || ma<-rad) return 0;

/*======================== Y-tests ========================*/
#define AXISTEST_Y02(a, b, fa, fb)			   \
	p0 = -a*v0.x + b*v0.z;		      	   \
	p2 = -a*v2.x + b*v2.z;	       	       	   \
    if(p0<p2) {mi=p0; ma=p2;} else {mi=p2; ma=p0;} \
	rad = fa * boxhalfsize.x + fb * boxhalfsize.z;   \
	if(mi>rad || ma<-rad) return 0;

#define AXISTEST_Y1(a, b, fa, fb)			   \
	p0 = -a*v0.x + b*v0.z;		      	   \
	p1 = -a*v1.x + b*v1.z;	     	       	   \
    if(p0<p1) {mi=p0; ma=p1;} else {mi=p1; ma=p0;} \
	rad = fa * boxhalfsize.x + fb * boxhalfsize.z;   \
	if(mi>rad || ma<-rad) return 0;

/*======================== Z-tests ========================*/

#define AXISTEST_Z12(a, b, fa, fb)			   \
	p1 = a*v1.x - b*v1.y;			           \
	p2 = a*v2.x - b*v2.y;			       	   \
    if(p2<p1) {mi=p2; ma=p1;} else {mi=p1; ma=p2;} \
	rad = fa * boxhalfsize.x + fb * boxhalfsize.y;   \
	if(mi>rad || ma<-rad) return 0;

#define AXISTEST_Z0(a, b, fa, fb)			   \
	p0 = a*v0.x - b*v0.y;				   \
	p1 = a*v1.x - b*v1.y;			           \
    if(p0<p1) {mi=p0; ma=p1;} else {mi=p1; ma=p0;} \
	rad = fa * boxhalfsize.x + fb * boxhalfsize.y;   \
	if(mi>rad || ma<-rad) return 0;

int triBoxOverlap(vec3 boxcenter, vec3 boxhalfsize, vec3 tv0, vec3 tv1, vec3 tv2)
{

  /*    use separating axis theorem to test overlap between triangle and box */
  /*    need to test for overlap in these directions: */
  /*    1) the {x,y,z}-directions (actually, since we use the AABB of the triangle */
  /*       we do not even need to test these) */
  /*    2) normal of the triangle */
  /*    3) crossproduct(edge from tri, {x,y,z}-directin) */
  /*       this gives 3x3=9 more tests */
   vec3 v0, v1, v2;
   float mi,ma,p0,p1,p2,rad,fex,fey,fez;		// -NJMP- "d" local variable removed
   vec3 normal, e0, e1, e2;

   /* This is the fastest branch on Sun */
   /* move everything so that the boxcenter is in (0,0,0) */
   v0 = tv0 - boxcenter;
   v1 = tv1 - boxcenter;
   v2 = tv2 - boxcenter;

   /* compute triangle edges */
   e0 = v1 - v0;	/* tri edge 0 */	
   e1 = v2 - v1;	/* tri edge 1 */
   e2 = v0 - v2;	/* tri edge 2 */

   /* Bullet 3:  */
   /*  test the 9 tests first (this was faster) */
   fex = abs(e0.x);
   fey = abs(e0.y);
   fez = abs(e0.z);
   AXISTEST_X01(e0.z, e0.y, fez, fey);
   AXISTEST_Y02(e0.z, e0.x, fez, fex);
   AXISTEST_Z12(e0.y, e0.x, fey, fex);

   fex = abs(e1.x);
   fey = abs(e1.y);
   fez = abs(e1.z);
   AXISTEST_X01(e1.z, e1.y, fez, fey);
   AXISTEST_Y02(e1.z, e1.x, fez, fex);
   AXISTEST_Z0(e1.y, e1.x, fey, fex);

   fex = abs(e2.x);
   fey = abs(e2.y);
   fez = abs(e2.z);
   AXISTEST_X2(e2.z, e2.y, fez, fey);
   AXISTEST_Y1(e2.z, e2.x, fez, fex);
   AXISTEST_Z12(e2.y, e2.x, fey, fex);

   /* Bullet 1: */
   /*  first test overlap in the {x,y,z}-directions */
   /*  find min, max of the triangle each direction, and test for overlap in */
   /*  that direction -- this is equivalent to testing a minimal AABB around */
   /*  the triangle against the AABB */

   /* test in X-direction */
   vec3 mi3 = min(v0, min(v1, v2));
   vec3 ma3 = max(v0, max(v1, v2));
   bool no_overlap = any(greaterThan(mi3, boxhalfsize)) || any(lessThan(ma3, -boxhalfsize));
   if (no_overlap)
	return 0;
   //if(mi3.x>boxhalfsize.x || ma3.x<-boxhalfsize.x) return 0;

   /* test in Y-direction */
   //if(mi3.y>boxhalfsize.y || ma3.y<-boxhalfsize.y) return 0;

   /* test in Z-direction */
   //if(mi3.z>boxhalfsize.z || ma3.z<-boxhalfsize.z) return 0;

   /* Bullet 2: */
   /*  test if the box intersects the plane of the triangle */
   /*  compute plane equation of triangle: normal*x+d=0 */
   normal = cross(e0, e1);
   // -NJMP- (line removed here)
   
   return planeBoxOverlap(normal,v0,boxhalfsize);
}
#endif

struct ObjectBox
{
	vec3 pos;
	vec3 size;
};

bool checkTriCellOverlap(vec3 p0, vec3 p1, vec3 p2, vec3 ccenter, vec3 csize)
{
#ifndef TRIBOX_VEC
	float triverts[3][3];
	triverts[0][0] = p0.x;
	triverts[0][1] = p0.y;
	triverts[0][2] = p0.z;

	triverts[1][0] = p1.x;
	triverts[1][1] = p1.y;
	triverts[1][2] = p1.z;
	
	triverts[2][0] = p2.x;
	triverts[2][1] = p2.y;
	triverts[2][2] = p2.z;

	float boxcenter[3];
	boxcenter[0] = ccenter.x;
	boxcenter[1] = ccenter.y;
	boxcenter[2] = ccenter.z;

	float boxhalfsize[3];
	boxhalfsize[0] = csize.x;
	boxhalfsize[1] = csize.y;
	boxhalfsize[2] = csize.z;
	return triBoxOverlap(boxcenter, boxhalfsize, triverts) != 0;
#else
	return triBoxOverlap(ccenter, csize, p0, p1, p2) != 0;
#endif
}






#define laneID gl_LocalInvocationIndex
//#define gridID (gl_LocalInvocationID.x + gl_LocalInvocationID.y * GRID_RES + gl_LocalInvocationID.z * GRID_RES * GRID_RES)
//#define gridID3 ivec3(gl_LocalInvocationID.x, gl_LocalInvocationID.y, gl_LocalInvocationID.z)

#define GRID_SIZE grid_size

#ifndef USE_LINKED_LISTS
void main()
{
	int grid_min_z;
    int grid_max_z; // exclusive
    int grid_min_y;
    int grid_max_y; // exclusive
	int grid_min_x;
	int grid_max_x; // exclusive

	uint workgroup_index = gl_WorkGroupID.x + gl_WorkGroupID.y * gl_NumWorkGroups.x + gl_WorkGroupID.z * gl_NumWorkGroups.x * gl_NumWorkGroups.y;

    // copy to not use memory reads
    vec3 grid_size = in_bbox_data.grid_size_raytrace;
    ivec3 grid_shift = in_bbox_data.grid_shift_raytrace;

    // actual pass

    int grid_x = int(gl_WorkGroupID.x) * GRID_RES_PARTITION_X;
	int grid_y = int(gl_WorkGroupID.y) * GRID_RES_PARTITION_Y;
	int grid_z = int(gl_WorkGroupID.z) * GRID_RES_PARTITION_Z;

	grid_min_x = grid_x;
    grid_max_x = grid_x + GRID_RES_PARTITION_X;

    {
        grid_min_y = grid_y;
        grid_max_y = grid_y + GRID_RES_PARTITION_Y;

        {   
		    grid_min_z = grid_z;
            grid_max_z = grid_z + GRID_RES_PARTITION_Z;

            uint cell_id_base = GRID_RES * GRID_RES * grid_z + GRID_RES * grid_y + grid_x;
			uint cell_id_local = GRID_RES * GRID_RES * gl_LocalInvocationID.z + GRID_RES * gl_LocalInvocationID.y + gl_LocalInvocationID.x;

			#ifdef PASS_2
			if (in_buckets.voxelization_group_sizes[workgroup_index] == 0)
			{
				in_buckets.sizes16[cell_id_base + cell_id_local] = uint16_t(0);
				return;
			}
			#endif

    
            #define build_cell_mi ivec3 cell_mi = ivec3(floor(mi / GRID_SIZE) + grid_shift);
            #define build_cell_ma ivec3 cell_ma = ivec3(ceil(ma / GRID_SIZE) + grid_shift);

		#ifdef PASS_1
	        grid_counters[gl_LocalInvocationID.z][gl_LocalInvocationID.y][gl_LocalInvocationID.x] = 0;

            // cells centered around 0,0,0
            for(int i = int(laneID); i < numFaces; i+= GRID_RES_PARTITION_X * GRID_RES_PARTITION_Y * GRID_RES_PARTITION_Z)
            {
		        vec3 p1 = rt_get_vertex(i * 3 + 0);
		        vec3 p2 = rt_get_vertex(i * 3 + 1);
		        vec3 p3 = rt_get_vertex(i * 3 + 2);

                // bbox
                vec3 mi = min(p1, min(p2, p3));
                vec3 ma = max(p1, max(p2, p3));

                build_cell_mi
                build_cell_ma

                if (cell_ma.x < grid_min_x || cell_ma.y < grid_min_y || cell_ma.z < grid_min_z || cell_mi.x > grid_max_x || cell_mi.y > grid_max_y || cell_mi.z > grid_max_z)
                    continue;

                cell_mi = clamp(cell_mi, ivec3(grid_min_x, grid_min_y, grid_min_z), ivec3(grid_max_x, grid_max_y, grid_max_z));
                cell_ma = clamp(cell_ma, ivec3(grid_min_x, grid_min_y, grid_min_z), ivec3(grid_max_x, grid_max_y, grid_max_z));

                for(int cz = cell_mi.z; cz < cell_ma.z; cz++)
                {
                    for(int cy = cell_mi.y; cy < cell_ma.y; cy++)
                    {
                        for(int cx = cell_mi.x; cx < cell_ma.x; cx++)
                        {
                            if (cx < grid_min_x || cy < grid_min_y || cz < grid_min_z || cx > grid_max_x || cy > grid_max_y || cz > grid_max_z)
                                continue;

                            vec3 cellPos = vec3(ivec3(cx, cy, cz) - grid_shift) * GRID_SIZE;
			                vec3 cellSize = vec3(GRID_SIZE);

                            if (!checkTriCellOverlap(p1, p2, p3, cellPos + cellSize * 0.5, cellSize * 0.5))
                                continue;

                            // #define gridID (gl_LocalInvocationID.x + gl_LocalInvocationID.y * GRID_RES + sliceID * GRID_RES * GRID_RES)
                            atomicAdd(grid_counters[cz - grid_z][cy - grid_y][cx - grid_x], 1);
                        }
                    }
                }
            }

	        memoryBarrierShared();
	        barrier();

	        // calculate bucket offsets.
            // in the first pass we just count. then store sizes + indices and also sum all the buckets for this group
            // so in second pass we can update indices in other groups.
            {   
                int cell_count = grid_counters[gl_LocalInvocationID.z][gl_LocalInvocationID.y][gl_LocalInvocationID.x];
		        in_buckets.sizes[cell_id_base + cell_id_local] = cell_count;
            }

            memoryBarrierShared();
	        barrier();

	        if (laneID == 0)
	        {
                int indices_offset = 0;
                for(uint gz = 0; gz < GRID_RES_PARTITION_Z; gz++)
                {
                    for(uint gy = 0; gy < GRID_RES_PARTITION_Y; gy++)
                    {
                        for(uint gx = 0; gx < GRID_RES_PARTITION_X; gx++)
						{
							int cell_count = grid_counters[gz][gy][gx];
							uint cell_id = gx + gy * GRID_RES + gz * GRID_RES * GRID_RES;
							in_buckets.offsets[cell_id_base + cell_id] = indices_offset;
							indices_offset += cell_count;
						}
					}
				}

		        in_buckets.total_intersections = 0;
                in_buckets.voxelization_group_sizes[workgroup_index] = indices_offset;
	        }

	        memoryBarrierShared();
	        barrier();

        #endif // PASS_1

        #ifdef PASS_2
            // Move data back into groupshared, update all offsets to include sizes from all previous groups
            
			{
				// first new base offset. This sucks a bit...

				uint indices_offset = 0;
				
				//memoryBarrierShared();
				//barrier();

				// parallel sum. this one is slower for the whole run.
				// probably because some threads only do sum just a few elements. maybe manual unrolling would work
				// but probably not worth it
				#if 0

				uint workgroups_num = (GRID_RES / GRID_RES_PARTITION_X) * (GRID_RES / GRID_RES_PARTITION_Y) * (GRID_RES / GRID_RES_PARTITION_Z);
				for(uint i = laneID; i < workgroups_num; i += GRID_RES_PARTITION_X * GRID_RES_PARTITION_Y * GRID_RES_PARTITION_Z)
				{
					if (i < workgroup_index)
						shared_voxelization_group_sizes[i] = in_buckets.voxelization_group_sizes[i];
					else
						shared_voxelization_group_sizes[i] = 0;
				}

				const uint steps = uint(log2(workgroups_num));
				
				// Reduce
				uint stride = workgroups_num / 2;	// will always be even
				for (uint step = 0; step < steps; step++)
				{
					for(uint i = laneID; i < stride; i+= GRID_RES_PARTITION_X * GRID_RES_PARTITION_Y * GRID_RES_PARTITION_Z)
					{
						shared_voxelization_group_sizes[i] = shared_voxelization_group_sizes[i] + shared_voxelization_group_sizes[i + stride];
					}

					stride >>= 1;

					memoryBarrierShared();
					barrier();
				}
				
				indices_offset = shared_voxelization_group_sizes[0];
				#else

				for(uint g = 0; g < workgroup_index; g++)
				{
					indices_offset += in_buckets.voxelization_group_sizes[g];
				}

				memoryBarrierShared();
				barrier();

				#endif

				// add this base offset to all of our offsets

				grid_counters[gl_LocalInvocationID.z][gl_LocalInvocationID.y][gl_LocalInvocationID.x] = int(in_buckets.offsets[cell_id_base + cell_id_local] + indices_offset);
				in_buckets.offsets[cell_id_base + cell_id_local] = grid_counters[gl_LocalInvocationID.z][gl_LocalInvocationID.y][gl_LocalInvocationID.x];
			}

	        // actually move triangle into the buckets

	        for(int i = int(laneID); i < numFaces; i+= GRID_RES_PARTITION_X * GRID_RES_PARTITION_Y * GRID_RES_PARTITION_Z)
	        {
		        vec3 p1 = rt_get_vertex(i * 3 + 0);
		        vec3 p2 = rt_get_vertex(i * 3 + 1);
		        vec3 p3 = rt_get_vertex(i * 3 + 2);

                // bbox
                vec3 mi = min(p1, min(p2, p3));
                vec3 ma = max(p1, max(p2, p3));

                build_cell_mi
                build_cell_ma
  
                 if (cell_ma.x < grid_min_x || cell_ma.y < grid_min_y || cell_ma.z < grid_min_z || cell_mi.x > grid_max_x || cell_mi.y > grid_max_y || cell_mi.z > grid_max_z)
                    continue;

                cell_mi = clamp(cell_mi, ivec3(grid_min_x, grid_min_y, grid_min_z), ivec3(grid_max_x, grid_max_y, grid_max_z));
                cell_ma = clamp(cell_ma, ivec3(grid_min_x, grid_min_y, grid_min_z), ivec3(grid_max_x, grid_max_y, grid_max_z));

                for(int cz = cell_mi.z; cz < cell_ma.z; cz++)
                {
                    for(int cy = cell_mi.y; cy < cell_ma.y; cy++)
                    {
                        for(int cx = cell_mi.x; cx < cell_ma.x; cx++)
                        {
                            if (cx < grid_min_x || cy < grid_min_y || cz < grid_min_z || cx > grid_max_x || cy > grid_max_y || cz > grid_max_z)
                                continue;

                            vec3 cellPos = vec3(ivec3(cx, cy, cz) - grid_shift) * GRID_SIZE;
			                vec3 cellSize = vec3(GRID_SIZE);

                            if (!checkTriCellOverlap(p1, p2, p3, cellPos + cellSize * 0.5, cellSize * 0.5))
                                continue;

                            int insert_pos = atomicAdd(grid_counters[cz - grid_z][cy - grid_y][cx - grid_x], 1);
                            atomicExchange(in_buckets.indices[insert_pos], i);
                        }
                    }
                }
	        }

	        memoryBarrier();
	        barrier();

            in_buckets.offsets16[cell_id_base + cell_id_local] = uint16_t(in_buckets.offsets[cell_id_base + cell_id_local]);
            in_buckets.sizes16[cell_id_base + cell_id_local] = uint16_t(in_buckets.sizes[cell_id_base + cell_id_local]);

		#endif // PASS_2
        }
    }

    memoryBarrier();
	barrier();
}




#else // -- linked lists version ---------------------------------------------





void main()
{
	int grid_min_z;
    int grid_max_z; // exclusive
    int grid_min_y;
    int grid_max_y; // exclusive
	int grid_min_x;
	int grid_max_x; // exclusive

	uint workgroup_index = gl_WorkGroupID.x + gl_WorkGroupID.y * gl_NumWorkGroups.x + gl_WorkGroupID.z * gl_NumWorkGroups.x * gl_NumWorkGroups.y;
    
    // we are responsible to initialize linked list before we attempt to use it. this is another pass of compute
    #ifdef PASS_1

    int grid_x = int(gl_WorkGroupID.x) * GRID_RES_PARTITION_X;
	int grid_y = int(gl_WorkGroupID.y) * GRID_RES_PARTITION_Y;
	int grid_z = int(gl_WorkGroupID.z) * GRID_RES_PARTITION_Z;

    // this condition should always be true
    uint cell_id_base = GRID_RES * GRID_RES * grid_z + GRID_RES * grid_y + grid_x;
	uint cell_id_local = GRID_RES * GRID_RES * gl_LocalInvocationID.z + GRID_RES * gl_LocalInvocationID.y + gl_LocalInvocationID.x;

    {
        // fill in histogram info from previous frame data
        uint size = in_buckets.sizes[cell_id_base + cell_id_local];
        uint bucket = (size > 0) ? uint(log2(float(size))) + 1 : 0;

		// TODO: make this conditional. it is really slowing us down... (0.8ms only for this)
        //atomicAdd(in_buckets.sizes_histogram[bucket], 1);

        link_list_init(int(cell_id_base + cell_id_local));
        in_faces_list_data.buffer_counter = GRID_RES * GRID_RES * GRID_RES;
        in_buckets.sizes[cell_id_base + cell_id_local] = 0;

		mark_high_level_cell_empty(int(grid_x + gl_LocalInvocationID.x), int(grid_y + gl_LocalInvocationID.y), int(grid_z + gl_LocalInvocationID.z));
    }
    
    memoryBarrierShared();
	barrier();

    #else

    // copy to not use memory reads
    vec3 grid_size = in_buckets.grid_size;
	vec3 cellSize = vec3(GRID_SIZE);

    // actual pass

    int grid_x = int(gl_WorkGroupID.x) * GRID_RES_PARTITION_X;
	int grid_y = int(gl_WorkGroupID.y) * GRID_RES_PARTITION_Y;
	int grid_z = int(gl_WorkGroupID.z) * GRID_RES_PARTITION_Z;

	grid_min_x = grid_x;
    grid_max_x = grid_x + GRID_RES_PARTITION_X;

    {
        grid_min_y = grid_y;
        grid_max_y = grid_y + GRID_RES_PARTITION_Y;

        {   
		    grid_min_z = grid_z;
            grid_max_z = grid_z + GRID_RES_PARTITION_Z;

            uint cell_id_base = GRID_RES * GRID_RES * grid_z + GRID_RES * grid_y + grid_x;
			uint cell_id_local = GRID_RES * GRID_RES * gl_LocalInvocationID.z + GRID_RES * gl_LocalInvocationID.y + gl_LocalInvocationID.x;

            #define build_cell_mi ivec3 cell_mi = ivec3(floor(mi / GRID_SIZE));
            #define build_cell_ma ivec3 cell_ma = ivec3(ceil(ma / GRID_SIZE));

            // cells centered around 0,0,0
            for(int i = int(laneID); i < numFaces; i+= GRID_RES_PARTITION_X * GRID_RES_PARTITION_Y * GRID_RES_PARTITION_Z)
            {
		        vec3 p1 = rt_get_vertex(i * 3 + 0) - in_bbox_data.bbox_raytrace_min.xyz;
		        vec3 p2 = rt_get_vertex(i * 3 + 1) - in_bbox_data.bbox_raytrace_min.xyz;
		        vec3 p3 = rt_get_vertex(i * 3 + 2) - in_bbox_data.bbox_raytrace_min.xyz;

                // bbox
                vec3 mi = min(p1, min(p2, p3));
                vec3 ma = max(p1, max(p2, p3));

                build_cell_mi
                build_cell_ma

                if (cell_ma.x < grid_min_x || cell_ma.y < grid_min_y || cell_ma.z < grid_min_z || cell_mi.x > grid_max_x || cell_mi.y > grid_max_y || cell_mi.z > grid_max_z)
                    continue;

                cell_mi = clamp(cell_mi, ivec3(grid_min_x, grid_min_y, grid_min_z), ivec3(grid_max_x, grid_max_y, grid_max_z));
                cell_ma = clamp(cell_ma, ivec3(grid_min_x, grid_min_y, grid_min_z), ivec3(grid_max_x, grid_max_y, grid_max_z));

                for(int cz = cell_mi.z; cz <= cell_ma.z; cz++)
                {
                    for(int cy = cell_mi.y; cy <= cell_ma.y; cy++)
                    {
                        for(int cx = cell_mi.x; cx <= cell_ma.x; cx++)
                        {
                            //if (cx < grid_min_x || cy < grid_min_y || cz < grid_min_z || cx > grid_max_x || cy > grid_max_y || cz > grid_max_z)
                            //    continue;

                            vec3 cellPos = vec3(ivec3(cx, cy, cz)) * GRID_SIZE;

                            if (!checkTriCellOverlap(p1, p2, p3, cellPos + cellSize * 0.5, cellSize * 0.5))
                                continue;

                            // #define gridID (gl_LocalInvocationID.x + gl_LocalInvocationID.y * GRID_RES + sliceID * GRID_RES * GRID_RES)
                            {
                                //uint list_index = (cx - grid_x) + (cy - grid_y) * GRID_RES + (cz - grid_z) * GRID_RES * GRID_RES;
                                uint list_index = cx + cy * GRID_RES + cz * GRID_RES * GRID_RES;
                                link_list_push_value(list_index, i);

								mark_high_level_cell_occupied(cx, cy, cz);
                            }
                        }
                    }
                }
            }

	        memoryBarrierShared();
	        barrier();
        }
    }

    #endif // PASS_1

    memoryBarrier();
	barrier();
}

#endif

