diff --git a/res/org/lwjgl/demo/opengl/raytracing/cubetrace/grass_high.png b/res/org/lwjgl/demo/opengl/raytracing/cubetrace/grass_high.png
new file mode 100644
index 00000000..0af4567e
Binary files /dev/null and b/res/org/lwjgl/demo/opengl/raytracing/cubetrace/grass_high.png differ
diff --git a/res/org/lwjgl/demo/opengl/raytracing/cubetrace/kdtreeseparate32grass.glsl b/res/org/lwjgl/demo/opengl/raytracing/cubetrace/kdtreeseparate32grass.glsl
new file mode 100644
index 00000000..688aaacf
--- /dev/null
+++ b/res/org/lwjgl/demo/opengl/raytracing/cubetrace/kdtreeseparate32grass.glsl
@@ -0,0 +1,178 @@
+/*
+ * Copyright LWJGL. All rights reserved.
+ * License terms: http://lwjgl.org/license.php
+ */
+#version 430 core
+#if GL_NV_gpu_shader5
+#extension GL_NV_gpu_shader5 : enable
+#elif GL_AMD_gpu_shader_int16
+#extension GL_AMD_gpu_shader_int16 : enable
+#endif
+
+#define NO_AXIS uint(-1u)
+#define NO_NODE uint(-1u)
+#define MAX_DESCEND 200u
+#define MAX_ROPES 100u
+#define SKY_COLOR vec3(0.96, 0.98, 1.0)
+
+#define SIDE_Y_NEG 2u
+#define SIDE_Z_NEG 4u
+
+layout(binding = 0, rgba8) writeonly uniform image2D framebufferImage;
+uniform vec3 cam[5];
+uniform ivec2 off = ivec2(0,0), cbwidth = ivec2(1,1);
+uniform int scale = 1;
+uniform uint startNodeIdx = 0u;
+layout(location = 0) uniform sampler2D tex;
+
+struct node {
+  uvec4 rightOrLeaf_SplitAxis_SplitPos_GlobalOffset;
+  uint ropes[6];
+};
+layout(std430, binding = 0) readonly restrict buffer Nodes {
+  node[] nodes;
+};
+
+struct nodegeom {
+  u8vec3 min, max;
+};
+layout(std430, binding = 1) readonly restrict buffer NodeGeoms {
+  nodegeom[] nodegeoms;
+};
+
+struct leafnode {
+  uvec2 voxels;
+};
+layout(std430, binding = 2) readonly restrict buffer LeafNodes {
+  leafnode[] leafnodes;
+};
+
+layout(std430, binding = 3) readonly restrict buffer Voxels {
+  u8vec4[] voxels;
+};
+
+const bool intersectGrass(const in vec3 origin, const in vec3 dir, const in vec3 bmin, const in vec3 bmax,
+                          inout float t, out vec2 uv) {
+  vec3 p0l0 = (bmin + bmax) * 0.5 - origin;
+  float tf = (p0l0.x + p0l0.z) / (dir.x + dir.z), tf2 = (p0l0.z - p0l0.x) / (dir.z - dir.x);
+  vec3 p1 = origin + tf * dir, p2 = origin + tf2 * dir;
+  tf = mix(t, tf, all(lessThanEqual(p1, bmax)) && all(greaterThanEqual(p1, bmin)) && tf >= 0.0 && tf < t);
+  tf2 = mix(t, tf2, all(lessThanEqual(p2, bmax)) && all(greaterThanEqual(p2, bmin)) && tf2 >= 0.0 && tf2 < t);
+  vec2 uv1 = vec2(p1.z - bmin.z, 1.0 - p1.y + bmin.y), uv2 = vec2(p2.z - bmin.z, 1.0 - p2.y + bmin.y);
+  bool hit1 = tf < t && texture(tex, uv1).a > 0.8, hit2 = tf2 < t && texture(tex, uv2).a > 0.8;
+  if (hit1 || hit2) {
+    uv = mix(uv2, uv1, hit1);
+    t = min(tf, tf2);
+    return true;
+  }
+  return false;
+}
+
+const bool intersectVoxels(const in vec3 origin, const in vec3 dir, const uint firstVoxel,
+                           const uint numVoxels, inout float t, out uint vindex, out vec2 uv) {
+  bool hit = false;
+  for (uint i = 0; i < numVoxels; i++) {
+    const u8vec4 v = voxels[i + firstVoxel];
+    const vec3 vx = vec3(v.xyz)*scale;
+    if (intersectGrass(origin, dir, vx, vx + vec3(scale), t, uv)) {
+      hit = true;
+      vindex = i + firstVoxel;
+    }
+  }
+  return hit;
+}
+
+const float exitSide(const in vec3 origin, const in vec3 invdir,
+                     const in uvec3 boxMin, const in uvec3 boxMax,
+                     out uint exitSide) {
+  const bvec3 lt = lessThan(invdir, vec3(0.0));
+  const vec3 tmax = (mix(vec3(boxMax)+vec3(1.0), vec3(boxMin), lt)*scale - origin) * invdir;
+  const ivec3 signs = ivec3(lt);
+  vec2 vals;
+  vals = mix(vec2(tmax.y, SIDE_Y_NEG + signs.y), vec2(tmax.x, signs.x), tmax.y > tmax.x);
+  vals = mix(vec2(tmax.z, SIDE_Z_NEG + signs.z), vals, tmax.z > vals.x);
+  exitSide = uint(vals.y);
+  return vals.x;
+}
+
+const vec2 intersectBox(const in vec3 origin, const in vec3 invdir, const in u8vec3 bmin, const in u8vec3 bmax) {
+  const bvec3 lt = lessThan(invdir, vec3(0.0));
+  const vec3 m1 = (vec3(bmin)*scale - origin) * invdir, m2 = ((vec3(bmax) + vec3(1.0))*scale - origin) * invdir;
+  const vec3 tmin = mix(m1, m2, lt), tmax = mix(m2, m1, lt);
+  return vec2(max(max(tmin.x, tmin.y), tmin.z), min(min(tmax.x, tmax.y), tmax.z));
+}
+
+struct scenehitinfo {
+  float t;
+  uint vindex;
+  uint descends, ropes;
+  vec2 uv;
+};
+const bool intersectScene(in uint nodeIdx, const in vec3 origin, const in vec3 dir,
+                          const in vec3 invdir, out scenehitinfo shinfo) {
+  node n = nodes[nodeIdx];
+  vec3 o = origin;
+  nodegeom ng = nodegeoms[nodeIdx];
+  vec2 lambda = intersectBox(o, invdir, ng.min, ng.max);
+  if (lambda.y < 0.0 || lambda.x > lambda.y)
+    return false;
+  lambda.x = max(lambda.x, 0.0);
+  shinfo.descends = 0u;
+  shinfo.ropes = 0u;
+  shinfo.t = 1.0/0.0;
+  shinfo.uv = vec2(0.0);
+  uvec3 nmin = ng.min, nmax = ng.max;
+  while (true) {
+    vec3 pEntry = dir * lambda.x + o;
+    while (n.rightOrLeaf_SplitAxis_SplitPos_GlobalOffset.y != NO_AXIS) {
+      if (float(n.rightOrLeaf_SplitAxis_SplitPos_GlobalOffset.z)*scale <= pEntry[n.rightOrLeaf_SplitAxis_SplitPos_GlobalOffset.y]) {
+        nodeIdx = n.rightOrLeaf_SplitAxis_SplitPos_GlobalOffset.x;
+        nmin[n.rightOrLeaf_SplitAxis_SplitPos_GlobalOffset.y] = n.rightOrLeaf_SplitAxis_SplitPos_GlobalOffset.z;
+      } else {
+        nodeIdx = nodeIdx + 1u;
+        nmax[n.rightOrLeaf_SplitAxis_SplitPos_GlobalOffset.y] = n.rightOrLeaf_SplitAxis_SplitPos_GlobalOffset.z - 1u;
+      }
+      n = nodes[nodeIdx];
+      if (shinfo.descends++ > MAX_DESCEND)
+        return false;
+    }
+    if (n.rightOrLeaf_SplitAxis_SplitPos_GlobalOffset.x != NO_NODE) {
+      leafnode ln = leafnodes[n.rightOrLeaf_SplitAxis_SplitPos_GlobalOffset.x];
+      if (intersectVoxels(o, dir, ln.voxels.x, ln.voxels.y, shinfo.t, shinfo.vindex, shinfo.uv))
+        return true;
+    }
+    uint exit;
+    lambda.x = exitSide(o, invdir, nmin, nmax, exit);
+    nodeIdx = n.ropes[exit];
+    if (nodeIdx == NO_NODE)
+      return false;
+    n = nodes[nodeIdx];
+    ng = nodegeoms[nodeIdx];
+    nmin = ng.min;
+    nmax = ng.max;
+    if (shinfo.ropes++ > MAX_ROPES)
+      return false;
+    o = origin;
+  }
+  return false;
+}
+
+layout (local_size_x = 8, local_size_y = 8) in;
+void main(void) {
+  const ivec2 pix = ivec2(gl_GlobalInvocationID.xy) * cbwidth + off;
+  const ivec2 size = imageSize(framebufferImage);
+  if (any(greaterThanEqual(pix, size)))
+    return;
+  const vec2 p = (vec2(pix) + vec2(0.5)) / vec2(size);
+  const vec3 dir = normalize(mix(mix(cam[1], cam[2], p.y), mix(cam[3], cam[4], p.y), p.x));
+  vec3 acc = vec3(0.0);
+  vec3 att = vec3(1.0);
+  scenehitinfo shinfo;
+  if (intersectScene(startNodeIdx, cam[0], dir, 1.0/dir, shinfo)) {
+    vec3 col = texture(tex, shinfo.uv).rgb;
+    acc += col;// * 4.6 * vec3(float(shinfo.ropes) * 1.4E-2, float(shinfo.descends) * 4E-3, 0.0);
+  } else {
+    acc += SKY_COLOR;
+  }
+  imageStore(framebufferImage, pix, vec4(acc, 1.0));
+}
diff --git a/res/org/lwjgl/demo/opengl/raytracing/cubetrace/quad.fs.glsl b/res/org/lwjgl/demo/opengl/raytracing/cubetrace/quad.fs.glsl
deleted file mode 100644
index 3232c3d1..00000000
--- a/res/org/lwjgl/demo/opengl/raytracing/cubetrace/quad.fs.glsl
+++ /dev/null
@@ -1,13 +0,0 @@
-/*
- * Copyright LWJGL. All rights reserved.
- * License terms: https://www.lwjgl.org/license
- */
-#version 330 core
-
-uniform sampler2D tex;
-in vec2 texcoord;
-out vec4 color;
-
-void main(void) {
-  color = texture(tex, texcoord);
-}
diff --git a/res/org/lwjgl/demo/opengl/raytracing/cubetrace/quad.vs.glsl b/res/org/lwjgl/demo/opengl/raytracing/cubetrace/quad.vs.glsl
deleted file mode 100644
index 9a8a0d23..00000000
--- a/res/org/lwjgl/demo/opengl/raytracing/cubetrace/quad.vs.glsl
+++ /dev/null
@@ -1,13 +0,0 @@
-/*
- * Copyright LWJGL. All rights reserved.
- * License terms: https://www.lwjgl.org/license
- */
-#version 330 core
-
-in vec2 vertex;
-out vec2 texcoord;
-
-void main(void) {
-  gl_Position = vec4(vertex, 0.0, 1.0);
-  texcoord = vertex * 0.5 + vec2(0.5, 0.5);
-}
diff --git a/res/org/lwjgl/demo/opengl/raytracing/cubetrace/raytracing.glsl b/res/org/lwjgl/demo/opengl/raytracing/cubetrace/raytracing.glsl
deleted file mode 100644
index 06a2ed96..00000000
--- a/res/org/lwjgl/demo/opengl/raytracing/cubetrace/raytracing.glsl
+++ /dev/null
@@ -1,146 +0,0 @@
-/*
- * Copyright LWJGL. All rights reserved.
- * License terms: https://www.lwjgl.org/license
- */
-#version 430 core
-
-#define NO_NODE -1
-#define ERROR_COLOR vec3(1.0, 0.0, 1.0)
-#define MAX_FOLLOWS 750
-
-layout(binding = 0, rgba8) writeonly uniform image2D framebufferImage;
-uniform vec3 eye, ray00, ray01, ray10, ray11;
-uniform ivec2 off, cbwidth;
-
-struct node {
-  vec3 min; uint left;
-  vec3 max; uint right;
-  uint parent, firstVoxel, numVoxels;
-};
-layout(std430, binding = 0) readonly buffer Nodes {
-  node[] nodes;
-};
-
-struct voxel {
-  uvec3 p; uint c;
-};
-layout(std430, binding = 1) readonly buffer Voxels {
-  voxel[] voxels;
-};
-
-vec3 decodeColor(uint col) {
-  return vec3(col & 0x3FFu, (col >> 10u) & 0x3FF, (col >> 20u) & 0x3FF) / 1024.0;
-}
-
-bool intersectBox(const in vec3 origin, const in vec3 invdir,
-                  const in vec3 boxMin, const in vec3 boxMax, inout float t) {
-  bvec3 lt = lessThan(invdir, vec3(0.0));
-  vec3 m1 = (boxMin - origin) * invdir, m2 = (boxMax - origin) * invdir;
-  vec3 tmin = mix(m1, m2, lt), tmax = mix(m2, m1, lt);
-  float mint = max(max(tmin.x, tmin.y), tmin.z);
-  float maxt = min(min(tmax.x, tmax.y), tmax.z);
-  bool cond = mint < t && maxt >= 0.0 && mint < maxt;
-  t = mix(t, mint, cond);
-  return cond;
-}
-
-bool intersectVoxels(const in vec3 origin, const in vec3 invdir, const uint firstVoxel,
-                     const uint numVoxels, inout float t, out voxel hitvoxel) {
-  bool hit = false;
-  for (uint i = 0; i < numVoxels; i++) {
-    const voxel v = voxels[i + firstVoxel];
-    if (intersectBox(origin, invdir, vec3(v.p), vec3(v.p) + vec3(1.0), t)) {
-      hit = true;
-      hitvoxel = v;
-    }
-  }
-  return hit;
-}
-
-vec3 normalForVoxel(const in vec3 hit, const in vec3 dir, const in uvec3 v) {
-  float dx = dir.x > 0.0 ? abs(hit.x - float(v.x)) : abs(hit.x - float(v.x+1));
-  float dy = dir.y > 0.0 ? abs(hit.y - float(v.y)) : abs(hit.y - float(v.y+1));
-  float dz = dir.z > 0.0 ? abs(hit.z - float(v.z)) : abs(hit.z - float(v.z+1));
-  if (dx < dy && dx < dz)
-    return dir.x < 0.0 ? vec3(1.0, 0.0, 0.0) : vec3(-1.0, 0.0, 0.0);
-  else if (dy < dz)
-    return dir.y < 0.0 ? vec3(0.0, 1.0, 0.0) : vec3(0.0, -1.0, 0.0);
-  else
-    return dir.z < 0.0 ? vec3(0.0, 0.0, 1.0) : vec3(0.0, 0.0, -1.0);
-}
-
-uint closestChild(const in vec3 origin, const in node n, inout uint leftRightStack) {
-  if (n.left == NO_NODE)
-    return NO_NODE;
-  const node left = nodes[n.left], right = nodes[n.right];
-  const vec3 midLeft = (left.min + left.max) * vec3(0.5);
-  const vec3 midRight = (right.min + right.max) * vec3(0.5);
-  const vec3 dl = origin - midLeft, dr = origin - midRight;
-  uint res;
-  if (dot(dl, dl) < dot(dr, dr)) {
-    leftRightStack <<= 1;
-    res = n.left;
-  } else {
-    leftRightStack = leftRightStack << 1 | 1;
-    res = n.right;
-  }
-  return res;
-}
-
-bool processNextFarChild(inout uint nearFarStack, inout uint leftRightStack,
-                         const in uint pidx, inout uint nextIdx) {
-  node parent = nodes[pidx];
-  while ((nearFarStack & 1u) == 1u) {
-    nearFarStack >>= 1;
-    leftRightStack >>= 1;
-    if (parent.parent == NO_NODE)
-      return false;
-    parent = nodes[parent.parent];
-  }
-  nextIdx = (leftRightStack & 1u) == 0u ? parent.right : parent.left;
-  nearFarStack |= 1;
-  return true;
-}
-
-vec3 trace(const in vec3 origin, const in vec3 dir, const in vec3 invdir) {
-  float nt = 1.0/0.0, bt = 1.0/0.0;
-  vec3 normal = vec3(0.0), col = vec3(1.0);
-  uint nextIdx = 0u, iterations = 0u, leftRightStack = 0u, nearFarStack = 0u;
-  while (true) {
-    if (iterations++ > MAX_FOLLOWS)
-      return ERROR_COLOR;
-    const node next = nodes[nextIdx];
-    if (!intersectBox(origin, invdir, next.min, next.max, bt)) {
-      if (!processNextFarChild(nearFarStack, leftRightStack, next.parent, nextIdx))
-        break;
-    } else {
-      if (next.numVoxels > 0) {
-        voxel hitvoxel;
-        if (intersectVoxels(origin, invdir, next.firstVoxel, next.numVoxels, nt, hitvoxel)) {
-          normal = normalForVoxel(origin + nt * dir, dir, hitvoxel.p);
-          col = decodeColor(hitvoxel.c);
-        }
-        if (!processNextFarChild(nearFarStack, leftRightStack, next.parent, nextIdx))
-          break;
-      } else {
-        nextIdx = closestChild(origin, next, leftRightStack);
-        nearFarStack <<= 1;
-      }
-      bt = nt;
-    }
-  }
-  return col * max(0.2, dot(normal, vec3(0.4, 0.82, 0.4)));
-}
-
-layout (local_size_x = 8, local_size_y = 8) in;
-
-void main(void) {
-  ivec2 pix = ivec2(gl_GlobalInvocationID.xy) * cbwidth + off;
-  ivec2 size = imageSize(framebufferImage);
-  if (any(greaterThanEqual(pix, size)))
-    return;
-  vec2 p = (vec2(pix) + vec2(0.5)) / vec2(size);
-  vec3 dir = normalize(mix(mix(ray00, ray01, p.y), mix(ray10, ray11, p.y), p.x));
-  vec3 color = trace(eye, dir, 1.0 / dir);
-  imageStore(framebufferImage, pix, vec4(color, 1.0));
-}
diff --git a/res/org/lwjgl/demo/opengl/raytracing/cubetrace/raytracing_16t.glsl b/res/org/lwjgl/demo/opengl/raytracing/cubetrace/raytracing_16t.glsl
deleted file mode 100644
index 8d4368bc..00000000
--- a/res/org/lwjgl/demo/opengl/raytracing/cubetrace/raytracing_16t.glsl
+++ /dev/null
@@ -1,152 +0,0 @@
-/*
- * Copyright LWJGL. All rights reserved.
- * License terms: https://www.lwjgl.org/license
- */
-#version 430 core
-#if GL_NV_gpu_shader5
-#extension GL_NV_gpu_shader5 : enable
-#elif GL_AMD_gpu_shader_int16
-#extension GL_AMD_gpu_shader_int16 : enable
-#endif
-
-#define NO_NODE -1u
-#define ERROR_COLOR vec3(1.0, 0.0, 1.0)
-#define MAX_ITERATIONS 750u
-
-layout(binding = 0, rgba8) writeonly uniform image2D framebufferImage;
-uniform vec3 eye, ray00, ray01, ray10, ray11;
-uniform ivec2 off, cbwidth;
-
-struct node {
-  u16vec3 min; uint8_t numVoxels;
-  u16vec3 max;
-  uint left, right, parent, firstVoxel;
-};
-layout(std430, binding = 0) readonly buffer Nodes {
-  node[] nodes;
-};
-
-struct voxel {
-  u16vec3 p; uint16_t c;
-};
-layout(std430, binding = 1) readonly buffer Voxels {
-  voxel[] voxels;
-};
-
-vec3 decodeColor(uint16_t col) {
-  return vec3(col & 0x1Fu, (col >> 5u) & 0x1Fu, (col >> 10u) & 0x1Fu) / 31.0;
-}
-
-bool intersectBox(const in vec3 origin, const in vec3 invdir,
-                  const in u16vec3 bmin, const in u16vec3 bmax, inout float t) {
-  bvec3 lt = lessThan(invdir, vec3(0.0));
-  vec3 m1 = (bmin - origin) * invdir, m2 = (bmax - origin) * invdir;
-  vec3 tmin = mix(m1, m2, lt), tmax = mix(m2, m1, lt);
-  float mint = max(max(tmin.x, tmin.y), tmin.z);
-  float maxt = min(min(tmax.x, tmax.y), tmax.z);
-  bool cond = mint < t && maxt >= 0.0 && mint < maxt;
-  t = mix(t, mint, cond);
-  return cond;
-}
-
-bool intersectVoxels(const in vec3 origin, const in vec3 invdir, const uint firstVoxel,
-                     const uint8_t numVoxels, inout float t, out voxel hitvoxel) {
-  bool hit = false;
-  for (uint i = 0; i < numVoxels; i++) {
-    const voxel v = voxels[i + firstVoxel];
-    if (intersectBox(origin, invdir, v.p, v.p+u16vec3(1), t)) {
-      hit = true;
-      hitvoxel = v;
-    }
-  }
-  return hit;
-}
-
-vec3 normalForVoxel(const in vec3 hit, const in vec3 dir, const in u16vec3 v) {
-  float dx = dir.x > 0.0 ? abs(hit.x - float(v.x)) : abs(hit.x - float(v.x+1u));
-  float dy = dir.y > 0.0 ? abs(hit.y - float(v.y)) : abs(hit.y - float(v.y+1u));
-  float dz = dir.z > 0.0 ? abs(hit.z - float(v.z)) : abs(hit.z - float(v.z+1u));
-  if (dx < dy && dx < dz)
-    return dir.x < 0.0 ? vec3(1.0, 0.0, 0.0) : vec3(-1.0, 0.0, 0.0);
-  else if (dy < dz)
-    return dir.y < 0.0 ? vec3(0.0, 1.0, 0.0) : vec3(0.0, -1.0, 0.0);
-  else
-    return dir.z < 0.0 ? vec3(0.0, 0.0, 1.0) : vec3(0.0, 0.0, -1.0);
-}
-
-uint closestChild(const in vec3 origin, const in node n, inout uint leftRightStack) {
-  if (n.left == NO_NODE)
-    return NO_NODE;
-  const node left = nodes[n.left], right = nodes[n.right];
-  const vec3 midLeft = (left.min + left.max) * vec3(0.5);
-  const vec3 midRight = (right.min + right.max) * vec3(0.5);
-  const vec3 dl = origin - midLeft, dr = origin - midRight;
-  uint res;
-  if (dot(dl, dl) < dot(dr, dr)) {
-    leftRightStack <<= 1;
-    res = n.left;
-  } else {
-    leftRightStack = leftRightStack << 1 | 1;
-    res = n.right;
-  }
-  return res;
-}
-
-bool processNextFarChild(inout uint nearFarStack, inout uint leftRightStack,
-                         const in uint pidx, inout uint nextIdx) {
-  node parent = nodes[pidx];
-  while ((nearFarStack & 1u) == 1u) {
-    nearFarStack >>= 1;
-    leftRightStack >>= 1;
-    if (parent.parent == NO_NODE)
-      return false;
-    parent = nodes[parent.parent];
-  }
-  nextIdx = (leftRightStack & 1u) == 0u ? parent.right : parent.left;
-  nearFarStack |= 1;
-  return true;
-}
-
-vec3 trace(const in vec3 origin, const in vec3 dir, const in vec3 invdir) {
-  float nt = 1.0/0.0, bt = 1.0/0.0;
-  vec3 normal = vec3(0.0), col = vec3(1.0);
-  uint nextIdx = 0u, iterations = 0u;
-  uint leftRightStack = 0u, nearFarStack = 0u;
-  while (true) {
-    if (iterations++ > MAX_ITERATIONS)
-      return ERROR_COLOR;
-    const node next = nodes[nextIdx];
-    if (!intersectBox(origin, invdir, next.min, next.max, bt)) {
-      if (!processNextFarChild(nearFarStack, leftRightStack, next.parent, nextIdx))
-        break;
-    } else {
-      if (next.numVoxels > 0u) {
-        voxel hitvoxel;
-        if (intersectVoxels(origin, invdir, next.firstVoxel, next.numVoxels, nt, hitvoxel)) {
-          normal = normalForVoxel(origin + nt * dir, dir, hitvoxel.p);
-          col = decodeColor(hitvoxel.c);
-        }
-        if (!processNextFarChild(nearFarStack, leftRightStack, next.parent, nextIdx))
-          break;
-      } else {
-        nextIdx = closestChild(origin, next, leftRightStack);
-        nearFarStack <<= 1;
-      }
-      bt = nt;
-    }
-  }
-  return col * max(0.2, dot(normal, vec3(0.4, 0.82, 0.4)));
-}
-
-layout (local_size_x = 8, local_size_y = 8) in;
-
-void main(void) {
-  ivec2 pix = ivec2(gl_GlobalInvocationID.xy) * cbwidth + off;
-  ivec2 size = imageSize(framebufferImage);
-  if (any(greaterThanEqual(pix, size)))
-    return;
-  vec2 p = (vec2(pix) + vec2(0.5)) / vec2(size);
-  vec3 dir = normalize(mix(mix(ray00, ray01, p.y), mix(ray10, ray11, p.y), p.x));
-  vec3 color = trace(eye, dir, 1.0 / dir);
-  imageStore(framebufferImage, pix, vec4(color, 1.0));
-}
diff --git a/src/org/lwjgl/demo/opengl/raytracing/CubeTrace.java b/src/org/lwjgl/demo/opengl/raytracing/CubeTrace.java
index 820eccce..6cbad4c4 100644
--- a/src/org/lwjgl/demo/opengl/raytracing/CubeTrace.java
+++ b/src/org/lwjgl/demo/opengl/raytracing/CubeTrace.java
@@ -1,196 +1,87 @@
-/*
- * Copyright LWJGL. All rights reserved.
- * License terms: https://www.lwjgl.org/license
+/*	
+ * Copyright LWJGL. All rights reserved.	
+ * License terms: https://www.lwjgl.org/license	
  */
 package org.lwjgl.demo.opengl.raytracing;
 
 import static org.lwjgl.glfw.GLFW.*;
 import static org.lwjgl.opengl.GL43C.*;
