Generated part exports: data sheet, VTK, 3MF and OpenVDB

/** The field on the node lattice of `grid` over `box` [x0,y0,z0,x1,y1,z1] (default: every stored brick), distance clamped to ±band. */
export function toVTI(grid, { box = null, maxNodes = 6e7 } = {}) {
  const h = grid.h;
  let lo, hi;
  if (box) { lo = [0, 1, 2].map((a) => Math.floor(box[a] / h)); hi = [0, 1, 2].map((a) => Math.ceil(box[a + 3] / h)); }
  else { const e = grid.extent(); if (!e) throw Object.assign(new Error('export: the field is empty'), { code: 'empty' }); lo = e.lo.map((b) => b * 8); hi = e.hi.map((b) => b * 8 + 7); }
  const n = [hi[0] - lo[0] + 1, hi[1] - lo[1] + 1, hi[2] - lo[2] + 1], total = n[0] * n[1] * n[2];
  if (total > maxNodes) throw Object.assign(new Error('export: too many nodes'), { code: 'toobig', vars: { n: total } });
  const dist = new Float32Array(total), mask = new Uint8Array(total);
  for (let k = 0; k < n[2]; k++) for (let j = 0; j < n[1]; j++) for (let i = 0; i < n[0]; i++) {
    const v = grid.at(lo[0] + i, lo[1] + j, lo[2] + k), c = i + n[0] * (j + n[1] * k);
    dist[c] = v; mask[c] = v < 0 ? 1 : 0;
  }
  const ext = `0 ${n[0] - 1} 0 ${n[1] - 1} 0 ${n[2] - 1}`;
  return `<?xml version="1.0"?>\n<VTKFile type="ImageData" version="1.0" byte_order="LittleEndian" header_type="UInt32">\n` +
    `  <ImageData WholeExtent="${ext}" Origin="${num(lo[0] * h)} ${num(lo[1] * h)} ${num(lo[2] * h)}" Spacing="${num(h)} ${num(h)} ${num(h)}">\n    <Piece Extent="${ext}">\n` +
    `      <PointData Scalars="distance">\n        <DataArray type="Float32" Name="distance" format="binary">${blob(dist)}</DataArray>\n        <DataArray type="UInt8" Name="inside" format="binary">${blob(mask)}</DataArray>\n      </PointData>\n      <CellData/>\n    </Piece>\n  </ImageData>\n</VTKFile>\n`;
}

/** Surface triangles as an UnstructuredGrid: the points, triangle cells (VTK type 5), and per cell the area and the unit normal. */
export function toVTU({ positions, indices }, { name = 'surface' } = {}) {
  const nt = indices.length / 3, nv = positions.length / 3;
  const P = Float32Array.from(positions), conn = Int32Array.from(indices), off = new Int32Array(nt), types = new Uint8Array(nt).fill(5), area = new Float32Array(nt), nrm = new Float32Array(nt * 3);
  for (let t = 0; t < nt; t++) {
    off[t] = 3 * (t + 1);
    const a = indices[t * 3] * 3, b = indices[t * 3 + 1] * 3, c = indices[t * 3 + 2] * 3;
    const ux = P[b] - P[a], uy = P[b + 1] - P[a + 1], uz = P[b + 2] - P[a + 2], vx = P[c] - P[a], vy = P[c + 1] - P[a + 1], vz = P[c + 2] - P[a + 2];
    const x = uy * vz - uz * vy, y = uz * vx - ux * vz, z = ux * vy - uy * vx, l = Math.sqrt(x * x + y * y + z * z);
    area[t] = l / 2; if (l > 0) { nrm[t * 3] = x / l; nrm[t * 3 + 1] = y / l; nrm[t * 3 + 2] = z / l; }
  }
  return `<?xml version="1.0"?>\n<VTKFile type="UnstructuredGrid" version="1.0" byte_order="LittleEndian" header_type="UInt32">\n  <UnstructuredGrid>\n` +
    `    <Piece NumberOfPoints="${nv}" NumberOfCells="${nt}">\n      <Points>\n        <DataArray type="Float32" Name="${esc(name)}" NumberOfComponents="3" format="binary">${blob(P)}</DataArray>\n      </Points>\n` +
    `      <Cells>\n        <DataArray type="Int32" Name="connectivity" format="binary">${blob(conn)}</DataArray>\n        <DataArray type="Int32" Name="offsets" format="binary">${blob(off)}</DataArray>\n        <DataArray type="UInt8" Name="types" format="binary">${blob(types)}</DataArray>\n      </Cells>\n` +
    `      <CellData Normals="normal">\n        <DataArray type="Float32" Name="area" format="binary">${blob(area)}</DataArray>\n        <DataArray type="Float32" Name="normal" NumberOfComponents="3" format="binary">${blob(nrm)}</DataArray>\n      </CellData>\n    </Piece>\n  </UnstructuredGrid>\n</VTKFile>\n`;
}

