normals calculation fix

[voxconv.git] / vox_data.d

module vox_data;

import iv.bclamp;
import iv.glbinds.utils;
import iv.cmdcongl;
import iv.strex;
import iv.vfs;
import iv.vfs.io;
import iv.vmath;

static import iv.timer;

import vox_3dbmp;

version = voxdata_debug;
//version = vox_check_invariants;


// ////////////////////////////////////////////////////////////////////////// //
static align(1) struct VoxXYZ16 {
align(1):
  ushort x = void, y = void, z = void;
}


// ////////////////////////////////////////////////////////////////////////// //
/*
  this is a struct that keeps info about an individual voxel.
  it keeps voxel color, and face visibility info.
 */
align(1) struct VoxPix {
align(1):
  ubyte b, g, r;
  ubyte cull;
  uint nextz; // voxel with the next z; 0 means "no more"
  ushort z; // z of the current voxel

  uint rgb () const pure nothrow @safe @nogc {
    pragma(inline, true);
    return 0xff000000U|b|(cast(uint)g<<8)|(cast(uint)r<<16);
  }

  uint rgbcull () const pure nothrow @safe @nogc {
    pragma(inline, true);
    return b|(cast(uint)g<<8)|(cast(uint)r<<16)|(cast(uint)cull<<24);
  }
}


// ////////////////////////////////////////////////////////////////////////// //
/*
  this keeps voxel "voxmap" (it's like a pixmap, but with voxels instead
  of pixels). the representation is slightly optimised, tho: it keeps only
  actually used voxels, in vertical "slabs". as most voxel models have only
  contour voxels stored, this greatly reduces memory consumption. quering
  the individial voxel data is slightly slower, tho, but for our case it is
  not really important.

  there are some methods to "optimise" the voxmap. for example, to fix voxel
  face visibility info, and to remove "internal" voxels (it is called "hollow fill").

  also note that some optimisation methods may leave empty voxels in the voxmap.
  they have `cull` field equal to zero.

  to build a voxmap you simply call `addVoxel()` method, optionally providing
  the face visibility info. but you can use `0x3f` as a visibility, and then ask
  `VoxelData` object to calculate the proprer "cull" values for you. you can also
  add "internal" voxels, and then ask the object to fix the vixmap for you, removing
  all the unnecessary data.
 */
struct VoxelData {
public:
  uint xsize = 0, ysize = 0, zsize = 0;
  float cx = 0.0f, cy = 0.0f, cz = 0.0f;

  VoxPix[] data; // [0] is never used
  // xsize*ysize array, offsets in `data`; 0 means "no data here"
  // slabs are sorted from bottom to top, and never intersects
  uint[] xyofs;
  uint freelist = 0;
  uint voxpixtotal = 0;

public:
  void save (VFile fl) {
    fl.rawWriteExact("K8VOXDT0");
    fl.writeNum!ushort(cast(ushort)xsize);
    fl.writeNum!ushort(cast(ushort)ysize);
    fl.writeNum!ushort(cast(ushort)zsize);
    fl.writeNum(cx);
    fl.writeNum(cy);
    fl.writeNum(cz);

    uint voxelCount = 0;
    foreach (int y; 0..ysize) {
      foreach (int x; 0..xsize) {
        for (uint dofs = getDOfs(x, y); dofs; dofs = data[dofs].nextz) {
          if (data[dofs].cull) ++voxelCount;
        }
      }
    }
    fl.writeNum(voxelCount);

    foreach (int y; 0..ysize) {
      foreach (int x; 0..xsize) {
        uint count = 0;
        for (uint dofs = getDOfs(x, y); dofs; dofs = data[dofs].nextz) {
          if (data[dofs].cull) ++count;
        }
        fl.writeNum!ushort(cast(ushort)count);
        for (uint dofs = getDOfs(x, y); dofs; dofs = data[dofs].nextz) {
          if (data[dofs].cull) {
            fl.writeNum(data[dofs].b);
            fl.writeNum(data[dofs].g);
            fl.writeNum(data[dofs].r);
            fl.writeNum(data[dofs].cull);
            fl.writeNum(data[dofs].z);
          }
        }
      }
    }
  }

