Freeform bodies: subdivision surface and cage edits

/** One Catmull–Clark step. `track` (optional): per face { root, uv:[[u,v]…] }, refined alongside (quads). */
export function subdivide(m, track = null) {
  const t = topo(m);
  const nV = m.P.length;
  const fp = m.F.map((f) => avg(f.map((v) => m.P[v])));
  const eIdx = new Map();
  const P = m.P.map(() => null);
  const edges = [...t.E.values()];
  edges.forEach((e, i) => {
    eIdx.set(e.key, nV + i);
    const mid = lerp(m.P[e.a], m.P[e.b], 0.5);
    let p;
    if (e.w >= 1) p = mid;
    else {
      const sm = avg([m.P[e.a], m.P[e.b], fp[e.f[0]], fp[e.f[1]]]);
      p = e.w > 0 ? lerp(sm, mid, e.w) : sm;
    }
    P.push(p);
  });
  for (const p of fp) P.push(p);
  for (let v = 0; v < nV; v++) {
    const pv = m.P[v];
    const n = t.ve[v].length;
    if (!n) { P[v] = [...pv]; continue; }
    const Q = avg(t.vf[v].map((fi) => fp[fi]));
    const R = avg(t.ve[v].map((k) => { const e = t.E.get(k); return lerp(m.P[e.a], m.P[e.b], 0.5); }));
    const smooth = comb([[1 / n, Q], [2 / n, R], [(n - 3) / n, pv]]);
    const s = sharpAt(t, v);
    const corner = isCorner(m, t, v, s);
    if (s.length < 2 && !corner) { P[v] = smooth; continue; }
    const sharpP = corner ? [...pv] : comb([[1 / 8, m.P[other(s[0], v)]], [6 / 8, pv], [1 / 8, m.P[other(s[1], v)]]]);
    const sh = m.C.has(v) ? Infinity : s.reduce((a, e) => a + e.w, 0) / s.length;
    P[v] = sh >= 1 ? sharpP : lerp(smooth, sharpP, sh);
  }
  const F = [], nt = track ? [] : null;
  m.F.forEach((f, fi) => {
    const fv = nV + edges.length + fi;
    const L = f.length;
    for (let k = 0; k < L; k++) {
      const c = f[k];
      F.push([c, eIdx.get(ek(c, f[(k + 1) % L])), fv, eIdx.get(ek(f[(k + L - 1) % L], c))]);
      if (track) {
        const U = track[fi].uv, mid = (a, b) => [(a[0] + b[0]) / 2, (a[1] + b[1]) / 2];
        nt.push({ root: track[fi].root, uv: [U[k], mid(U[k], U[(k + 1) % L]), [U.reduce((s, q) => s + q[0], 0) / L, U.reduce((s, q) => s + q[1], 0) / L], mid(U[(k + L - 1) % L], U[k])] });
      }
    }
  });
  const W = new Map();
  for (const e of edges) {
    if (e.f.length < 2) continue;
    const w = m.W.get(e.key);
    if (!(w > 0)) continue;
    const w2 = w === Infinity ? Infinity : w - 1;
    if (w2 > 0) { const mi = eIdx.get(e.key); W.set(ek(e.a, mi), w2); W.set(ek(mi, e.b), w2); }
  }
  return { P, F, W, C: new Set(m.C), track: nt };
}

/** Exact limit position of vertex v of an all-quad mesh (remaining creases treated as sharp). */
export function limitPoint(m, t, v) {
  const s = sharpAt(t, v);
  const pv = m.P[v];
  if (isCorner(m, t, v, s)) return [...pv];
  if (s.length === 2) return comb([[1 / 6, m.P[other(s[0], v)]], [4 / 6, pv], [1 / 6, m.P[other(s[1], v)]]]);
  const n = t.ve[v].length;
  const terms = [[n * n, pv]];
  for (const k of t.ve[v]) terms.push([4, m.P[other(t.E.get(k), v)]]);
  for (const fi of t.vf[v]) { const f = m.F[fi]; const p = f.indexOf(v); terms.push([1, m.P[f[(p + 2) % 4]]]); }
  return comb(terms.map(([w, p]) => [w / (n * (n + 5)), p]));
}