+import static org.lwjgl.opengl.NVDrawTexture.*;
+import static org.lwjgl.opengl.NVFillRectangle.*;
+import static org.lwjgl.stb.STBImage.*;
 import static org.lwjgl.system.MemoryUtil.*;
 
-import java.io.IOException;
+import java.io.*;
 import java.lang.Math;
 import java.nio.*;
 import java.util.*;
 
 import org.joml.*;
-import org.lwjgl.demo.opengl.raytracing.CubeTrace.IBVHMortonTree.*;
 import org.lwjgl.demo.opengl.util.*;
 import org.lwjgl.glfw.*;
 import org.lwjgl.opengl.*;
 import org.lwjgl.system.*;
 
-/**
- * Tracing Minecraft-like axis-aligned cube worlds using a BVH tree with morton
- * code partitioning and bitstack traversal.
- * <p>
- * Axis-aligned cubes are really a natural fit for BVH trees, since they align
- * very well with the BVH subdivisions and a terrain with cubes at regular
- * positions can be subdivided optimally in log2(N) iterations.
- * 
- * @author Kai Burjack
- */
 public class CubeTrace {
-    /**
-     * Bounding Volume Hierarchy for integer lattices using morton code
-     * partitioning.
-     *
-     * @author Kai Burjack
-     */
-    public static class IBVHMortonTree {
-        public static class Voxel implements Comparable<Voxel> {
-            public int x, y, z;
-            public long morton;
-            public float r, g, b;
-
-            public Voxel(int x, int y, int z, float r, float g, float b) {
-                this.x = x;
-                this.y = y;
-                this.z = z;
-                this.r = r;
-                this.g = g;
-                this.b = b;
-            }
-
-            @Override
-            public int compareTo(Voxel o) {
-                return Long.compare(morton, o.morton);
-            }
-        }
-
-        public static final int MAX_POINTS_IN_NODE = 8;
-        public IBVHMortonTree parent;
-        public IBVHMortonTree left;
-        public IBVHMortonTree right;
-        public List<Voxel> voxels;
-        public int index;
-        public int first, last = -1;
-        public int minX = Integer.MAX_VALUE, minY = Integer.MAX_VALUE, minZ = Integer.MAX_VALUE;
-        public int maxX = Integer.MIN_VALUE, maxY = Integer.MIN_VALUE, maxZ = Integer.MIN_VALUE;
-
-        IBVHMortonTree(IBVHMortonTree parent) {
-            this.parent = parent;
-        }
-
-        IBVHMortonTree(IBVHMortonTree parent, List<Voxel> voxels, int first, int last) {
-            this.parent = parent;
-            this.first = first;
-            this.last = last;
-            this.voxels = voxels;
-            voxels.forEach(v -> {
-                this.minX = this.minX < v.x ? this.minX : v.x;
-                this.minY = this.minY < v.y ? this.minY : v.y;
-                this.minZ = this.minZ < v.z ? this.minZ : v.z;
-                this.maxX = this.maxX > v.x ? this.maxX : v.x;
-                this.maxY = this.maxY > v.y ? this.maxY : v.y;
-                this.maxZ = this.maxZ > v.z ? this.maxZ : v.z;
-            });
-        }
-
-        private static long expandBits(long v) {
-            v &= 0x00000000001fffffL;
-            v = (v | v << 32L) & 0x001f00000000ffffL;
-            v = (v | v << 16L) & 0x001f0000ff0000ffL;
-            v = (v | v << 8L) & 0x010f00f00f00f00fL;
-            v = (v | v << 4L) & 0x10c30c30c30c30c3L;
-            v = (v | v << 2L) & 0x1249249249249249L;
-            return v;
-        }
-
-        private static long morton3d(int x, int y, int z) {
-            return (expandBits(y) << 2L) + (expandBits(z) << 1L) + expandBits(x);
-        }
-
-        private static int findSplit(List<Voxel> sortedMortonCodes, int first, int last) {
-            long firstCode = sortedMortonCodes.get(first).morton;
-            long lastCode = sortedMortonCodes.get(last).morton;
-            if (firstCode == lastCode)
-                return first + last >> 1;
-            int commonPrefix = Long.numberOfLeadingZeros(firstCode ^ lastCode);
-            int split = first;
-            int step = last - first;
-            do {
-                step = step + 1 >> 1;
-                int newSplit = split + step;
-                if (newSplit < last) {
-                    long splitCode = sortedMortonCodes.get(newSplit).morton;
-                    int splitPrefix = Long.numberOfLeadingZeros(firstCode ^ splitCode);
-                    if (splitPrefix > commonPrefix)
-                        split = newSplit;
-                }
-            } while (step > 1);
-            return split;
-        }
-
-        public static IBVHMortonTree build(List<Voxel> voxels) {
-            voxels.forEach(v -> v.morton = morton3d(v.x, v.y, v.z));
-            Collections.sort(voxels);
-            return build(null, voxels, 0, voxels.size() - 1);
-        }
-
-        private static IBVHMortonTree build(IBVHMortonTree parent, List<Voxel> mortonSortedTriangles, int first,
-                int last) {
-            if (first > last - MAX_POINTS_IN_NODE)
-                return new IBVHMortonTree(parent, mortonSortedTriangles.subList(first, last + 1), first, last);
-            int split = findSplit(mortonSortedTriangles, first, last);
-            IBVHMortonTree tree = new IBVHMortonTree(parent);
-            tree.left = build(tree, mortonSortedTriangles, first, split);
-            tree.right = build(tree, mortonSortedTriangles, split + 1, last);
-            tree.minX = java.lang.Math.min(tree.left.minX, tree.right.minX);
-            tree.minY = java.lang.Math.min(tree.left.minY, tree.right.minY);
-            tree.minZ = java.lang.Math.min(tree.left.minZ, tree.right.minZ);
-            tree.maxX = java.lang.Math.max(tree.left.maxX, tree.right.maxX);
-            tree.maxY = java.lang.Math.max(tree.left.maxY, tree.right.maxY);
-            tree.maxZ = java.lang.Math.max(tree.left.maxZ, tree.right.maxZ);
-            return tree;
-        }
-    }
-
     private long window;
     private int width = 1920;
     private int height = 1080;
     private boolean resetFramebuffer;
     private int cbWidth = 3;
-    private byte[] samplePattern = samplePattern();
     private int cbPixel;
-    private int levelWidth = 2048;
-    private int levelDepth = 1024;
-    private int levelHeight = 64;
+    private byte[] samplePattern = samplePattern();
+    private int scale = 1;
+    private int levelWidth = 256 / scale;
+    private int levelDepth = 256 / scale;
+    private int levelHeight = 128 / scale;
     private int pttex;
     private int vao;
-    private int computeProgram;
+    private int tex;
+    private int finalGatherProgram;
     private int quadProgram;
-    private int sampler;
-    private int eyeUniform;
-    private int offUniform, cbwidthUniform;
-    private int ray00Uniform, ray10Uniform, ray01Uniform, ray11Uniform;
-    private int framebufferImageBinding;
+    private int camUniform;
+    private KDTreei.Node<KDTreei.Voxel> cameraNode;
+    private int startNodeIdxUniform;
+    private int prevVPUniform;
+    private int texUniform;
+    private int offUniform, cbwidthUniform, scaleUniform;
     private int nodesSsbo;
-    private int nodesSsboBinding;
+    private int nodeGeomsSsbo;
+    private int leafNodesSsbo;
     private int voxelsSsbo;
-    private int voxelsSsboBinding;
-    private int workGroupSizeX;
-    private int workGroupSizeY;
+    private int finalGatherTimeUniform;
+    private int finalGatherWorkGroupSizeX;
+    private int finalGatherWorkGroupSizeY;
     private float mouseX, mouseY;
-    private boolean mouseDown;
+    private boolean leftMouseDown;
     private boolean[] keydown = new boolean[GLFW.GLFW_KEY_LAST + 1];
     private Matrix4d projMatrix = new Matrix4d();
     private Matrix4d viewMatrix = new Matrix4d();
+    private Matrix4d viewProjMatrix = new Matrix4d();
+    private Matrix4d prevViewProjMatrix = new Matrix4d();
     private Matrix4d invViewProjMatrix = new Matrix4d();
-    private Vector3d tmpVector = new Vector3d();
-    private Vector3d cameraPosition = new Vector3d(levelWidth * 0.9f, levelHeight * 1.4f, levelDepth * 0.5f);
-    private Vector3d cameraLookAt = new Vector3d(levelWidth * 0.5f, levelHeight * 0.5f, levelDepth * 0.5f);
-    private Vector3d cameraUp = new Vector3d(0.0f, 1.0f, 0.0f);
+    private Vector3d tmpVector = new Vector3d(), tmpVector2 = new Vector3d();
+    private Vector3d cameraPosition = new Vector3d(levelWidth * 0.3f * scale, levelHeight * 0.4f * scale,
+            levelDepth * 0.5f * scale);
+    private Vector3d cameraLookDir = new Vector3d(0, -0.1f, 1).normalize();
+    private Quaterniond quaternion = new Quaterniond();
+    private KDTreei<KDTreei.Voxel> root;
     private GLFWErrorCallback errCallback;
     private GLFWKeyCallback keyCallback;
     private GLFWFramebufferSizeCallback fbCallback;
     private GLFWCursorPosCallback cpCallback;
     private GLFWMouseButtonCallback mbCallback;
     private Callback debugProc;
-    private long lastTime = System.nanoTime();
+    private GLCapabilities caps;
+    private boolean GL_NV_draw_texture, GL_NV_fill_rectangle;
     private int frame = 0;
-    private float avgTime = 0.0f;
-    private boolean hasShortsInShader;
+    private int[] viewport = { 0, 0, width, height };
+    private List<KDTreei.Voxel> voxels;
+
+    {
+        cameraPosition.set(2.033E+2f, 4.511E+0f, 1.232E+2f);
+        cameraLookDir.set(-1.504E-1f, 4.513E-2f, 9.876E-1f);
+    }
 
     private void init() throws IOException {
         glfwSetErrorCallback(errCallback = GLFWErrorCallback.createPrint(System.err));
@@ -202,13 +93,12 @@ private void init() throws IOException {
         glfwWindowHint(GLFW_CONTEXT_VERSION_MINOR, 3);
         glfwWindowHint(GLFW_VISIBLE, GLFW_FALSE);
         glfwWindowHint(GLFW_RESIZABLE, GLFW_TRUE);
-        window = glfwCreateWindow(width, height, "Tracing axis-aligned cubes with BVH", NULL, NULL);
+//        long monitor = glfwGetMonitors().get(1);
+//        GLFWVidMode vidmode = glfwGetVideoMode(monitor);
+//        window = glfwCreateWindow(vidmode.width(), vidmode.height(), "Raytracing water", monitor, NULL);
+        window = glfwCreateWindow(width, height, "Raytracing water", NULL, NULL);
         if (window == NULL)
             throw new AssertionError("Failed to create the GLFW window");
-        System.out.println("Press WSAD, LCTRL, SPACE to move around in the scene.");
-        System.out.println("Press Q/E to roll left/right.");
-        System.out.println("Hold down left shift to move faster.");
-        System.out.println("Move the mouse to look around.");
         glfwSetKeyCallback(window, keyCallback = new GLFWKeyCallback() {
             public void invoke(long window, int key, int scancode, int action, int mods) {
                 if (key == GLFW_KEY_ESCAPE && action != GLFW_RELEASE)
@@ -216,11 +106,18 @@ public void invoke(long window, int key, int scancode, int action, int mods) {
                 if (key == -1)
                     return;
                 keydown[key] = action == GLFW_PRESS || action == GLFW_REPEAT;
+                if (action == GLFW_PRESS && key == GLFW_KEY_R) {
+                    frame = 0;
+                } else if (action == GLFW_PRESS && key == GLFW_KEY_C) {
+                    System.out.println(cameraPosition);
+                    System.out.println(cameraLookDir);
+                }
             }
         });
         glfwSetFramebufferSizeCallback(window, fbCallback = new GLFWFramebufferSizeCallback() {
             public void invoke(long window, int width, int height) {
-                if (width > 0 && height > 0 && (CubeTrace.this.width != width || CubeTrace.this.height != height)) {
+                if (width > 0 && height > 0 && (CubeTrace.this.width != width
+                        || CubeTrace.this.height != height)) {
                     CubeTrace.this.width = width;
                     CubeTrace.this.height = height;
                     CubeTrace.this.resetFramebuffer = true;
@@ -229,12 +126,12 @@ public void invoke(long window, int width, int height) {
         });
         glfwSetCursorPosCallback(window, cpCallback = new GLFWCursorPosCallback() {
             public void invoke(long window, double x, double y) {
-                if (mouseDown) {
+                if (leftMouseDown) {
                     float deltaX = (float) x - CubeTrace.this.mouseX;
                     float deltaY = (float) y - CubeTrace.this.mouseY;
-                    CubeTrace.this.viewMatrix.rotateLocalY(deltaX * 0.002f);
-                    CubeTrace.this.viewMatrix.rotateLocalX(deltaY * 0.002f);
-                    CubeTrace.this.viewMatrix.normalize3x3();
+                    cameraLookDir
+                            .rotate(quaternion.identity().rotateAxis(-deltaY * 0.002f, viewMatrix.positiveX(tmpVector))
+                                    .rotateAxis(-deltaX * 0.002f, viewMatrix.positiveY(tmpVector)));
                 }
                 CubeTrace.this.mouseX = (float) x;
                 CubeTrace.this.mouseY = (float) y;
@@ -242,14 +139,14 @@ public void invoke(long window, double x, double y) {
         });
         glfwSetMouseButtonCallback(window, mbCallback = new GLFWMouseButtonCallback() {
             public void invoke(long window, int button, int action, int mods) {
-                if (action == GLFW_PRESS)
-                    CubeTrace.this.mouseDown = true;
-                else if (action == GLFW_RELEASE)
-                    CubeTrace.this.mouseDown = false;
+                if (action == GLFW_PRESS && button == GLFW_MOUSE_BUTTON_1)
+                    leftMouseDown = true;
+                else if (action == GLFW_RELEASE && button == GLFW_MOUSE_BUTTON_1)
+                    leftMouseDown = false;
             }
         });
-        GLFWVidMode vidmode = glfwGetVideoMode(glfwGetPrimaryMonitor());
-        glfwSetWindowPos(window, (vidmode.width() - width) / 2, (vidmode.height() - height) / 2);
+//        GLFWVidMode vidmode = glfwGetVideoMode(glfwGetPrimaryMonitor());
+//        glfwSetWindowPos(window, (vidmode.width() - width) / 2, (vidmode.height() - height) / 2);
         glfwMakeContextCurrent(window);
         glfwSwapInterval(0);
         try (MemoryStack frame = MemoryStack.stackPush()) {
@@ -258,28 +155,53 @@ else if (action == GLFW_RELEASE)
             width = framebufferSize.get(0);
             height = framebufferSize.get(1);
         }
-        GLCapabilities caps = GL.createCapabilities();
-        hasShortsInShader = caps.GL_NV_gpu_shader5 || caps.GL_AMD_gpu_shader_int16;
-        // debugProc = GLUtil.setupDebugMessageCallback();
-        viewMatrix.setLookAt(cameraPosition, cameraLookAt, cameraUp);
-        System.out.println("Building terrain...");
+        caps = GL.createCapabilities();
+        GL_NV_draw_texture = caps.GL_NV_draw_texture;
+        GL_NV_fill_rectangle = caps.GL_NV_fill_rectangle;
+        debugProc = GLUtil.setupDebugMessageCallback();
+        loadTextures();
         createFramebufferTextures();
-        createSampler();
-        quadFullScreenVao();
-        createComputeProgram();
-        initComputeProgram();
-        createQuadProgram();
-        initQuadProgram();
-        createSceneSSBOs(buildTerrainVoxels());
+        createFinalGatherProgram();
+        if (!GL_NV_draw_texture) {
+            createFullScreenVao();
+            createQuadProgram();
+            initQuadProgram();
+        }
+        System.out.println("Building terrain...");
+        voxels = buildTerrainVoxels();
+        createSceneSSBOs(voxels);
         glfwShowWindow(window);
     }
 
+    private void loadTextures() throws IOException {
+        try (MemoryStack frame = MemoryStack.stackPush()) {
+            IntBuffer width = frame.mallocInt(1);
+            IntBuffer height = frame.mallocInt(1);
+            IntBuffer components = frame.mallocInt(1);
+            tex = glGenTextures();
+            glBindTexture(GL_TEXTURE_2D, tex);
+            glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_MIN_FILTER, GL_LINEAR_MIPMAP_NEAREST);
+            glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_MAG_FILTER, GL_NEAREST);
+            glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_WRAP_S, GL_CLAMP_TO_EDGE);
+            glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_WRAP_T, GL_CLAMP_TO_EDGE);
+            ByteBuffer data;
+            data = stbi_load_from_memory(DemoUtils.ioResourceToByteBuffer("org/lwjgl/demo/opengl/raytracing/cubetrace/grass_high.png", 1024), width,
+                    height, components, 4);
+            int w = width.get(0), h = height.get(0);
+            glTexStorage2D(GL_TEXTURE_2D, 6, GL_RGBA8, w, h);
+            glTexSubImage2D(GL_TEXTURE_2D, 0, 0, 0, w, h, GL_RGBA, GL_UNSIGNED_BYTE, data);
+            stbi_image_free(data);
+            glGenerateMipmap(GL_TEXTURE_2D);
+            glBindTexture(GL_TEXTURE_2D, 0);
+        }
+    }
+
     private byte[] samplePattern() {
         switch (cbWidth) {
         case 1:
             return new byte[] { 0, 0 };
         case 2:
-        /* Use custom sample patterns for 2x2 up to 4x4 pixel groups */
+            /* Use custom sample patterns for 2x2 up to 4x4 pixel groups */
             return new byte[] { 0, 0, 1, 1, 1, 0, 0, 1 };
         case 3:
             return new byte[] { 0, 0, 1, 2, 2, 1, 0, 1, 2, 0, 1, 1, 2, 2, 0, 2, 1, 0 };
@@ -291,144 +213,165 @@ private byte[] samplePattern() {
         }
     }
 
-    private List<Voxel> buildTerrainVoxels() {
+    private static int idx(int x, int y, int z, int width, int depth) {
+        return x + z * width + y * width * depth;
+    }
+
+    private List<KDTreei.Voxel> buildTerrainVoxels() {
         int width = levelWidth, height = levelHeight, depth = levelDepth;
-        float xzScale = 0.01353f, yScale = 0.0322f;
+        float xzScale = 0.01543f * scale;
         byte[] field = new byte[width * depth * height];
-        Thread[] threads = new Thread[Runtime.getRuntime().availableProcessors()];
-        for (int i = 0; i < threads.length; i++) {
-            final int ti = i;
-            threads[i] = new Thread(() -> {
-                int depthSlice = depth / threads.length;
-                for (int z = ti * depthSlice; z < (ti + 1) * depthSlice; z++)
-                    for (int y = 0; y < height; y++)
-                        for (int x = 0; x < width; x++) {
-                            float v = SimplexNoise.noise(x * xzScale, z * xzScale, y * yScale);
-                            if (y == 0 || v > 0.15f)
-                                field[x + y * width + z * width * height] = ~0;
-                        }
-            });
-            threads[i].start();
-        }
-        for (int i = 0; i < threads.length; i++) {
-            try {
-                threads[i].join();
-            } catch (InterruptedException e) {
+        for (int z = 0; z < depth; z++) {
+            for (int x = 0; x < width; x++) {
+                float dx = (x - width * 0.5f) / width * 2, dz = (z - depth * 0.5f) / depth * 2;
+                float dist = (float) Math.min(1.0, Math.sqrt(dx * dx + dz * dz));
+                float y = height * (SimplexNoise.noise(x * xzScale, z * xzScale) * 0.5f + 0.5f)
+                        * ((float) Math.pow(1.0 - dist * 0.8, 5));
+                y = Math.min(y, height);
+                for (int y0 = 0; y0 <= y; y0++) {
+                    field[idx(x, y0, z, width, depth)] = ~0;
+                }
             }
         }
         /* Remove voxels that have neighbors at all sides */
-        int removed = 0;
-        for (int z = 1; z < depth - 1; z++)
-            for (int y = 1; y < height - 1; y++)
-                for (int x = 1; x < width - 1; x++) {
-                    int idx = x + y * width + z * width * height;
-                    int left = (x - 1) + y * width + z * width * height;
-                    int right = (x + 1) + y * width + z * width * height;
-                    int up = x + (y + 1) * width + z * width * height;
-                    int down = x + (y - 1) * width + z * width * height;
-                    int front = x + y * width + (z - 1) * width * height;
-                    int back = x + y * width + (z + 1) * width * height;
-                    if ((field[idx] & 1) == 1 && (field[left] & 1) == 1 && (field[right] & 1) == 1
-                            && (field[up] & 1) == 1 && (field[down] & 1) == 1 && (field[front] & 1) == 1
-                            && (field[back] & 1) == 1) {
-                        field[idx] = 3; // <- remember that this once was a filled voxel!
-                        removed++;
+        for (int z = 0; z < depth; z++) {
+            for (int y = 0; y < height; y++) {
+                for (int x = 0; x < width; x++) {
+                    int idx = idx(x, y, z, width, depth);
+                    int left = idx(x - 1, y, z, width, depth);
+                    int right = idx(x + 1, y, z, width, depth);
+                    int up = idx(x, y + 1, z, width, depth);
+                    int down = idx(x, y - 1, z, width, depth);
+                    int front = idx(x, y, z - 1, width, depth);
+                    int back = idx(x, y, z + 1, width, depth);
+                    if ((field[idx] & 1) == 1 && (x == 0 || (field[left] & 1) == 1)
+                            && (x == width - 1 || (field[right] & 1) == 1) && (y == height - 1 || (field[up] & 1) == 1)
+                            && (y == 0 || (field[down] & 1) == 1) && (z == 0 || (field[front] & 1) == 1)
+                            && (z == depth - 1 || (field[back] & 1) == 1)) {
+                        field[idx] = 1; // <- remember that this once was a filled voxel!
                     }
                 }
-        List<Voxel> voxels = new ArrayList<>();
-        for (int z = 0; z < depth; z++)
-            for (int y = 0; y < height; y++)
-                for (int x = 0; x < width; x++) {
-                    int idx = x + y * width + z * width * height;
-                    float r = 0.9f * (SimplexNoise.noise(x * 0.012f, y * 0.04f, z * 0.012f) * 0.5f + 0.5f);
-                    float g = r;
-                    float b = g;
+            }
+        }
+        List<KDTreei.Voxel> voxels = new ArrayList<>();
+        for (short z = 0; z < depth; z++) {
+            for (short y = 0; y < height; y++) {
+                for (short x = 0; x < width; x++) {
+                    int idx = idx(x, y, z, width, depth);
+                    int left = x > 0 && (field[idx(x - 1, y, z, width, depth)] & 1) == 1 ? 1 : 0;
+                    int right = x < width - 1 && (field[idx(x + 1, y, z, width, depth)] & 1) == 1 ? 1 : 0;
+                    int down = y > 0 && (field[idx(x, y - 1, z, width, depth)] & 1) == 1 ? 1 : 0;
+                    int up = y < height - 1 && (field[idx(x, y + 1, z, width, depth)] & 1) == 1 ? 1 : 0;
+                    int back = z > 0 && (field[idx(x, y, z - 1, width, depth)] & 1) == 1 ? 1 : 0;
+                    int front = z < depth - 1 && (field[idx(x, y, z + 1, width, depth)] & 1) == 1 ? 1 : 0;
+                    int sidesFlag = left | (right << 1) | (down << 2) | (up << 3) | (back << 4) | (front << 5);
                     if (field[idx] == ~0)
-                        voxels.add(new Voxel(x, y, z, r, g, b));
+                        voxels.add(new KDTreei.Voxel((byte) x, (byte) y, (byte) z, (byte) ((sidesFlag << 1) | 1)));
                 }
-        System.out.println("Removed voxels: " + removed);
+            }
+        }
         System.out.println("Retained voxels: " + voxels.size());
         return voxels;
     }
 
-    private static List<IBVHMortonTree> allocate(IBVHMortonTree node) {
-        List<IBVHMortonTree> linearNodes = new ArrayList<>();
-        Queue<IBVHMortonTree> nodes = new ArrayDeque<>();
-        int index = 0;
+    private static List<KDTreei.Node<KDTreei.Voxel>> allocate(KDTreei.Node<KDTreei.Voxel> node) {
+        List<KDTreei.Node<KDTreei.Voxel>> linearNodes = new ArrayList<>();
+        LinkedList<KDTreei.Node<KDTreei.Voxel>> nodes = new LinkedList<>();
+        int index = 0, leafIndex = 0;
         nodes.add(node);
         while (!nodes.isEmpty()) {
-            IBVHMortonTree n = nodes.poll();
+            KDTreei.Node<KDTreei.Voxel> n = nodes.removeFirst();
             linearNodes.add(n);
             n.index = index++;
             if (n.left != null) {
-                nodes.add(n.left);
-                nodes.add(n.right);
+                nodes.addFirst(n.right);
+                nodes.addFirst(n.left);
+            } else {
+                if (!n.voxels.isEmpty())
+                    n.leafIndex = leafIndex++;
+                else
+                    n.leafIndex = -1;
+                n.voxels.forEach(v -> {
+                    v.nindex = n.index;
+                });
             }
         }
         return linearNodes;
     }
 
-    private void createSceneSSBOs(List<Voxel> voxels) {
-        System.out.println("Building BVH...");
-        IBVHMortonTree root = IBVHMortonTree.build(voxels);
-        System.out.println("Writing BVH to buffers...");
-        DynamicByteBuffer voxelsBuffer = new DynamicByteBuffer(voxels.size() * 4 * 4);
-        if (hasShortsInShader) {
-            voxels.forEach(v -> {
-                int r = (int) (v.r * 32.0f);
-                int g = (int) (v.g * 32.0f);
-                int b = (int) (v.b * 32.0f);
-                voxelsBuffer.putShort(v.x).putShort(v.y).putShort(v.z)
-                        .putShort((r & 0x1F) | ((g & 0x1F) << 5) | ((b & 0x1F) << 10));
-            });
-        } else {
-            voxels.forEach(v -> {
-                int r = (int) (v.r * 1024.0f);
-                int g = (int) (v.g * 1024.0f);
-                int b = (int) (v.b * 1024.0f);
-                voxelsBuffer.putInt(v.x).putInt(v.y).putInt(v.z)
-                        .putInt((r & 0x3FF) | ((g & 0x3FF) << 10) | ((b & 0x3FF) << 20));
-            });
+    private void createSceneSSBOs(List<KDTreei.Voxel> voxels) {
+        System.out.println("Building kd tree...");
+        for (int i = 0; i < 1; i++) {
+            long time1 = System.nanoTime();
+            root = KDTreei.build(voxels, 17);
+            long time2 = System.nanoTime();
+            System.out.println("KD tree build took: " + (time2 - time1) * 1E-6f + " ms.");
         }
-        voxelsBuffer.flip();
-        System.out.println("Voxels SSBO size: " + voxelsBuffer.remaining() / 1024 / 1024 + " MB");
-        DynamicByteBuffer nodesBuffer = new DynamicByteBuffer(8192);
-        bhvToBuffers(root, nodesBuffer);
-        nodesBuffer.flip();
-        System.out.println("BVH SSBO size: " + nodesBuffer.remaining() / 1024 / 1024 + " MB");
-        this.nodesSsbo = glGenBuffers();
+        System.out.println("Writing kd and voxels to buffers...");
+        DynamicByteBuffer voxelsBuffer = new DynamicByteBuffer();
+        DynamicByteBuffer nodesBuffer = new DynamicByteBuffer();
+        DynamicByteBuffer nodeGeomsBuffer = new DynamicByteBuffer();
+        DynamicByteBuffer leafNodesBuffer = new DynamicByteBuffer();
+        long time1 = System.nanoTime();
+        kdTreeToBuffers(root, 0, 0, nodesBuffer, nodeGeomsBuffer, leafNodesBuffer, voxelsBuffer);
+        long time2 = System.nanoTime();
+        System.out.println("kd tree and voxels to buffers took: " + (time2 - time1) * 1E-6f + " ms.");
+        System.out.println("Voxels SSBO size: " + String.format("%.2f", voxelsBuffer.pos / 1024.0 / 1024.0) + " MB");
+        System.out.println("Nodes SSBO size: " + String.format("%.2f", nodesBuffer.pos / 1024.0 / 1024.0) + " MB");
+        System.out.println(
+                "Nodes Geometry SSBO size: " + String.format("%.2f", nodeGeomsBuffer.pos / 1024.0 / 1024.0) + " MB");
+        System.out.println(
+                "Leaf Nodes SSBO size: " + String.format("%.2f", leafNodesBuffer.pos / 1024.0 / 1024.0) + " MB");
+        nodesSsbo = glGenBuffers();
         glBindBuffer(GL_ARRAY_BUFFER, nodesSsbo);
-        glBufferData(GL_ARRAY_BUFFER, nodesBuffer.bb, GL_STATIC_DRAW);
+        nglBufferData(GL_ARRAY_BUFFER, nodesBuffer.pos, nodesBuffer.addr, GL_STATIC_DRAW);
         nodesBuffer.free();
-        glBindBuffer(GL_ARRAY_BUFFER, 0);
-        this.voxelsSsbo = glGenBuffers();
+        nodeGeomsSsbo = glGenBuffers();
+        glBindBuffer(GL_ARRAY_BUFFER, nodeGeomsSsbo);
+        nglBufferData(GL_ARRAY_BUFFER, nodeGeomsBuffer.pos, nodeGeomsBuffer.addr, GL_STATIC_DRAW);
+        nodeGeomsBuffer.free();
+        leafNodesSsbo = glGenBuffers();
+        glBindBuffer(GL_ARRAY_BUFFER, leafNodesSsbo);
+        nglBufferData(GL_ARRAY_BUFFER, leafNodesBuffer.pos, leafNodesBuffer.addr, GL_STATIC_DRAW);
+        leafNodesBuffer.free();
+        voxelsSsbo = glGenBuffers();
         glBindBuffer(GL_ARRAY_BUFFER, voxelsSsbo);
-        glBufferData(GL_ARRAY_BUFFER, voxelsBuffer.bb, GL_STATIC_DRAW);
+        nglBufferData(GL_ARRAY_BUFFER, voxelsBuffer.pos, voxelsBuffer.addr, GL_STATIC_DRAW);
         voxelsBuffer.free();
         glBindBuffer(GL_ARRAY_BUFFER, 0);
     }
 
-    private void bhvToBuffers(IBVHMortonTree root, DynamicByteBuffer nodesBuffer) {
-        if (hasShortsInShader) {
-            for (IBVHMortonTree n : allocate(root)) {
-                nodesBuffer.putShort(n.minX).putShort(n.minY).putShort(n.minZ).putByte(n.last - n.first + 1).putByte(0);
-                nodesBuffer.putShort(n.maxX + 1).putShort(n.maxY + 1).putShort(n.maxZ + 1).putShort(0);
-                nodesBuffer.putInt(n.left != null ? n.left.index : -1).putInt(n.right != null ? n.right.index : -1);
-                nodesBuffer.putInt(n.parent != null ? n.parent.index : -1).putInt(n.first);
-            }
-        } else {
-            for (IBVHMortonTree n : allocate(root)) {
-                nodesBuffer.putFloat(n.minX).putFloat(n.minY).putFloat(n.minZ)
-                        .putInt(n.left != null ? n.left.index : -1);
-                nodesBuffer.putFloat(n.maxX + 1).putFloat(n.maxY + 1).putFloat(n.maxZ + 1)
-                        .putInt(n.right != null ? n.right.index : -1);
-                nodesBuffer.putInt(n.parent != null ? n.parent.index : -1).putInt(n.first).putInt(n.last - n.first + 1)
-                        .putInt(-1);
+    private void kdTreeToBuffers(KDTreei<KDTreei.Voxel> root, int nodeIndexOffset, int voxelIndexOffset,
+            DynamicByteBuffer nodesBuffer, DynamicByteBuffer nodeGeomsBuffer, DynamicByteBuffer leafNodesBuffer,
+            DynamicByteBuffer voxelsBuffer) {
+        int first = 0;
+        List<KDTreei.Node<KDTreei.Voxel>> nodes = allocate(root.root);
+        System.out.println("Num nodes: " + nodes.size());
+        for (KDTreei.Node<KDTreei.Voxel> n : nodes) {
+            int numVoxels = 0;
+            if (n.left == null) {
+                numVoxels = n.voxels.size();
+                n.voxels.forEach(v -> {
+                    voxelsBuffer.putByte(v.x).putByte(v.y).putByte(v.z).putByte(v.paletteIndex);
+                });
+                if (n.leafIndex != -1)
+                    leafNodesBuffer.putInt(first).putInt(numVoxels);
             }
+            nodeGeomsBuffer.putByte(n.boundingBox.minX).putByte(n.boundingBox.minY).putByte(n.boundingBox.minZ)
+                    .putByte(0);
+            nodeGeomsBuffer.putByte(n.boundingBox.maxX - 1).putByte(n.boundingBox.maxY - 1)
+                    .putByte(n.boundingBox.maxZ - 1).putByte(0);
+            nodesBuffer.putInt(n.right != null ? n.right.index + nodeIndexOffset : n.leafIndex);
+            nodesBuffer.putInt(n.splitAxis).putInt(n.splitPos);
+            nodesBuffer.putByte(0).putByte(0).putByte(0).putByte(0);
+            for (int i = 0; i < 6; i++)
+                nodesBuffer.putInt(n.ropes[i] != null ? n.ropes[i].index : -1);
+            nodesBuffer.putInt(0).putInt(0);
+            first += numVoxels;
         }
     }
 
-    private void quadFullScreenVao() {
+    private void createFullScreenVao() {
         this.vao = glGenVertexArrays();
         int vbo = glGenBuffers();
         glBindVertexArray(vao);
@@ -438,9 +381,11 @@ private void quadFullScreenVao() {
         fv.put(-1.0f).put(-1.0f);
         fv.put(1.0f).put(-1.0f);
         fv.put(1.0f).put(1.0f);
-        fv.put(1.0f).put(1.0f);
-        fv.put(-1.0f).put(1.0f);
-        fv.put(-1.0f).put(-1.0f);
+        if (!GL_NV_fill_rectangle) {
+            fv.put(1.0f).put(1.0f);
+            fv.put(-1.0f).put(1.0f);
+            fv.put(-1.0f).put(-1.0f);
+        }
         glBufferData(GL_ARRAY_BUFFER, bb, GL_STATIC_DRAW);
         MemoryUtil.memFree(bb);
         glEnableVertexAttribArray(0);
@@ -451,10 +396,8 @@ private void quadFullScreenVao() {
 
     private void createQuadProgram() throws IOException {
         int program = glCreateProgram();
-        int vshader = DemoUtils.createShader("org/lwjgl/demo/opengl/raytracing/cubetrace/quad.vs.glsl",
-                GL_VERTEX_SHADER);
-        int fshader = DemoUtils.createShader("org/lwjgl/demo/opengl/raytracing/cubetrace/quad.fs.glsl",
-                GL_FRAGMENT_SHADER);
+        int vshader = DemoUtils.createShader("cubetrace/experiments/quad.vs.glsl", GL_VERTEX_SHADER);
+        int fshader = DemoUtils.createShader("cubetrace/experiments/quad.fs.glsl", GL_FRAGMENT_SHADER);
         glAttachShader(program, vshader);
         glAttachShader(program, fshader);
         glBindAttribLocation(program, 0, "vertex");
@@ -462,76 +405,58 @@ private void createQuadProgram() throws IOException {
         glLinkProgram(program);
         int linked = glGetProgrami(program, GL_LINK_STATUS);
         String programLog = glGetProgramInfoLog(program);
-        if (programLog.trim().length() > 0)
+        if (programLog != null && programLog.trim().length() > 0)
             System.err.println(programLog);
         if (linked == 0)
             throw new AssertionError("Could not link program");
         this.quadProgram = program;
     }
 
-    private void createComputeProgram() throws IOException {
+    private void initQuadProgram() {
+        glUseProgram(quadProgram);
+        int texUniform = glGetUniformLocation(quadProgram, "tex");
+        glUniform1i(texUniform, 0);
+        glUseProgram(0);
+    }
+
+    private void createFinalGatherProgram() throws IOException {
         int program = glCreateProgram();
-        int cshader = DemoUtils.createShader(
-                "org/lwjgl/demo/opengl/raytracing/cubetrace/raytracing" + (hasShortsInShader ? "_16t" : "") + ".glsl",
-                GL_COMPUTE_SHADER);
+        int cshader = DemoUtils.createShader("org/lwjgl/demo/opengl/raytracing/cubetrace/kdtreeseparate32grass.glsl", GL_COMPUTE_SHADER);
         glAttachShader(program, cshader);
         glLinkProgram(program);
         int linked = glGetProgrami(program, GL_LINK_STATUS);
         String programLog = glGetProgramInfoLog(program);
-        if (programLog.trim().length() > 0)
+        if (programLog != null && programLog.trim().length() > 0)
             System.err.println(programLog);
         if (linked == 0)
             throw new AssertionError("Could not link program");
-        this.computeProgram = program;
-    }
-
-    private void initQuadProgram() {
-        glUseProgram(quadProgram);
-        int texUniform = glGetUniformLocation(quadProgram, "tex");
-        glUniform1i(texUniform, 0);
-        glUseProgram(0);
-    }
-
-    private void initComputeProgram() {
-        glUseProgram(computeProgram);
+        this.finalGatherProgram = program;
+        glUseProgram(finalGatherProgram);
         try (MemoryStack frame = MemoryStack.stackPush()) {
             IntBuffer workGroupSize = frame.mallocInt(3);
-            glGetProgramiv(computeProgram, GL_COMPUTE_WORK_GROUP_SIZE, workGroupSize);
-            workGroupSizeX = workGroupSize.get(0);
-            workGroupSizeY = workGroupSize.get(1);
-            cbwidthUniform = glGetUniformLocation(computeProgram, "cbwidth");
-            offUniform = glGetUniformLocation(computeProgram, "off");
-            eyeUniform = glGetUniformLocation(computeProgram, "eye");
-            ray00Uniform = glGetUniformLocation(computeProgram, "ray00");
-            ray10Uniform = glGetUniformLocation(computeProgram, "ray10");
-            ray01Uniform = glGetUniformLocation(computeProgram, "ray01");
-            ray11Uniform = glGetUniformLocation(computeProgram, "ray11");
-            IntBuffer props = frame.mallocInt(1);
-            IntBuffer params = frame.mallocInt(1);
-            props.put(0, GL_BUFFER_BINDING);
-            int nodesResourceIndex = glGetProgramResourceIndex(computeProgram, GL_SHADER_STORAGE_BLOCK, "Nodes");
-            glGetProgramResourceiv(computeProgram, GL_SHADER_STORAGE_BLOCK, nodesResourceIndex, props, null, params);
-            nodesSsboBinding = params.get(0);
-            int voxelsResourceIndex = glGetProgramResourceIndex(computeProgram, GL_SHADER_STORAGE_BLOCK, "Voxels");
-            glGetProgramResourceiv(computeProgram, GL_SHADER_STORAGE_BLOCK, voxelsResourceIndex, props, null, params);
-            voxelsSsboBinding = params.get(0);
-            int loc = glGetUniformLocation(computeProgram, "framebufferImage");
-            glGetUniformiv(computeProgram, loc, params);
-            framebufferImageBinding = params.get(0);
+            glGetProgramiv(finalGatherProgram, GL_COMPUTE_WORK_GROUP_SIZE, workGroupSize);
+            finalGatherWorkGroupSizeX = workGroupSize.get(0);
+            finalGatherWorkGroupSizeY = workGroupSize.get(1);
         }
+        cbwidthUniform = glGetUniformLocation(finalGatherProgram, "cbwidth");
+        scaleUniform = glGetUniformLocation(finalGatherProgram, "scale");
+        texUniform = glGetUniformLocation(finalGatherProgram, "tex");
+        glUniform1i(texUniform, 0);
+        offUniform = glGetUniformLocation(finalGatherProgram, "off");
+        camUniform = glGetUniformLocation(finalGatherProgram, "cam");
+        startNodeIdxUniform = glGetUniformLocation(finalGatherProgram, "startNodeIdx");
+        prevVPUniform = glGetUniformLocation(finalGatherProgram, "prevVP");
+        finalGatherTimeUniform = glGetUniformLocation(finalGatherProgram, "time");
+        glUniform1i(scaleUniform, scale);
         glUseProgram(0);
     }
 
     private void createFramebufferTextures() {
         pttex = glGenTextures();
         glBindTexture(GL_TEXTURE_2D, pttex);
-        glTexStorage2D(GL_TEXTURE_2D, 1, GL_RGBA32F, width, height);
-    }
-
-    private void createSampler() {
-        this.sampler = glGenSamplers();
-        glSamplerParameteri(this.sampler, GL_TEXTURE_MIN_FILTER, GL_NEAREST);
-        glSamplerParameteri(this.sampler, GL_TEXTURE_MAG_FILTER, GL_NEAREST);
+        glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_MIN_FILTER, GL_NEAREST);
+        glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_MAG_FILTER, GL_NEAREST);
+        glTexStorage2D(GL_TEXTURE_2D, 1, GL_RGBA8, width, height);
     }
 
     private void resizeFramebufferTexture() {
@@ -540,104 +465,115 @@ private void resizeFramebufferTexture() {
     }
 
     private void update(float dt) {
-        float factor = 20.0f;
+        if (resetFramebuffer) {
+            resizeFramebufferTexture();
+            resetFramebuffer = false;
+        }
+        viewport[2] = width;
+        viewport[3] = height;
+        float factor = 10.0f;
         if (keydown[GLFW_KEY_LEFT_SHIFT])
-            factor = 50.0f;
+            factor = 20.0f;
         if (keydown[GLFW_KEY_W])
-            viewMatrix.translateLocal(0, 0, factor * dt);
+            cameraPosition.sub(viewMatrix.positiveZ(tmpVector).mul(dt * factor));
         if (keydown[GLFW_KEY_S])
-            viewMatrix.translateLocal(0, 0, -factor * dt);
+            cameraPosition.add(viewMatrix.positiveZ(tmpVector).mul(dt * factor));
         if (keydown[GLFW_KEY_A])
-            viewMatrix.translateLocal(factor * dt, 0, 0);
+            cameraPosition.sub(viewMatrix.positiveX(tmpVector).mul(dt * factor));
         if (keydown[GLFW_KEY_D])
-            viewMatrix.translateLocal(-factor * dt, 0, 0);
-        if (keydown[GLFW_KEY_Q])
-            viewMatrix.rotateLocalZ(factor * 0.02f * -dt);
-        if (keydown[GLFW_KEY_E])
-            viewMatrix.rotateLocalZ(factor * 0.02f * dt);
+            cameraPosition.add(viewMatrix.positiveX(tmpVector).mul(dt * factor));
         if (keydown[GLFW_KEY_LEFT_CONTROL])
-            viewMatrix.translateLocal(0, factor * dt, 0);
+            cameraPosition.sub(tmpVector.set(0, 1, 0).mul(dt * factor));
         if (keydown[GLFW_KEY_SPACE])
-            viewMatrix.translateLocal(0, -factor * dt, 0);
-        projMatrix.setPerspective((float) Math.toRadians(60.0f), (float) width / height, 0.01f, 100.0f);
+            cameraPosition.add(tmpVector.set(0, 1, 0).mul(dt * factor));
+        projMatrix.setPerspective((float) Math.toRadians(60.0f), (float) width / height, 1f, 1000.0f);
+        viewMatrix.setLookAt(cameraPosition, tmpVector.set(cameraPosition).add(cameraLookDir), tmpVector2.set(0, 1, 0));
+        viewProjMatrix.set(projMatrix).mul(viewMatrix).invert(invViewProjMatrix);
+        /* Determine camera node */
+        cameraNode = root.findNode(cameraPosition);
     }
 
-    private void trace() {
-        glUseProgram(computeProgram);
-        if (resetFramebuffer) {
-            resizeFramebufferTexture();
-            resetFramebuffer = false;
-        }
-        projMatrix.invertPerspectiveView(viewMatrix, invViewProjMatrix);
+    private void trace(float elapsedSeconds) {
+        glUseProgram(finalGatherProgram);
         viewMatrix.originAffine(cameraPosition);
+        glUniform1ui(startNodeIdxUniform, cameraNode.index);
+        glUniform1f(finalGatherTimeUniform, elapsedSeconds);
         glUniform2i(cbwidthUniform, cbWidth, cbWidth);
         glUniform2i(offUniform, samplePattern[cbPixel << 1], samplePattern[(cbPixel << 1) + 1]);
         cbPixel = (cbPixel + 1) % (cbWidth * cbWidth);
-        glUniform3f(eyeUniform, (float) cameraPosition.x, (float) cameraPosition.y, (float) cameraPosition.z);
-        invViewProjMatrix.transformProject(tmpVector.set(-1, -1, 0)).sub(cameraPosition);
-        glUniform3f(ray00Uniform, (float) tmpVector.x, (float) tmpVector.y, (float) tmpVector.z);
-        invViewProjMatrix.transformProject(tmpVector.set(-1, 1, 0)).sub(cameraPosition);
-        glUniform3f(ray01Uniform, (float) tmpVector.x, (float) tmpVector.y, (float) tmpVector.z);
-        invViewProjMatrix.transformProject(tmpVector.set(1, -1, 0)).sub(cameraPosition);
-        glUniform3f(ray10Uniform, (float) tmpVector.x, (float) tmpVector.y, (float) tmpVector.z);
-        invViewProjMatrix.transformProject(tmpVector.set(1, 1, 0)).sub(cameraPosition);
-        glUniform3f(ray11Uniform, (float) tmpVector.x, (float) tmpVector.y, (float) tmpVector.z);
-        glBindImageTexture(framebufferImageBinding, pttex, 0, false, 0, GL_WRITE_ONLY, GL_RGBA32F);
-        glBindBufferBase(GL_SHADER_STORAGE_BUFFER, nodesSsboBinding, nodesSsbo);
-        glBindBufferBase(GL_SHADER_STORAGE_BUFFER, voxelsSsboBinding, voxelsSsbo);
-        int numGroupsX = (int) Math.ceil((double)width / workGroupSizeX);
-        int numGroupsY = (int) Math.ceil((double)height / workGroupSizeY);
-        int q0 = ARBOcclusionQuery.glGenQueriesARB();
-        int q1 = ARBOcclusionQuery.glGenQueriesARB();
-        ARBTimerQuery.glQueryCounter(q0, ARBTimerQuery.GL_TIMESTAMP);
-        glDispatchCompute(numGroupsX, numGroupsY, 1);
-        ARBTimerQuery.glQueryCounter(q1, ARBTimerQuery.GL_TIMESTAMP);
-        while (ARBOcclusionQuery.glGetQueryObjectiARB(q0, ARBOcclusionQuery.GL_QUERY_RESULT_AVAILABLE_ARB) == 0
-                || ARBOcclusionQuery.glGetQueryObjectiARB(q1, ARBOcclusionQuery.GL_QUERY_RESULT_AVAILABLE_ARB) == 0)
-            Thread.yield();
-        long time1 = ARBTimerQuery.glGetQueryObjectui64(q0, ARBOcclusionQuery.GL_QUERY_RESULT_ARB);
-        long time2 = ARBTimerQuery.glGetQueryObjectui64(q1, ARBOcclusionQuery.GL_QUERY_RESULT_ARB);
-        float factor = (float) frame / (frame + 1);
-        avgTime = (1.0f - factor) * avgTime + factor * (time2 - time1);
-        frame++;
-        if (System.nanoTime() - lastTime >= 1E9) {
-            System.err.println(avgTime * 1E-6 + " ms.");
-            lastTime = System.nanoTime();
-            frame = 0;
+        glUniform3f(camUniform, (float) cameraPosition.x, (float) cameraPosition.y, (float) cameraPosition.z);
+        invViewProjMatrix.transformProject(tmpVector.set(-1, -1, -1)).sub(cameraPosition);
+        glUniform3f(camUniform + 1, (float) tmpVector.x, (float) tmpVector.y, (float) tmpVector.z);
+        invViewProjMatrix.transformProject(tmpVector.set(-1, 1, -1)).sub(cameraPosition);
+        glUniform3f(camUniform + 2, (float) tmpVector.x, (float) tmpVector.y, (float) tmpVector.z);
+        invViewProjMatrix.transformProject(tmpVector.set(1, -1, -1)).sub(cameraPosition);
+        glUniform3f(camUniform + 3, (float) tmpVector.x, (float) tmpVector.y, (float) tmpVector.z);
+        invViewProjMatrix.transformProject(tmpVector.set(1, 1, -1)).sub(cameraPosition);
+        glUniform3f(camUniform + 4, (float) tmpVector.x, (float) tmpVector.y, (float) tmpVector.z);
+        try (MemoryStack frame = MemoryStack.stackPush()) {
+            glUniformMatrix4fv(prevVPUniform, false, prevViewProjMatrix.get(frame.mallocFloat(16)));
         }
-        ARBOcclusionQuery.glDeleteQueriesARB(q0);
-        ARBOcclusionQuery.glDeleteQueriesARB(q1);
+        glBindImageTexture(0, pttex, 0, false, 0, GL_WRITE_ONLY, GL_RGBA8);
+        glBindBufferBase(GL_SHADER_STORAGE_BUFFER, 0, nodesSsbo);
+        glBindBufferBase(GL_SHADER_STORAGE_BUFFER, 1, nodeGeomsSsbo);
+        glBindBufferBase(GL_SHADER_STORAGE_BUFFER, 2, leafNodesSsbo);
+        glBindBufferBase(GL_SHADER_STORAGE_BUFFER, 3, voxelsSsbo);
+        glBindTexture(GL_TEXTURE_2D, tex);
+        int worksizeX = (int) Math.ceil((double) width / cbWidth / finalGatherWorkGroupSizeX);
+        int worksizeY = (int) Math.ceil((double) height / cbWidth / finalGatherWorkGroupSizeY);
+        glDispatchCompute(worksizeX, worksizeY, 1);
         glMemoryBarrier(GL_SHADER_IMAGE_ACCESS_BARRIER_BIT);
-        glBindBufferBase(GL_SHADER_STORAGE_BUFFER, nodesSsboBinding, 0);
-        glBindBufferBase(GL_SHADER_STORAGE_BUFFER, voxelsSsboBinding, 0);
-        glBindImageTexture(framebufferImageBinding, 0, 0, false, 0, GL_WRITE_ONLY, GL_RGBA32F);
+        glBindTexture(GL_TEXTURE_2D, 0);
+        glBindBufferBase(GL_SHADER_STORAGE_BUFFER, 0, 0);
+        glBindBufferBase(GL_SHADER_STORAGE_BUFFER, 1, 0);
+        glBindBufferBase(GL_SHADER_STORAGE_BUFFER, 2, 0);
+        glBindBufferBase(GL_SHADER_STORAGE_BUFFER, 3, 0);
+        glBindImageTexture(0, 0, 0, false, 0, GL_WRITE_ONLY, GL_RGBA8);
         glUseProgram(0);
     }
 
     private void present() {
-        glUseProgram(quadProgram);
-        glBindVertexArray(vao);
-        glBindTexture(GL_TEXTURE_2D, pttex);
-        glBindSampler(0, this.sampler);
-        glDrawArrays(GL_TRIANGLES, 0, 6);
-        glBindSampler(0, 0);
-        glBindTexture(GL_TEXTURE_2D, 0);
-        glBindVertexArray(0);
-        glUseProgram(0);
+        if (GL_NV_draw_texture) {
+            glDrawTextureNV(pttex, 0, 0, 0, width, height, 0, 0, 0, 1, 1);
+        } else {
+            glUseProgram(quadProgram);
+            glBindVertexArray(vao);
+            glBindTexture(GL_TEXTURE_2D, pttex);
+            if (GL_NV_fill_rectangle) {
+                glPolygonMode(GL_FRONT_AND_BACK, GL_FILL_RECTANGLE_NV);
+                glDrawArrays(GL_TRIANGLES, 0, 3);
+                glPolygonMode(GL_FRONT_AND_BACK, GL_FILL);
+            } else {
+                glDrawArrays(GL_TRIANGLES, 0, 6);
+            }
+            glBindTexture(GL_TEXTURE_2D, 0);
+            glBindVertexArray(0);
+            glUseProgram(0);
+        }
     }
 
+    private int averageOver = 400;
+
     private void loop() {
-        float lastTime = System.nanoTime();
+        long lastTime = System.nanoTime();
+        long lastAvg = System.nanoTime();
         while (!glfwWindowShouldClose(window)) {
-            float thisTime = System.nanoTime();
-            float dt = (thisTime - lastTime) / 1E9f;
-            lastTime = thisTime;
             glfwPollEvents();
+            long thisTime = System.nanoTime();
+            float dt = (thisTime - lastTime) * 1E-9f;
+            if ((frame % averageOver) == 0) {
+                String text = String.format("%.3f", (thisTime - lastAvg) * 1E-6f / averageOver) + " ms.";
+                glfwSetWindowTitle(window, text);
+                System.out.println(text);
+                lastAvg = thisTime;
+            }
+            lastTime = thisTime;
+            update((float) dt);
             glViewport(0, 0, width, height);
-            update(dt);
-            trace();
+            trace(thisTime * 1E-9f);
             present();
             glfwSwapBuffers(window);
+            frame++;
         }
     }
 
diff --git a/src/org/lwjgl/demo/opengl/raytracing/DemoSsboTrianglesStacklessKdTree.java b/src/org/lwjgl/demo/opengl/raytracing/DemoSsboTrianglesStacklessKdTree.java
index 23ef4f61..ba4b74e6 100644
--- a/src/org/lwjgl/demo/opengl/raytracing/DemoSsboTrianglesStacklessKdTree.java
+++ b/src/org/lwjgl/demo/opengl/raytracing/DemoSsboTrianglesStacklessKdTree.java
@@ -337,16 +337,14 @@ void createSceneSSBO() {
         DynamicByteBuffer nodesBuffer = new DynamicByteBuffer();
         DynamicByteBuffer trianglesBuffer = new DynamicByteBuffer();
         kdTreeToBuffers(kdtree, nodesBuffer, trianglesBuffer);
-        nodesBuffer.flip();
-        trianglesBuffer.flip();
 
         this.nodesSsbo = glGenBuffers();
         glBindBuffer(GL_ARRAY_BUFFER, nodesSsbo);
-        glBufferData(GL_ARRAY_BUFFER, nodesBuffer.bb, GL_STATIC_DRAW);
+        nglBufferData(GL_ARRAY_BUFFER, nodesBuffer.pos, nodesBuffer.addr, GL_STATIC_DRAW);
         glBindBuffer(GL_ARRAY_BUFFER, 0);
         this.trianglesSsbo = glGenBuffers();
         glBindBuffer(GL_ARRAY_BUFFER, trianglesSsbo);
-        glBufferData(GL_ARRAY_BUFFER, trianglesBuffer.bb, GL_STATIC_DRAW);
+        nglBufferData(GL_ARRAY_BUFFER, trianglesBuffer.pos, trianglesBuffer.addr, GL_STATIC_DRAW);
         glBindBuffer(GL_ARRAY_BUFFER, 0);
     }
 
diff --git a/src/org/lwjgl/demo/opengl/raytracing/HybridDemoSsboTriangles.java b/src/org/lwjgl/demo/opengl/raytracing/HybridDemoSsboTriangles.java
index 9173b230..17bc9e5c 100644
--- a/src/org/lwjgl/demo/opengl/raytracing/HybridDemoSsboTriangles.java
+++ b/src/org/lwjgl/demo/opengl/raytracing/HybridDemoSsboTriangles.java
@@ -298,8 +298,7 @@ private void createSceneSSBO() {
             objects.add(obj);
         }
         Std430Writer.write(objects, GPUObject.class, objectsBuffer);
-        objectsBuffer.flip();
-        glBufferData(GL_ARRAY_BUFFER, objectsBuffer.bb, GL_STATIC_DRAW);
+        nglBufferData(GL_ARRAY_BUFFER, objectsBuffer.pos, objectsBuffer.addr, GL_STATIC_DRAW);
         glBindBuffer(GL_ARRAY_BUFFER, 0);
     }
 
diff --git a/src/org/lwjgl/demo/opengl/raytracing/tutorial/Tutorial6.java b/src/org/lwjgl/demo/opengl/raytracing/tutorial/Tutorial6.java
index 7e4f3b02..706fafca 100644
--- a/src/org/lwjgl/demo/opengl/raytracing/tutorial/Tutorial6.java
+++ b/src/org/lwjgl/demo/opengl/raytracing/tutorial/Tutorial6.java
@@ -502,18 +502,16 @@ private void createSceneSSBOs() {
         DynamicByteBuffer nodesBuffer = new DynamicByteBuffer();
         DynamicByteBuffer trianglesBuffer = new DynamicByteBuffer();
         bhvToBuffers(root, nodesBuffer, trianglesBuffer);
-        nodesBuffer.flip();
-        trianglesBuffer.flip();
         /*
          * And finally we upload the two buffers to the SSBOs.
          */
         this.nodesSsbo = glGenBuffers();
         glBindBuffer(GL_ARRAY_BUFFER, nodesSsbo);
-        glBufferData(GL_ARRAY_BUFFER, nodesBuffer.bb, GL_STATIC_DRAW);
+        nglBufferData(GL_ARRAY_BUFFER, nodesBuffer.pos, nodesBuffer.addr, GL_STATIC_DRAW);
         glBindBuffer(GL_ARRAY_BUFFER, 0);
         this.trianglesSsbo = glGenBuffers();
         glBindBuffer(GL_ARRAY_BUFFER, trianglesSsbo);
-        glBufferData(GL_ARRAY_BUFFER, trianglesBuffer.bb, GL_STATIC_DRAW);
+        nglBufferData(GL_ARRAY_BUFFER, trianglesBuffer.pos, trianglesBuffer.addr, GL_STATIC_DRAW);
         glBindBuffer(GL_ARRAY_BUFFER, 0);
     }
 
diff --git a/src/org/lwjgl/demo/opengl/raytracing/tutorial/Tutorial6_2.java b/src/org/lwjgl/demo/opengl/raytracing/tutorial/Tutorial6_2.java
index 34723a9b..28c26745 100644
--- a/src/org/lwjgl/demo/opengl/raytracing/tutorial/Tutorial6_2.java
+++ b/src/org/lwjgl/demo/opengl/raytracing/tutorial/Tutorial6_2.java
@@ -497,18 +497,16 @@ private void createSceneSSBOs() {
         DynamicByteBuffer nodesBuffer = new DynamicByteBuffer();
         DynamicByteBuffer trianglesBuffer = new DynamicByteBuffer();
         bhvToBuffers(root, nodesBuffer, trianglesBuffer);
-        nodesBuffer.flip();
-        trianglesBuffer.flip();
         /*
          * And finally we upload the two buffers to the SSBOs.
          */
         this.nodesSsbo = glGenBuffers();
         glBindBuffer(GL_ARRAY_BUFFER, nodesSsbo);
-        glBufferData(GL_ARRAY_BUFFER, nodesBuffer.bb, GL_STATIC_DRAW);
+        nglBufferData(GL_ARRAY_BUFFER, nodesBuffer.pos, nodesBuffer.addr, GL_STATIC_DRAW);
         glBindBuffer(GL_ARRAY_BUFFER, 0);
         this.trianglesSsbo = glGenBuffers();
         glBindBuffer(GL_ARRAY_BUFFER, trianglesSsbo);
-        glBufferData(GL_ARRAY_BUFFER, trianglesBuffer.bb, GL_STATIC_DRAW);
+        nglBufferData(GL_ARRAY_BUFFER, trianglesBuffer.pos, trianglesBuffer.addr, GL_STATIC_DRAW);
         glBindBuffer(GL_ARRAY_BUFFER, 0);
     }
 
diff --git a/src/org/lwjgl/demo/opengl/raytracing/tutorial/Tutorial7.java b/src/org/lwjgl/demo/opengl/raytracing/tutorial/Tutorial7.java
index b905956d..f43dacbc 100644
--- a/src/org/lwjgl/demo/opengl/raytracing/tutorial/Tutorial7.java
+++ b/src/org/lwjgl/demo/opengl/raytracing/tutorial/Tutorial7.java
@@ -495,16 +495,14 @@ private void createSceneSSBOs() {
         DynamicByteBuffer nodesBuffer = new DynamicByteBuffer();
         DynamicByteBuffer trianglesBuffer = new DynamicByteBuffer();
         kdTreeToBuffers(kdtree, nodesBuffer, trianglesBuffer);
-        nodesBuffer.flip();
-        trianglesBuffer.flip();
 
         this.nodesSsbo = glGenBuffers();
         glBindBuffer(GL_ARRAY_BUFFER, nodesSsbo);
-        glBufferData(GL_ARRAY_BUFFER, nodesBuffer.bb, GL_STATIC_DRAW);
+        nglBufferData(GL_ARRAY_BUFFER, nodesBuffer.pos, nodesBuffer.addr, GL_STATIC_DRAW);
         glBindBuffer(GL_ARRAY_BUFFER, 0);
         this.trianglesSsbo = glGenBuffers();
         glBindBuffer(GL_ARRAY_BUFFER, trianglesSsbo);
-        glBufferData(GL_ARRAY_BUFFER, trianglesBuffer.bb, GL_STATIC_DRAW);
+        nglBufferData(GL_ARRAY_BUFFER, trianglesBuffer.pos, trianglesBuffer.addr, GL_STATIC_DRAW);
         glBindBuffer(GL_ARRAY_BUFFER, 0);
     }
 
diff --git a/src/org/lwjgl/demo/opengl/util/Boundable.java b/src/org/lwjgl/demo/opengl/util/Boundable.java
new file mode 100644
index 00000000..a9af51f6
--- /dev/null
+++ b/src/org/lwjgl/demo/opengl/util/Boundable.java
@@ -0,0 +1,9 @@
+package org.lwjgl.demo.opengl.util;
+
+public interface Boundable<T> {
+    int min(int axis);
+    int max(int axis);
+    void extend(int axis);
+    T splitLeft(int splitAxis, int splitPos);
+    T splitRight(int splitAxis, int splitPos);
+}
diff --git a/src/org/lwjgl/demo/opengl/util/DynamicByteBuffer.java b/src/org/lwjgl/demo/opengl/util/DynamicByteBuffer.java
index 184e91b0..abe2529e 100644
--- a/src/org/lwjgl/demo/opengl/util/DynamicByteBuffer.java
+++ b/src/org/lwjgl/demo/opengl/util/DynamicByteBuffer.java
@@ -1,12 +1,8 @@
-/*
- * Copyright LWJGL. All rights reserved.
- * License terms: https://www.lwjgl.org/license
- */
 package org.lwjgl.demo.opengl.util;
 
-import java.nio.ByteBuffer;
+import static org.lwjgl.system.MemoryUtil.*;
 
-import org.lwjgl.system.*;
+import java.nio.*;
 
 /**
  * Dynamically growable {@link ByteBuffer}.
@@ -15,70 +11,76 @@
  */
 public class DynamicByteBuffer {
 
-    public ByteBuffer bb;
+    public long addr;
+    public int pos;
+    public int cap;
 
     public DynamicByteBuffer() {
         this(8192);
     }
+
     public DynamicByteBuffer(int initialSize) {
-        bb = MemoryUtil.memAlloc(initialSize);
+        addr = nmemAlloc(initialSize);
+        cap = initialSize;
     }
 
     private void grow() {
-        ByteBuffer newbb = MemoryUtil.memRealloc(bb, (int) (bb.capacity() * 1.5));
-        bb = newbb;
+        int newCap = (int) (cap * 1.5f);
+        long newAddr = nmemRealloc(addr, newCap);
+        cap = newCap;
+        addr = newAddr;
     }
 
     public void free() {
-        MemoryUtil.memFree(bb);
+        nmemFree(addr);
     }
 
     public DynamicByteBuffer putFloat(float v) {
-        if (bb.remaining() < 4) {
+        if (cap - pos < 4)
             grow();
-        }
-        bb.putFloat(v);
+        memPutFloat(addr + pos, v);
+        pos += 4;
         return this;
     }
 
     public DynamicByteBuffer putLong(long v) {
-        if (bb.remaining() < 8) {
+        if (cap - pos < 8)
             grow();
-        }
-        bb.putLong(v);
+        memPutLong(addr + pos, v);
+        pos += 8;
         return this;
     }
 
     public DynamicByteBuffer putInt(int v) {
-        if (bb.remaining() < 4) {
+        if (cap - pos < 4)
             grow();
-        }
-        bb.putInt(v);
+        memPutInt(addr + pos, v);
+        pos += 4;
         return this;
     }
 
     public DynamicByteBuffer putShort(int v) {
-        if (bb.remaining() < 2) {
+        if (v > 1 << 16)
+            throw new IllegalArgumentException();
+        if (cap - pos < 2)
             grow();
-        }
-        bb.putShort((short) (v & 0xFFFF));
+        memPutShort(addr + pos, (short) v);
+        pos += 2;
         return this;
     }
 
     public DynamicByteBuffer putByte(int v) {
-        if (bb.remaining() < 1) {
+        if (v > 255)
+            throw new IllegalArgumentException();
+        if (cap - pos < 2)
             grow();
-        }
-        bb.put((byte) (v & 0xFF));
+        memPutByte(addr + pos, (byte) (v & 0xFF));
+        pos++;
         return this;
     }
 
-    public void flip() {
-        bb.flip();
-    }
-
     public int remaining() {
-        return bb.remaining();
+        return (int) (cap - pos);
     }
 
 }
diff --git a/src/org/lwjgl/demo/opengl/util/KDTreei.java b/src/org/lwjgl/demo/opengl/util/KDTreei.java
new file mode 100644
index 00000000..250c5af4
--- /dev/null
+++ b/src/org/lwjgl/demo/opengl/util/KDTreei.java
@@ -0,0 +1,538 @@
+package org.lwjgl.demo.opengl.util;
+
+import java.lang.Math;
+import java.nio.*;
+import java.util.*;
+import org.joml.*;
+
+public class KDTreei<T extends Boundable<T>> {
+
+    public Node<T> root;
+    private int maxVoxelCount = 2;
+    private float nodeIntersectCosts = 1.0f;
+    private float voxelIntersectCosts = 1.0f;
+
+    private class IntervalBoundary {
+        int type;
+        int pos;
+
+        IntervalBoundary(int t, int p) {
+            type = t;
+            pos = p;
+        }
+
+        int compareTo(IntervalBoundary sib) {
+            return Integer.compare(pos, sib.pos);
+        }
+    }
+
+    public static class Box implements Boundable<Box> {
+        public int minX, minY, minZ;
+        public int maxX, maxY, maxZ;
+
+        public Box() {
+        }
+
+        public Box(int minX, int minY, int minZ, int maxX, int maxY, int maxZ) {
+            this.minX = minX;
+            this.minY = minY;
+            this.minZ = minZ;
+            this.maxX = maxX;
+            this.maxY = maxY;
+            this.maxZ = maxZ;
+        }
+
+        public Box(Box bbox) {
+            this.minX = bbox.minX;
+            this.minY = bbox.minY;
+            this.minZ = bbox.minZ;
+            this.maxX = bbox.maxX;
+            this.maxY = bbox.maxY;
+            this.maxZ = bbox.maxZ;
+        }
+
+        public float diagonal() {
+            float dx = maxX - minX;
+            float dy = maxY - minY;
+            float dz = maxZ - minZ;
+            return (float) Math.sqrt(dx * dx + dy * dy + dz * dz);
+        }
+
+        public int min(int axis) {
+            switch (axis) {
+            case 0:
+                return minX;
+            case 1:
+                return minY;
+            case 2:
+                return minZ;
+            default:
+                throw new IllegalArgumentException();
+            }
+        }
+
+        public int max(int axis) {
+            switch (axis) {
+            case 0:
+                return maxX;
+            case 1:
+                return maxY;
+            case 2:
+                return maxZ;
+            default:
+                throw new IllegalArgumentException();
+            }
+        }
+
+        public void setMax(int splitAxis, int split) {
+            switch (splitAxis) {
+            case 0:
+                maxX = split;
+                break;
+            case 1:
+                maxY = split;
+                break;
+            case 2:
+                maxZ = split;
+                break;
+            }
+        }
+
+        public void setMin(int splitAxis, int split) {
+            switch (splitAxis) {
+            case 0:
+                minX = split;
+                break;
+            case 1:
+                minY = split;
+                break;
+            case 2:
+                minZ = split;
+                break;
+            }
+        }
+
+        public boolean intersectsWithBox(Boundable<?> vx) {
+            return maxX >= vx.min(0) && maxY >= vx.min(1) && maxZ >= vx.min(2) && minX <= vx.max(0) && minY <= vx.max(1)
+                    && minZ <= vx.max(2);
+        }
+
+        @Override
+        public void extend(int axis) {
+            throw new UnsupportedOperationException();
+        }
+
+        @Override
+        public Box splitLeft(int splitAxis, int splitPos) {
+            throw new UnsupportedOperationException();
+        }
+
+        @Override
+        public Box splitRight(int splitAxis, int splitPos) {
+            throw new UnsupportedOperationException();
+        }
+    }
+
+    public static class Voxel implements Boundable<Voxel>, Comparable<Voxel> {
+        public byte x, y, z;
+        public byte ex, ey, ez;
+        public int morton;
+        public byte paletteIndex;
+        public int nindex;
+
+        public Voxel() {
+        }
+
+        public Voxel(byte x, byte y, byte z, byte paletteIndex) {
+            this.x = x;
+            this.y = y;
+            this.z = z;
+            this.paletteIndex = paletteIndex;
+        }
+
+        public Voxel(byte x, byte y, byte z, byte ex, byte ey, byte ez, byte paletteIndex) {
+            this.x = x;
+            this.y = y;
+            this.z = z;
+            this.ex = ex;
+            this.ey = ey;
+            this.ez = ez;
+            this.paletteIndex = paletteIndex;
+        }
+
+        public int min(int axis) {
+            switch (axis) {
+            case 0:
+                return x & 0xFF;
+            case 1:
+                return y & 0xFF;
+            case 2:
+                return z & 0xFF;
+            default:
+                throw new IllegalArgumentException();
+            }
+        }
+
+        public int max(int axis) {
+            switch (axis) {
+            case 0:
+                return (x & 0xFF) + 1 + (ex & 0xFF);
+            case 1:
+                return (y & 0xFF) + 1 + (ey & 0xFF);
+            case 2:
+                return (z & 0xFF) + 1 + (ez & 0xFF);
+            default:
+                throw new IllegalArgumentException();
+            }
+        }
+
+        @Override
+        public void extend(int axis) {
+            switch (axis) {
+            case 0:
+                this.ex++;
+                break;
+            case 1:
+                this.ey++;
+                break;
+            case 2:
+                this.ez++;
+                break;
+            default:
+                throw new IllegalArgumentException();
+            }
+        }
+
+        @Override
+        public Voxel splitLeft(int axis, int pos) {
+            switch (axis) {
+            case 0:
+                return new Voxel(x, y, z, (byte) (pos - (x & 0xFF) - 1), ey, ez, paletteIndex);
+            case 1:
+                return new Voxel(x, y, z, ex, (byte) (pos - (y & 0xFF) - 1), ez, paletteIndex);
+            case 2:
+                return new Voxel(x, y, z, ex, ey, (byte) (pos - (z & 0xFF) - 1), paletteIndex);
+            default:
+                throw new IllegalArgumentException();
+            }
+        }
+
+        @Override
+        public Voxel splitRight(int axis, int pos) {
+            switch (axis) {
+            case 0:
+                return new Voxel((byte) pos, y, z, (byte) ((ex & 0xFF) - (pos - (x & 0xFF))), ey, ez, paletteIndex);
+            case 1:
+                return new Voxel(x, (byte) pos, z, ex, (byte) ((ey & 0xFF) - (pos - (y & 0xFF))), ez, paletteIndex);
+            case 2:
+                return new Voxel(x, y, (byte) pos, ex, ey, (byte) ((ez & 0xFF) - (pos - (z & 0xFF))), paletteIndex);
+            default:
+                throw new IllegalArgumentException();
+            }
+        }
+
+        @Override
+        public int compareTo(Voxel o) {
+            return Long.compare(morton, o.morton);
+        }
+
+        @Override
+        public String toString() {
+            return "[" + x + ", " + y + ", " + z + "]";
+        }
+    }
+
+    public static class Node<B> {
+        private static final int SIDE_X_POS = 0;
+        private static final int SIDE_X_NEG = 1;
+        private static final int SIDE_Y_POS = 2;
+        private static final int SIDE_Y_NEG = 3;
+        private static final int SIDE_Z_POS = 4;
+        private static final int SIDE_Z_NEG = 5;
+
+        public int splitAxis = -1;
+        public int splitPos;
+        public Box boundingBox;
+        public Node<B> left;
+        public Node<B> right;
+        public List<B> voxels = new ArrayList<>();
+        public Node<B>[] ropes;
+        public int index;
+        public int leafIndex;
+        public int first;
+        public int count;
+
+        int isParallelTo(int side) {
+            switch (splitAxis) {
+            case 0:
+                return side == SIDE_X_NEG ? -1 : side == SIDE_X_POS ? +1 : 0;
+            case 1:
+                return side == SIDE_Y_NEG ? -1 : side == SIDE_Y_POS ? +1 : 0;
+            case 2:
+                return side == SIDE_Z_NEG ? -1 : side == SIDE_Z_POS ? +1 : 0;
+            default:
+                throw new IllegalArgumentException();
+            }
+        }
+
+        public final boolean isLeafNode() {
+            return splitAxis == -1;
+        }
+
+        @SuppressWarnings("unchecked")
+        protected void processNode(Node<B>[] ropes) {
+            if (isLeafNode()) {
+                this.ropes = ropes;
+            } else {
+                int sideLeft;
+                int sideRight;
+                if (splitAxis == 0) {
+                    sideLeft = SIDE_X_NEG;
+                    sideRight = SIDE_X_POS;
+                } else if (splitAxis == 1) {
+                    sideLeft = SIDE_Y_NEG;
+                    sideRight = SIDE_Y_POS;
+                } else if (splitAxis == 2) {
+                    sideLeft = SIDE_Z_NEG;
+                    sideRight = SIDE_Z_POS;
+                } else {
+                    throw new AssertionError();
+                }
+                this.left.ropes = new Node[6];
+                System.arraycopy(ropes, 0, this.left.ropes, 0, 6);
+                this.left.ropes[sideRight] = this.right;
+                this.left.processNode(this.left.ropes);
+                this.right.ropes = new Node[6];
+                System.arraycopy(ropes, 0, this.right.ropes, 0, 6);
+                this.right.ropes[sideLeft] = this.left;
+                this.right.processNode(this.right.ropes);
+            }
+        }
+
+        protected void optimizeRopes() {
+            for (int i = 0; i < 6; i++) {
+                ropes[i] = optimizeRope(ropes[i], i);
+            }
+            if (left != null)
+                left.optimizeRopes();
+            if (right != null)
+                right.optimizeRopes();
+        }
+
+        protected Node<B> optimizeRope(Node<B> rope, int side) {
+            if (rope == null) {
+                return rope;
+            }
+            Node<B> r = rope;
+            while (!r.isLeafNode()) {
+                int parallelSide = r.isParallelTo(side);
+                if (parallelSide == +1) {
+                    r = r.left;
+                } else if (parallelSide == -1) {
+                    r = r.right;
+                } else {
+                    if (r.splitPos < boundingBox.min(r.splitAxis)) {
+                        r = r.right;
+                    } else if (r.splitPos > boundingBox.max(r.splitAxis)) {
+                        r = r.left;
+                    } else {
+                        break;
+                    }
+                }
+            }
+            return r;
+        }
+
+        public void frustumCull(Matrix4d viewProjection, Vector3d p, int depth, int maxDepth, int maxNodes,
+                PriorityQueue<Node<B>> nodes) {
+            Box b = boundingBox;
+            if (nodes.size() > maxNodes)
+                return;
+            if (!viewProjection.testAab(b.minX, b.minY, b.minZ, b.maxX, b.maxY, b.maxZ))
+                return;
+            if (left == null || depth >= maxDepth) {
+                if (this.count > 0)
+                    nodes.add(this);
+            } else if (left != null) {
+                if (p.get(splitAxis) < splitPos) {
+                    left.frustumCull(viewProjection, p, depth + 1, maxDepth, maxNodes, nodes);
+                    right.frustumCull(viewProjection, p, depth + 1, maxDepth, maxNodes, nodes);
+                } else {
+                    right.frustumCull(viewProjection, p, depth + 1, maxDepth, maxNodes, nodes);
+                    left.frustumCull(viewProjection, p, depth + 1, maxDepth, maxNodes, nodes);
+                }
+            }
+        }
+
+        public Node<B> findNode(Vector3d cameraPosition) {
+            Vector3d p = cameraPosition;
+            Box b = boundingBox;
+            if (p.x < b.minX || p.x > b.maxX || p.y < b.minY || p.y > b.maxY || p.z < b.minZ || p.z > b.maxZ
+                    || left == null)
+                return this;
+            if (cameraPosition.get(splitAxis) < splitPos)
+                return left.findNode(cameraPosition);
+            return right.findNode(cameraPosition);
+        }
+    }
+
+    @SuppressWarnings("unchecked")
+    public static <T extends Boundable<T>> KDTreei<T> build(List<T> voxels, int maxDepth) {
+        return build(voxels, new Node[6], maxDepth);
+    }
+
+    public static <T extends Boundable<T>> KDTreei<T> build(List<T> voxels, Node<T>[] neighbors, int maxDepth) {
+        Box b = new Box(Integer.MAX_VALUE, Integer.MAX_VALUE, Integer.MAX_VALUE, Integer.MIN_VALUE, Integer.MIN_VALUE,
+                Integer.MIN_VALUE);
+        for (T v : voxels) {
+            b.minX = b.minX < v.min(0) ? b.minX : v.min(0);
+            b.minY = b.minY < v.min(1) ? b.minY : v.min(1);
+            b.minZ = b.minZ < v.min(2) ? b.minZ : v.min(2);
+            b.maxX = b.maxX > v.max(0) ? b.maxX : v.max(0);
+            b.maxY = b.maxY > v.max(1) ? b.maxY : v.max(1);
+            b.maxZ = b.maxZ > v.max(2) ? b.maxZ : v.max(2);
+        }
+        KDTreei<T> root = new KDTreei<T>();
+        root.buildTree(voxels, b, neighbors, maxDepth);
+        return root;
+    }
+
+    public void frustumCull(Matrix4d viewProjection, Vector3d p, int maxDepth, int maxNodes,
+            PriorityQueue<Node<T>> nodes) {
+        root.frustumCull(viewProjection, p, 0, maxDepth, maxNodes, nodes);
+    }
+
+    public Node<T> findNode(Vector3d cameraPosition) {
+        return root.findNode(cameraPosition);
+    }
+
+    public void level(int level, ShortBuffer nodes) {
+        level(root, level, nodes);
+    }
+
+    private void level(Node<T> n, int level, ShortBuffer nodes) {
+        if (level == 0 || n.left == null) {
+            nodes.put((short) n.index);
+        } else {
+            level(n.left, level - 1, nodes);
+            level(n.right, level - 1, nodes);
+        }
+    }
+
+    public void buildTree(List<T> list, Box bbox, Node<T>[] neighbors, int maxDepth) {
+        if (root != null)
+            root = null;
+        root = new Node<T>();
+        root.voxels = list;
+        root.boundingBox = bbox;
+        buildTree(root, 0, 0, maxDepth);
+        root.processNode(root.ropes = neighbors);
+        root.optimizeRopes();
+    }
+
+    private void buildTree(Node<T> node, int axis, int depth, int maxDepth) {
+        if (node.voxels.size() > maxVoxelCount && depth < maxDepth) {
+            node.splitAxis = axis;
+            node.splitPos = findSplitPlane(node);
+        } else {
+            node.splitAxis = -1;
+        }
+        if (node.splitAxis == -1) {
+            return;
+        }
+        if (node.voxels.size() > maxVoxelCount) {
+            node.left = new Node<T>();
+            node.right = new Node<T>();
+            node.left.boundingBox = new Box(node.boundingBox);
+            node.left.boundingBox.setMax(node.splitAxis, node.splitPos);
+            node.right.boundingBox = new Box(node.boundingBox);
+            node.right.boundingBox.setMin(node.splitAxis, node.splitPos);
+            node.voxels.forEach(vx -> {
+                if (vx.min(node.splitAxis) >= node.splitPos) {
+                    node.right.voxels.add(vx);
+                } else if (vx.max(node.splitAxis) <= node.splitPos) {
+                    node.left.voxels.add(vx);
+                } else {
+                    T left = vx.splitLeft(node.splitAxis, node.splitPos);
+                    T right = vx.splitRight(node.splitAxis, node.splitPos);
+                    if (left.max(node.splitAxis) > node.splitPos)
+                        throw new AssertionError();
+                    if (right.min(node.splitAxis) < node.splitPos)
+                        throw new AssertionError();
+                    node.left.voxels.add(left);
+                    node.right.voxels.add(right);
+                }
+            });
+            node.voxels.clear();
+            int nextAxis = (axis + 1) % 3;
+            buildTree(node.left, nextAxis, depth + 1, maxDepth);
+            buildTree(node.right, nextAxis, depth + 1, maxDepth);
+        }
+    }
+
+    private int findSplitPlane(Node<T> node) {
+        if (node == null)
+            return -1;
+        Box bb = node.boundingBox;
+        int nPrims = node.voxels.size();
+        int xw = bb.maxX - bb.minX;
+        int yw = bb.maxY - bb.minY;
+        int zw = bb.maxZ - bb.minZ;
+        int box_width;
+        int ax;
+        if (xw > yw && xw > zw) {
+            ax = 0;
+            box_width = xw;
+        } else if (yw > zw) {
+            ax = 1;
+            box_width = yw;
+        } else {
+            ax = 2;
+            box_width = zw;
+        }
+        float inv_box_width = 1.0f / box_width;
+        List<IntervalBoundary> intervals = new ArrayList<IntervalBoundary>();
+        final int count = node.voxels.size();
+        final int divisor = (int) Math.ceil(count / 100.0);
+        nPrims /= divisor;
+        for (int i = 0; i < count; i += divisor) {
+            T vx = node.voxels.get(i);
+            if (!bb.intersectsWithBox(vx)) {
+                throw new IllegalStateException("!!! KDTree.findSplitPlane: no intersection of boxes");
+            }
+            intervals.add(new IntervalBoundary(0, vx.min(ax)));
+            intervals.add(new IntervalBoundary(1, vx.max(ax)));
+        }
+        Collections.sort(intervals, (sib1, sib2) -> sib1.compareTo(sib2));
+        int done_intervals = 0;
+        int open_intervals = 0;
+        float alpha;
+        int minid = 0;
+        float mincost = Float.MAX_VALUE;
+        for (int i = 0; i < intervals.size(); i++) {
+            IntervalBoundary in = intervals.get(i);
+            if (in.type == 1) {
+                open_intervals--;
+                done_intervals++;
+            }
+            alpha = (in.pos - bb.min(ax)) * inv_box_width;
+            float cost = voxelIntersectCosts + nodeIntersectCosts
+                    * ((done_intervals + open_intervals) * alpha + (nPrims - done_intervals) * (1.0f - alpha));
+            if (cost < mincost) {
+                minid = i;
+                mincost = cost;
+            }
+            if (in.type == 0) {
+                open_intervals++;
+            }
+        }
+        int splitPlane = intervals.get(minid).pos;
+        if (splitPlane == bb.min(ax) || splitPlane == bb.max(ax)) {
+            node.splitAxis = -1;
+            return -1;
+        }
+        node.splitAxis = ax;
+        return intervals.get(minid).pos;
+    }
+
+}