  void load (VFile fl) {
    char[8] sign = void;
    fl.rawReadExact(sign[]);
    if (sign[] != "K8VOXDT0") throw new Exception("invalid voxel data");

    xsize = fl.readNum!ushort;
    ysize = fl.readNum!ushort;
    zsize = fl.readNum!ushort;
    cx = fl.readNum!float;
    cy = fl.readNum!float;
    cz = fl.readNum!float;

    const uint voxelCount = fl.readNum!uint;
    delete xyofs;
    xyofs = new uint[xsize*ysize];
    delete data;
    data = new VoxPix[voxelCount+1]; // [0] is unused

    freelist = 0;
    voxpixtotal = voxelCount;

    uint cpos = 1;
    foreach (int y; 0..ysize) {
      foreach (int x; 0..xsize) {
        uint count = fl.readNum!ushort;
        if (count == 0) {
          xyofs[cast(uint)y*xsize+cast(uint)x] = 0;
          continue;
        }
        xyofs[cast(uint)y*xsize+cast(uint)x] = cpos;
        int prevz = -1;
        while (count--) {
          VoxPix *vx = &data[cpos];
          vx.b = fl.readNum!ubyte;
          vx.g = fl.readNum!ubyte;
          vx.r = fl.readNum!ubyte;
          vx.cull = fl.readNum!ubyte;
          vx.nextz = -1;
          if (!vx.cull) throw new Exception("invalid voxel data (cull)");
          vx.z = fl.readNum!ushort;
          if (prevz >= 0) {
            if (data[prevz].z >= vx.z) throw new Exception("invalid voxel data (z)");
            data[prevz].nextz = cast(int)cpos;
          }
          prevz = cast(int)cpos;
          ++cpos;
        }
      }
    }
    assert(cpos == data.length);
  }

private:
  uint allocVox () nothrow @safe {
    assert(data.length);
    ++voxpixtotal;
    if (!freelist) {
      if (data.length >= 0x3fffffffU) assert(0, "too many voxels");
      const uint lastel = cast(uint)data.length;
      data.length += 1024;
      freelist = cast(uint)data.length-1;
      while (freelist >= lastel) {
        data[freelist].nextz = freelist+1;
        --freelist;
      }
      freelist = lastel;
      data[$-1].nextz = 0;
    }
    const uint res = freelist;
    freelist = data[res].nextz;
    return res;
  }

public:
  void clear () {
    delete data;
    data = null;
    delete xyofs;
    xyofs = null;
    xsize = ysize = zsize = 0;
    cx = cy = cz = 0.0f;
    freelist = 0;
    voxpixtotal = 0;
  }

  void setSize (uint xs, uint ys, uint zs) {
    clear();
    if (!xs || !ys || !zs) return;
    xsize = xs;
    ysize = ys;
    zsize = zs;
    xyofs = new uint[xsize*ysize];
    xyofs[] = 0;
    data.length = 1; // data[0] is never used
  }

  uint getDOfs (int x, int y) const nothrow @trusted @nogc {
    pragma(inline, true);
    if (x < 0 || y < 0) return 0;
    if (cast(uint)x >= xsize || cast(uint)y >= ysize) return 0;
    return xyofs.ptr[cast(uint)y*xsize+cast(uint)x];
  }

  // high byte is cull info
  // returns 0 if there is no such voxel
  uint voxofs (int x, int y, int z) const nothrow @safe @nogc {
    uint dofs = getDOfs(x, y);
    while (dofs) {
      if (data[dofs].z == cast(ushort)z) return dofs;
      if (data[dofs].z > cast(ushort)z) return 0;
      dofs = data[dofs].nextz;
    }
    return 0;
  }