/** Unit normal of the limit surface at an interior vertex (Halstead's tangent masks), or null. */
export function tangentNormal(m, t, v) {
  const r = ring(m, t, v);
  if (!r) return null;
  const n = r.length;
  const c = Math.cos(2 * Math.PI / n);
  const A = 1 + c + Math.cos(Math.PI / n) * Math.sqrt(2 * (9 + c));
  let t1 = [0, 0, 0], t2 = [0, 0, 0];
  r.forEach(({ e, f }, j) => {
    const a = 2 * Math.PI * j / n, b = 2 * Math.PI * (j + 1) / n;
    t1 = V.add(t1, V.add(V.mul(m.P[e], A * Math.cos(a)), V.mul(m.P[f], Math.cos(a) + Math.cos(b))));
    t2 = V.add(t2, V.add(V.mul(m.P[e], A * Math.sin(a)), V.mul(m.P[f], Math.sin(a) + Math.sin(b))));
  });
  const nn = V.cross(t1, t2);
  const l = V.len(nn);
  return l > 1e-300 ? V.mul(nn, 1 / l) : null;
}

/** Bézier net (4×4, [i][j]) of quad fi, and whether it is exact (all corners regular). */
export function quadPatch(m, t, fi, limits = null) {
  const f = m.F[fi];
  if (f.length !== 4) fail('form.err.quad_patch');
  const G = [0, 1, 2, 3].map(() => [null, null, null, null]);
  G[1][1] = m.P[f[0]]; G[2][1] = m.P[f[1]]; G[2][2] = m.P[f[2]]; G[1][2] = m.P[f[3]];
  [G[1][0], G[2][0]] = across(m, t, fi, 0);
  [G[3][1], G[3][2]] = across(m, t, fi, 1);
  [G[2][3], G[1][3]] = across(m, t, fi, 2);
  [G[0][2], G[0][1]] = across(m, t, fi, 3);
  const b = [0, 1, 2, 3].map(() => [null, null, null, null]);
  for (let i = 0; i < 4; i++) for (let j = 0; j < 4; j++) {
    if ((i === 0 || i === 3) && (j === 0 || j === 3)) continue;
    const s = [0, 0, 0];
    for (let a = 0; a < 4; a++) {
      if (!MB[i][a]) continue;
      for (let c = 0; c < 4; c++) {
        const w = MB[i][a] * MB[j][c];
        if (!w) continue;
        const p = G[a][c];
        s[0] += w * p[0]; s[1] += w * p[1]; s[2] += w * p[2];
      }
    }
    b[i][j] = s;
  }
  const lim = (v) => (limits ? limits[v] : limitPoint(m, t, v));
  const CORN = [[0, 0, 1, 0, 0, 1], [3, 0, 2, 0, 3, 1], [3, 3, 3, 2, 2, 3], [0, 3, 0, 2, 1, 3]];   // corner, its two edge neighbours
  let exact = true;
  f.forEach((v, k) => {
    const [ci, cj, ai, aj, bi, bj] = CORN[k];
    b[ci][cj] = lim(v);
    if (regularCorner(m, t, fi, v)) return;
    exact = false;
    const s = sharpAt(t, v);
    if (s.length >= 2 || isCorner(m, t, v, s)) return;   // crease/corner: the curves along sharp edges are already exact
    const n = tangentNormal(m, t, v);
    if (!n) return;
    for (const [i, j] of [[ai, aj], [bi, bj]]) b[i][j] = V.sub(b[i][j], V.mul(n, V.dot(V.sub(b[i][j], b[ci][cj]), n)));
  });
  return { b, exact };
}

