Surface display meshes, curvature map and zebra stripes

/** Display meshes of the surface bodies left after a rebuild: [{ feature, positions, indices }]. */
export function sheetMeshes(ctx) {
  const out = [];
  if (!ctx.sheets) return out;
  const { oc, T } = ctx;
  for (const [id, sh] of ctx.sheets) {
    try {
      const bb = ctx.bbox(sh.shape);
      const diag = Math.max(1e-6, V.len(V.sub(bb.max, bb.min)));
      const md = T(oc.ReplicadMeshExtractor.extract(sh.shape, diag * 2e-4, 2 * Math.PI / 64, false));
      const mem = oc.wasmMemory.buffer;
      const positions = new Float32Array(mem, md.getVerticesPtr(), md.getVerticesSize()).slice();
      const indices = Uint32Array.from(new Int32Array(mem, md.getTrianglesPtr(), md.getTrianglesSize()));
      oc.BRepTools.Clean(sh.shape, true);
      let area = 0;
      for (const f of explore(oc, T, sh.shape, 'FACE')) { const g = T(new oc.GProp_GProps()); oc.BRepGProp.SurfaceProperties(f, g, false, false); area += g.Mass(); }
      out.push({ feature: id, positions, indices, area, faces: [...sh.names.entries()].length });
    } catch { /* a sheet that cannot be meshed is left out of the view */ }
  }
  return out;
}

/**
 * Discrete curvature per vertex of a closed or open triangle mesh { positions, indices }:
 *   gauss = (2π − Σ angles) / A_mixed   (π − Σ on boundary vertices, so open edges read near 0)
 *   mean  = |Σ (cot α + cot β)(p_i − p_j)| / (4 A_mixed), signed by the vertex normal
 * (Meyer, Desbrun, Schröder, Barr 2003, with the barycentric area as A_mixed). On a sphere of radius r:
 * gauss → 1/r², mean → 1/r as the mesh is refined.
 */
export function curvature(mesh) {
  const P = mesh.positions, I = mesh.indices;
  const n = P.length / 3;
  const angle = new Float64Array(n), area = new Float64Array(n);
  const lap = new Float64Array(n * 3), nrm = new Float64Array(n * 3);
  const p = (i) => [P[i * 3], P[i * 3 + 1], P[i * 3 + 2]];
  const edgeCount = new Map();
  for (let t = 0; t < I.length; t += 3) {
    const ids = [I[t], I[t + 1], I[t + 2]];
    const q = ids.map(p);
    const cr = V.cross(V.sub(q[1], q[0]), V.sub(q[2], q[0]));
    const A = V.len(cr) / 2;
    if (!(A > 0)) continue;
    for (let k = 0; k < 3; k++) {
      const i = ids[k], j = ids[(k + 1) % 3], l = ids[(k + 2) % 3];
      const a = V.sub(q[(k + 1) % 3], q[k]), b = V.sub(q[(k + 2) % 3], q[k]);
      angle[i] += Math.atan2(V.len(V.cross(a, b)), V.dot(a, b));
      area[i] += A / 3;
      for (let c = 0; c < 3; c++) nrm[i * 3 + c] += cr[c];
      // cot of the angle at l, opposite edge (i, j)
      const u = V.sub(q[k], q[(k + 2) % 3]), w = V.sub(q[(k + 1) % 3], q[(k + 2) % 3]);
      const cot = V.dot(u, w) / V.len(V.cross(u, w));
      const d = V.sub(q[k], q[(k + 1) % 3]);
      for (let c = 0; c < 3; c++) { lap[i * 3 + c] += cot * d[c]; lap[j * 3 + c] -= cot * d[c]; }
      const key = i < j ? i + ',' + j : j + ',' + i;
      edgeCount.set(key, (edgeCount.get(key) || 0) + 1);
      void l;
    }
  }
  const boundary = new Uint8Array(n);
  for (const [k, c] of edgeCount) if (c === 1) for (const v of k.split(',')) boundary[+v] = 1;
  const gauss = new Float64Array(n), mean = new Float64Array(n);
  for (let i = 0; i < n; i++) {
    if (!(area[i] > 0)) continue;
    gauss[i] = ((boundary[i] ? Math.PI : 2 * Math.PI) - angle[i]) / area[i];
    const L = [lap[i * 3], lap[i * 3 + 1], lap[i * 3 + 2]];
    const N = [nrm[i * 3], nrm[i * 3 + 1], nrm[i * 3 + 2]];
    mean[i] = boundary[i] ? 0 : (V.dot(L, N) >= 0 ? 1 : -1) * V.len(L) / (4 * area[i]);
  }
  return { gauss, mean, boundary };
}

/** Zebra stripe value (0 or 1) of a unit normal seen along `view`: stripes of a vertical reflected light pattern. */
export function zebra(n, view = [0, 0, -1], stripes = 12) {
  const d = V.norm(view);
  const r = V.sub(d, V.mul(n, 2 * V.dot(d, n)));   // reflected view ray
  const a = Math.atan2(r[0], r[2]);
  return Math.floor((a / Math.PI + 1) * stripes / 2) % 2;
}