NURBS curves and surfaces: points, derivatives, normals and curvature

export function point(c, u) {
  u = clampU(c, u);
  const n = c.cvs.length - 1, p = c.deg, k = findSpan(n, p, u, c.knots), N = basisFuns(k, u, p, c.knots);
  let x = 0, y = 0, z = 0, w = 0;
  for (let j = 0; j <= p; j++) {
    const i = k - p + j, q = c.cvs[i], wi = c.w ? c.w[i] : 1, f = N[j] * wi;
    x += f * q[0]; y += f * q[1]; z += f * q[2]; w += f;
  }
  return [x / w, y / w, z / w];
}

/** [C(u), C'(u), …, C^(d)(u)]; derivatives above the degree are zero for a polynomial curve. */
export function derivs(c, u, d = 2) {
  u = clampU(c, u);
  const n = c.cvs.length - 1, p = c.deg, k = findSpan(n, p, u, c.knots);
  const du = Math.min(d, p), nd = dersBasisFuns(k, u, p, du, c.knots);
  const A = [], W = [];
  for (let r = 0; r <= d; r++) {
    let x = 0, y = 0, z = 0, w = 0;
    if (r <= du) for (let j = 0; j <= p; j++) {
      const i = k - p + j, q = c.cvs[i], wi = c.w ? c.w[i] : 1, f = nd[r][j] * wi;
      x += f * q[0]; y += f * q[1]; z += f * q[2]; w += f;
    }
    A.push([x, y, z]); W.push(w);
  }
  const CK = [];
  for (let r = 0; r <= d; r++) {
    let v = A[r].slice();
    for (let i = 1; i <= r; i++) { const f = binomial(r, i) * W[i]; v = [v[0] - f * CK[r - i][0], v[1] - f * CK[r - i][1], v[2] - f * CK[r - i][2]]; }
    CK.push(V.mul(v, 1 / W[0]));
  }
  return CK;
}

export function curvature(c, u) {
  const [, d1, d2] = derivs(c, u, 2), l = V.len(d1);
  return l < 1e-300 ? 0 : V.len(V.cross(d1, d2)) / (l * l * l);
}

/** κN: the curvature vector (points to the centre of curvature, length = curvature). */
export function curvatureVector(c, u) {
  const [, d1, d2] = derivs(c, u, 2), l2 = V.dot(d1, d1);
  if (l2 < 1e-300) return [0, 0, 0];
  const t = V.mul(d1, 1 / Math.sqrt(l2));
  return V.mul(V.sub(d2, V.mul(t, V.dot(d2, t))), 1 / l2);
}

export function torsion(c, u) {
  const [, d1, d2, d3] = derivs(c, u, 3), x = V.cross(d1, d2), l = V.dot(x, x);
  return l < 1e-300 ? 0 : V.dot(x, d3) / l;
}

export function fromH(deg, knots, Pw) {
  const rational = Pw.some((q) => Math.abs(q[3] - 1) > 1e-14);
  return {
    deg, knots: knots.slice(),
    cvs: Pw.map((q) => [q[0] / q[3], q[1] / q[3], q[2] / q[3]]),
    w: rational ? Pw.map((q) => q[3]) : null,
  };
}

/** A7.1: a circular arc from angle a0 to a1 (radians, counter-clockwise about xAxis × yAxis), as a rational quadratic. */
export function arc(center, xAxis, yAxis, r, a0, a1) {
  let th = a1 - a0;
  if (!(Math.abs(th) > 1e-12)) fail('nurbs.err.arc');
  if (Math.abs(th) > 2 * Math.PI) th = Math.sign(th) * 2 * Math.PI;
  const X = V.norm(xAxis), Y = V.norm(yAxis);
  const narcs = Math.abs(th) <= Math.PI / 2 ? 1 : Math.abs(th) <= Math.PI ? 2 : Math.abs(th) <= 1.5 * Math.PI ? 3 : 4;
  const dth = th / narcs, w1 = Math.cos(dth / 2);
  const at = (ang) => V.add(center, V.add(V.mul(X, r * Math.cos(ang)), V.mul(Y, r * Math.sin(ang))));
  const tan = (ang) => V.add(V.mul(X, -Math.sin(ang)), V.mul(Y, Math.cos(ang)));
  const cvs = [at(a0)], w = [1];
  let ang = a0, P0 = at(a0), T0 = tan(a0);
  for (let i = 1; i <= narcs; i++) {
    ang += dth;
    const P2 = at(ang), T2 = tan(ang);
    // the middle point: where the end tangents meet (they meet at distance r·tan(dth/2) from each end)
    const P1 = V.add(P0, V.mul(T0, r * Math.tan(dth / 2)));
    cvs.push(P1, P2); w.push(w1, 1);
    P0 = P2; T0 = T2;
  }
  const knots = [0, 0, 0];
  for (let i = 1; i < narcs; i++) knots.push(i / narcs, i / narcs);
  knots.push(1, 1, 1);
  return { deg: 2, knots, cvs, w };
}

