Turning a face mesh back into an extrusion

/** A closed triangle mesh (flat xyz positions, flat indices) that is exactly a straight prism: two parallel caps and sides along their normal.
 *  → { origin, x, y, n, depth, outer: [[u, v]…] (counter-clockwise), holes: [[[u, v]…]] } or null for anything else (it is then left alone). */
export function meshToPrism(positions, indices, { tol = 1e-6 } = {}) {
  const nv = Math.floor(positions.length / 3), nt = Math.floor(indices.length / 3);
  if (nv < 6 || nt < 8 || nt > MAX_TRIANGLES) return null;
  const V = (i) => [positions[3 * i], positions[3 * i + 1], positions[3 * i + 2]];
  let lo = [Infinity, Infinity, Infinity], hi = [-Infinity, -Infinity, -Infinity];
  for (let i = 0; i < nv; i++) for (let k = 0; k < 3; k++) { const v = positions[3 * i + k]; if (!Number.isFinite(v)) return null; if (v < lo[k]) lo[k] = v; if (v > hi[k]) hi[k] = v; }
  const diag = Math.hypot(hi[0] - lo[0], hi[1] - lo[1], hi[2] - lo[2]);
  if (!(diag > 0)) return null;
  const eps = tol * diag;
  // closed and consistent: every undirected edge in exactly two triangles, once each way
  const dir = new Map();
  const tri = [];
  const buckets = new Map();
  for (let t = 0; t < nt; t++) {
    const i = indices[3 * t], j = indices[3 * t + 1], k = indices[3 * t + 2];
    if (!(i >= 0 && i < nv && j >= 0 && j < nv && k >= 0 && k < nv)) return null;
    tri.push([i, j, k]);
    for (const [a, b] of [[i, j], [j, k], [k, i]]) { const key = a * nv + b; if (dir.has(key)) return null; dir.set(key, 1); }
    const c = cross(sub(V(j), V(i)), sub(V(k), V(i))), a = len(c);
    if (a > 0) {
      const u = [c[0] / a, c[1] / a, c[2] / a];
      const key = u.map((q) => Math.round(q * 1e4) || 0).join(',');   // one bucket per direction
      const b = buckets.get(key) || { n: u, area: 0 };
      b.area += a; buckets.set(key, b);
    }
  }
  for (const key of dir.keys()) { const a = Math.floor(key / nv), b = key % nv; if (!dir.has(b * nv + a)) return null; }
  // the axis is the normal of the caps: try the directions with the most area first (a prism's sides can outweigh its caps)
  const cands = [...buckets.values()].sort((p, q) => q.area - p.area).slice(0, 8);
  for (const cd of cands) {
    const r = prismAlong(positions, indices, cd.n, { nv, tri, V, eps, diag });
    if (r) return r;
  }
  return null;
}

function prismAlong(positions, indices, n, { nv, tri, V, eps }) {
  // layers along n
  const h = new Float64Array(nv);
  let h0 = Infinity, h1 = -Infinity;
  for (let i = 0; i < nv; i++) { h[i] = dot(V(i), n); if (h[i] < h0) h0 = h[i]; if (h[i] > h1) h1 = h[i]; }
  const depth = h1 - h0;
  if (!(depth > eps)) return null;
  const layer = new Int8Array(nv);
  for (let i = 0; i < nv; i++) { if (Math.abs(h[i] - h0) <= eps) layer[i] = 0; else if (Math.abs(h[i] - h1) <= eps) layer[i] = 1; else return null; }
  // frame: x ⟂ n from the world axis least along n; y = n × x
  const ax = Math.abs(n[0]) <= Math.abs(n[1]) && Math.abs(n[0]) <= Math.abs(n[2]) ? [1, 0, 0] : Math.abs(n[1]) <= Math.abs(n[2]) ? [0, 1, 0] : [0, 0, 1];
  let x = cross(ax, n); const lx = len(x); x = [x[0] / lx, x[1] / lx, x[2] / lx];
  const y = cross(n, x);
  const origin = [n[0] * h0, n[1] * h0, n[2] * h0];
  const uv = (i) => { const d = sub(V(i), origin); return [dot(d, x), dot(d, y)]; };
  // caps and sides
  const bottom = [], top = [];
  for (const [i, j, k] of tri) {
    const L = layer[i] + layer[j] + layer[k];
    const c = cross(sub(V(j), V(i)), sub(V(k), V(i))), a = len(c);
    if (a <= 0) continue;
    const nd = dot(c, n) / a;
    if (L === 0) { if (nd > -1 + 1e-6) return null; bottom.push([i, j, k]); }
    else if (L === 3) { if (nd < 1 - 1e-6) return null; top.push([i, j, k]); }
    else if (Math.abs(nd) > 1e-6) return null;
  }
  if (!bottom.length || !top.length) return null;
  const boundary = (T) => {
    const cnt = new Map();
    for (const [i, j, k] of T) for (const [a, b] of [[i, j], [j, k], [k, i]]) cnt.set(a * nv + b, true);
    const edges = [];
    for (const key of cnt.keys()) { const a = Math.floor(key / nv), b = key % nv; if (!cnt.has(b * nv + a)) edges.push([a, b]); }
    return loops(edges);
  };
  const lb = boundary(bottom), lt = boundary(top);
  if (!lb || !lt || !lb.length || lb.length !== lt.length) return null;
  const pts = (l) => l.map(uv);
  const key = (p) => Math.round(p[0] / eps) + ',' + Math.round(p[1] / eps);
  // the top outline, seen from below, must be the bottom outline (same set of edges)
  const edgeSet = (ls) => { const s = new Set(); for (const l of ls) { const p = pts(l); for (let i = 0; i < p.length; i++) { const a = key(p[i]), b = key(p[(i + 1) % p.length]); s.add(a < b ? a + '|' + b : b + '|' + a); } } return s; };
  const eb = edgeSet(lb), et = edgeSet(lt);
  if (eb.size !== et.size) return null;
  for (const e of eb) if (!et.has(e)) return null;
  // outer = the loop of largest area; the rest must be inside it
  let ps = lb.map(pts);
  ps = ps.map((p) => (area2(p) < 0 ? p.slice().reverse() : p)).map((p) => simplify(p, eps * eps));
  ps.sort((a, b) => area2(b) - area2(a));
  const outer = ps[0], holes = ps.slice(1);
  if (outer.length < 3 || holes.some((p) => p.length < 3 || !inside(outer, p[0]))) return null;
  return { origin, x, y, n, depth, outer, holes: holes.map((p) => p.slice().reverse()) };   // holes run clockwise
}

