NURBS knot insertion and removal, degree change, split and join

export function insertKnot(c, u, r = 1) {
  const p = c.deg, U = c.knots, n = c.cvs.length - 1;
  const [a, b] = domain(c);
  // inside the domain by more than rounding (a knot within 1e-12 of an end would be read as that end's knot)
  if (!(u > a + 1e-12 * (b - a) && u < b - 1e-12 * (b - a))) fail('nurbs.err.knot_outside');
  const s = multiplicity(U, u);
  if (r <= 0) return { deg: p, knots: U.slice(), cvs: c.cvs.map((q) => q.slice()), w: c.w ? c.w.slice() : null };
  if (r + s > p) fail('nurbs.err.knot_mult');
  const k = findSpan(n, p, u, U), Pw = hcvs(c), mp = n + p + 1, nq = n + r;
  const UQ = new Array(mp + r + 1), Qw = new Array(nq + 1);
  for (let i = 0; i <= k; i++) UQ[i] = U[i];
  for (let i = 1; i <= r; i++) UQ[k + i] = u;
  for (let i = k + 1; i <= mp; i++) UQ[i + r] = U[i];
  for (let i = 0; i <= k - p; i++) Qw[i] = Pw[i];
  for (let i = k - s; i <= n; i++) Qw[i + r] = Pw[i];
  const Rw = [];
  for (let i = 0; i <= p - s; i++) Rw[i] = Pw[k - p + i].slice();
  let L = 0;
  for (let j = 1; j <= r; j++) {
    L = k - p + j;
    for (let i = 0; i <= p - j - s; i++) {
      const al = (u - U[L + i]) / (U[i + k + 1] - U[L + i]);
      Rw[i] = V.comb(Rw[i + 1], al, Rw[i], 1 - al);
    }
    Qw[L] = Rw[0];
    Qw[k + r - j - s] = Rw[p - j - s];
  }
  for (let i = L + 1; i < k - s; i++) Qw[i] = Rw[i - L];
  return fromH(p, UQ, Qw);
}

/** Homogeneous-space tolerance for a geometric tolerance (The NURBS Book § 5.4, eq. 5.30). */
function hTol(c, tol) {
  if (!c.w) return tol;
  const wmin = Math.min(...c.w);
  let pmax = 0;
  for (const q of c.cvs) pmax = Math.max(pmax, V.len(q));
  return tol * wmin / (1 + pmax);
}

/** Remove up to `num` copies of the interior knot u while the curve moves no more than tol. */
export function removeKnot(c, u, num = 1, tol = 1e-9) {
  const p = c.deg, n = c.cvs.length - 1, ord = p + 1;
  const U = c.knots.slice(), Pw = hcvs(c), m = U.length - 1;
  const s = multiplicity(U, u);
  const [a, b] = domain(c);
  if (s === 0 || !(u > a && u < b)) return { curve: c, removed: 0 };
  let r = -1;
  for (let i = 0; i <= m; i++) if (Math.abs(U[i] - u) <= 1e-12 * (b - a)) r = i;
  const TOL = hTol(c, tol);
  const fout = (2 * r - s - p) >> 1;
  let first = r - p, last = r - s, t = 0;
  const temp = [];
  for (t = 0; t < Math.min(num, s); t++) {
    const off = first - 1;
    temp[0] = Pw[off]; temp[last + 1 - off] = Pw[last + 1];
    let i = first, j = last, ii = 1, jj = last - off, remflag = false;
    while (j - i > t) {
      const alfi = (u - U[i]) / (U[i + ord + t] - U[i]);
      const alfj = (u - U[j - t]) / (U[j + ord] - U[j - t]);
      temp[ii] = V.mul(V.sub(Pw[i], V.mul(temp[ii - 1], 1 - alfi)), 1 / alfi);
      temp[jj] = V.mul(V.sub(Pw[j], V.mul(temp[jj + 1], alfj)), 1 / (1 - alfj));
      i++; ii++; j--; jj--;
    }
    if (j - i < t) { if (V.dist(temp[ii - 1], temp[jj + 1]) <= TOL) remflag = true; }
    else {
      const alfi = (u - U[i]) / (U[i + ord + t] - U[i]);
      if (V.dist(Pw[i], V.comb(temp[ii + t + 1], alfi, temp[ii - 1], 1 - alfi)) <= TOL) remflag = true;
    }
    if (!remflag) break;
    i = first; j = last;
    while (j - i > t) { Pw[i] = temp[i - off]; Pw[j] = temp[j - off]; i++; j--; }
    first--; last++;
  }
  if (t === 0) return { curve: c, removed: 0 };
  for (let k = r + 1; k <= m; k++) U[k - t] = U[k];
  U.length = m + 1 - t;
  let j = fout, i = j;
  for (let k = 1; k < t; k++) { if (k % 2 === 1) i++; else j--; }
  for (let k = i + 1; k <= n; k++) { Pw[j] = Pw[k]; j++; }
  Pw.length = n + 1 - t;
  return { curve: fromH(p, U, Pw), removed: t };
}