/** Point of a bicubic Bézier patch (b[i][j], i along u). */
export function evalPatch(b, u, v) {
  const bu = bern(u), bv = bern(v);
  const s = [0, 0, 0];
  for (let i = 0; i < 4; i++) for (let j = 0; j < 4; j++) { const w = bu[i] * bv[j], p = b[i][j]; s[0] += w * p[0]; s[1] += w * p[1]; s[2] += w * p[2]; }
  return s;
}

/** Unit normal (∂u × ∂v) of a patch. */
export function patchNormal(b, u, v) {
  const bu = bern(u), bv = bern(v), du = bernD(u), dv = bernD(v);
  const su = [0, 0, 0], sv = [0, 0, 0];
  for (let i = 0; i < 4; i++) for (let j = 0; j < 4; j++) {
    const p = b[i][j];
    for (let k = 0; k < 3; k++) { su[k] += du[i] * bv[j] * p[k]; sv[k] += bu[i] * dv[j] * p[k]; }
  }
  return V.norm(V.cross(su, sv));
}

/** The refinement level a conversion needs: all quads and no semi-sharp creases left; `level` may ask for more (≤ 3). */
export function conversionLevel(form, level = null) {
  const wmax = Math.max(0, ...(form.creases || []).map((c) => (c[2] >= SHARP ? 0 : c[2])));
  const nonQuad = form.faces.some((f) => f.length !== 4);
  let need = Math.max(nonQuad ? 1 : 0, Math.ceil(wmax - 1e-12));
  if (level === null) {
    const m = meshOf(form), t = topo(m);
    if (!need && m.F.some((f, fi) => f.some((v) => !regularCorner(m, t, fi, v)))) need = 1;   // isolate extraordinary points
    return Math.min(4, need);
  }
  return Math.min(4, Math.max(need, Math.round(level)));
}

/**
 * How well the patches meet: largest gap and normal angle (degrees) across shared edges — overall and
 * between exact patches only — sampled at `n` points per edge.
 */
export function patchReport(pr, n = 5) {
  const { mesh: m, topo: t, patches: ps } = pr;
  const SIDE = [(s) => [s, 0], (s) => [1, s], (s) => [1 - s, 1], (s) => [0, 1 - s]];
  let gap = 0, deg = 0, degExact = 0;
  for (const e of t.E.values()) {
    if (e.f.length !== 2) continue;
    const [f0, f1] = e.f;
    const k0 = m.F[f0].findIndex((v, i) => ek(v, m.F[f0][(i + 1) % 4]) === e.key);
    const k1 = m.F[f1].findIndex((v, i) => ek(v, m.F[f1][(i + 1) % 4]) === e.key);
    const sharp = e.w > 0;
    for (let i = 0; i <= n; i++) {
      const s = i / n;
      const [u0, v0] = SIDE[k0](s), [u1, v1] = SIDE[k1](1 - s);
      const p0 = evalPatch(ps[f0].b, u0, v0), p1 = evalPatch(ps[f1].b, u1, v1);
      gap = Math.max(gap, V.len(V.sub(p0, p1)));
      if (sharp) continue;
      const n0 = patchNormal(ps[f0].b, u0, v0), n1 = patchNormal(ps[f1].b, u1, v1);
      const a = Math.acos(Math.min(1, Math.max(-1, V.dot(n0, n1)))) * 180 / Math.PI;
      if (!Number.isFinite(a)) continue;
      deg = Math.max(deg, a);
      if (ps[f0].exact && ps[f1].exact) degExact = Math.max(degExact, a);
    }
  }
  return { gap, deg, degExact };
}

/**
 * Move cage vertices so their limit points land on a target surface (`proj` = closest point on it):
 * p ← p + (proj(L) − L), a few rounds. The limit points then lie on it; between them the surface is
 * the subdivision surface, not the analytic one.
 */