/** The tree: root → 32³ internal (level 2) → 16³ internal (level 1) → 8³ leaves (our bricks), every list in OpenVDB's order. */
function buildTree(grid) {
  const bg = Math.fround(grid.bw);
  const l2 = new Map();                                  // key "x,y,z" of the level-2 origin → { o:[x,y,z], kids: Map(offset → level-1) }
  let voxels = 0;
  const lo = [Infinity, Infinity, Infinity], hi = [-Infinity, -Infinity, -Infinity];
  for (const key of grid.sortedKeys()) {
    const b = grid.bricks.get(key), [bx, by, bz] = brickOf(key);
    const o2 = [bx >> 9 << 12, by >> 9 << 12, bz >> 9 << 12];                     // 512 bricks = 4096 voxels per level-2 node
    const k2 = o2.join(',');
    let n2 = l2.get(k2); if (!n2) { n2 = { o: o2, kids: new Map() }; l2.set(k2, n2); }
    const off2 = ((bx >> 4) & 31) << 10 | ((by >> 4) & 31) << 5 | ((bz >> 4) & 31);
    let n1 = n2.kids.get(off2); if (!n1) { n1 = { leaves: new Map(), inside: new Set() }; n2.kids.set(off2, n1); }
    const off1 = (bx & 15) << 8 | (by & 15) << 4 | (bz & 15);
    if (b === INSIDE) { n1.inside.add(off1); continue; }
    // reorder x-fastest → z-fastest, and mark the narrow band active
    const vals = new Float32Array(BN), mask = newMask(BN);
    for (let z = 0; z < 8; z++) for (let y = 0; y < 8; y++) for (let x = 0; x < 8; x++) {
      const v = b[x | y << 3 | z << 6], n = x << 6 | y << 3 | z;
      vals[n] = v;
      if (v < bg && v > -bg) {
        setBit(mask, n); voxels++;
        const p = [bx * 8 + x, by * 8 + y, bz * 8 + z];
        for (let a = 0; a < 3; a++) { if (p[a] < lo[a]) lo[a] = p[a]; if (p[a] > hi[a]) hi[a] = p[a]; }
      }
    }
    n1.leaves.set(off1, { o: [bx * 8, by * 8, bz * 8], vals, mask });
  }
  // the root's table is a std::map on Coord: x, then y, then z
  const roots = [...l2.values()].sort((a, b) => a.o[0] - b.o[0] || a.o[1] - b.o[1] || a.o[2] - b.o[2]);
  return { bg, roots, voxels, bbox: voxels ? { lo, hi } : null };
}

/**
 * io::writeCompressedValues with COMPRESS_ACTIVE_MASK, as Compression.h decides it: look at up to three distinct inactive values
 * (skipping child slots), pick the metadata byte, then write the inactive values it needs, the selection mask when there are two,
 * and the active values (or every value when there are more than two inactive ones). Float32 comparisons are exact, as there.
 */
function writeValues(out, vals, count, valueMask, childMask, bg) {
  const f = Math.fround, bgf = f(bg), mbg = f(-bg);
  const inactive = [bgf, bgf];
  let unique = 0;
  for (let i = 0; i < count && unique < 3; i++) {
    if (isOn(valueMask, i) || (childMask && isOn(childMask, i))) continue;
    const v = f(vals[i]);
    const seen = (unique > 0 && v === inactive[0]) || (unique > 1 && v === inactive[1]);
    if (!seen) { if (unique < 2) inactive[unique] = v; unique++; }
  }
  let meta = NO_MASK_OR_INACTIVE_VALS;
  if (unique === 1) {
    if (inactive[0] !== bgf) meta = inactive[0] === mbg ? NO_MASK_AND_MINUS_BG : NO_MASK_AND_ONE_INACTIVE_VAL;
  } else if (unique === 2) {
    if (inactive[0] !== bgf && inactive[1] !== bgf) meta = MASK_AND_TWO_INACTIVE_VALS;
    else if (inactive[1] === bgf) meta = inactive[0] === mbg ? MASK_AND_NO_INACTIVE_VALS : MASK_AND_ONE_INACTIVE_VAL;
    else if (inactive[0] === bgf) { meta = inactive[1] === mbg ? MASK_AND_NO_INACTIVE_VALS : MASK_AND_ONE_INACTIVE_VAL; const s = inactive[0]; inactive[0] = inactive[1]; inactive[1] = s; }
  } else if (unique > 2) meta = NO_MASK_AND_ALL_VALS;
  out.u8(meta);
  if (meta === NO_MASK_AND_ONE_INACTIVE_VAL || meta === MASK_AND_ONE_INACTIVE_VAL || meta === MASK_AND_TWO_INACTIVE_VALS) {
    out.f32(inactive[0]);
    if (meta === MASK_AND_TWO_INACTIVE_VALS) out.f32(inactive[1]);
  }
  if (meta === NO_MASK_AND_ALL_VALS) { out.f32s(vals, count); return; }
  const active = [];
  if (meta === NO_MASK_OR_INACTIVE_VALS || meta === NO_MASK_AND_MINUS_BG || meta === NO_MASK_AND_ONE_INACTIVE_VAL) {
    for (let i = 0; i < count; i++) if (isOn(valueMask, i)) active.push(vals[i]);
  } else {
    const sel = newMask(count);
    for (let i = 0; i < count; i++) {
      if (isOn(valueMask, i)) active.push(vals[i]);
      else if (f(vals[i]) === inactive[1]) setBit(sel, i);
    }
    out.mask(sel);
  }
  out.f32s(active, active.length);
}