                                                       ic_y0 + ic_y1 + ic_y2,
                                                       ic_y0 + ic_y1 + ic_y2 + ic_y3);
           }
-          if (vsum) {
+          if (indexed && vsum) {
             out_vertex_level0_ptr[32 * z] = make_uint4(vc_y0,
                                                        vc_y0 + vc_y1,
                                                        vc_y0 + vc_y1 + vc_y2,

                                                                           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,


namespace {

  // 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
