Fitted curves and surfaces: sweeps, networks, Coons patches, offsets, rebuilds and extensions

export function fitFixedKnots(Q, t, p, U, { fixEnds = true, weights = null } = {}) {
  const n = U.length - p - 2, m = Q.length - 1;
  const P = new Array(n + 1);
  if (fixEnds) { P[0] = Q[0].slice(0, 3); P[n] = Q[m].slice(0, 3); }
  const lo = fixEnds ? 1 : 0, hi = fixEnds ? n - 1 : n, k = hi - lo + 1;
  if (k <= 0) return { deg: p, knots: U.slice(), cvs: P.map((q, i) => q || Q[Math.round(i * m / n)].slice(0, 3)), w: null };
  const NtN = Array.from({ length: k }, () => new Array(k).fill(0)), NtR = Array.from({ length: k }, () => [0, 0, 0]);
  for (let i = 0; i <= m; i++) {
    if (fixEnds && (i === 0 || i === m)) continue;
    const row = rowAt(p, U, n, t[i]), wt = weights ? weights[i] : 1;
    let r = Q[i].slice(0, 3);
    if (fixEnds) r = V.sub(V.sub(r, V.mul(P[0], row[0])), V.mul(P[n], row[n]));
    for (let a = lo; a <= hi; a++) {
      if (!row[a]) continue;
      for (let b = lo; b <= hi; b++) if (row[b]) NtN[a - lo][b - lo] += wt * row[a] * row[b];
      for (let d = 0; d < 3; d++) NtR[a - lo][d] += wt * row[a] * r[d];
    }
  }
  // a control point no sample reaches would leave the system singular: tie it lightly to its neighbours (smoothing)
  for (let a = 0; a < k; a++) if (NtN[a][a] < 1e-12) {
    NtN[a][a] += 2e-6;
    if (a > 0) NtN[a][a - 1] -= 1e-6;
    if (a + 1 < k) NtN[a][a + 1] -= 1e-6;
  }
  const X = V.solve(NtN, NtR);
  if (!X) fail('nurbs.err.fit_singular');
  for (let a = lo; a <= hi; a++) P[a] = X[a - lo];
  return { deg: p, knots: U.slice(), cvs: P, w: null };
}

function lsqMany(lists, t, p, U) {
  // least squares with fixed knots and ends, many lists of points (any dimension) sharing one matrix
  const n = U.length - p - 2, m = t.length - 1, dim = lists[0][0].length;
  if (n === m) {
    const A = t.map((u) => rowAt(p, U, n, u)), X = V.solve(A, t.map((_, i) => lists.flatMap((L) => L[i])));
    if (!X) fail('nurbs.err.fit_singular');
    return lists.map((_, k) => X.map((r) => r.slice(k * dim, k * dim + dim)));
  }
  const k = n - 1, rows = [];
  for (let i = 1; i < m; i++) rows.push(rowAt(p, U, n, t[i]));
  const NtN = Array.from({ length: k }, (_, a) => Array.from({ length: k }, (_, b) => { let s = 0; for (const r of rows) s += r[a + 1] * r[b + 1]; return s; }));
  for (let a = 0; a < k; a++) if (NtN[a][a] < 1e-12) NtN[a][a] = 1e-12;
  const R = Array.from({ length: k }, () => new Array(lists.length * dim).fill(0));
  rows.forEach((r, ii) => {
    const i = ii + 1;
    lists.forEach((L, li) => {
      for (let d = 0; d < dim; d++) {
        const rv = L[i][d] - r[0] * L[0][d] - r[n] * L[m][d];
        for (let a = 0; a < k; a++) if (r[a + 1]) R[a][li * dim + d] += r[a + 1] * rv;
      }
    });
  });
  const X = V.solve(NtN, R);
  if (!X) fail('nurbs.err.fit_singular');
  return lists.map((L, li) => [L[0].slice(), ...X.map((r) => r.slice(li * dim, li * dim + dim)), L[m].slice()]);
}