  // high byte is cull info
  // returns 0 if there is no such voxel
  uint query (int x, int y, int z) const nothrow @safe @nogc {
    immutable uint dofs = voxofs(x, y, z);
    if (!dofs) return 0;
    if (!data[dofs].cull) return 0;
    return data[dofs].rgbcull();
  }

  VoxPix *queryVP (int x, int y, int z) nothrow @trusted @nogc {
    pragma(inline, true);
    immutable uint dofs = voxofs(x, y, z);
    return (dofs ? &data[dofs] : null);
  }

  // high byte is cull info
  // returns 0 if there is no such voxel
  ubyte queryCull (int x, int y, int z) const nothrow @safe @nogc {
    immutable uint dofs = voxofs(x, y, z);
    return (dofs ? data[dofs].cull : 0);
  }

  void removeVoxel (int x, int y, int z) nothrow @safe @nogc {
    uint dofs = getDOfs(x, y);
    uint prevdofs = 0;
    while (dofs) {
      if (data[dofs].z == cast(ushort)z) {
        // remove this voxel
        if (prevdofs) {
          data[prevdofs].nextz = data[dofs].nextz;
        } else {
          xyofs[cast(uint)y*xsize+cast(uint)x] = data[dofs].nextz;
        }
        data[dofs].nextz = freelist;
        freelist = dofs;
        --voxpixtotal;
        return;
      }
      if (data[dofs].z > cast(ushort)z) return;
      prevdofs = dofs;
      dofs = data[dofs].nextz;
    }
  }

  void addVoxel (int x, int y, int z, uint rgb, ubyte cull) nothrow @safe {
    cull &= 0x3f;
    if (!cull) { removeVoxel(x, y, z); return; }
    if (x < 0 || y < 0 || z < 0) return;
    if (cast(uint)x >= xsize || cast(uint)y >= ysize || cast(uint)z >= zsize) return;
    uint dofs = getDOfs(x, y);
    uint prevdofs = 0;
    while (dofs) {
      if (data[dofs].z == cast(ushort)z) {
        // replace this voxel
        data[dofs].b = rgb&0xff;
        data[dofs].g = (rgb>>8)&0xff;
        data[dofs].r = (rgb>>16)&0xff;
        data[dofs].cull = cull;
        return;
      }
      if (data[dofs].z > cast(ushort)z) break;
      prevdofs = dofs;
      dofs = data[dofs].nextz;
    }
    // insert before dofs
    immutable uint vidx = allocVox();
    data[vidx].b = rgb&0xff;
    data[vidx].g = (rgb>>8)&0xff;
    data[vidx].r = (rgb>>16)&0xff;
    data[vidx].cull = cull;
    data[vidx].z = cast(ushort)z;
    data[vidx].nextz = dofs;
    if (prevdofs) {
      assert(data[prevdofs].nextz == dofs);
      data[prevdofs].nextz = vidx;
    } else {
      xyofs[cast(uint)y*xsize+cast(uint)x] = vidx;
    }
  }

  // only for existing voxels; won't remove empty voxels
  void setVoxelCull (int x, int y, int z, ubyte cull) nothrow @safe @nogc {
    VoxPix *vp = queryVP(x, y, z);
    if (vp) vp.cull = cast(ubyte)(cull&0x3f);
  }

  version(vox_check_invariants)
  void checkInvariants () const nothrow /*@safe*/ @trusted /*@nogc*/ {
    version(voxdata_debug) conwriteln("checking invariants...");
    uint voxcount = 0;
    foreach (uint y; 0..ysize) {
      foreach (uint x; 0..xsize) {
        uint dofs = getDOfs(x, y);
        if (!dofs) continue;
        ++voxcount;
        ushort prevz = data[dofs].z;
        dofs = data[dofs].nextz;
        while (dofs) {
          ++voxcount;
          assert(prevz < data[dofs].z, "broken voxel data Z invariant");
          prevz = data[dofs].z;
          dofs = data[dofs].nextz;
        }
      }
    }
    assert(voxcount == voxpixtotal, "invalid number of voxels");
  }