export function fitTo(form, proj, rounds = 400) {
  let f = clone(form);
  for (let r = 0; r < rounds; r++) {
    const m = subdivide(meshOf(f));
    const t = topo(m);
    let worst = 0;
    f.verts = f.verts.map((p, v) => { const L = limitPoint(m, t, v); const d = V.sub(proj(L), L); worst = Math.max(worst, V.len(d)); return V.add(p, d); });
    if (worst < 1e-12 * sizeOf(f)) break;
  }
  return f;
}

/** A torus about Y: n around the axis, m around the tube; R major, r minor. Every vertex is regular. */
export function torus(n = 8, m = 6, R = 50, r = 20) {
  n = Math.max(3, Math.min(64, Math.round(n))); m = Math.max(3, Math.min(64, Math.round(m)));
  const verts = [], faces = [];
  for (let i = 0; i < n; i++) for (let j = 0; j < m; j++) {
    const a = 2 * Math.PI * i / n, b = 2 * Math.PI * j / m;
    verts.push([(R + r * Math.cos(b)) * Math.cos(a), r * Math.sin(b), (R + r * Math.cos(b)) * Math.sin(a)]);
  }
  const id = (i, j) => (i % n) * m + (j % m);
  for (let i = 0; i < n; i++) for (let j = 0; j < m; j++) faces.push([id(i, j), id(i, j + 1), id(i + 1, j + 1), id(i + 1, j)]);
  return fromQuads(verts, faces);
}

/** Mirror partner of every vertex (itself on the plane, −1 when none). */
export function partners(form) {
  if (!form.sym) return form.verts.map(() => -1);
  const tol = 1e-6 * sizeOf(form);
  const cell = (p) => p.map((x) => Math.round(x / (tol * 10))).join(',');
  const grid = new Map();
  form.verts.forEach((p, i) => { const k = cell(p); if (!grid.has(k)) grid.set(k, []); grid.get(k).push(i); });
  return form.verts.map((p) => {
    const q = mirrorP(form.sym, p);
    const c = q.map((x) => Math.round(x / (tol * 10)));
    for (let dx = -1; dx <= 1; dx++) for (let dy = -1; dy <= 1; dy++) for (let dz = -1; dz <= 1; dz++) {
      for (const j of grid.get([c[0] + dx, c[1] + dy, c[2] + dz].join(',')) || []) if (V.len(V.sub(form.verts[j], q)) <= tol) return j;
    }
    return -1;
  });
}

/** Rebuild the mirrored half from the kept half (stable order: kept faces first). */
export function resym(form) {
  if (!form.sym) return form;
  const { sym } = form;
  const a = AXES[sym.axis];
  const tol = 1e-6 * sizeOf(form);
  const sd = (p) => (p[a] - sym.at) * sym.side;
  const keepF = form.faces.filter((f) => sd(avg(f.map((v) => form.verts[v]))) > tol);
  if (!keepF.length) fail('form.err.sym_empty');
  const used = [...new Set(keepF.flat())].sort((x, y) => x - y);
  const map = new Map(used.map((v, i) => [v, i]));
  const verts = used.map((v) => { const p = [...form.verts[v]]; if (sd(p) < tol) p[a] = sym.at; return p; });
  const mirrorOf = new Map();
  verts.forEach((p, i) => {
    if (Math.abs(p[a] - sym.at) <= tol) { mirrorOf.set(i, i); return; }
    mirrorOf.set(i, verts.length); verts.push(mirrorP(sym, p));
  });
  const faces = keepF.map((f) => f.map((v) => map.get(v)));
  for (const f of keepF.map((g) => g.map((v) => map.get(v)))) faces.push(f.map((v) => mirrorOf.get(v)).reverse());
  const creases = [];
  for (const [x, y, w] of form.creases || []) {
    if (!map.has(x) || !map.has(y)) continue;
    const X = map.get(x), Y = map.get(y);
    creases.push([X, Y, w]);
    const MX = mirrorOf.get(X), MY = mirrorOf.get(Y);
    if (MX !== X || MY !== Y) creases.push([MX, MY, w]);
  }
  const corners = [];
  for (const v of form.corners || []) if (map.has(v)) { const X = map.get(v); corners.push(X); if (mirrorOf.get(X) !== X) corners.push(mirrorOf.get(X)); }
  return checked({ ...emptyForm(), verts, faces, creases: dedupeCreases(creases), corners: [...new Set(corners)], sym: { ...sym } });
}