export function fitTo(f, p, tol, { maxCvs = 400, samples = 6, start = null } = {}) {
  let count = Math.max(p + 1, start || p + 1);
  let best = null;
  for (;;) {
    const S = Math.max(count * samples, 24), t = [], Q = [];
    for (let i = 0; i <= S; i++) { t.push(i / S); Q.push(f(i / S)); }
    const c = fitFixedKnots(Q, t, p, clampedUniform(count, p));
    let max = 0;
    for (let i = 0; i < S; i++) { const s = (i + 0.5) / S, d = V.dist(C.point(c, s), f(s)); if (d > max) max = d; }
    best = { curve: c, deviation: max };
    if (max <= tol || count >= maxCvs) return best;
    count = Math.min(maxCvs, Math.ceil(count * 1.5) + 1);
  }
}

export function fitSurfaceTo(f, du, dv, tol, { maxCvs = 60, samples = 4, start = null } = {}) {
  let nu = Math.max(du + 1, start ? start[0] : du + 1), nv = Math.max(dv + 1, start ? start[1] : dv + 1);
  const evalS = (s, u, v) => {
    // inline tensor evaluation (polynomial surface) to keep fit.js free of surface.js
    const ku = findSpan(s.cvs.length - 1, s.degU, u, s.knotsU), kv = findSpan(s.cvs[0].length - 1, s.degV, v, s.knotsV);
    const Nu = basisFuns(ku, u, s.degU, s.knotsU), Nv = basisFuns(kv, v, s.degV, s.knotsV), o = [0, 0, 0];
    for (let a = 0; a <= s.degU; a++) for (let b = 0; b <= s.degV; b++) { const q = s.cvs[ku - s.degU + a][kv - s.degV + b], g = Nu[a] * Nv[b]; o[0] += g * q[0]; o[1] += g * q[1]; o[2] += g * q[2]; }
    return o;
  };
  for (;;) {
    const Su = Math.max(nu * samples, 12), Sv = Math.max(nv * samples, 12);
    const tu = Array.from({ length: Su + 1 }, (_, i) => i / Su), tv = Array.from({ length: Sv + 1 }, (_, j) => j / Sv);
    const Q = tu.map((u) => tv.map((v) => f(u, v)));
    const s = fitSurfaceGrid(Q, tu, tv, du, dv, clampedUniform(nu, du), clampedUniform(nv, dv));
    let maxU = 0, maxV = 0;
    for (let i = 0; i < Su; i++) for (let j = 0; j <= Sv; j += Math.max(1, Sv >> 3)) { const u = (i + 0.5) / Su; maxU = Math.max(maxU, V.dist(evalS(s, u, tv[j]), f(u, tv[j]))); }
    for (let j = 0; j < Sv; j++) for (let i = 0; i <= Su; i += Math.max(1, Su >> 3)) { const v = (j + 0.5) / Sv; maxV = Math.max(maxV, V.dist(evalS(s, tu[i], v), f(tu[i], v))); }
    let mid = 0;
    for (let i = 0; i < Su; i += Math.max(1, Su >> 3)) for (let j = 0; j < Sv; j += Math.max(1, Sv >> 3)) { const u = (i + 0.5) / Su, v = (j + 0.5) / Sv; mid = Math.max(mid, V.dist(evalS(s, u, v), f(u, v))); }
    const dev = Math.max(maxU, maxV, mid);
    if (dev <= tol || (nu >= maxCvs && nv >= maxCvs)) return { surface: s, deviation: dev };
    const nu0 = nu, nv0 = nv;
    if (maxU >= maxV * 0.5 && nu < maxCvs) nu = Math.min(maxCvs, Math.ceil(nu * 1.5) + 1);
    if (maxV >= maxU * 0.5 && nv < maxCvs) nv = Math.min(maxCvs, Math.ceil(nv * 1.5) + 1);
    if (nu === nu0 && nv === nv0) { if (nu < maxCvs) nu = Math.min(maxCvs, Math.ceil(nu * 1.5) + 1); else nv = Math.min(maxCvs, Math.ceil(nv * 1.5) + 1); }
  }
}

