 
       unsigned sum = 0;
       bool zmask = cell.z < cells.z;
-      if (xmask && zmask) {
-
-        uint6 next = fetch(ptr, field_end_d, field_row_stride, threshold);
+      if (zmask) {
+        uint6 next = fetch(ptr, field_end_d, field_row_stride, threshold, thread);
         ptr += field_slice_stride;
-
-        uint5 cases = mergeY(mergeZ(prev, next));
-        prev = next;
-
-        uint5 counts{
-          cell.y + 0u < cells.y ? index_count[cases.e0] : 0u,
-          cell.y + 1u < cells.y ? index_count[cases.e1] : 0u,
-          cell.y + 2u < cells.y ? index_count[cases.e2] : 0u,
-          cell.y + 3u < cells.y ? index_count[cases.e3] : 0u,
-          cell.y + 4u < cells.y ? index_count[cases.e4] : 0u,
-        };
-
-        sum = counts.e0 + counts.e1 + counts.e2 + counts.e3 + counts.e4;
-
-        if (sum) {
-          // 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_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 (xmask && thread < 31) {
+
+
+          uint5 cases = mergeY(mergeZ(prev, next));
+          prev = next;
+
+          uint5 counts{
+            cell.y + 0u < cells.y ? index_count[cases.e0] : 0u,
+            cell.y + 1u < cells.y ? index_count[cases.e1] : 0u,
+            cell.y + 2u < cells.y ? index_count[cases.e2] : 0u,
+            cell.y + 3u < cells.y ? index_count[cases.e3] : 0u,
+            cell.y + 4u < cells.y ? index_count[cases.e4] : 0u,
+          };
+
+          sum = counts.e0 + counts.e1 + counts.e2 + counts.e3 + counts.e4;
+
+          if (sum) {
+            // 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_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);
+          }
         }
       }
       cell.z++;

  };

  template<typename T>
  __device__ uint6 fetch(const T* ptr, const T* end, const size_t stride, const float t)
  {
    unsigned t0 = ((ptr < end) && (ptr[0] < t) ? 1 : 0) | ((ptr <= end) && (ptr[1] < t) ? 2 : 0); ptr += stride;
    unsigned t1 = ((ptr < end) && (ptr[0] < t) ? 1 : 0) | ((ptr <= end) && (ptr[1] < t) ? 2 : 0); ptr += stride;
    unsigned t2 = ((ptr < end) && (ptr[0] < t) ? 1 : 0) | ((ptr <= end) && (ptr[1] < t) ? 2 : 0); ptr += stride;
    unsigned t3 = ((ptr < end) && (ptr[0] < t) ? 1 : 0) | ((ptr <= end) && (ptr[1] < t) ? 2 : 0); ptr += stride;
    unsigned t4 = ((ptr < end) && (ptr[0] < t) ? 1 : 0) | ((ptr <= end) && (ptr[1] < t) ? 2 : 0); ptr += stride;
    unsigned t5 = ((ptr < end) && (ptr[0] < t) ? 1 : 0) | ((ptr <= end) && (ptr[1] < t) ? 2 : 0); ptr += stride;
    return { t0, t1, t2, t3, t4, t5 };
  }


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