  void removeEmptyVoxels () nothrow /*@safe*/ @trusted /*@nogc*/ {
    version(voxdata_debug) conwriteln("removing empty voxels...");
    uint count = 0;
    foreach (uint y; 0..ysize) {
      foreach (uint x; 0..xsize) {
        uint dofs = getDOfs(x, y);
        if (!dofs) continue;
        uint prevdofs = 0;
        while (dofs) {
          if (!data[dofs].cull) {
            // remove it
            const uint ndofs = data[dofs].nextz;
            if (prevdofs) {
              data[prevdofs].nextz = ndofs;
            } else {
              xyofs[cast(uint)y*xsize+cast(uint)x] = ndofs;
            }
            data[dofs].nextz = freelist;
            freelist = dofs;
            --voxpixtotal;
            dofs = ndofs;
            ++count;
          } else {
            prevdofs = dofs;
            dofs = data[dofs].nextz;
          }
        }
      }
    }
    if (count) conwriteln("removed ", count, " empty voxel", (count != 1 ? "s" : ""));
  }

  static immutable int[3][6] cullofs = [
    [ 1, 0, 0], // left
    [-1, 0, 0], // right
    [ 0,-1, 0], // near
    [ 0, 1, 0], // far
    [ 0, 0, 1], // top
    [ 0, 0,-1], // bottom
  ];

  static ubyte cullmask (uint cidx) pure nothrow @safe @nogc {
    pragma(inline, true);
    return cast(ubyte)(1U<<cidx);
  }

  // opposite mask
  static ubyte cullopmask (uint cidx) pure nothrow @safe @nogc {
    pragma(inline, true);
    return cast(ubyte)(1U<<(cidx^1));
  }

  // remove inside voxels, leaving only contour
  void removeInsideFaces () nothrow /*@safe*/ @trusted {
    version(voxdata_debug) conwriteln("removing inside voxels...");
    foreach (int y; 0..ysize) {
      foreach (int x; 0..xsize) {
        for (uint dofs = getDOfs(x, y); dofs; dofs = data[dofs].nextz) {
          if (!data[dofs].cull) continue;
          // check
          const int z = cast(int)data[dofs].z;
          foreach (uint cidx; 0..6) {
            // go in this dir, removing the corresponding voxel side
            immutable ubyte cmask = cullmask(cidx);
            immutable ubyte opmask = cullopmask(cidx);
            immutable ubyte checkmask = cmask|opmask;
            immutable int dx = cullofs[cidx][0];
            immutable int dy = cullofs[cidx][1];
            immutable int dz = cullofs[cidx][2];
            int vx = x, vy = y, vz = z;
            uint myofs = dofs;
            while (myofs && (data[myofs].cull&cmask)) {
              immutable int sx = vx+dx;
              immutable int sy = vy+dy;
              immutable int sz = vz+dz;
              immutable uint sofs = voxofs(sx, sy, sz);
              if (!sofs) break;
              if (!(data[sofs].cull&checkmask)) break;
              // fix culls
              data[myofs].cull ^= cmask;
              data[sofs].cull &= cast(ubyte)(~cast(uint)opmask);
              vx = sx;
              vy = sy;
              vz = sz;
              myofs = sofs;
            }
          }
        }
      }
    }
  }

  // if we have ANY voxel at the corresponding side, don't render that face
  // return number of fixed voxels
  uint fixFaceVisibility () nothrow @safe @nogc {
    uint count = 0;
    foreach (int y; 0..ysize) {
      foreach (int x; 0..xsize) {
        for (uint dofs = getDOfs(x, y); dofs; dofs = data[dofs].nextz) {
          const ubyte ocull = data[dofs].cull;
          if (!ocull) continue;
          const int z = cast(int)data[dofs].z;
          // if we have ANY voxel at the corresponding side, don't render that face
          foreach (uint cidx; 0..6) {
            const ubyte cmask = cullmask(cidx);
            if (data[dofs].cull&cmask) {
              if (queryCull(x+cullofs[cidx][0], y+cullofs[cidx][1], z+cullofs[cidx][2])) {
                data[dofs].cull ^= cmask; // reset bit
              }
            }
          }
          count += (data[dofs].cull != ocull);
        }
      }
    }
    return count;
  }

