// bound
grid_aabb[0][0] = grid_aabb[0][1] = grid_aabb[0][2] = FLT_MAX;
grid_aabb[1][0] = grid_aabb[1][1] = grid_aabb[1][2] = FLT_MIN;
for (std::size_t i = 0; i < number_of_vertices; ++i) {
       if (vertices[i][0] < grid_aabb[0][0]) grid_aabb[0][0] = vertices[i][0];
  else if (vertices[i][0] > grid_aabb[1][0]) grid_aabb[1][0] = vertices[i][0];
       if (vertices[i][1] < grid_aabb[0][1]) grid_aabb[0][1] = vertices[i][1];
  else if (vertices[i][1] > grid_aabb[1][1]) grid_aabb[1][1] = vertices[i][1];
       if (vertices[i][2] < grid_aabb[0][2]) grid_aabb[0][2] = vertices[i][2];
  else if (vertices[i][2] > grid_aabb[1][2]) grid_aabb[1][2] = vertices[i][2];
}

// compute grid size
float grid_density = 4.0;
float grid_volume = (grid_aabb[1][0] - grid_aabb[0][0]) * (grid_aabb[1][1] - grid_aabb[0][1])
  * (grid_aabb[1][2] - grid_aabb[0][2]);
float cell_density = std::pow((grid_density * number_of_triangles) / grid_volume, 1.0f / 3.0f);
grid_size[0] = std::size_t(((grid_aabb[1][0] - grid_aabb[0][0]) * cell_density) + 0.5f);
grid_size[1] = std::size_t(((grid_aabb[1][1] - grid_aabb[0][1]) * cell_density) + 0.5f);
grid_size[2] = std::size_t(((grid_aabb[1][2] - grid_aabb[0][2]) * cell_density) + 0.5f);
if (grid_size[0] == 0) grid_size[0] = 1;
if (grid_size[1] == 0) grid_size[1] = 1;
if (grid_size[2] == 0) grid_size[2] = 1;

// allocate and clear grid
std::size_t grid_cells_size = (grid_size[0] * grid_size[1] * grid_size[2]) + 1;
grid_cells = (index_type*) RTKernel::kernel_malloc(grid_cells_size * sizeof(index_type), 32);
std::memset((void*) grid_cells, 0, grid_cells_size * sizeof(index_type));

// count
for (std::size_t i = 0; i < number_of_triangles; ++i) {

  float aabb_triangle[2][3];
  aabb_triangle[0][0] = aabb_triangle[1][0] = triangles[i][0][0];
  aabb_triangle[0][1] = aabb_triangle[1][1] = triangles[i][0][1];
  aabb_triangle[0][2] = aabb_triangle[1][2] = triangles[i][0][2];
       if (triangles[i][1][0] < aabb_triangle[0][0]) aabb_triangle[0][0] = triangles[i][1][0];
  else if (triangles[i][1][0] > aabb_triangle[1][0]) aabb_triangle[1][0] = triangles[i][1][0];
       if (triangles[i][1][1] < aabb_triangle[0][1]) aabb_triangle[0][1] = triangles[i][1][1];
  else if (triangles[i][1][1] > aabb_triangle[1][1]) aabb_triangle[1][1] = triangles[i][1][1];
       if (triangles[i][1][2] < aabb_triangle[0][2]) aabb_triangle[0][2] = triangles[i][1][2];
  else if (triangles[i][1][2] > aabb_triangle[1][2]) aabb_triangle[1][2] = triangles[i][1][2];
       if (triangles[i][2][0] < aabb_triangle[0][0]) aabb_triangle[0][0] = triangles[i][2][0];
  else if (triangles[i][2][0] > aabb_triangle[1][0]) aabb_triangle[1][0] = triangles[i][2][0];
       if (triangles[i][2][1] < aabb_triangle[0][1]) aabb_triangle[0][1] = triangles[i][2][1];
  else if (triangles[i][2][1] > aabb_triangle[1][1]) aabb_triangle[1][1] = triangles[i][2][1];
       if (triangles[i][2][2] < aabb_triangle[0][2]) aabb_triangle[0][2] = triangles[i][2][2];
  else if (triangles[i][2][2] > aabb_triangle[1][2]) aabb_triangle[1][2] = triangles[i][2][2];

  std::size_t min[3];
  min[0] = ((aabb_triangle[0][0] - grid_aabb[0][0]) / (grid_aabb[1][0] - grid_aabb[0][0])) * grid_size[0];
  min[1] = ((aabb_triangle[0][1] - grid_aabb[0][1]) / (grid_aabb[1][1] - grid_aabb[0][1])) * grid_size[1];
  min[2] = ((aabb_triangle[0][2] - grid_aabb[0][2]) / (grid_aabb[1][2] - grid_aabb[0][2])) * grid_size[2];
  if (min[0] == grid_size[0]) min[0] = grid_size[0] - 1;
  if (min[1] == grid_size[1]) min[1] = grid_size[1] - 1;
  if (min[2] == grid_size[2]) min[2] = grid_size[2] - 1;

  std::size_t max[3];
  max[0] = ((aabb_triangle[1][0] - grid_aabb[0][0]) / (grid_aabb[1][0] - grid_aabb[0][0])) * grid_size[0];
  max[1] = ((aabb_triangle[1][1] - grid_aabb[0][1]) / (grid_aabb[1][1] - grid_aabb[0][1])) * grid_size[1];
  max[2] = ((aabb_triangle[1][2] - grid_aabb[0][2]) / (grid_aabb[1][2] - grid_aabb[0][2])) * grid_size[2];
  if (max[0] == grid_size[0]) max[0] = grid_size[0] - 1;
  if (max[1] == grid_size[1]) max[1] = grid_size[1] - 1;
  if (max[2] == grid_size[2]) max[2] = grid_size[2] - 1;

  for (std::size_t k = min[2]; k <= max[2]; ++k) {
    for (std::size_t j = min[1]; j <= max[1]; ++j) {
      for (std::size_t i = min[0]; i <= max[0]; ++i) {
        ++grid_cells[(k * (grid_size[0] * grid_size[1])) + (j * grid_size[0]) + i];
      }
    }
  }
}