export function rebuild(c, count, p = c.deg, { samples = null } = {}) {
  const [a, b] = C.domain(c), S = Math.max(samples || 0, count * 4, 32);
  const Q = [], t = [];
  for (let i = 0; i <= S; i++) { const u = a + (b - a) * i / S; Q.push(C.point(c, u)); t.push(i / S); }
  const deg = Math.min(p, count - 1);
  const out = approximate(Q, deg, count, { params: t, knots: clampedUniform(count, deg) });
  // measure between the fitting samples
  const M = [];
  for (let i = 0; i < S; i++) M.push(C.point(c, a + (b - a) * (i + 0.5) / S));
  return { curve: out, deviation: deviation(out, M) };
}

export function deviation(c, Q) {
  let max = 0, sum = 0, at = 0;
  Q.forEach((q, i) => { const d = C.closestParam(c, q).d; sum += d; if (d > max) { max = d; at = i; } });
  return { max, mean: Q.length ? sum / Q.length : 0, at };
}

/** Frames along a curve at parameters ts: rotation-minimising (double reflection) from the start's frame, or fixed. */
export function frames(rail, ts, { frame = 'rmf', up = null } = {}) {
  const fine = [];
  const [a, b] = C.domain(rail), N = Math.max(64, ts.length * 8);
  for (let i = 0; i <= N; i++) fine.push(a + (b - a) * i / N);
  const all = [...new Set([...fine, ...ts])].sort((x, y) => x - y);
  const T0 = C.tangent(rail, all[0]);
  let r = up ? V.norm(V.sub(up, V.mul(T0, V.dot(up, T0)))) : null;
  if (!r || V.len(r) < 1e-9) {
    const kv = C.curvatureVector(rail, all[0]);
    r = V.len(kv) > 1e-9 ? V.norm(kv) : V.perp(T0);
  }
  const out = new Map();
  let x = C.point(rail, all[0]), t = T0;
  out.set(all[0], { o: x, t, r, s: V.cross(t, r) });
  for (let i = 1; i < all.length; i++) {
    const x1 = C.point(rail, all[i]), t1 = C.tangent(rail, all[i]);
    if (frame === 'fixed') r = V.norm(V.sub(r, V.mul(t1, V.dot(r, t1))));
    else {
      const v1 = V.sub(x1, x), c1 = V.dot(v1, v1);
      if (c1 > 1e-300) {
        const rL = V.sub(r, V.mul(v1, 2 / c1 * V.dot(v1, r))), tL = V.sub(t, V.mul(v1, 2 / c1 * V.dot(v1, t)));
        const v2 = V.sub(t1, tL), c2 = V.dot(v2, v2);
        r = c2 > 1e-300 ? V.sub(rL, V.mul(v2, 2 / c2 * V.dot(v2, rL))) : rL;
      }
      r = V.norm(V.sub(r, V.mul(t1, V.dot(r, t1))));
    }
    x = x1; t = t1;
    out.set(all[i], { o: x1, t: t1, r, s: V.cross(t1, r) });
  }
  return ts.map((u) => out.get(u));
}

function placeProfile(profile, f0, f, scale = 1, twist = 0) {
  // profile coordinates in the start frame (r, s, t), turned by `twist` about t, scaled, set in frame f
  const c = Math.cos(twist), sn = Math.sin(twist);
  return {
    deg: profile.deg, knots: profile.knots.slice(), w: profile.w ? profile.w.slice() : null,
    cvs: profile.cvs.map((p) => {
      const d = V.sub(p, f0.o);
      let a = V.dot(d, f0.r), b = V.dot(d, f0.s);
      const z = V.dot(d, f0.t);
      [a, b] = [a * c - b * sn, a * sn + b * c];
      return V.add(f.o, V.add(V.add(V.mul(f.r, a * scale), V.mul(f.s, b * scale)), V.mul(f.t, z * scale)));
    }),
  };
}

