                                                           const size_t                    field_slice_stride,
                                                           const float                     threshold,
                                                           const uint3                     chunks,
-                                                          const uint3                     cell_max_index)
+                                                          const uint3                     grid_max_index)
   {
     const uint32_t warp = threadIdx.x / 32;
     const uint32_t thread = threadIdx.x % 32;

    float r1 = __fmaf_rn(float(1 << 16), t5, __fmaf_rn(float(1 << 8), t4, t3));
    return make_float2(r0, r1);
  }

  __device__ __forceinline__ uint32_t piercingAxesFromCase(uint32_t c)
  {

    const FieldType* ptr = field_d
                         + cell.z * field_slice_stride
                         + cell.y * field_row_stride
                         + min(grid_max_index.x, cell.x);
    float2 prev = fetch2(ptr, field_row_stride, cell.y, grid_max_index.y, threshold);
    ptr += field_slice_stride;

    for (unsigned z = 0; z < 5; z++) {

      unsigned isum = 0;
      unsigned vsum = 0;
      bool zmask = cell.z <= grid_max_index.z;
      if (zmask) {
        float2 next = fetch2(ptr, field_row_stride, cell.y, grid_max_index.y, threshold);
        ptr += cell.z + 1 < grid_max_index.z ? field_slice_stride : 0;

        float2 t0 = make_float2(__fmaf_rn(float(1 << 4), next.x, prev.x),
                                __fmaf_rn(float(1 << 4), next.y, prev.y));

        prev = next;

        float2 tt = make_float2(__shfl_down_sync(0xffffffff, t0.x, 1),
                                __shfl_down_sync(0xffffffff, t0.y, 1));

        if (cell.x <= grid_max_index.x && thread < 31) {
          uint32_t g0 = static_cast<uint32_t>(__fmaf_rn(2.f, tt.x, t0.x));
          uint32_t g1 = static_cast<uint32_t>(__fmaf_rn(2.f, tt.y, t0.y));
          uint32_t s0 = __byte_perm(g0, g1, (4 << 12) | (2 << 8) | (1 << 4) | (0 << 0));
          g0 = g0 | (s0 >> 6);
          g1 = g1 | (g1 >> 6);

          uint32_t case_y0 = __byte_perm(g0, 0, (4 << 12) | (4 << 8) | (4 << 4) | (0 << 0));
          uint32_t case_y1 = __byte_perm(g0, 0, (4 << 12) | (4 << 8) | (4 << 4) | (1 << 0));
          uint32_t case_y2 = __byte_perm(g0, 0, (4 << 12) | (4 << 8) | (4 << 4) | (2 << 0));
          uint32_t case_y3 = __byte_perm(g1, 0, (4 << 12) | (4 << 8) | (4 << 4) | (0 << 0));
          uint32_t case_y4 = __byte_perm(g1, 0, (4 << 12) | (4 << 8) | (4 << 4) | (1 << 0));

          uint32_t ic_y0 = (cell.z < grid_max_index.z) && (cell.y + 0u) <= grid_max_index.y ? index_count[case_y0] : 0u;
          uint32_t ic_y1 = (cell.z < grid_max_index.z) && (cell.y + 1u) <= grid_max_index.y ? index_count[case_y1] : 0u;
          uint32_t ic_y2 = (cell.z < grid_max_index.z) && (cell.y + 2u) <= grid_max_index.y ? index_count[case_y2] : 0u;
          uint32_t ic_y3 = (cell.z < grid_max_index.z) && (cell.y + 3u) <= grid_max_index.y ? index_count[case_y3] : 0u;
          uint32_t ic_y4 = (cell.z < grid_max_index.z) && (cell.y + 4u) <= grid_max_index.y ? index_count[case_y4] : 0u;

          uint32_t vc_y0 = axisCountFromCase(case_y0);
          uint32_t vc_y1 = axisCountFromCase(case_y1);
          uint32_t vc_y2 = axisCountFromCase(case_y2);
          uint32_t vc_y3 = axisCountFromCase(case_y3);
          uint32_t vc_y4 = axisCountFromCase(case_y4);

          vsum = vc_y0 + vc_y1 + vc_y2 + vc_y3 + vc_y4;
          isum = ic_y0 + ic_y1 + ic_y2 + ic_y3 + ic_y4;

          if (isum | vsum) {
            // 8-bit arithmetic should suffice here
            out_index_level0_d[32 * 5 * blockIdx.x + 32 * z + thread] = make_uint4(vc_y0,
                                                                                   vc_y0 + vc_y1,
                                                                                   vc_y0 + vc_y1 + vc_y2,
                                                                                   vc_y0 + vc_y1 + vc_y2 + vc_y3);

            index_cases_d[5 * (32 * 5 * blockIdx.x + 32 * z + thread) + 0] = case_y0;
            index_cases_d[5 * (32 * 5 * blockIdx.x + 32 * z + thread) + 1] = case_y1;
            index_cases_d[5 * (32 * 5 * blockIdx.x + 32 * z + thread) + 2] = case_y2;
            index_cases_d[5 * (32 * 5 * blockIdx.x + 32 * z + thread) + 3] = case_y3;
            index_cases_d[5 * (32 * 5 * blockIdx.x + 32 * z + thread) + 4] = case_y4;
            out_index_level0_d[32 * 5 * blockIdx.x + 32 * z + thread] = make_uint4(ic_y0,
                                                                                   ic_y0 + ic_y1,
                                                                                   ic_y0 + ic_y1 + ic_y2,
                                                                                   ic_y0 + ic_y1 + ic_y2 + ic_y3);
          }
          if (vsum) {
            out_vertex_level0_d[32 * 5 * blockIdx.x + 32 * z + thread] = make_uint4(vc_y0,
                                                                                    vc_y0 + vc_y1,
                                                                                    vc_y0 + vc_y1 + vc_y2,
                                                                                    vc_y0 + vc_y1 + vc_y2 + vc_y3);

          }
          
        }
      }
      cell.z++;
      out_vertex_sideband_d[32 * (5 * blockIdx.x + z) + threadIdx.x] = vsum;
      out_index_sideband_d[32 * (5 * blockIdx.x + z) + threadIdx.x] = isum;
    }
  }