// accumulate 
for (std::size_t i = 1; i <= (grid_size[0] * grid_size[1] * grid_size[2]); ++i) {
  grid_cells[i] += grid_cells[i - 1];
}
grid_lists_size = grid_cells[grid_size[0] * grid_size[1] * grid_size[2]];
grid_lists = (index_type*) RTKernel::kernel_malloc(grid_lists_size * sizeof(index_type), 32);

// insert
for (std::size_t j = 0; j < number_of_triangles; ++j) {
  std::size_t i = number_of_triangles - 1 - j;

  float aabb_triangle[2][3];
  aabb_triangle[0][0] = aabb_triangle[1][0] = triangles[i][0][0];
  aabb_triangle[0][1] = aabb_triangle[1][1] = triangles[i][0][1];
  aabb_triangle[0][2] = aabb_triangle[1][2] = triangles[i][0][2];
       if (triangles[i][1][0] < aabb_triangle[0][0]) aabb_triangle[0][0] = triangles[i][1][0];
  else if (triangles[i][1][0] > aabb_triangle[1][0]) aabb_triangle[1][0] = triangles[i][1][0];
       if (triangles[i][1][1] < aabb_triangle[0][1]) aabb_triangle[0][1] = triangles[i][1][1];
  else if (triangles[i][1][1] > aabb_triangle[1][1]) aabb_triangle[1][1] = triangles[i][1][1];
       if (triangles[i][1][2] < aabb_triangle[0][2]) aabb_triangle[0][2] = triangles[i][1][2];
  else if (triangles[i][1][2] > aabb_triangle[1][2]) aabb_triangle[1][2] = triangles[i][1][2];
       if (triangles[i][2][0] < aabb_triangle[0][0]) aabb_triangle[0][0] = triangles[i][2][0];
  else if (triangles[i][2][0] > aabb_triangle[1][0]) aabb_triangle[1][0] = triangles[i][2][0];
       if (triangles[i][2][1] < aabb_triangle[0][1]) aabb_triangle[0][1] = triangles[i][2][1];
  else if (triangles[i][2][1] > aabb_triangle[1][1]) aabb_triangle[1][1] = triangles[i][2][1];
       if (triangles[i][2][2] < aabb_triangle[0][2]) aabb_triangle[0][2] = triangles[i][2][2];
  else if (triangles[i][2][2] > aabb_triangle[1][2]) aabb_triangle[1][2] = triangles[i][2][2];

  std::size_t min[3];
  min[0] = ((aabb_triangle[0][0] - grid_aabb[0][0]) / (grid_aabb[1][0] - grid_aabb[0][0])) * grid_size[0];
  min[1] = ((aabb_triangle[0][1] - grid_aabb[0][1]) / (grid_aabb[1][1] - grid_aabb[0][1])) * grid_size[1];
  min[2] = ((aabb_triangle[0][2] - grid_aabb[0][2]) / (grid_aabb[1][2] - grid_aabb[0][2])) * grid_size[2];
  if (min[0] == grid_size[0]) min[0] = grid_size[0] - 1;
  if (min[1] == grid_size[1]) min[1] = grid_size[1] - 1;
  if (min[2] == grid_size[2]) min[2] = grid_size[2] - 1;

  std::size_t max[3];
  max[0] = ((aabb_triangle[1][0] - grid_aabb[0][0]) / (grid_aabb[1][0] - grid_aabb[0][0])) * grid_size[0];
  max[1] = ((aabb_triangle[1][1] - grid_aabb[0][1]) / (grid_aabb[1][1] - grid_aabb[0][1])) * grid_size[1];
  max[2] = ((aabb_triangle[1][2] - grid_aabb[0][2]) / (grid_aabb[1][2] - grid_aabb[0][2])) * grid_size[2];
  if (max[0] == grid_size[0]) max[0] = grid_size[0] - 1;
  if (max[1] == grid_size[1]) max[1] = grid_size[1] - 1;
  if (max[2] == grid_size[2]) max[2] = grid_size[2] - 1;

  index_type triangle_index = i;
  for (std::size_t k = min[2]; k <= max[2]; ++k) {
    for (std::size_t j = min[1]; j <= max[1]; ++j) {
      for (std::size_t i = min[0]; i <= max[0]; ++i) {
        grid_lists[--grid_cells[(k * (grid_size[0] * grid_size[1])) + (j * grid_size[0]) + i]] = triangle_index;
      }
    }
  }
}