/**
 * Transform the selection by a 4×4 column-major matrix. `soft` > 0 pulls the vertices within that
 * distance along with a smooth falloff. With symmetry on, the mirror side follows (mirrored
 * transform) and seam vertices stay on the plane.
 */
export function transform(form, sel, M4, { soft = 0 } = {}) {
  const f = clone(form);
  const ids = selVerts(form, sel);
  if (!ids.length) return f;
  const w = new Map(ids.map((v) => [v, 1]));
  const pt = form.sym ? partners(form) : null;
  const src = new Set(ids);
  if (pt) for (const v of ids) if (pt[v] >= 0) src.add(pt[v]);
  if (soft > 0) {
    const S = [...src].map((v) => form.verts[v]);
    form.verts.forEach((p, v) => { if (w.has(v)) return; let d = Infinity; for (const q of S) d = Math.min(d, V.len(V.sub(p, q))); const x = falloff(d, soft); if (x > 0) w.set(v, x); });
  }
  let MR = null, sideSel = 1, a = 0;
  if (form.sym) {
    a = AXES[form.sym.axis];
    const R = V.reflection(V.mul([0, 0, 0].map((_, i) => (i === a ? 1 : 0)), form.sym.at), [0, 0, 0].map((_, i) => (i === a ? 1 : 0)));
    MR = V.matMul(R, V.matMul(M4, R));
    const s = ids.reduce((acc, v) => acc + (form.verts[v][a] - form.sym.at), 0);
    sideSel = s < 0 ? -1 : 1;
    for (const v of ids) if (pt[v] >= 0 && pt[v] !== v && !w.has(pt[v])) w.set(pt[v], 1);
  }
  const tol = 1e-6 * sizeOf(form);
  for (const [v, wt] of w) {
    const p = form.verts[v];
    let M = M4;
    if (MR) { const d = p[a] - form.sym.at; if (Math.abs(d) > tol && Math.sign(d) !== sideSel) M = MR; }
    const q = V.applyMat(M, p);
    let r = lerp(p, q, wt);
    if (pt && pt[v] === v) r[a] = form.sym.at;
    f.verts[v] = r;
  }
  return f;
}

/** Extrude faces as one region: the region moves `dist` along its averaged normals and side walls join it. */
export function extrudeFaces(form, faceIds, dist = 0) {
  return symEdit(form, { faces: faceIds }, (f, sel) => {
    const set = new Set(sel.faces.filter((i) => f.faces[i]));
    if (!set.size) fail('form.err.pick_extrude');
    const nrm = new Map();
    for (const fi of set) {
      const g = f.faces[fi];
      const n = newell(g.map((v) => f.verts[v]));
      for (const v of g) nrm.set(v, V.add(nrm.get(v) || [0, 0, 0], n));
    }
    const dup = new Map();
    for (const [v, n] of nrm) {
      let d = V.mul(V.norm(n), dist);
      if (f.sym) { const a = AXES[f.sym.axis]; if (Math.abs(f.verts[v][a] - f.sym.at) <= 1e-6 * sizeOf(f)) d[a] = 0; }
      dup.set(v, f.verts.length); f.verts.push(V.add(f.verts[v], d));
    }
    const dirUse = new Map();
    for (const fi of set) { const g = f.faces[fi]; g.forEach((v, k) => dirUse.set(ek(v, g[(k + 1) % g.length]), (dirUse.get(ek(v, g[(k + 1) % g.length])) || 0) + 1)); }
    const sides = [];
    for (const fi of set) {
      const g = f.faces[fi];
      g.forEach((v, k) => {
        const w = g[(k + 1) % g.length];
        if (dirUse.get(ek(v, w)) === 1) sides.push([v, w, dup.get(w), dup.get(v)]);
      });
      f.faces[fi] = g.map((v) => dup.get(v));
    }
    // with symmetry, a wall standing on the seam is inside the body: leave it out
    const seam = (v) => f.sym && Math.abs(f.verts[v][AXES[f.sym.axis]] - f.sym.at) <= 1e-6 * sizeOf(f);
    for (const s of sides) if (!(seam(s[0]) && seam(s[1]))) f.faces.push(s);
    // creases on the region's edges move up with it
    f.creases = f.creases.map(([a, b, w]) => (dirUse.has(ek(a, b)) ? [dup.get(a), dup.get(b), w] : [a, b, w]));
    return compact(f);
  });
}