  void create3DBitmap (ref Vox3DBitmap bmp) {
    bmp.setSize(xsize, ysize, zsize);
    foreach (int y; 0..ysize) {
      foreach (int x; 0..xsize) {
        for (uint dofs = getDOfs(x, y); dofs; dofs = data[dofs].nextz) {
          if (data[dofs].cull) bmp.setPixel(x, y, cast(int)data[dofs].z);
        }
      }
    }
  }

  // this fills everything outside of the voxel, and
  // then resets culling bits for all invisible faces
  // i don't care about memory yet
  uint hollowFill () @trusted {
    Vox3DBitmap bmp;
    scope(exit) { bmp.clear(); }
    bmp.setSize(xsize+2, ysize+2, zsize+2);

    VoxXYZ16[] stack;
    scope(exit) delete stack;
    uint stackpos;

    stack.length = 32768;
    stack.assumeSafeAppend;
    stackpos = 0;
    assert(xsize <= cast(uint)stack.length);

    // this is definitely empty
    VoxXYZ16 xyz;
    xyz.x = xyz.y = xyz.z = 0;
    bmp.setPixel(cast(int)xyz.x, cast(int)xyz.y, cast(int)xyz.z);
    stack[stackpos++] = xyz;

    auto tm = iv.timer.Timer(true);
    while (stackpos) {
      xyz = stack[--stackpos];
      for (uint dd = 0; dd < 6; ++dd) {
        const int nx = cast(int)xyz.x+cullofs[dd][0];
        const int ny = cast(int)xyz.y+cullofs[dd][1];
        const int nz = cast(int)xyz.z+cullofs[dd][2];
        if (bmp.setPixel(nx, ny, nz)) continue;
        if (queryCull(nx-1, ny-1, nz-1)) continue;
        if (stackpos == cast(uint)stack.length) {
          stack.length += 32768;
          stack.assumeSafeAppend;
        }
        stack[stackpos++] = VoxXYZ16(cast(ushort)nx, cast(ushort)ny, cast(ushort)nz);
      }
    }
    tm.stop();
    version(voxdata_debug) conwriteln("*** flooded in ", tm.toString(), " ", stack.length, " stack items used.");

    // unmark contour voxels
    // this is required for proper face removing
    foreach (int y; 0..ysize) {
      foreach (int x; 0..xsize) {
        for (uint dofs = getDOfs(x, y); dofs; dofs = data[dofs].nextz) {
          if (!data[dofs].cull) continue;
          const int z = cast(int)data[dofs].z;
          bmp.resetPixel(x+1, y+1, z+1);
        }
      }
    }

    // now check it
    uint changed = 0;
    foreach (int y; 0..ysize) {
      foreach (int x; 0..xsize) {
        for (uint dofs = getDOfs(x, y); dofs; dofs = data[dofs].nextz) {
          immutable ubyte omask = data[dofs].cull;
          if (!omask) continue;
          data[dofs].cull = 0x3f;
          // check
          const int z = cast(int)data[dofs].z;
          foreach (uint cidx; 0..6) {
            immutable ubyte cmask = cullmask(cidx);
            if (!(data[dofs].cull&cmask)) continue;
            const int nx = x+cullofs[cidx][0];
            const int ny = y+cullofs[cidx][1];
            const int nz = z+cullofs[cidx][2];
            if (bmp.getPixel(nx+1, ny+1, nz+1)) continue;
            // reset this cull bit
            data[dofs].cull ^= cmask;
          }
          changed += (omask != data[dofs].cull);
        }
      }
    }
    return changed;
  }

