                                                           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;

    const FieldType* ptr = field_d
                         + cell.z * field_slice_stride
                         + cell.y * field_row_stride
                         + min(cell_max_index.x, cell.x);
    float2 prev = fetch(ptr, field_end_d, field_row_stride, threshold, thread);
    ptr += field_slice_stride;

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

      unsigned isum = 0;
      unsigned vsum = 0;
      bool zmask = cell.z <= cell_max_index.z;
      if (zmask) {
        float2 next = fetch(ptr, field_end_d, field_row_stride, threshold, thread);
        ptr += field_slice_stride;

      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 <= cell_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);

          uint5 cases = {
            __byte_perm(g0, 0, (4 << 12) | (4 << 8) | (4 << 4) | (0 << 0)),
            __byte_perm(g0, 0, (4 << 12) | (4 << 8) | (4 << 4) | (1 << 0)),
            __byte_perm(g0, 0, (4 << 12) | (4 << 8) | (4 << 4) | (2 << 0)),
            __byte_perm(g1, 0, (4 << 12) | (4 << 8) | (4 << 4) | (0 << 0)),
            __byte_perm(g1, 0, (4 << 12) | (4 << 8) | (4 << 4) | (1 << 0))
          };

          uint32_t vc0 = axisCountFromCase(cases.e0);
          uint32_t vc1 = axisCountFromCase(cases.e1);
          uint32_t vc2 = axisCountFromCase(cases.e2);
          uint32_t vc3 = axisCountFromCase(cases.e3);
          uint32_t vc4 = axisCountFromCase(cases.e4);

          vsum = vc0 + vc1 + vc2 + vc3 + vc4;
          if (vsum) {
            out_index_level0_d[32 * 5 * blockIdx.x + 32 * y + thread] = make_uint4(vc0,
                                                                                   vc0 + vc1,
                                                                                   vc0 + vc1 + vc2,
                                                                                   vc0 + vc1 + vc2 + vc3);
          }

          uint5 counts{
            cell.y + 0u <= cell_max_index.y ? index_count[cases.e0] : 0u,
            cell.y + 1u <= cell_max_index.y ? index_count[cases.e1] : 0u,
            cell.y + 2u <= cell_max_index.y ? index_count[cases.e2] : 0u,
            cell.y + 3u <= cell_max_index.y ? index_count[cases.e3] : 0u,
            cell.y + 4u <= cell_max_index.y ? index_count[cases.e4] : 0u,
          };

          isum = counts.e0 + counts.e1 + counts.e2 + counts.e3 + counts.e4;
          if (isum) {
            // MC cases and HP base level is only visited if sum is nonzero
            index_cases_d[5 * (32 * 5 * blockIdx.x + 32 * y + thread) + 0] = cases.e0;
            index_cases_d[5 * (32 * 5 * blockIdx.x + 32 * y + thread) + 1] = cases.e1;
            index_cases_d[5 * (32 * 5 * blockIdx.x + 32 * y + thread) + 2] = cases.e2;
            index_cases_d[5 * (32 * 5 * blockIdx.x + 32 * y + thread) + 3] = cases.e3;
            index_cases_d[5 * (32 * 5 * blockIdx.x + 32 * y + thread) + 4] = cases.e4;
            out_index_level0_d[32 * 5 * blockIdx.x + 32 * y + thread] = make_uint4(counts.e0,
                                                                             counts.e0 + counts.e1,
                                                                             counts.e0 + counts.e1 + counts.e2,
                                                                             counts.e0 + counts.e1 + counts.e2 + counts.e3);
          }
          if (vsum) {
            out_vertex_level0_d[32 * 5 * blockIdx.x + 32 * y + thread] = make_uint4(vc0,
                                                                                    vc0 + vc1,
                                                                                    vc0 + vc1 + vc2,
                                                                                    vc0 + vc1 + vc2 + vc3);

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


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

  template<typename T>
  __device__ float2 fetch2(const T* ptr, const size_t stride, const uint32_t y, const uint32_t max_y, const float t)
  {
    float t0 = (ptr[0] < t) ? 1.f : 0.f; ptr += (y + 1 <= max_y) ? stride : 0;
    float t1 = (ptr[0] < t) ? 1.f : 0.f; ptr += (y + 2 <= max_y) ? stride : 0;
    float t2 = (ptr[0] < t) ? 1.f : 0.f; ptr += (y + 3 <= max_y) ? stride : 0;
    float t3 = (ptr[0] < t) ? 1.f : 0.f; ptr += (y + 4 <= max_y) ? stride : 0;
    float t4 = (ptr[0] < t) ? 1.f : 0.f; ptr += (y + 5 <= max_y) ? stride : 0;
    float t5 = (ptr[0] < t) ? 1.f : 0.f;

    float r0 = __fmaf_rn(float(1 << 16), t2, __fmaf_rn(float(1 << 8), t1, t0));
    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)
  {