function newell(pts) {
  const n = [0, 0, 0];
  for (let i = 0; i < pts.length; i++) { const p = pts[i], q = pts[(i + 1) % pts.length]; n[0] += (p[1] - q[1]) * (p[2] + q[2]); n[1] += (p[2] - q[2]) * (p[0] + q[0]); n[2] += (p[0] - q[0]) * (p[1] + q[1]); }
  return n;
}

/** Flatten vertices onto a plane (default: the least-squares plane through them). */
export function flatten(form, sel, pl = null) {
  const ids = selVerts(form, sel);
  if (ids.length < 3 && !pl) fail('form.err.pick_flatten');
  let o, n;
  if (pl) { o = pl.o; n = V.norm(pl.n); } else { o = avg(ids.map((v) => form.verts[v])); n = bestNormal(ids.map((v) => V.sub(form.verts[v], o))); }
  const f = clone(form);
  const pt = form.sym ? partners(form) : null;
  const todo = new Set(ids);
  if (pt) for (const v of ids) if (pt[v] >= 0) todo.add(pt[v]);
  for (const v of todo) {
    const isMirror = pt && !ids.includes(v);
    const oo = isMirror ? mirrorP(form.sym, o) : o;
    const nn = isMirror ? mirrorP({ ...form.sym, at: 0 }, n) : n;
    const p = form.verts[v];
    f.verts[v] = V.sub(p, V.mul(nn, V.dot(V.sub(p, oo), nn)));
  }
  return f;
}

/** Smallest-variance direction of centred points (Jacobi eigen-solve of the 3×3 covariance). */
export function bestNormal(pts) {
  const A = [[0, 0, 0], [0, 0, 0], [0, 0, 0]];
  for (const p of pts) for (let i = 0; i < 3; i++) for (let j = 0; j < 3; j++) A[i][j] += p[i] * p[j];
  const Vm = [[1, 0, 0], [0, 1, 0], [0, 0, 1]];
  for (let sweep = 0; sweep < 50; sweep++) {
    let off = 0;
    for (let p = 0; p < 3; p++) for (let q = p + 1; q < 3; q++) off += A[p][q] ** 2;
    if (off < 1e-30) break;
    for (let p = 0; p < 3; p++) for (let q = p + 1; q < 3; q++) {
      if (Math.abs(A[p][q]) < 1e-300) continue;
      const th = (A[q][q] - A[p][p]) / (2 * A[p][q]);
      const tt = Math.sign(th || 1) / (Math.abs(th) + Math.sqrt(th * th + 1));
      const c = 1 / Math.sqrt(tt * tt + 1), s = tt * c;
      for (let k = 0; k < 3; k++) { const akp = A[k][p], akq = A[k][q]; A[k][p] = c * akp - s * akq; A[k][q] = s * akp + c * akq; }
      for (let k = 0; k < 3; k++) { const apk = A[p][k], aqk = A[q][k]; A[p][k] = c * apk - s * aqk; A[q][k] = s * apk + c * aqk; }
      for (let k = 0; k < 3; k++) { const vkp = Vm[k][p], vkq = Vm[k][q]; Vm[k][p] = c * vkp - s * vkq; Vm[k][q] = s * vkp + c * vkq; }
    }
  }
  let mi = 0;
  for (let i = 1; i < 3; i++) if (A[i][i] < A[mi][mi]) mi = i;
  return V.norm([Vm[0][mi], Vm[1][mi], Vm[2][mi]]);
}