function sweepBy(place, { tol, maxSections, degree = 3 }) {
  // place(t) → the section at rail parameter t ∈ [0, 1]; sections added until the mid-sections are within tol
  let n = 4, last = null;
  for (;;) {
    const ts = Array.from({ length: n + 1 }, (_, i) => i / n);
    const secs = place(ts);
    const { surface } = loft(secs, { degree, params: ts });
    const mids = ts.slice(0, -1).map((x, i) => (x + ts[i + 1]) / 2), ref = place(mids);
    let dev = 0;
    mids.forEach((v, k) => {
      const sec = ref[k], [a, b] = C.domain(sec);
      for (let i = 0; i <= 12; i++) { const u = a + (b - a) * i / 12; dev = Math.max(dev, V.dist(S.point(surface, (u - a) / (b - a), v), C.point(sec, u))); }
    });
    last = { surface, report: { exact: false, deviation: dev, sections: n + 1, tol } };
    if (dev <= tol || n >= maxSections) return last;
    n = Math.min(maxSections, n * 2);
  }
}

const frameAt = (t) => {
  const A = C.point(r1, t), Bp = C.point(r2, t), x = V.sub(Bp, A), d = V.len(x);
  if (d < 1e-12) fail('nurbs.err.rails_meet');
  const X = V.mul(x, 1 / d);
  let Z = V.add(C.tangent(r1, t), C.tangent(r2, t));
  Z = V.sub(Z, V.mul(X, V.dot(Z, X)));
  if (V.len(Z) < 1e-12) Z = V.perp(X);
  Z = V.norm(Z);
  return { o: A, X, Y: V.cross(Z, X), Z, d };
};

/** Closest pair of parameters between two curves (polylines first, then Newton on |A(s) − B(t)|²). */
export function curveCurve(A, Bc) {
  const [a0, a1] = C.domain(A), [b0, b1] = C.domain(Bc);
  const K = 96;
  let best = { s: a0, t: b0, d: Infinity };
  for (let i = 0; i <= K; i++) {
    const s = a0 + (a1 - a0) * i / K, p = C.point(A, s), q = C.closestParam(Bc, p);
    if (q.d < best.d) best = { s, t: q.u, d: q.d };
  }
  let { s, t } = best;
  for (let it = 0; it < 40; it++) {
    const da = C.derivs(A, s, 2), db = C.derivs(Bc, t, 2), r = V.sub(da[0], db[0]);
    const f = V.dot(r, da[1]), g = -V.dot(r, db[1]);
    const J00 = V.dot(da[1], da[1]) + V.dot(r, da[2]), J01 = -V.dot(da[1], db[1]), J11 = V.dot(db[1], db[1]) - V.dot(r, db[2]);
    const det = J00 * J11 - J01 * J01;
    if (Math.abs(det) < 1e-300) break;
    const ds = -(J11 * f - J01 * g) / det, dt = -(J00 * g - J01 * f) / det;
    const sn = Math.min(a1, Math.max(a0, s + ds)), tn = Math.min(b1, Math.max(b0, t + dt));
    const moved = Math.abs(sn - s) + Math.abs(tn - t);
    s = sn; t = tn;
    if (moved < 1e-15) break;
  }
  const d = V.dist(C.point(A, s), C.point(Bc, t));
  return d <= best.d ? { s, t, d } : best;
}

function pchip(x, y) {
  const n = x.length, h = [], dlt = [], m = new Array(n).fill(0);
  for (let i = 0; i + 1 < n; i++) { h.push(x[i + 1] - x[i]); dlt.push((y[i + 1] - y[i]) / h[i]); }
  if (n === 2) { m[0] = m[1] = dlt[0]; }
  else {
    for (let i = 1; i + 1 < n; i++) {
      if (dlt[i - 1] * dlt[i] <= 0) m[i] = 0;
      else { const w1 = 2 * h[i] + h[i - 1], w2 = h[i] + 2 * h[i - 1]; m[i] = (w1 + w2) / (w1 / dlt[i - 1] + w2 / dlt[i]); }
    }
    m[0] = ((2 * h[0] + h[1]) * dlt[0] - h[0] * dlt[1]) / (h[0] + h[1]);
    if (Math.sign(m[0]) !== Math.sign(dlt[0])) m[0] = 0; else if (Math.sign(dlt[0]) !== Math.sign(dlt[1]) && Math.abs(m[0]) > 3 * Math.abs(dlt[0])) m[0] = 3 * dlt[0];
    const k = n - 1;
    m[k] = ((2 * h[k - 1] + h[k - 2]) * dlt[k - 1] - h[k - 1] * dlt[k - 2]) / (h[k - 1] + h[k - 2]);
    if (Math.sign(m[k]) !== Math.sign(dlt[k - 1])) m[k] = 0; else if (Math.sign(dlt[k - 1]) !== Math.sign(dlt[k - 2]) && Math.abs(m[k]) > 3 * Math.abs(dlt[k - 1])) m[k] = 3 * dlt[k - 1];
  }
  return (t) => {
    let i = 0;
    while (i + 2 < n && t > x[i + 1]) i++;
    const hh = h[i], s = (t - x[i]) / hh, s2 = s * s, s3 = s2 * s;
    return (2 * s3 - 3 * s2 + 1) * y[i] + (s3 - 2 * s2 + s) * hh * m[i] + (-2 * s3 + 3 * s2) * y[i + 1] + (s3 - s2) * hh * m[i + 1];
  };
}