  void optimise (bool doHollowFill) @trusted {
    version(vox_check_invariants) checkInvariants();
    version(voxdata_debug) conwriteln("optimising mesh with ", voxpixtotal, " individual voxels...");
    if (doHollowFill) {
      version(voxdata_debug) conwriteln("optimising voxel culling...");
      uint count = hollowFill();
      if (count) conwriteln("fixed ", count, " voxel", (count != 1 ? "s" : ""));
      //count = fixFaceVisibility();
      //if (count) conwriteln("fixed ", count, " voxel", (count != 1 ? "s" : ""));
    } else {
      removeInsideFaces();
      version(voxdata_debug) conwriteln("optimising voxel culling...");
      uint count = fixFaceVisibility();
      if (count) conwriteln("fixed ", count, " voxel", (count != 1 ? "s" : ""));
    }
    removeEmptyVoxels();
    version(vox_check_invariants) checkInvariants();
    version(voxdata_debug) conwriteln("final optimised mesh contains ", voxpixtotal, " individual voxels...");
  }
}


// ////////////////////////////////////////////////////////////////////////// //
/*
  voxel data optimised for queries

  each slab is kept in this format:
    ushort zlo, zhi
    ushort runcount
  then run index follows:
    ushort z0
    ushort z1+1
    ushort ofs (from this z)
    ushort reserved
  index always ends with zhi+1, zhi+1
  then (b,g,r,cull) array
 */
struct VoxelDataSmall {
public:
  uint xsize = 0, ysize = 0, zsize = 0;
  float cx = 0.0f, cy = 0.0f, cz = 0.0f;

  ubyte[] data;
  // xsize*ysize array, offsets in `data`; 0 means "no data here"
  // slabs are sorted from bottom to top, and never intersects
  uint[] xyofs;

public:
  void save (VFile fl) {
    fl.rawWriteExact("K8VOXDC0");
    fl.writeNum!ushort(cast(ushort)xsize);
    fl.writeNum!ushort(cast(ushort)ysize);
    fl.writeNum!ushort(cast(ushort)zsize);
    fl.writeNum(cx);
    fl.writeNum(cy);
    fl.writeNum(cz);
    fl.rawWriteExact(xyofs[]);
    fl.writeNum!uint(cast(uint)data.length);
    fl.rawWriteExact(data[]);
  }


  void load (VFile fl) {
    char[8] sign = void;
    fl.rawReadExact(sign[]);
    if (sign[] != "K8VOXDC0") throw new Exception("invalid compressed voxel data");
    xsize = fl.readNum!ushort;
    ysize = fl.readNum!ushort;
    zsize = fl.readNum!ushort;
    cx = fl.readNum!float;
    cy = fl.readNum!float;
    cz = fl.readNum!float;
    delete xyofs;
    xyofs = new uint[xsize*ysize];
    fl.rawReadExact(xyofs[]);
    const uint dsz = fl.readNum!uint;
    delete data;
    data = new ubyte[dsz];
    fl.rawReadExact(data[]);
  }

private:
  void appendByte (ubyte v) @trusted {
    pragma(inline, true);
    data.assumeSafeAppend;
    data ~= v;
  }

