   // and create a full set of 32 cells, but that was slower. Also tried to add
   // extra end-fetch for the last thread to get the missing samnple for a full
   // set of 32 cells, but that was slower as well.
-  template<class FieldType>
-  __global__ __launch_bounds__(32) void reduceBase(uint8_t* __restrict__           index_cases_d,
-                                                   uint4* __restrict__             out_index_level0_d,
-                                                   uint32_t* __restrict__          out_index_sideband_d,
-                                                   const uint8_t* __restrict__     index_count,
-                                                   const FieldType* __restrict__   field_d,
-                                                   const FieldType* __restrict__   field_end_d,
-                                                   const size_t                    field_row_stride,
-                                                   const size_t                    field_slice_stride,
-                                                   const float                     threshold,
-                                                   const uint3                     chunks,
-                                                   const uint3                     cells)
-  {
-    const uint32_t warp = threadIdx.x / 32;
-    const uint32_t thread = threadIdx.x % 32;
-    const uint32_t chunk_ix = blockIdx.x + warp;
-
-    // Figure out which chunk we are in.
-    // FIXME: Try to use grid extents to mach chunk grid to avoid all these modulo/divisions.
-    uint3 chunk = make_uint3(chunk_ix % chunks.x,
-                             (chunk_ix / chunks.x) % chunks.y,
-                             (chunk_ix / chunks.x) / chunks.y);
-    uint3 cell = make_uint3(31 * chunk.x + thread,
-                            5 * chunk.y,
-                            5 * chunk.z);
-
-    // One warp processes a 32 * 5 * 5 chunk, outputs 800 mc cases to base level and 800/5=160 sums to sideband.
-    const FieldType* ptr = field_d
-      + cell.z * field_slice_stride
-      + cell.y * field_row_stride
-      + cell.x;
-
-    bool xmask = cell.x < cells.x;
-
-    // First fetch initial z-slice, and then fetch 5 additional slices, to
-    // produce 5 slices of cells.
-    float2 prev = fetch(ptr, field_end_d, field_row_stride, threshold, thread);
-    ptr += field_slice_stride;
-    for (unsigned z = 0; z < 5; z++) {
-      unsigned isum = 0;
-      bool zmask = cell.z < cells.z;
-      if (zmask) {
-        float2 next = fetch(ptr, field_end_d, field_row_stride, threshold, thread);
-        ptr += field_slice_stride;
-
-        // Merge bits from previous and current slice.
-        float2 t0 = make_float2(__fmaf_rn(float(1 << 4), next.x, prev.x),
-                                __fmaf_rn(float(1 << 4), next.y, prev.y));
-        prev = next;
-
-        // Fetch case from x+1 using warp-shuffle. Last thread has garbage data,
-        // but that result is masked out.
-        float2 tt = make_float2(__shfl_down_sync(0xffffffff, t0.x, 1),
-                                __shfl_down_sync(0xffffffff, t0.y, 1));
-
-        if (xmask && thread < 31) {
-          // Merge in results from the warp shuffle. 
-          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));
-
-          // Fetch y+1 cases and merge along Y
-          uint32_t s0 = __byte_perm(g0, g1, (4 << 12) | (2 << 8) | (1 << 4) | (0 << 0));
-          g0 = g0 | (s0 >> 6);
-          g1 = g1 | (g1 >> 6);
-
-          // Extract cases for the 5-element y-column of cases that each thread has
-          uint8_t case_y0 = __byte_perm(g0, 0, (4 << 12) | (4 << 8) | (4 << 4) | (0 << 0));
-          uint8_t case_y1 = __byte_perm(g0, 0, (4 << 12) | (4 << 8) | (4 << 4) | (1 << 0));
-          uint8_t case_y2 = __byte_perm(g0, 0, (4 << 12) | (4 << 8) | (4 << 4) | (2 << 0));
-          uint8_t case_y3 = __byte_perm(g1, 0, (4 << 12) | (4 << 8) | (4 << 4) | (0 << 0));
-          uint8_t case_y4 = __byte_perm(g1, 0, (4 << 12) | (4 << 8) | (4 << 4) | (1 << 0));
-
-          // Look up vertex count. Set count for cells outside of domain to zero.
-          uint8_t count_y0 = cell.y + 0u < cells.y ? index_count[case_y0] : 0u;
-          uint8_t count_y1 = cell.y + 1u < cells.y ? index_count[case_y1] : 0u;
-          uint8_t count_y2 = cell.y + 2u < cells.y ? index_count[case_y2] : 0u;
-          uint8_t count_y3 = cell.y + 3u < cells.y ? index_count[case_y3] : 0u;
-          uint8_t count_y4 = cell.y + 4u < cells.y ? index_count[case_y4] : 0u;
-
-          // Find sum that goes to sideband. If this is zero, there is no
-          // geometry produced in this cell, and the HP baselevel nor the
-          // array of MC cases will be touched, so no point in storing it
-          // This gives an significant speedup.
-          isum = count_y0 + count_y1 + count_y2 + count_y3 + count_y4;
-          if (isum) {
-            // MC cases and HP base level is only visited if sum is nonzero
-            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(count_y0,
-                                                                                   count_y0 + count_y1,
-                                                                                   count_y0 + count_y1 + count_y2,
-                                                                                   count_y0 + count_y1 + count_y2 + count_y3);
-          }
-        }
-      }
-      cell.z++;
-      // Store sideband data.
-      out_index_sideband_d[32 * (5 * blockIdx.x + z) + threadIdx.x] = isum;
-    }
-  }
-
-  // Indexed baselevel.
-  // ------------------                                                   
-
-  // Fetch 32[x] x 5[y] samples from the scalar field, used for the indexed
-  // baselevel buildup. A bit more care testing for valid Y so samples outside
-  // will have a GL_CLAMP-ish behaviour, which we rely on to not produce
-  // edge intersections (which produces vertices) outside of the domiain. 
-  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;
-
-    // Simd-in-float3, packing 6 uint8 bitsets into 2 float32.
-    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);
-  }
-
-  // Given a MC case (8-bits), create a mask for the (0,0,0)-(1,0,0), (0,0,0)-(0,1,0),
-  // and (0,0,0)-(0,0,1) axes where a bit is set if the bits of the edge end-points in
-  // the MC case are different (and thus edge pierces the iso-surface).
-  //
-  // This apporachs works by broadcasting the bit of (0,0,0) to all bits, and xor this
-  // with the case, and mask out the relevant bits.
-  __device__ __forceinline__ uint32_t piercingAxesFromCase(uint32_t c)
-  {
-    uint32_t a;
-    asm("{\n"
-        "  .reg .u32 t1, t2, t3;\n"
-        "  bfe.s32 t1, %1, 0, 1;\n"     // Sign-extend bit 0 (0,0,0) over whole word
-        "  xor.b32 t2, %1, t1;\n"       // XOR with (0,0,0) to see which corners change sign wrt (0,0,0).
-        "  and.b32 %0, t2, 0b10110;\n"  // Mask out (1,0,0), (0,1,0) and (0,0,1).
-        "}"
-        : "=r"(a) : "r"(c));
-    return a;
-  }
-
-  // Count the number of edge-intersections of edges belonging to that cell (0-3).
-  // We just count the number of set bits in a bitmask.
-  __device__ __forceinline__ uint32_t axisCountFromCase(const uint32_t c)
-  {
-    uint32_t a = piercingAxesFromCase(c);
-    uint32_t n;
-    asm("popc.b32 %0, %1;": "=r"(n) : "r"(a));         // Count number of 1's (should be 0..3).
-    return n;
-  }
-
-  // Fetch samples from scalar-field, deduce MC cases for cells, and write out
-  // first level ot the HistoPyramid and sideband.
-  //
-  // See comments for reduceBase
   //
   // In addition to calculating the MC cell case, it also calculates the number
   // of unique vertices belong to that cell (0-3). This grid is one cell-layer

                                                                           ctx->index_sidebands[0],                    // block writes 5*32 uint32's
                                                                           ctx->tables->index_count,
                                                                           field_d,
                                                                           field_d + size_t(field_size.x) * field_size.y * field_size.z,
                                                                           size_t(field_size.x),
                                                                           size_t(field_size.x) * field_size.y,
                                                                           threshold,

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

        if ((indexed ? cell.x <= grid_max_index.x : cell.x < grid_max_index.x) && thread < 31) {
          // Merge in results from the warp shuffle. 
          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));
