Smooth shading and texture placement

/**
 * Smooth shading (0.4 M18): faces joined by soft edges are one smooth surface. → Map `${fid}|${v}` → the unit normal
 * (definition coordinates) at that corner: the area-weighted average over the faces of its surface that meet at v.
 * Faces with no soft edge are left out (they shade flat).
 */
export function smoothNormals(def) {
  const out = new Map();
  if (![...def.edges.values()].some((e) => e.soft)) return out;
  const ef = edgeFaces(def), up = new Map();
  const find = (x) => { while (up.get(x) !== x) { up.set(x, up.get(up.get(x))); x = up.get(x); } return x; };
  for (const fid of def.faces.keys()) up.set(fid, fid);
  const joined = new Set();
  for (const e of def.edges.values()) {
    if (!e.soft) continue;
    const fs = ef.get(key(e.a, e.b)) || [];
    for (let i = 1; i < fs.length; i++) { up.set(find(fs[i]), find(fs[0])); joined.add(fs[0]); joined.add(fs[i]); }
  }
  const acc = new Map();
  for (const fid of joined) {
    const f = def.faces.get(fid), n = faceNormal(def, f), a = polyArea3(f.loops[0].map((v) => def.verts.get(v))) || 1e-12, r = find(fid);
    for (const loop of f.loops) for (const v of loop) { const k = r + '|' + v; acc.set(k, add(acc.get(k) || [0, 0, 0], mul(n, a))); }
  }
  for (const fid of joined) {
    const f = def.faces.get(fid), r = find(fid), fn = faceNormal(def, f);
    for (const loop of f.loops) for (const v of loop) {
      const s = acc.get(r + '|' + v), L = len(s);
      out.set(fid + '|' + v, L > 1e-12 ? mul(s, 1 / L) : fn);
    }
  }
  return out;
}

/**
 * The image paint that puts texture coordinates `uvs` at the points `pts` of a flat face (normal n): the affine map
 * u = U·p + u0 with U·n = 0, solved from the first three points that make a triangle. → paint | null (degenerate).
 */
export function imagePaintFit(id, pts, uvs, n, avg = null) {
  for (let i = 0; i + 2 < pts.length; i++) for (let j = i + 1; j + 1 < pts.length; j++) for (let k = j + 1; k < pts.length; k++) {
    const A = [[...pts[i], 1], [...pts[j], 1], [...pts[k], 1], [n[0], n[1], n[2], 0]];
    const solve = (b) => {   // 4×4 Gaussian elimination with partial pivoting
      const M = A.map((r, q) => [...r, b[q]]);
      for (let c = 0; c < 4; c++) {
        let piv = c; for (let r = c + 1; r < 4; r++) if (Math.abs(M[r][c]) > Math.abs(M[piv][c])) piv = r;
        if (Math.abs(M[piv][c]) < 1e-12) return null;
        [M[c], M[piv]] = [M[piv], M[c]];
        for (let r = 0; r < 4; r++) if (r !== c) { const f = M[r][c] / M[c][c]; for (let q = c; q < 5; q++) M[r][q] -= f * M[c][q]; }
      }
      return M.map((r, q) => r[4] / r[q]);
    };
    const U = solve([uvs[i][0], uvs[j][0], uvs[k][0], 0]), V = solve([uvs[i][1], uvs[j][1], uvs[k][1], 0]);
    if (U && V && [...U, ...V].every(Number.isFinite)) return imagePaint(id, U, V, avg);
  }
  return null;
}

const r7 = (x) => { const v = Number(x.toPrecision(7)); return Object.is(v, -0) ? 0 : v; };