Voxel grid and lattice struts

export function toVoxelGrid(grid, { bbox = null } = {}) {
  const h = grid.h;
  let lo, hi;
  if (bbox) { lo = [0, 1, 2].map((a) => Math.floor(bbox[a] / h)); hi = [0, 1, 2].map((a) => Math.ceil(bbox[a + 3] / h)); }
  else {
    const e = grid.extent();
    if (!e) return { h, nx: 0, ny: 0, nz: 0, origin: [0, 0, 0], solid: new Uint8Array(0), count: 0, oddRows: 0 };
    lo = e.lo.map((b) => b * 8); hi = e.hi.map((b) => b * 8 + 8);
  }
  const nx = hi[0] - lo[0], ny = hi[1] - lo[1], nz = hi[2] - lo[2];
  const solid = new Uint8Array(nx * ny * nz);
  let count = 0;
  const bw = grid.bw;
  for (const key of grid.sortedKeys()) {
    const [bx, by, bz] = brickOf(key), b = grid.bricks.get(key);
    for (let k = 0; k < 8; k++) for (let j = 0; j < 8; j++) for (let i = 0; i < 8; i++) {
      const I = bx * 8 + i, J = by * 8 + j, K = bz * 8 + k;
      if (I < lo[0] || I >= hi[0] || J < lo[1] || J >= hi[1] || K < lo[2] || K >= hi[2]) continue;
      let s = 0;
      if (b === INSIDE) s = -8 * bw;
      else for (let c = 0; c < 8; c++) s += grid.at(I + (c & 1), J + (c >> 1 & 1), K + (c >> 2 & 1));
      if (s < 0) { solid[(I - lo[0]) + nx * ((J - lo[1]) + ny * (K - lo[2]))] = 1; count++; }
    }
  }
  // cubes whose own brick is absent but whose corners reach into a stored one are handled by the stored brick's own loop only
  // when they start inside it; the rest are outside by construction (every node of an absent brick is outside the band).
  return { h, nx, ny, nz, origin: [lo[0] * h, lo[1] * h, lo[2] * h], solid, count, oddRows: 0 };
}

/**
 * Struts of a lattice filling the box [x0,y0,z0,x1,y1,z1] with cubic cells of side `cell` (whole cells, starting at the box's
 * low corner; the last row may stick out past the box by less than a cell). Shared struts are listed once. Radius r.
 * → Float64Array(8·n), row-major beams.
 */
export function latticeBeams(kind, box, cell, r) {
  if (!unit[kind]) throw new FieldError('lattice_kind', 'field: unknown lattice ' + kind + ' (' + LATTICE_KINDS.join(', ') + ')');
  if (!(cell > 0)) throw new FieldError('cell', 'field: lattice cell must be positive');
  const nx = Math.max(1, Math.ceil((box[3] - box[0]) / cell - 1e-9)), ny = Math.max(1, Math.ceil((box[4] - box[1]) / cell - 1e-9)), nz = Math.max(1, Math.ceil((box[5] - box[2]) / cell - 1e-9));
  if (nx * ny * nz > 4e6) throw new FieldError('cells', 'field: lattice would have more than 4 million cells');
  const cellBeams = unit[kind](), seen = new Set(), out = [];
  const q = (v) => Math.round(v * 1e6);
  for (let i = 0; i < nx; i++) for (let j = 0; j < ny; j++) for (let k = 0; k < nz; k++) {
    for (const [a, b] of cellBeams) {
      const A = [q(i + a[0]), q(j + a[1]), q(k + a[2])], B = [q(i + b[0]), q(j + b[1]), q(k + b[2])];
      const lo = A[0] < B[0] || (A[0] === B[0] && (A[1] < B[1] || (A[1] === B[1] && A[2] <= B[2]))) ? [A, B] : [B, A];
      const key = lo[0].join(',') + '>' + lo[1].join(',');
      if (seen.has(key)) continue;
      seen.add(key);
      out.push(box[0] + (i + a[0]) * cell, box[1] + (j + a[1]) * cell, box[2] + (k + a[2]) * cell, box[0] + (i + b[0]) * cell, box[1] + (j + b[1]) * cell, box[2] + (k + b[2]) * cell, r, r);
    }
  }
  return Float64Array.from(out);
}