Makes a curve that passes through every given point, with optional end directions.
How the precision is kept
Parameters come from chord lengths, knots from averaging them, and the control points from one dense linear system solved by Gaussian elimination with partial pivoting in 64-bit. The curve passes through the points up to the rounding of that solve.
Calculations under this policy
5 calculations in 2 source files. The code is read from the released app when this page loads. It is the calculation people run, not a copy kept for this page.
A parameter between 0 and 1 for each point to fit.
dᵢ = dᵢ₋₁ + |Qᵢ − Qᵢ₋₁| (or its square root, or 1); tᵢ = dᵢ ÷ d(last)
A running sum of distances divided by the total. If the total is zero the parameters are evenly spaced.
nurbs/fit.js
lines 32 to 43
ModelCraft 0.8.1
export function params(Q, kind = 'chord') {
const n = Q.length;
if (n < 2) return [0];
const d = [0];
for (let i = 1; i < n; i++) {
const l = V.dist(Q[i], Q[i - 1]);
d.push(d[i - 1] + (kind === 'uniform' ? 1 : kind === 'centripetal' ? Math.sqrt(l) : l));
}
const total = d[n - 1];
if (!(total > 0)) return Q.map((_, i) => i / (n - 1));
return d.map((x) => x / total);
}
Knots by averaging
64-bit
A clamped knot vector for a given number of control points from the points' parameters.
equal counts: U(j + p) = (t(j) + … + t(j + p − 1)) ÷ p; more points than control points: U = (1 − α)·t(i − 1) + α·t(i), i = floor(j·d), α = j·d − i
Each interior knot is a mean of p parameters, or a blend of two neighbors when there are more points than control points. End knots are 0 and 1 exactly.
nurbs/fit.js
lines 45 to 57
ModelCraft 0.8.1
export function averagedKnots(t, p, count) {
const U = new Array(p + 1).fill(0);
const inner = count - p - 1;
// A9.1 / eq. 9.8 generalised to a parameter list whose length may exceed `count` (approximation) or equal it
if (t.length === count) {
for (let j = 1; j <= inner; j++) { let s = 0; for (let i = j; i < j + p; i++) s += t[i]; U.push(s / p); }
} else {
const d = t.length / (count - p);
for (let j = 1; j <= inner; j++) { const i = Math.floor(j * d), a = j * d - i; U.push((1 - a) * t[i - 1] + a * t[i]); }
}
for (let i = 0; i <= p; i++) U.push(1);
return U;
}
Interpolating curve
64-bit
The curve of a given degree through every point.
Σⱼ Nⱼ(tᵢ)·Pⱼ = Qᵢ for every point; Σⱼ N′ⱼ(t)·Pⱼ = D at an end with a given derivative
One row of basis values per point, and one row of basis derivatives per end derivative, solved for the control points by vec.solve. The degree is lowered when there are too few points.
nurbs/fit.js
lines 66 to 84
ModelCraft 0.8.1
export function interpolate(Q, p, { kind = 'chord', params: t0 = null, d0 = null, d1 = null } = {}) {
const m = Q.length;
if (m < 2) fail('nurbs.err.cv_count', { min: 2 });
const t = t0 ? t0.slice() : params(Q, kind);
const count = m + (d0 ? 1 : 0) + (d1 ? 1 : 0);
const deg = Math.min(p, count - 1);
const tk = [...(d0 ? [t[0]] : []), ...t, ...(d1 ? [t[m - 1]] : [])];
const U = averagedKnots(tk, deg, count), n = count - 1;
const A = [], B = [];
Q.forEach((q, i) => {
A.push(rowAt(deg, U, n, t[i]));
B.push([q[0], q[1], q.length > 2 ? q[2] : 0]);
if (i === 0 && d0) { A.push(rowAt(deg, U, n, t[0], 1)); B.push(d0.slice(0, 3)); }
if (i === m - 1 && d1) { A.push(rowAt(deg, U, n, t[m - 1], 1)); B.push(d1.slice(0, 3)); }
});
const P = V.solve(A, B);
if (!P) fail('nurbs.err.fit_singular');
return { deg, knots: U, cvs: P, w: null };
}
Interpolating many lists at once
64-bit
Control points for several point lists that share parameters, used for the columns of a loft or a curve network.
The same matrix as for one curve with all lists as right-hand sides, in any dimension, so homogeneous points of rational sections are interpolated as they are.
nurbs/fit.js
lines 157 to 171
ModelCraft 0.8.1
export function interpolateMany(lists, t, p, { d0 = null, d1 = null } = {}) {
const m = t.length, count = m + (d0 ? 1 : 0) + (d1 ? 1 : 0), deg = Math.min(p, count - 1);
const tk = [...(d0 ? [t[0]] : []), ...t, ...(d1 ? [t[m - 1]] : [])];
const U = averagedKnots(tk, deg, count), n = count - 1, dim = lists[0][0].length;
const A = [], B = [];
for (let i = 0; i < m; i++) {
A.push(rowAt(deg, U, n, t[i])); B.push(lists.flatMap((L) => L[i]));
if (i === 0 && d0) { A.push(rowAt(deg, U, n, t[0], 1)); B.push(lists.flatMap((_, k) => d0[k])); }
if (i === m - 1 && d1) { A.push(rowAt(deg, U, n, t[m - 1], 1)); B.push(lists.flatMap((_, k) => d1[k])); }
}
const X = V.solve(A, B);
if (!X) fail('nurbs.err.fit_singular');
const rows = lists.map((_, k) => X.map((r) => r.slice(k * dim, k * dim + dim)));
return { deg, knots: U, rows };
}
Gaussian elimination with partial pivoting, then back-substitution
Works on a copy in plain arrays. The largest entry of each column is chosen as the pivot; if it is smaller than 1e-14 in absolute value the system is reported singular. The threshold does not scale with the size of the matrix entries.
nurbs/vec.js
lines 37 to 63
ModelCraft 0.8.1
/** Solve A x = B (A n×n, B n×k, both arrays of rows) by Gaussian elimination with partial pivoting; null when singular. */
export function solve(A, B) {
const n = A.length, k = B[0].length;
const M = A.map((r, i) => r.slice().concat(B[i].slice()));
for (let c = 0; c < n; c++) {
let piv = c, best = Math.abs(M[c][c]);
for (let r = c + 1; r < n; r++) { const v = Math.abs(M[r][c]); if (v > best) { best = v; piv = r; } }
if (best < 1e-14) return null;
if (piv !== c) { const t = M[c]; M[c] = M[piv]; M[piv] = t; }
const row = M[c], inv = 1 / row[c];
for (let r = c + 1; r < n; r++) {
const f = M[r][c] * inv;
if (f === 0) continue;
const rr = M[r];
for (let j = c; j < n + k; j++) rr[j] -= f * row[j];
}
}
const X = Array.from({ length: n }, () => new Array(k).fill(0));
for (let r = n - 1; r >= 0; r--) {
for (let j = 0; j < k; j++) {
let s = M[r][n + j];
for (let c = r + 1; c < n; c++) s -= M[r][c] * X[c][j];
X[r][j] = s / M[r][r];
}
}
return X;
}
In detail
Method
The NURBS Book § 9.2: chord-length or centripetal parameters (eq. 9.5, 9.6), averaged knots (eq. 9.8), global interpolation A9.1 with end derivatives (§ 9.2.2), as cited in the file's header.
Number types
64-bit numbers in plain arrays.
Tolerances
The solve reports a singular system when a pivot is smaller than 1e-14 in absolute value.
Limits
The pivot threshold is absolute, not scaled to the matrix. There is no iterative refinement and no estimate of the system's conditioning.
Approach: 64-bit, direct
Every value is a 64-bit number from input to result. Formulas are closed form or a direct solve, and every comparison uses a stated tolerance.