/** Add one face through points (3+), reusing existing vertices within `tol` (builds a surface click by click). */
export function addFace(form, pts, { tol = 1e-6 } = {}) {
  if (!Array.isArray(pts) || pts.length < 3 || pts.length > 64) fail('form.err.face_points');
  const f = clone(form);
  const ids = pts.map((p) => {
    if (typeof p === 'number') { if (!f.verts[p]) fail('form.err.face_point_gone'); return p; }
    const i = f.verts.findIndex((q) => V.len(V.sub(q, p)) <= tol);
    if (i >= 0) return i;
    f.verts.push([...p]);
    return f.verts.length - 1;
  });
  if (new Set(ids).size !== ids.length) fail('form.err.face_dup');
  f.faces.push(ids);
  // match the neighbours' orientation: reverse when it shares an edge in the same direction
  const used = new Set();
  form.faces.forEach((g) => g.forEach((v, k) => used.add(v + '>' + g[(k + 1) % g.length])));
  if (ids.some((v, k) => used.has(v + '>' + ids[(k + 1) % ids.length]))) f.faces[f.faces.length - 1] = ids.reverse();
  return checked(f);
}

/** Bridge two faces or two open loops (A, B: { face } | { edge }) of equal size with `segments` rows of quads. */
export function bridge(form, A, B, { segments = 1 } = {}) {
  const sel = { faces: [A.face, B.face].filter((x) => x !== undefined), edges: [A.edge, B.edge].filter((x) => x !== undefined) };
  return symEdit(form, sel, (f, s) => {
    let fi = 0, ei = 0;
    const take = (X) => (X.face !== undefined ? { face: s.faces[fi++] } : { edge: s.edges[ei++] });
    const a = take(A), b = take(B);
    const loopOf = (X) => (X.face !== undefined ? (f.faces[X.face] || fail('form.err.bridge_gone')) : boundaryLoop(f, X.edge));
    const la = [...loopOf(a)], lb = [...loopOf(b)];
    if (la.length !== lb.length) fail('form.err.bridge_count', { a: la.length, b: lb.length });
    if (la.some((v) => lb.includes(v))) fail('form.err.bridge_share');
    const n = la.length;
    let best = 0, bd = Infinity;
    for (let o = 0; o < n; o++) { let d = 0; for (let i = 0; i < n; i++) d += V.len(V.sub(f.verts[la[i]], f.verts[lb[((o - i) % n + n) % n]])) ** 2; if (d < bd) { bd = d; best = o; } }
    const gone = new Set([a.face, b.face].filter((x) => x !== undefined));
    f.faces = f.faces.filter((_, i) => !gone.has(i));
    const segs = Math.max(1, Math.min(32, Math.round(segments)));
    const rows = [la];
    for (let r = 1; r < segs; r++) rows.push(la.map((v, i) => { f.verts.push(lerp(f.verts[v], f.verts[lb[((best - i) % n + n) % n]], r / segs)); return f.verts.length - 1; }));
    rows.push(la.map((_, i) => lb[((best - i) % n + n) % n]));
    for (let r = 0; r < segs; r++) for (let i = 0; i < n; i++) {
      const i1 = (i + 1) % n;
      f.faces.push([rows[r][i], rows[r][i1], rows[r + 1][i1], rows[r + 1][i]]);
    }
    return compact(f);
  });
}