       }
       printf("CUDA %d\n", index_count);
       if (index_count) {
-        extractP<<<(index_count + 255) / 256, 256>>>(output, index_count, pyramid);
+        extractP<<<(index_count + 255) / 256, 256>>>(output, index_count, pyramid, chunks, scale);
       }
     }
   }

                            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

  }



  __device__ __forceinline__ uint32_t processHP5Item(uint32_t& key, uint32_t offset, const uint4 item)
  {
    if (key < item.x) {
    }
    else if (key < item.y) {
      key -= item.x;
      offset += 1;
    }
    else if (key < item.z) {
      key -= item.y;
      offset += 2;
    }
    else if (key < item.w) {
      key -= item.z;
      offset += 3;
    }
    else {
      key -= item.w;
      offset += 4;
    }
    return offset;
  }

  __global__ void extractP(float* __restrict__        output,
                           const uint32_t             output_count,
                           const uint4* __restrict__  pyramid,
                           const uint2                chunks,
                           const float3               scale)
  {
    const uint32_t* level_offset = (const uint32_t*)(pyramid);
    const uint32_t level_count = level_offset[31];

    uint32_t ix = blockDim.x * blockIdx.x + threadIdx.x;
    if (ix < output_count) {
      uint32_t key = ix;
      // Traverse apex
      uint32_t offset = 0;
      offset = processHP5Item(key, 0, pyramid[9]);
      offset = processHP5Item(key, 5*offset, pyramid[10 + offset]);
      offset = processHP5Item(key, 5*offset, pyramid[15 + offset]);
      for (unsigned l = level_count - 4; l < level_count; l--) {
        offset = processHP5Item(key, 5 * offset, pyramid[level_offset[l] + offset]);
      }

      uint32_t q = offset / 800;

      uint3 chunk = make_uint3(q % chunks.x,
                               (q / chunks.x) % chunks.y,
                               (q / chunks.x) / chunks.y);

      uint32_t r = offset % 800;
      uint3 cell = make_uint3(32 * chunk.x + ((r/5) % 32),
                              5 * chunk.y + ((r % 5)),
                              5 * chunk.z + ((r / 5) / 32));


      output[3 * ix + 0] = scale.x * cell.x;
      output[3 * ix + 1] = scale.y * cell.y;
      output[3 * ix + 2] = scale.z * cell.z;
    }
  }


  __global__ void runExtractionP(float* __restrict__        output,
                                 const uint32_t             output_count,
                                 const uint4* __restrict__  pyramid,
                                 const uint2                chunks,
                                 const float3               scale)
  {
    if (threadIdx.x == 0) {
      uint32_t index_count = min(output_count, pyramid[8].x);