export function polyline(pts) {
  if (pts.length < 2) fail('nurbs.err.cv_count', { min: 2 });
  const L = [0];
  for (let i = 1; i < pts.length; i++) L.push(L[i - 1] + Math.max(1e-12, V.dist(pts[i], pts[i - 1])));
  const knots = [0, ...L.map((x) => x / L[L.length - 1]), 1];
  return { deg: 1, knots, cvs: pts.map((p) => [p[0], p[1], p.length > 2 ? p[2] : 0]), w: null };
}

export function transform(c, m) {
  return { deg: c.deg, knots: c.knots.slice(), cvs: c.cvs.map((q) => V.xform(m, q)), w: c.w ? c.w.slice() : null };
}

export function point(s, u, v) {
  u = clamp(u, domainU(s)); v = clamp(v, domainV(s));
  const nu = s.cvs.length - 1, nv = s.cvs[0].length - 1;
  const ku = findSpan(nu, s.degU, u, s.knotsU), kv = findSpan(nv, s.degV, v, s.knotsV);
  const Nu = basisFuns(ku, u, s.degU, s.knotsU), Nv = basisFuns(kv, v, s.degV, s.knotsV);
  let x = 0, y = 0, z = 0, w = 0;
  for (let a = 0; a <= s.degU; a++) {
    const i = ku - s.degU + a;
    for (let b = 0; b <= s.degV; b++) {
      const j = kv - s.degV + b, q = s.cvs[i][j], wij = s.w ? s.w[i][j] : 1, f = Nu[a] * Nv[b] * wij;
      x += f * q[0]; y += f * q[1]; z += f * q[2]; w += f;
    }
  }
  return [x / w, y / w, z / w];
}

export function derivs(s, u, v, d = 2) {
  u = clamp(u, domainU(s)); v = clamp(v, domainV(s));
  const pu = s.degU, pv = s.degV, nu = s.cvs.length - 1, nv = s.cvs[0].length - 1;
  const du = Math.min(d, pu), dv = Math.min(d, pv);
  const ku = findSpan(nu, pu, u, s.knotsU), kv = findSpan(nv, pv, v, s.knotsV);
  const Nu = dersBasisFuns(ku, u, pu, du, s.knotsU), Nv = dersBasisFuns(kv, v, pv, dv, s.knotsV);
  const A = [], W = [];
  for (let k = 0; k <= d; k++) {
    A.push([]); W.push([]);
    for (let l = 0; l <= d - k; l++) {
      let x = 0, y = 0, z = 0, w = 0;
      if (k <= du && l <= dv) for (let a = 0; a <= pu; a++) {
        const i = ku - pu + a;
        for (let b = 0; b <= pv; b++) {
          const j = kv - pv + b, q = s.cvs[i][j], wij = s.w ? s.w[i][j] : 1, f = Nu[k][a] * Nv[l][b] * wij;
          x += f * q[0]; y += f * q[1]; z += f * q[2]; w += f;
        }
      }
      A[k].push([x, y, z]); W[k].push(w);
    }
  }
  const S = [];
  for (let k = 0; k <= d; k++) {
    S.push([]);
    for (let l = 0; l <= d - k; l++) {
      let p = A[k][l].slice();
      for (let j = 1; j <= l; j++) p = V.sub(p, V.mul(S[k][l - j], binomial(l, j) * W[0][j]));
      for (let i = 1; i <= k; i++) {
        p = V.sub(p, V.mul(S[k - i][l], binomial(k, i) * W[i][0]));
        let p2 = [0, 0, 0];
        for (let j = 1; j <= l; j++) p2 = V.add(p2, V.mul(S[k - i][l - j], binomial(l, j) * W[i][j]));
        p = V.sub(p, V.mul(p2, binomial(k, i)));
      }
      S[k].push(V.mul(p, 1 / W[0][0]));
    }
  }
  return S;
}

