Printability checks on a voxel solid

// ---- squared Euclidean distance transform in voxel units -------------------------------------------------------------------
function edt1d(f, n, d, v, z) {                     // lower envelope of parabolas rooted at the samples f[q]
  let k = 0;
  v[0] = 0; z[0] = -INF; z[1] = INF;
  for (let q = 1; q < n; q++) {
    let s;
    for (;;) {
      const p = v[k];
      s = ((f[q] + q * q) - (f[p] + p * p)) / (2 * q - 2 * p);
      if (s <= z[k] && k > 0) k--; else break;
    }
    if (s <= z[k]) { v[0] = q; z[0] = -INF; z[1] = INF; k = 0; continue; }
    k++; v[k] = q; z[k] = s; z[k + 1] = INF;
  }
  k = 0;
  for (let q = 0; q < n; q++) {
    while (z[k + 1] < q) k++;
    const dq = q - v[k];
    d[q] = dq * dq + f[v[k]];
  }
}

/** Connected components (6-neighbours) of the voxels where mask[i] = 1 → { label: Int32Array (−1 outside), sizes: number[], sum: [[x,y,z sums]] }. */
export function components(mask, nx, ny, nz) {
  const n = mask.length, label = new Int32Array(n).fill(-1), sizes = [], sums = [], stack = new Int32Array(n);
  const nxy = nx * ny;
  for (let s = 0; s < n; s++) {
    if (!mask[s] || label[s] >= 0) continue;
    const id = sizes.length;
    let sp = 0, cnt = 0, sx = 0, sy = 0, sz = 0;
    stack[sp++] = s; label[s] = id;
    while (sp) {
      const c = stack[--sp], i = c % nx, j = Math.floor(c / nx) % ny, k = Math.floor(c / nxy);
      cnt++; sx += i; sy += j; sz += k;
      if (i > 0 && mask[c - 1] && label[c - 1] < 0) { label[c - 1] = id; stack[sp++] = c - 1; }
      if (i < nx - 1 && mask[c + 1] && label[c + 1] < 0) { label[c + 1] = id; stack[sp++] = c + 1; }
      if (j > 0 && mask[c - nx] && label[c - nx] < 0) { label[c - nx] = id; stack[sp++] = c - nx; }
      if (j < ny - 1 && mask[c + nx] && label[c + nx] < 0) { label[c + nx] = id; stack[sp++] = c + nx; }
      if (k > 0 && mask[c - nxy] && label[c - nxy] < 0) { label[c - nxy] = id; stack[sp++] = c - nxy; }
      if (k < nz - 1 && mask[c + nxy] && label[c + nxy] < 0) { label[c + nxy] = id; stack[sp++] = c + nxy; }
    }
    sizes.push(cnt); sums.push([sx, sy, sz]);
  }
  return { label, sizes, sums };
}

const sol = components(P.a, P.nx, P.ny, P.nz);
if (sol.sizes.length > 1) {
  const main = sol.sizes.indexOf(Math.max(...sol.sizes));
  sol.sizes.forEach((n, c) => { if (c !== main) out.islands.push({ volume: n * h ** 3, at: at(sol.sums[c].map((s) => s / n)) }); });
}
// ---- empty space: reached from the outside or sealed
const empty = new Uint8Array(N);
for (let i = 0; i < N; i++) empty[i] = P.a[i] ? 0 : 1;
const emp = components(empty, P.nx, P.ny, P.nz), outsideId = emp.label[0];
emp.sizes.forEach((n, c) => { if (c !== outsideId) out.sealed.push({ volume: n * h ** 3, at: at(emp.sums[c].map((s) => s / n)) }); });

const mm = (x) => +x.toFixed(3);

/**
 * Downward-facing surface that is steeper than the process prints without support, with the build direction +z. A face needs support when
 * its overhang angle from vertical (asin of the downward part of its normal) is above `angle`; the faces on the lowest plane lie on the bed and need none.
 * mesh = { positions, indices } in the printer's frame (z up). → { area, share (of the total area), at: [centroids of the 20 biggest], count }
 */
export function overhang(mesh, angle) {
  const P = mesh.positions, I = mesh.indices, lim = Math.sin(angle * Math.PI / 180);
  let zmin = Infinity;
  for (let i = 2; i < P.length; i += 3) if (P[i] < zmin) zmin = P[i];
  let area = 0, total = 0, count = 0;
  const hits = [];
  for (let t = 0; t < I.length; t += 3) {
    const a = I[t] * 3, b = I[t + 1] * 3, c = I[t + 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 nx = uy * vz - uz * vy, ny = uz * vx - ux * vz, nz = ux * vy - uy * vx, len = Math.sqrt(nx * nx + ny * ny + nz * nz);
    if (!(len > 0)) continue;
    const ar = len / 2;
    total += ar;
    if (-nz / len > lim + 1e-9 && Math.max(P[a + 2], P[b + 2], P[c + 2]) > zmin + 1e-6 * (1 + Math.abs(zmin))) {
      // a face right on the lowest plane is the bed face
      if (Math.abs(P[a + 2] - zmin) < 1e-6 && Math.abs(P[b + 2] - zmin) < 1e-6 && Math.abs(P[c + 2] - zmin) < 1e-6) continue;
      area += ar; count++;
      hits.push([ar, (P[a] + P[b] + P[c]) / 3, (P[a + 1] + P[b + 1] + P[c + 1]) / 3, (P[a + 2] + P[b + 2] + P[c + 2]) / 3]);
    }
  }
  hits.sort((p, q) => q[0] - p[0]);
  return { area, total, share: total > 0 ? area / total : 0, count, at: hits.slice(0, 20).map((h) => [h[1], h[2], h[3]]) };
}