  void appendShort (ushort v) @trusted {
    pragma(inline, true);
    data.assumeSafeAppend;
    data ~= cast(ubyte)v;
    data.assumeSafeAppend;
    data ~= cast(ubyte)(v>>8);
  }

private:
  private uint createSlab (ref VoxelData vox, uint dofs0) {
    while (dofs0 && !vox.data[dofs0].cull) dofs0 = vox.data[dofs0].nextz;
    if (!dofs0) return 0;
    // calculate zlo and zhi, and count runs
    ushort runcount = 0;
    ushort z0 = vox.data[dofs0].z;
    ushort z1 = z0;
    ushort nxz = cast(ushort)(z0-1);
    for (uint dofs = dofs0; dofs; dofs = vox.data[dofs].nextz) {
      if (!vox.data[dofs].cull) continue;
      z1 = vox.data[dofs].z;
      if (z1 != nxz) ++runcount;
      nxz = cast(ushort)(z1+1);
    }
    assert(runcount);
    if (data.length == 0) appendByte(0); // unused
    const uint startofs = cast(uint)data.length;
    // zlo
    appendShort(z0);
    // zhi
    appendShort(z1);
    // runcount
    appendShort(runcount);
    // run index (will be filled later)
    uint idxofs = cast(uint)data.length;
    for (uint f = 0; f < runcount; ++f) {
      appendShort(0); // z0
      appendShort(0); // z1
      appendShort(0); // offset
      appendShort(0); // reserved
    }
    // last index item
    appendShort(cast(ushort)(z1+1));
    appendShort(cast(ushort)(z1+1));
    appendShort(0); // offset
    appendShort(0); // reserved
    nxz = cast(ushort)(z0-1);
    ushort lastz = ushort.max;
    // put runs
    for (uint dofs = dofs0; dofs; dofs = vox.data[dofs].nextz) {
      if (!vox.data[dofs].cull) continue;
      z1 = vox.data[dofs].z;
      if (z1 != nxz) {
        // new run
        // fix prev run?
        if (lastz != ushort.max) {
          data[idxofs-6] = cast(ubyte)lastz;
          data[idxofs-5] = cast(ubyte)(lastz>>8);
        }
        // set run info
        // calculate offset
        const uint rofs = cast(uint)data.length-idxofs;
        assert(rofs <= 0xffffU);
        // z0
        data[idxofs++] = cast(ubyte)z1;
        data[idxofs++] = cast(ubyte)(z1>>8);
        // skip z1
        idxofs += 2;
        // offset
        data[idxofs++] = cast(ubyte)rofs;
        data[idxofs++] = cast(ubyte)(rofs>>8);
        // skip reserved
        idxofs += 2;
      }
      lastz = nxz = cast(ushort)(z1+1);
      // b, g, r, cull
      appendByte(vox.data[dofs].b);
      appendByte(vox.data[dofs].g);
      appendByte(vox.data[dofs].r);
      appendByte(vox.data[dofs].cull);
    }
    // fix prev run?
    assert(lastz != ushort.max);
    data[idxofs-6] = cast(ubyte)lastz;
    data[idxofs-5] = cast(ubyte)(lastz>>8);
    return startofs;
  }


  void checkValidity (ref VoxelData vox) {
    version(none) {
      foreach (uint y; 0..ysize) {
        foreach (uint x; 0..xsize) {
          for (uint dofs = vox.getDOfs(x, y); dofs; dofs = vox.data[dofs].nextz) {
            if (!vox.data[dofs].cull) continue;
            const uint vd = queryVox(x, y, vox.data[dofs].z);
            if (vd != vox.data[dofs].rgbcull()) assert(0);
          }
        }
      }
    } else {
      foreach (uint y; 0..ysize) {
        foreach (uint x; 0..xsize) {
          foreach (uint z; 0..zsize) {
            const uint vd = vox.query(x, y, z);
            if (vd != queryVox(x, y, z)) assert(0);
          }
        }
      }
    }
  }

public:
  void checkAccessTimes (ref VoxelData vox) {
    uint vvn = 0;
    auto tm = iv.timer.Timer(true);
    foreach (uint count; 0..3) {
      foreach (uint y; 0..ysize) {
        foreach (uint x; 0..xsize) {
          foreach (uint z; 0..zsize) {
            vvn += vox.query(x, y, z);
          }
        }
      }
    }
    tm.stop();
    conwriteln("linked list data query time: ", tm.toString(), " (", vvn, ")");

    vvn = 0;
    tm.restart();
    foreach (uint count; 0..3) {
      foreach (uint y; 0..ysize) {
        foreach (uint x; 0..xsize) {
          foreach (uint z; 0..zsize) {
            vvn += queryVox(x, y, z);
          }
        }
      }
    }
    tm.stop();
    conwriteln("sorted array data query time: ", tm.toString(), " (", vvn, ")");
  }

public:
  void clear () {
    delete data; data = null;
    delete xyofs; xyofs = null;
    xsize = ysize = zsize = 0;
    cx = cy = cz = 0.0f;
  }