/** Unit normal Su × Sv. Where it vanishes (a pole or a collapsed edge) the normal is the limit from inside the domain. */
export function normal(s, u, v) {
  const D = derivs(s, u, v, 1), l1 = V.len(D[1][0]), l2 = V.len(D[0][1]);
  let n = V.cross(D[1][0], D[0][1]);
  // a derivative a billion times shorter than the other is rounding noise on a collapsed edge, not a direction
  const degenerate = l1 <= 1e-9 * l2 || l2 <= 1e-9 * l1;
  if (!degenerate && V.len(n) > 1e-10 * l1 * l2) return V.norm(n);
  const [ua, ub] = domainU(s), [va, vb] = domainV(s);
  // a collapsed edge: near v = v0 with Sv → 0, Sv ≈ (v − v0)·Svu… — the limit of Su × Sv is ±Su × Suv (and Suv × Sv when Su
  // vanishes), the sign set by which side of the domain the point is on
  const D2 = derivs(s, u, v, 2), lu = V.len(D2[1][0]), lv = V.len(D2[0][1]);
  if (lv <= 1e-9 * Math.max(lu, 1e-300) && lu > 1e-300) {
    const m = V.cross(D2[1][0], D2[1][1]), sg = u - ua <= ub - u ? 1 : -1;
    if (V.len(m) > 1e-12 * lu * V.len(D2[1][1])) return V.norm(V.mul(m, sg));
  }
  if (lu <= 1e-9 * Math.max(lv, 1e-300) && lv > 1e-300) {
    const m = V.cross(D2[1][1], D2[0][1]), sg = v - va <= vb - v ? 1 : -1;
    if (V.len(m) > 1e-12 * lv * V.len(D2[1][1])) return V.norm(V.mul(m, sg));
  }
  for (const h of [1e-7, 1e-5, 1e-3]) {
    const uu = u + (u - ua < (ub - ua) / 2 ? 1 : -1) * h * (ub - ua), vv = v + (v - va < (vb - va) / 2 ? 1 : -1) * h * (vb - va);
    const E = derivs(s, uu, vv, 1);
    n = V.cross(E[1][0], E[0][1]);
    if (V.len(n) > 1e-300) return V.norm(n);
  }
  return [0, 0, 1];
}

/** On a collapsed edge (one first derivative vanishes) the curvatures are the limits from just inside the domain. */
function inside(s, u, v) {
  const D = derivs(s, u, v, 1), l1 = V.len(D[1][0]), l2 = V.len(D[0][1]);
  if (!(l1 <= 1e-6 * l2 || l2 <= 1e-6 * l1)) return [u, v];
  const [ua, ub] = domainU(s), [va, vb] = domainV(s), h = 1e-5;
  const uu = l1 <= 1e-6 * l2 ? u : u + (u - ua <= ub - u ? 1 : -1) * h * (ub - ua);
  const vv = l2 <= 1e-6 * l1 ? v : v + (v - va <= vb - v ? 1 : -1) * h * (vb - va);
  return [uu, vv];
}

export function curvature(s, u, v) {
  const f = forms(s, u, v), det = f.E * f.G - f.F * f.F;
  if (!(det > 1e-300)) return { K: 0, H: 0, k1: 0, k2: 0, d1: f.Su, d2: f.Sv };
  const K = (f.L * f.N - f.M * f.M) / det;
  const H = (f.E * f.N - 2 * f.F * f.M + f.G * f.L) / (2 * det);
  const disc = Math.sqrt(Math.max(0, H * H - K));
  const k1 = H + disc, k2 = H - disc;
  // principal direction for k: (L − kE) a + (M − kF) b = 0  (or the second row when the first vanishes)
  const dirFor = (k) => {
    let a = -(f.M - k * f.F), b = f.L - k * f.E;
    if (Math.abs(a) + Math.abs(b) < 1e-12 * (Math.abs(f.L) + Math.abs(f.M) + Math.abs(k) * (f.E + f.G) + 1e-300)) { a = -(f.N - k * f.G); b = f.M - k * f.F; }
    if (Math.abs(a) + Math.abs(b) < 1e-300) return V.norm(f.Su);
    return V.norm(V.add(V.mul(f.Su, a), V.mul(f.Sv, b)));
  };
  let d1 = dirFor(k1), d2 = dirFor(k2);
  if (disc < 1e-12 * (Math.abs(H) + 1e-300)) { d1 = V.norm(f.Su); d2 = V.norm(V.cross(f.n, d1)); }
  return { K, H, k1, k2, d1, d2 };
}

