   }
 
 
+
+  __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,
-                           uint32_t                   output_count,
-                           const uint4* __restrict__  pyramid)
-  {
+                           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;
-    float t = (10.f * ix) / output_count;
-
-    output[3 * ix + 0] = 0.5f + 0.5f * cos(t);
-    output[3 * ix + 1] = 0.5f + 0.5f * cos(3.2*t);
-    output[3 * ix + 2] = 0.5f + 0.5f * cos(2.7*t);
+    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,
-                                 uint32_t                   output_count,
-                                 const uint4* __restrict__  pyramid)
+                                 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);

                              sb[5 * threadIdx.x + 1],
                              sb[5 * threadIdx.x + 2],
                              sb[5 * threadIdx.x + 3]);
        hp1_d[offset1] = hp;
        sb1_d[offset1] = hp.x + hp.y + hp.z + hp.w + sb[5 * threadIdx.x + 4];
      }
    }

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