  void createFrom (ref VoxelData vox) {
    clear();
    //vox.removeEmptyVoxels(); // just in case
    xsize = vox.xsize;
    ysize = vox.ysize;
    zsize = vox.zsize;
    xyofs.length = xsize*ysize;
    cx = vox.cx;
    cy = vox.cy;
    cz = vox.cz;
    foreach (uint y; 0..ysize) {
      foreach (uint x; 0..xsize) {
        const uint dofs = createSlab(vox, vox.getDOfs(x, y));
        xyofs[cast(uint)y*xsize+cast(uint)x] = dofs;
      }
    }
    conwriteln("*** created compressed voxel data; ", data.length, " bytes for ",
               xsize, "x", ysize, "x", zsize);
    checkValidity(vox);
  }


  uint queryVox (int x, int y, int z) nothrow @trusted @nogc {
    //pragma(inline, true);
    if (x < 0 || y < 0 || z < 0) return 0;
    if (cast(uint)x >= xsize || cast(uint)y >= ysize) return 0;
    uint dofs = xyofs.ptr[cast(uint)y*xsize+cast(uint)x];
    if (!dofs) return 0;
    const(ushort) *dptr = cast(const(ushort)*)(data.ptr+dofs);
      //conwriteln("z=", z, "; zlo=", dptr[0], "; zhi=", dptr[1], "; runcount=", dptr[2]);
    if (cast(ushort)z < *dptr++) return 0;
    if (cast(ushort)z > *dptr++) return 0;
    uint runcount = *dptr++;
    if (runcount <= 4) {
      // there is no reason to perform binary search here
      while (cast(ushort)z > *dptr) dptr += 4;
        //conwriteln("  rz0=", dptr[0], "; rz1=", dptr[1], "; rofs=", dptr[2]);
      if (z == *dptr) {
        const(uint) *dv = cast(const(uint)*)((cast(const(ubyte)*)dptr)+dptr[2]);
        return *dv;
      } else {
        dptr -= 4;
        const ushort cz = *dptr;
          //conwriteln("  cz=", cz, "; cz1=", dptr[1], "; rofs=", dptr[2]);
        assert(cz < z);
        if (cast(ushort)z >= dptr[1]) return 0; // no such voxel
        const(uint) *dv = cast(const(uint)*)((cast(const(ubyte)*)dptr)+dptr[2]);
        return *(dv+z-cz);
      }
    } else {
      //{ import core.stdc.stdio : printf; printf("runcount=%u\n", cast(uint)runcount); }
      //assert(0);
      // perform binary search
      uint lo = 0, hi = runcount-1;
      for (;;) {
        uint mid = (lo+hi)>>1;
        const(ushort) *dp = dptr+(mid<<2);
        if (cast(ushort)z >= dp[0] && cast(ushort)z < dp[1]) {
          const(uint) *dv = cast(const(uint)*)((cast(const(ubyte)*)dp)+dp[2]);
          return *(dv+z-*dp);
        }
        if (cast(ushort)z < dp[0]) {
          if (mid == lo) break;
          hi = mid-1;
        } else {
          if (mid == hi) { lo = hi; break; }
          lo = mid+1;
        }
      }
      const(ushort) *dp = dptr+(lo<<2);
      //{ import core.stdc.stdio : printf; printf("000: z=%u; lo=%u; cz0=%u; cz1=%u\n", cast(uint)z, cast(uint)lo, cast(uint)dp[0], cast(uint)dp[1]); }
      while (cast(ushort)z >= dp[1]) dp += 4;
      //{ import core.stdc.stdio : printf; printf("001:   lo=%u; cz0=%u; cz1=%u\n", cast(uint)z, cast(uint)lo, cast(uint)dp[0], cast(uint)dp[1]); }
      if (cast(ushort)z < dp[0]) return 0;
      const(uint) *dv = cast(const(uint)*)((cast(const(ubyte)*)dp)+dp[2]);
      return *(dv+z-*dp);
    }
  }
}