function simplify(pts, tol) {
  let out = pts.slice();
  for (let changed = true; changed && out.length > 3;) {
    changed = false;
    for (let i = 0; i < out.length && out.length > 3; i++) {
      const a = out[(i + out.length - 1) % out.length], b = out[i], c = out[(i + 1) % out.length];
      const cr = (b[0] - a[0]) * (c[1] - b[1]) - (b[1] - a[1]) * (c[0] - b[0]);
      if (Math.abs(cr) <= tol * Math.max(1, Math.hypot(c[0] - a[0], c[1] - a[1]))) { out.splice(i, 1); changed = true; i--; }
    }
  }
  return out;
}

/** Volume of a closed triangle mesh (divergence theorem) — for tests and checks. */
export function meshVolume(positions, indices) {
  let v = 0;
  for (let t = 0; t + 2 < indices.length; t += 3) {
    const a = indices[t] * 3, b = indices[t + 1] * 3, c = indices[t + 2] * 3;
    const A = [positions[a], positions[a + 1], positions[a + 2]], B = [positions[b], positions[b + 1], positions[b + 2]], C = [positions[c], positions[c + 1], positions[c + 2]];
    v += dot(A, cross(B, C)) / 6;
  }
  return v;
}

/** Area of the prism's profile (outer minus holes) — depth × this is its volume. */
export const prismArea = (p) => Math.abs(area2(p.outer)) - p.holes.reduce((s, h) => s + Math.abs(area2(h)), 0);

const real = (v) => { const s = String(+Number(v).toFixed(6)); return /[.eE]/.test(s) ? s : s + '.'; };

/** Write a prism as IFC lines through `add(text) → '#id'` (as step-ifc.js's writer has). Returns the IfcExtrudedAreaSolid reference. */
export function prismLines(add, prism) {
  const pt2 = (p) => add('IFCCARTESIANPOINT((' + real(p[0]) + ',' + real(p[1]) + '))');
  const closed = (pts) => { const ids = pts.map(pt2); return add('IFCPOLYLINE((' + [...ids, ids[0]].join(',') + '))'); };
  const outer = closed(prism.outer);
  const prof = prism.holes.length
    ? add("IFCARBITRARYPROFILEDEFWITHVOIDS(.AREA.,$," + outer + ',(' + prism.holes.map(closed).join(',') + '))')
    : add('IFCARBITRARYCLOSEDPROFILEDEF(.AREA.,$,' + outer + ')');
  const o = add('IFCCARTESIANPOINT((' + prism.origin.map(real).join(',') + '))');
  const zd = add('IFCDIRECTION((' + prism.n.map(real).join(',') + '))');
  const xd = add('IFCDIRECTION((' + prism.x.map(real).join(',') + '))');
  const pos = add('IFCAXIS2PLACEMENT3D(' + o + ',' + zd + ',' + xd + ')');
  const up = add('IFCDIRECTION((0.,0.,1.))');
  return add('IFCEXTRUDEDAREASOLID(' + prof + ',' + pos + ',' + up + ',' + real(prism.depth) + ')');
}