export function network(profiles, guides, { tol = 0.01, fitTol = null } = {}) {
  if (profiles.length < 2 || guides.length < 2) fail('nurbs.err.network_size');
  fitTol = fitTol == null ? tol / 10 : fitTol;
  let P = profiles.map((c) => C.reparam(c, 0, 1)), G = guides.map((c) => C.reparam(c, 0, 1));
  // order: guides along the first profile, profiles along the first guide; curves running backwards turned round
  const orderBy = (list, ref) => list.map((c, i) => ({ c, i, at: curveCurve(ref, c) })).sort((x, y) => x.at.s - y.at.s);
  G = orderBy(G, P[0]).map((x) => x.c);
  P = orderBy(P, G[0]).map((x) => x.c);
  const X = P.map((p) => G.map((g) => curveCurve(p, g)));
  let gap = 0;
  for (const row of X) for (const x of row) gap = Math.max(gap, x.d);
  if (gap > tol) fail('nurbs.err.network_gap', { gap });
  P = P.map((p, k) => (X[k][X[k].length - 1].s < X[k][0].s ? C.reverse(p) : p));
  G = G.map((g, j) => (X[X.length - 1][j].t < X[0][j].t ? C.reverse(g) : g));
  const Y = P.map((p) => G.map((g) => curveCurve(p, g)));
  const ownU = Y.map((row) => row.map((x) => x.s)), ownV = G.map((_, j) => Y.map((row) => row[j].t));
  const ubar = G.map((_, j) => ownU.reduce((s, r) => s + r[j], 0) / P.length), vbar = P.map((_, k) => ownV.reduce((s, r) => s + r[k], 0) / G.length);
  // the outer curves bound the surface: ū and v̄ run 0…1 between the first and last curve of the other family
  const norm01 = (a) => a.map((x) => (x - a[0]) / (a[a.length - 1] - a[0]));
  const ub = norm01(ubar), vb = norm01(vbar);
  const trimTo = (c, a, b) => C.reparam(C.segment(c, Math.max(a, 0), Math.min(b, 1)), 0, 1);
  const Pt = P.map((p, k) => trimTo(p, ownU[k][0], ownU[k][G.length - 1]));
  const Gt = G.map((g, j) => trimTo(g, ownV[j][0], ownV[j][P.length - 1]));
  const ownU2 = Pt.map((_, k) => norm01(ownU[k])), ownV2 = Gt.map((_, j) => norm01(ownV[j]));
  const fp = family(Pt, ownU2, ub, fitTol), fg = family(Gt, ownV2, vb, fitTol);
  const Lu = loft(fp.curves, { params: vb, degree: Math.min(3, P.length - 1) }).surface;
  const Lv = S.swap(loft(fg.curves, { params: ub, degree: Math.min(3, G.length - 1) }).surface);
  // T: tensor interpolation of the crossing points (the average of the two curves' points there)
  const Q = ub.map((u, j) => vb.map((v, k) => V.mul(V.add(C.point(fp.curves[k], u), C.point(fg.curves[j], v)), 0.5)));
  const du = Math.min(3, G.length - 1), dv = Math.min(3, P.length - 1);
  const rowsV = F.interpolateMany(Q, vb, dv);
  const cols = rowsV.rows[0].map((_, jj) => rowsV.rows.map((r) => r[jj]));
  const colsU = F.interpolateMany(cols, ub, du);
  const T = { degU: colsU.deg, degV: rowsV.deg, knotsU: colsU.knots, knotsV: rowsV.knots, cvs: colsU.rows[0].map((_, i) => colsU.rows.map((col) => col[i])), w: null };
  const [a, b, c] = makeSame([Lu, Lv, T]);
  const surface = sumNets(a, b, c);
  // measured: every curve of the network against the surface
  let dev = 0;
  fp.curves.forEach((cv, k) => { for (let i = 0; i <= 24; i++) dev = Math.max(dev, V.dist(S.point(surface, i / 24, vb[k]), C.point(Pt[k], ownU2[k].length ? pchip(ub, ownU2[k])(i / 24) : i / 24))); });
  fg.curves.forEach((cv, j) => { for (let i = 0; i <= 24; i++) dev = Math.max(dev, V.dist(S.point(surface, ub[j], i / 24), C.point(Gt[j], pchip(vb, ownV2[j])(i / 24)))); });
  return { surface, report: { exact: fp.exact && fg.exact, deviation: dev, gap, profiles: P.length, guides: G.length, uParams: ub, vParams: vb } };
}

