Freeform mesh objects: quantities and placement on the ground

/** STL (binary or ASCII) → mesh, vertices merged within 1e-6 of a unit. `scale` = mm per file unit. */
export function readStl(bytes, { scale = 1 } = {}) {
  const u8 = bytes instanceof Uint8Array ? bytes : new Uint8Array(bytes);
  const pos = [], idx = [], seen = new Map();
  const add = (x, y, z) => { const k = `${Math.round(x * 1e6)},${Math.round(y * 1e6)},${Math.round(z * 1e6)}`; let i = seen.get(k); if (i === undefined) { i = pos.length / 3; seen.set(k, i); pos.push(x * scale, y * scale, z * scale); } idx.push(i); };
  const dv = new DataView(u8.buffer, u8.byteOffset, u8.byteLength);
  const nb = u8.length >= 84 ? dv.getUint32(80, true) : -1;
  if (nb >= 0 && 84 + nb * 50 === u8.length) {
    if (nb > MAX_TRIS) throw new Error('too-big');
    for (let t = 0; t < nb; t++) { const o = 84 + t * 50 + 12; for (let v = 0; v < 3; v++) add(dv.getFloat32(o + v * 12, true), dv.getFloat32(o + v * 12 + 4, true), dv.getFloat32(o + v * 12 + 8, true)); }
  } else {
    const text = new TextDecoder().decode(u8);
    const re = /vertex\s+([-+0-9.eE]+)\s+([-+0-9.eE]+)\s+([-+0-9.eE]+)/g;
    let m;
    while ((m = re.exec(text))) { add(+m[1], +m[2], +m[3]); if (idx.length / 3 > MAX_TRIS) throw new Error('too-big'); }
    if (idx.length % 3) idx.length -= idx.length % 3;
  }
  const mesh = checkMesh({ positions: pos, indices: idx });
  if (!mesh) throw new Error('not-mesh');
  return mesh;
}

/** Volume (divergence theorem: Σ v0·(v1×v2)/6, a closed outward mesh), surface area, plan footprint area, height — in m. */
export function meshQuantities(mesh) {
  const P = mesh.positions, I = mesh.indices;
  let vol = 0, area = 0, zmin = Infinity, zmax = -Infinity;
  const edges = new Map();
  for (let k = 0; k < I.length; k += 3) {
    const a = I[k] * 3, b = I[k + 1] * 3, c = I[k + 2] * 3;
    const ax = P[a], ay = P[a + 1], az = P[a + 2], bx = P[b], by = P[b + 1], bz = P[b + 2], cx = P[c], cy = P[c + 1], cz = P[c + 2];
    vol += (ax * (by * cz - bz * cy) - ay * (bx * cz - bz * cx) + az * (bx * cy - by * cx)) / 6;
    const ux = bx - ax, uy = by - ay, uz = bz - az, vx = cx - ax, vy = cy - ay, vz = cz - az;
    area += Math.hypot(uy * vz - uz * vy, uz * vx - ux * vz, ux * vy - uy * vx) / 2;
    for (const [p, q] of [[I[k], I[k + 1]], [I[k + 1], I[k + 2]], [I[k + 2], I[k]]]) { const key = p < q ? p + ',' + q : q + ',' + p; edges.set(key, (edges.get(key) || 0) + 1); }
  }
  for (let i = 2; i < P.length; i += 3) { zmin = Math.min(zmin, P[i]); zmax = Math.max(zmax, P[i]); }
  const closed = [...edges.values()].every((n) => n === 2);
  const hull = convexHull(planPoints(mesh));
  return { volume: closed ? Math.abs(vol) / 1e9 : null, area: area / 1e6, footprint: Math.abs(polyArea(hull)) / 1e6, height: (zmax - zmin) / 1000, zmin: zmin / 1000, closed };
}

export function convexHull(pts) {
  const p = [...pts].sort((a, b) => a.x - b.x || a.y - b.y), cross = (o, a, b) => (a.x - o.x) * (b.y - o.y) - (a.y - o.y) * (b.x - o.x);
  if (p.length < 3) return p;
  const lo = [], hi = [];
  for (const q of p) { while (lo.length >= 2 && cross(lo[lo.length - 2], lo[lo.length - 1], q) <= 0) lo.pop(); lo.push(q); }
  for (let i = p.length - 1; i >= 0; i--) { const q = p[i]; while (hi.length >= 2 && cross(hi[hi.length - 2], hi[hi.length - 1], q) <= 0) hi.pop(); hi.push(q); }
  return [...lo.slice(0, -1), ...hi.slice(0, -1)];
}

/** mm (mesh) → project coordinates, turned by at.rot about the insertion point. */
function toWorld(rec, units) {
  const k = 1 / (1000 * unitM(units)), a = (rec.at.rot * Math.PI) / 180, c = Math.cos(a), s = Math.sin(a);
  return (x, y) => ({ x: rec.at.x + (x * c - y * s) * k, y: rec.at.y + (x * s + y * c) * k });
}

/** Where it sits: follow = the lowest ground under the footprint (corners, edge midpoints, centre); cut = baseZ (or the mean). */
export function placed(rec, surface, units) {
  const fp = footprint(rec, units), um = unitM(units);
  const zs = [];
  if (surface) {
    for (let i = 0; i < fp.length; i++) { const p = fp[i], q = fp[(i + 1) % fp.length]; for (const r of [p, { x: (p.x + q.x) / 2, y: (p.y + q.y) / 2 }]) { const z = surface.elevationAt(r.x, r.y, 1e-6); if (z !== null) zs.push(z); } }
    const cz = surface.elevationAt(rec.at.x, rec.at.y, 1e-6); if (cz !== null) zs.push(cz);
  }
  const ground = zs.length ? { min: Math.min(...zs), mean: zs.reduce((a, b) => a + b, 0) / zs.length } : { min: 0, mean: 0 };
  const zmin = Math.min(...rec.mesh.positions.filter((_, i) => i % 3 === 2)) / 1000 / um;
  const base = rec.base === 'cut' ? rec.baseZ ?? ground.mean : ground.min;
  const w = toWorld(rec, units), P = rec.mesh.positions, positions = [];
  for (let i = 0; i < P.length; i += 3) { const q = w(P[i], P[i + 1]); positions.push(q.x, q.y, base + P[i + 2] / 1000 / um - zmin); }
  // below the existing ground (follow): sampled on a grid over the footprint — an estimate, as the spike says
  let below = 0;
  if (surface && rec.base === 'follow') {
    let x0 = Infinity, y0 = Infinity, x1 = -Infinity, y1 = -Infinity;
    for (const p of fp) { x0 = Math.min(x0, p.x); y0 = Math.min(y0, p.y); x1 = Math.max(x1, p.x); y1 = Math.max(y1, p.y); }
    const cell = Math.max((x1 - x0) / 40, (y1 - y0) / 40, 1e-6);
    for (let x = x0 + cell / 2; x < x1; x += cell) for (let y = y0 + cell / 2; y < y1; y += cell) { if (!inside(x, y, fp)) continue; const g = surface.elevationAt(x, y, 1e-6); if (g !== null && g > base) below += (g - base) * cell * cell; }
    below *= um ** 3;
  }
  return { base, positions, footprint: fp, belowGround: below };
}