/** The curve as its Bézier pieces (each a clamped curve of one span on that span's own parameter interval). */
export function decompose(c) {
  const p = c.deg;
  let full = c;
  for (const k of distinctKnots(c.knots).slice(1, -1)) if (k.m < p) full = insertKnot(full, k.u, p - k.m);
  const Pw = hcvs(full), ks = distinctKnots(full.knots), out = [];
  for (let s = 0; s + 1 < ks.length; s++) {
    const a = ks[s].u, b = ks[s + 1].u;
    const seg = Pw.slice(s * p, s * p + p + 1);
    out.push(fromH(p, [...new Array(p + 1).fill(a), ...new Array(p + 1).fill(b)], seg));
  }
  return out;
}

export function elevate(c, t = 1, { removeKnots = true } = {}) {
  if (!Number.isInteger(t) || t < 0) fail('nurbs.err.degree');
  if (t === 0) return reparam(c, ...domain(c));
  const p = c.deg, ph = p + t;
  if (ph > MAX_DEG) fail('nurbs.err.degree', { max: MAX_DEG });
  const orig = distinctKnots(c.knots);
  const pieces = decompose(c).map((b) => {
    const P = hcvs(b), Q = [];
    for (let i = 0; i <= ph; i++) {
      let q = [0, 0, 0, 0];
      for (let j = Math.max(0, i - t); j <= Math.min(p, i); j++) q = V.add(q, V.mul(P[j], binomial(p, j) * binomial(t, i - j) / binomial(ph, i)));
      Q.push(q);
    }
    return Q;
  });
  let out = assemble(ph, pieces, orig.map((k) => k.u));
  if (!removeKnots) return out;
  const tol = 1e-10 * curveSize(c);
  for (const k of orig.slice(1, -1)) out = removeKnot(out, k.u, p - k.m, tol).curve;
  return out;
}

/** One degree lower, if the curve moves no more than tol; the error is a bound (control points of the re-elevated pieces). */
export function reduce(c, tol = 1e-6) {
  const p = c.deg;
  if (p < 2) return null;
  const orig = distinctKnots(c.knots);
  let err = 0;
  const pieces = decompose(c).map((b) => {
    const P = hcvs(b), q = p - 1, L = new Array(p), R = new Array(p);
    L[0] = P[0];
    for (let i = 1; i < p; i++) { const al = i / p; L[i] = V.mul(V.sub(P[i], V.mul(L[i - 1], al)), 1 / (1 - al)); }
    R[q] = P[p];
    for (let i = q - 1; i >= 0; i--) { const al = (i + 1) / p; R[i] = V.mul(V.sub(P[i + 1], V.mul(R[i + 1], 1 - al)), 1 / al); }
    const Q = [];
    for (let i = 0; i <= q; i++) { const lam = q ? i / q : 0.5; Q.push(V.comb(L[i], 1 - lam, R[i], lam)); }
    Q[0] = P[0]; Q[q] = P[p];
    // bound: raise Q back to degree p and compare control points (convex hull: the curves differ by no more than that)
    for (let i = 0; i <= p; i++) {
      const a = i / p, e = i === 0 ? Q[0] : i === p ? Q[q] : V.comb(Q[i - 1], a, Q[i], 1 - a);
      err = Math.max(err, V.dist(e.slice(0, 3).map((x) => x / e[3]), P[i].slice(0, 3).map((x) => x / P[i][3])));
    }
    return Q;
  });
  if (err > tol) return null;
  let out = assemble(p - 1, pieces, orig.map((k) => k.u));
  const rest = Math.max(0, tol - err);
  for (const k of orig.slice(1, -1)) out = removeKnot(out, k.u, p - k.m, rest).curve;
  return { curve: out, error: err };
}

export function split(c, u) {
  const [a, b] = domain(c);
  if (!(u > a + 1e-12 * (b - a) && u < b - 1e-12 * (b - a))) fail('nurbs.err.knot_outside');
  const p = c.deg, s = multiplicity(c.knots, u);
  const full = s < p ? insertKnot(c, u, p - s) : c;
  const Pw = hcvs(full), U = full.knots;
  let k = -1;
  for (let i = 0; i < U.length; i++) if (Math.abs(U[i] - u) <= 1e-12 * (b - a)) { k = i; break; }
  // knots U[k … k+p-1] are u (p copies); the left curve ends with p+1 copies, the right one starts with them
  const left = fromH(p, [...U.slice(0, k + p), u], Pw.slice(0, k));
  const right = fromH(p, [u, ...U.slice(k)], Pw.slice(k - 1));
  return [left, right];
}