export function coons(c0, c1, d0, d1, { tol = 1e-6 } = {}) {
  let [a, b, l, r] = [c0, c1, d0, d1].map((c) => C.reparam(c, 0, 1));
  // orient: a and b run from l to r; l and r run from a to b
  const near = (p, q) => V.dist(p, q);
  const E = (c) => C.endsOf(c);
  if (near(E(a).start, E(l).start) > near(E(a).end, E(l).start) && near(E(a).end, E(l).start) < near(E(a).start, E(l).end)) a = C.reverse(a);
  if (near(E(l).start, E(a).start) > near(E(l).end, E(a).start)) l = C.reverse(l);
  if (near(E(r).start, E(a).end) > near(E(r).end, E(a).end)) r = C.reverse(r);
  if (near(E(b).start, E(l).end) > near(E(b).end, E(l).end)) b = C.reverse(b);
  const P00 = E(a).start, P10 = E(a).end, P01 = E(b).start, P11 = E(b).end;
  const gap = Math.max(near(P00, E(l).start), near(P10, E(r).start), near(P01, E(l).end), near(P11, E(r).end));
  const size = Math.max(near(P00, P11), near(P10, P01), 1e-9);
  if (gap > Math.max(tol, 1e-9 * size) * 1e3) fail('nurbs.err.corners', { gap });
  let dev = 0;
  [a, b, l, r] = [a, b, l, r].map((c) => { const q = polynomial(c, tol); dev = Math.max(dev, q.deviation); return q.curve; });
  const Ru = ruled(a, b).surface, Rv = S.swap(ruled(l, r).surface);
  const T = planarPatch([P00, P10, P11, P01]);
  const [x, y, z] = makeSame([Ru, Rv, T]);
  return { surface: sumNets(x, y, z), report: { exact: dev === 0, deviation: dev, cornerGap: gap } };
}

export function offset(s, d, { tol = 0.01, maxCvs = 60 } = {}) {
  const s01 = S.reparam(s);
  const r = F.fitSurfaceTo((u, v) => V.add(S.point(s01, u, v), V.mul(S.normal(s01, u, v), d)), Math.max(3, s.degU), Math.max(3, s.degV), tol, { maxCvs, start: [s.cvs.length, s.cvs[0].length] });
  return { surface: r.surface, report: { exact: false, deviation: r.deviation, tol } };
}

export function rebuild(s, nu, nv, du = 3, dv = 3) {
  const s01 = S.reparam(s), Su = Math.max(nu * 4, 16), Sv = Math.max(nv * 4, 16);
  const tu = Array.from({ length: Su + 1 }, (_, i) => i / Su), tv = Array.from({ length: Sv + 1 }, (_, j) => j / Sv);
  const out = F.fitSurfaceGrid(tu.map((u) => tv.map((v) => S.point(s01, u, v))), tu, tv, Math.min(du, nu - 1), Math.min(dv, nv - 1), B.clampedUniform(nu, Math.min(du, nu - 1)), B.clampedUniform(nv, Math.min(dv, nv - 1)));
  let dev = 0;
  for (let i = 0; i < Su; i++) for (let j = 0; j < Sv; j++) { const u = (i + 0.5) / Su, v = (j + 0.5) / Sv; dev = Math.max(dev, V.dist(S.point(out, u, v), S.point(s01, u, v))); }
  return { surface: out, report: { exact: false, deviation: dev } };
}