export function normalCurvature(s, u, v, t) {
  const f = forms(s, u, v), det = f.E * f.G - f.F * f.F;
  if (!(det > 1e-300)) return 0;
  const p = V.dot(t, f.Su), q = V.dot(t, f.Sv);
  const a = (f.G * p - f.F * q) / det, b = (f.E * q - f.F * p) / det;
  const I = f.E * a * a + 2 * f.F * a * b + f.G * b * b;
  return I > 1e-300 ? (f.L * a * a + 2 * f.M * a * b + f.N * b * b) / I : 0;
}

export function isoCurve(s, dir, t) {
  const Pw = hgrid(s);
  if (dir === 'u') {
    // u fixed: blend the rows i by N_i(u); the result runs in v
    const nu = s.cvs.length - 1, u = clamp(t, domainU(s)), k = findSpan(nu, s.degU, u, s.knotsU), N = basisFuns(k, u, s.degU, s.knotsU);
    const row = Pw[0].map((_, j) => { let q = [0, 0, 0, 0]; for (let a = 0; a <= s.degU; a++) q = V.add(q, V.mul(Pw[k - s.degU + a][j], N[a])); return q; });
    return C.fromH(s.degV, s.knotsV, row);
  }
  const nv = s.cvs[0].length - 1, v = clamp(t, domainV(s)), k = findSpan(nv, s.degV, v, s.knotsV), N = basisFuns(k, v, s.degV, s.knotsV);
  const col = Pw.map((r) => { let q = [0, 0, 0, 0]; for (let b = 0; b <= s.degV; b++) q = V.add(q, V.mul(r[k - s.degV + b], N[b])); return q; });
  return C.fromH(s.degU, s.knotsU, col);
}

/** Angle between two vectors in radians (0…π); 0 when either is zero. */
export function angle(a, b) {
  const la = len(a), lb = len(b);
  if (la < 1e-300 || lb < 1e-300) return 0;
  // atan2 of |a×b| and a·b stays accurate for tiny angles (acos does not)
  const c = a.length === 3 ? len(cross(a, b)) : Math.sqrt(Math.max(0, la * la * lb * lb - dot(a, b) ** 2));
  return Math.atan2(c, dot(a, b));
}

/** Apply a 4×4 affine matrix (row-major, 16 numbers; the last row is ignored) to a point or (w = 0) a direction. */
export const xform = (m, p, w = 1) => [
  m[0] * p[0] + m[1] * p[1] + m[2] * p[2] + m[3] * w,
  m[4] * p[0] + m[5] * p[1] + m[6] * p[2] + m[7] * w,
  m[8] * p[0] + m[9] * p[1] + m[10] * p[2] + m[11] * w,
];

export function maps(s, { tol = 0.01, maxPerSpan = 64 } = {}) {
  const t = S.tessellate(s, { tol, maxPerSpan, curvature: true });
  const kmax = new Float64Array(t.K.length);
  for (let i = 0; i < kmax.length; i++) { const H = t.H[i], K = t.K[i], d = Math.sqrt(Math.max(0, H * H - K)); kmax[i] = Math.max(Math.abs(H + d), Math.abs(H - d)); }
  return { ...t, kmax, rangeK: ranges(t.K), rangeH: ranges(t.H), rangeKmax: ranges(kmax) };
}

export function comb(c, { samples = 64, scale = 1 } = {}) {
  const [a, b] = C.domain(c), teeth = [];
  let max = 0;
  for (let i = 0; i <= samples; i++) {
    const u = a + (b - a) * i / samples, p = C.point(c, u), kv = C.curvatureVector(c, u), k = V.len(kv);
    teeth.push({ u, p, q: V.sub(p, V.mul(kv, scale)), k });
    if (k > max) max = k;
  }
  return { teeth, max };
}

let max = 0;
for (const o of this.surfaces()) for (const k of o.geometry.userData.K || []) if (Number.isFinite(k)) max = Math.max(max, Math.abs(k));
this.kRange = { min: -max, max };