export function reparam(c, a = 0, b = 1) {
  const [ca, cb] = domain(c), f = (b - a) / (cb - ca);
  const knots = c.knots.map((x) => a + (x - ca) * f);
  for (let i = 0; i <= c.deg; i++) { knots[i] = a; knots[knots.length - 1 - i] = b; }
  return { deg: c.deg, knots, cvs: c.cvs.map((q) => q.slice()), w: c.w ? c.w.slice() : null };
}

export function reverse(c) {
  const [a, b] = domain(c);
  return { deg: c.deg, knots: c.knots.map((x) => a + b - x).reverse(), cvs: c.cvs.map((q) => q.slice()).reverse(), w: c.w ? c.w.slice().reverse() : null };
}

/** Join curves end to start (each next one turned round when its far end is the near one) into one curve, C0 at the joints. */
export function join(curves, tol = 1e-6) {
  if (!curves.length) fail('nurbs.err.cv_count', { min: 1 });
  const list = [curves[0]];
  for (let i = 1; i < curves.length; i++) {
    const prevEnd = endsOf(list[i - 1]).end, e = endsOf(curves[i]);
    if (V.dist(prevEnd, e.start) <= tol) list.push(curves[i]);
    else if (V.dist(prevEnd, e.end) <= tol) list.push(reverse(curves[i]));
    else if (i === 1 && V.dist(endsOf(list[0]).start, e.start) <= tol) { list[0] = reverse(list[0]); list.push(curves[i]); }
    else if (i === 1 && V.dist(endsOf(list[0]).start, e.end) <= tol) { list[0] = reverse(list[0]); list.push(reverse(curves[i])); }
    else fail('nurbs.err.join_gap', { index: i });
  }
  const p = Math.max(...list.map((c) => c.deg));
  const lens = list.map((c) => Math.max(1e-12, arcLength(c))), total = lens.reduce((x, y) => x + y, 0);
  let at = 0, knots = null, Pw = null;
  list.forEach((c0, i) => {
    let c = elevate(c0, p - c0.deg);
    c = reparam(c, at / total, (at + lens[i]) / total);
    at += lens[i];
    let H = hcvs(c);
    if (!Pw) { knots = c.knots.slice(); Pw = H; return; }
    // scale this piece's homogeneous points so its first weight equals the last weight so far (same curve, continuous weights)
    const f = Pw[Pw.length - 1][3] / H[0][3];
    H = H.map((q) => V.mul(q, f));
    Pw[Pw.length - 1] = V.mul(V.add(Pw[Pw.length - 1], H[0]), 0.5);
    Pw.push(...H.slice(1));
    knots = [...knots.slice(0, knots.length - 1), ...c.knots.slice(p + 1)];
  });
  knots[knots.length - 1] = 1;
  for (let i = 0; i <= p; i++) knots[knots.length - 1 - i] = 1;
  return fromH(p, knots, Pw);
}

function fromRows(s, dir, rows) {
  const H = rows.map((c) => C.hcvs(c)), K = rows[0].knots;
  for (const c of rows) if (c.knots.length !== K.length || c.knots.some((x, i) => Math.abs(x - K[i]) > 1e-12 * (1 + Math.abs(K[i])))) return null;
  if (dir === 'v') return fromHGrid(s.degU, rows[0].deg, s.knotsU, K, H);
  const Pw = H[0].map((_, i) => H.map((col) => col[i]));
  return fromHGrid(rows[0].deg, s.degV, K, s.knotsV, Pw);
}

export function elevate(s, dir, t = 1) {
  checkDir(dir);
  const rows = rowsOf(s, dir);
  return fromRows(s, dir, rows.map((c) => C.elevate(c, t))) || fromRows(s, dir, rows.map((c) => C.elevate(c, t, { removeKnots: false })));
}

function mapKnots(K, p, n, a, b) {
  const ka = K[p], kb = K[n], f = (b - a) / (kb - ka), out = K.map((x) => a + (x - ka) * f);
  for (let i = 0; i <= p; i++) { out[i] = a; out[out.length - 1 - i] = b; }
  return out;
}

function clamped(n, p) { const U = new Array(p + 1).fill(0); for (let i = 1; i < n - p; i++) U.push(i / (n - p)); for (let i = 0; i <= p; i++) U.push(1); return U; }

const at = (c, f) => { const [a, b] = C.domain(c); return a + (b - a) * f; };