/** Carry the surface on past an edge by `length` (mm, measured along the cross direction at the edge's middle). Exact: the
 *  last knot span's own polynomial is continued (de Casteljau past the end), so the joint is as smooth as the span itself;
 *  the joint knot is left at multiplicity p (the shape there is the same polynomial on both sides). */
export function extend(s, edge, length) {
  if (!(length > 0)) fail('nurbs.err.length');
  S.edgeParams(s, edge);
  const dir = edge[0], flip = edge[1] === '0';
  let w = s;
  if (dir === 'v') w = S.swap(w);
  if (flip) w = S.reverse(w, 'u');
  w = S.reparam(w);
  const p = w.degU, ks = B.distinctKnots(w.knotsU), last = ks[ks.length - 2].u;
  const piece = last > 0 ? S.split(w, 'u', last)[1] : w;
  const Hp = S.hgrid(piece);
  const build = (du) => {
    const t = (1 + du - last) / (1 - last), sCut = (1 - last) / (1 + du - last);
    const extRows = Hp[0].map((_, j) => {
      // Bézier points of the span's polynomial on [last, 1 + du]: the left chain of de Casteljau at t > 1 …
      let q = Hp.map((r) => r[j].slice());
      const left = [q[0]];
      for (let r = 1; r <= p; r++) { const nx = []; for (let i = 0; i + 1 < q.length; i++) nx.push(V.comb(q[i], 1 - t, q[i + 1], t)); q = nx; left.push(q[0]); }
      // … split at sCut (the old end, u = 1): the right part is the extension on [1, 1 + du]
      let z = left;
      const right = [];
      for (let r = 1; r <= p; r++) { right.push(z[z.length - 1]); const nx = []; for (let i = 0; i + 1 < z.length; i++) nx.push(V.comb(z[i], 1 - sCut, z[i + 1], sCut)); z = nx; }
      right.push(z[0]);
      return right.reverse();
    });
    const knotsU = w.knotsU.slice(0, w.knotsU.length - 1);
    for (let i = 0; i <= p; i++) knotsU.push(1 + du);
    const net = S.hgrid(w);
    for (let i = 1; i <= p; i++) net.push(extRows.map((row) => row[i]));
    const ext = S.fromHGrid(p, w.degV, [...new Array(p + 1).fill(0), ...new Array(p + 1).fill(1)], w.knotsV, extRows[0].map((_, i) => extRows.map((row) => row[i])));
    return { out: S.fromHGrid(p, w.degV, knotsU.map((x) => x / (1 + du)), w.knotsV, net), ext };
  };
  // the added length is measured along the middle iso-curve across the edge; secant steps on du until it is `length`
  const lenOf = (du) => C.arcLength(S.isoCurve(build(du).ext, 'v', 0.5));
  let d0 = length / (V.len(S.derivs(w, 1, 0.5, 1)[1][0]) || 1), l0 = lenOf(d0), d1 = d0 * length / Math.max(l0, 1e-300), l1 = lenOf(d1);
  for (let it = 0; it < 30 && Math.abs(l1 - length) > 1e-10 * length; it++) {
    const dn = l1 === l0 ? d1 * length / l1 : d1 + (length - l1) * (d1 - d0) / (l1 - l0);
    d0 = d1; l0 = l1; d1 = dn > 0 ? dn : d1 / 2; l1 = lenOf(d1);
  }
  let out = build(d1).out;
  if (flip) out = S.reverse(out, 'u');
  if (dir === 'v') out = S.swap(out);
  return { surface: out, report: { exact: true, length: l1, joint: flip ? d1 / (1 + d1) : 1 / (1 + d1) } };
}