Sketch splines: evaluation, fitting and exact pieces
Evaluates sketch splines and their derivatives, fits a spline through points, cuts an exact piece out of one, and gives the area a closed spline encloses.
How the precision is kept
Cubic B-splines with 64-bit knots and poles. Basis functions use the Cox–de Boor recurrence, which adds only non-negative terms. A fit spline is one dense linear solve with partial pivoting. A piece of a spline is cut by knot insertion, which keeps the shape up to rounding, and the enclosed area uses a 3-point Gauss rule that is exact for a cubic.
Calculations under this policy
13 calculations in 1 source file. 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 binary search that compares the parameter with the knots directly, with no arithmetic on them, so the decision is exact. The shift in (lo + hi + 1) >> 1 acts on array indexes.
sketch-spline.js
lines 27 to 34
PartCraft 0.8.1
/** Span k with U[k] ≤ u < U[k+1] among p ≤ k ≤ m − 1 (the last non-empty span at the end). */
export function findSpan(U, p, m, u) {
if (u >= U[m]) { let k = m - 1; while (k > p && U[k] >= U[k + 1]) k--; return k; }
if (u <= U[p]) { let k = p; while (k < m - 1 && U[k + 1] <= u) k++; return k; }
let lo = p, hi = m - 1;
while (lo < hi) { const mid = (lo + hi + 1) >> 1; if (U[mid] <= u) lo = mid; else hi = mid - 1; }
return lo;
}
Basis functions
64-bit
The values of the B-spline basis functions that are not zero at a parameter.
Within the span every term is a product of non-negative numbers, so nothing cancels. A zero denominator (a repeated knot) gives a zero term instead of a division by zero.
sketch-spline.js
lines 36 to 50
PartCraft 0.8.1
export function basisFuns(k, u, p, U) {
const N = [1], left = [], right = [];
for (let j = 1; j <= p; j++) {
left[j] = u - U[k + 1 - j]; right[j] = U[k + j] - u;
let saved = 0;
for (let r = 0; r < j; r++) {
const den = right[r + 1] + left[j - r];
const t = den === 0 ? 0 : N[r] / den;
N[r] = saved + right[r + 1] * t;
saved = left[j - r] * t;
}
N[j] = saved;
}
return N;
}
Basis function derivatives
64-bit
The basis functions and their derivatives up to a given order at a parameter.
The derivative coefficients are differences divided by knot differences, so unlike the values they can cancel; knot spans that are very uneven magnify this. Zero denominators are replaced by zero terms.
sketch-spline.js
lines 52 to 91
PartCraft 0.8.1
/** ders[i][j]: i-th derivative of basis function k − p + j at u (i ≤ n). */
export function dersBasisFuns(k, u, p, n, U) {
const ndu = [[1]], left = [], right = [];
for (let j = 1; j <= p; j++) {
ndu[j] = [];
left[j] = u - U[k + 1 - j]; right[j] = U[k + j] - u;
let saved = 0;
for (let r = 0; r < j; r++) {
ndu[j][r] = right[r + 1] + left[j - r];
const t = ndu[j][r] === 0 ? 0 : ndu[r][j - 1] / ndu[j][r];
ndu[r][j] = saved + right[r + 1] * t;
saved = left[j - r] * t;
}
ndu[j][j] = saved;
}
const ders = [];
for (let i = 0; i <= n; i++) ders.push(new Array(p + 1).fill(0));
for (let j = 0; j <= p; j++) ders[0][j] = ndu[j][p];
for (let r = 0; r <= p; r++) {
let s1 = 0, s2 = 1;
const a = [new Array(p + 1).fill(0), new Array(p + 1).fill(0)];
a[0][0] = 1;
for (let kk = 1; kk <= n; kk++) {
let d = 0;
const rk = r - kk, pk = p - kk;
if (r >= kk) { a[s2][0] = ndu[pk + 1][rk] === 0 ? 0 : a[s1][0] / ndu[pk + 1][rk]; d = a[s2][0] * ndu[rk][pk]; }
const j1 = rk >= -1 ? 1 : -rk, j2 = r - 1 <= pk ? kk - 1 : p - r;
for (let j = j1; j <= j2; j++) {
a[s2][j] = ndu[pk + 1][rk + j] === 0 ? 0 : (a[s1][j] - a[s1][j - 1]) / ndu[pk + 1][rk + j];
d += a[s2][j] * ndu[rk + j][pk];
}
if (r <= pk) { a[s2][kk] = ndu[pk + 1][r] === 0 ? 0 : -a[s1][kk - 1] / ndu[pk + 1][r]; d += a[s2][kk] * ndu[r][pk]; }
ders[kk][r] = d;
[s1, s2] = [s2, s1];
}
}
let f = p;
for (let kk = 1; kk <= n; kk++) { for (let j = 0; j <= p; j++) ders[kk][j] *= f; f *= p - kk; }
return ders;
}
Point on a spline
64-bit
The point of a spline at a parameter.
C(u) = Σⱼ Nⱼ(u) × Pⱼ (degree + 1 terms)
A sum of four weighted poles for a cubic, with weights that add up to 1, in 64-bit.
sketch-spline.js
lines 97 to 104
PartCraft 0.8.1
export function pointAt(c, u) {
const { deg: p, knots: U, poles: P } = c;
const k = findSpan(U, p, P.length, u);
const N = basisFuns(k, u, p, U);
let x = 0, y = 0;
for (let j = 0; j <= p; j++) { const q = P[k - p + j]; x += N[j] * q[0]; y += N[j] * q[1]; }
return [x, y];
}
Derivatives of a spline
64-bit
The point, tangent and higher derivatives of a spline at a parameter.
C⁽ⁱ⁾(u) = Σⱼ Nⱼ⁽ⁱ⁾(u) × Pⱼ
Derivative weights add up to zero, so for poles far from the origin the sum subtracts large numbers; the error is of the order of the rounding of the pole coordinates times the weights.
sketch-spline.js
lines 106 to 118
PartCraft 0.8.1
/** [C(u), C′(u), C″(u), …] up to the n-th derivative. */
export function derivsAt(c, u, n = 1) {
const { deg: p, knots: U, poles: P } = c;
const k = findSpan(U, p, P.length, u);
const D = dersBasisFuns(k, u, p, Math.min(n, p), U);
const out = [];
for (let i = 0; i <= n; i++) {
let x = 0, y = 0;
if (i <= p) for (let j = 0; j <= p; j++) { const q = P[k - p + j]; x += D[i][j] * q[0]; y += D[i][j] * q[1]; }
out.push([x, y]);
}
return out;
}
Chord-length parameters
64-bit
Gives each fit point a parameter equal to the summed straight distances up to it.
tᵢ = tᵢ₋₁ + |Qᵢ − Qᵢ₋₁|
A running 64-bit sum of Math.hypot distances. Two fit points at the same place are refused, because the knots would repeat.
sketch-spline.js
lines 167 to 177
PartCraft 0.8.1
function chordParams(Q, closed) {
const t = [0];
const n = Q.length;
for (let i = 1; i < n + (closed ? 1 : 0); i++) {
const a = Q[i - 1], b = Q[i % n];
const d = Math.hypot(b[0] - a[0], b[1] - a[1]);
if (!(d > 0)) throw new Error(tr('model.sketch_spline.fit_points_coincide'));
t.push(t[i - 1] + d);
}
return t;
}
Dense linear solve
64-bit
Solves the linear system that gives a fit spline its poles.
A · X = B by Gauss–Jordan elimination with partial pivoting
The largest entry of each column is taken as pivot and every other row is cleared against it. A pivot of 1e-300 or less stops with a message. There is no scaling of rows and no refinement of the result.
sketch-spline.js
lines 179 to 194
PartCraft 0.8.1
/** Dense solve A·X = B (B: columns), Gaussian elimination with partial pivoting. */
function solveDense(A, B) {
const n = A.length, M = A.map((r, i) => [...r, ...B[i]]), w = B[0].length;
for (let c = 0; c < n; c++) {
let piv = c;
for (let r = c + 1; r < n; r++) if (Math.abs(M[r][c]) > Math.abs(M[piv][c])) piv = r;
if (!(Math.abs(M[piv][c]) > 1e-300)) throw new Error(tr('model.sketch_spline.cannot_pass'));
[M[c], M[piv]] = [M[piv], M[c]];
for (let r = 0; r < n; r++) {
if (r === c || M[r][c] === 0) continue;
const f = M[r][c] / M[c][c];
for (let k = c; k < n + w; k++) M[r][k] -= f * M[c][k];
}
}
return M.map((r, i) => r.slice(n).map((v) => v / M[i][i]));
}
Spline through fit points
64-bit
The cubic spline that passes through the given points, open with natural or given end directions, or closed and periodic.
rows: C(tᵢ) = Qᵢ; ends: C″ = 0 or C′ = (cos a, sin a); solve for the poles
For an open spline the first and last poles are overwritten with the end fit points after the solve, so the curve ends on them exactly. The interior fit points are met to the accuracy of the linear solve.
sketch-spline.js
lines 196 to 232
PartCraft 0.8.1
/**
* Interpolating cubic through Q (chord-length parameters, a knot at every fit point).
* Open: natural ends, or a unit tangent at angle a0 / a1. Closed: periodic.
*/
export function interpFit(Q, { closed = false, a0 = null, a1 = null } = {}) {
if (closed) {
const n = Q.length;
if (n < 3) throw new Error(tr('model.sketch_spline.closed_three'));
const tau = chordParams(Q, true);
const flat = periodicCurve(Q, tau);
const A = Q.map(() => new Array(n).fill(0));
for (let i = 0; i < n; i++) {
const k = i + DEG, N = basisFuns(k, tau[i], DEG, flat.knots);
for (let j = 0; j <= DEG; j++) A[i][(k - DEG + j) % n] += N[j];
}
const P = solveDense(A, Q.map((q) => [q[0], q[1]]));
return periodicCurve(P, tau);
}
const n1 = Q.length;
if (n1 < 2) throw new Error(tr('model.sketch_spline.open_two'));
const t = chordParams(Q, false);
const U = [t[0], t[0], t[0], t[0], ...t.slice(1, -1), t[n1 - 1], t[n1 - 1], t[n1 - 1], t[n1 - 1]];
const m = n1 + 2, p = DEG;
const A = [], B = [];
const row = (u, d, rhs) => {
const k = findSpan(U, p, m, u), D = dersBasisFuns(k, u, p, d, U), r = new Array(m).fill(0);
for (let j = 0; j <= p; j++) r[k - p + j] = D[d][j];
A.push(r); B.push(rhs);
};
const end = (u, ang) => (Number.isFinite(ang) ? row(u, 1, [Math.cos(ang), Math.sin(ang)]) : row(u, 2, [0, 0]));
row(t[0], 0, Q[0].slice()); end(t[0], a0);
for (let i = 1; i < n1 - 1; i++) row(t[i], 0, Q[i].slice());
end(t[n1 - 1], a1); row(t[n1 - 1], 0, Q[n1 - 1].slice());
const P = solveDense(A, B);
P[0] = Q[0].slice(); P[m - 1] = Q[n1 - 1].slice(); // exactly through the end points (shared profile vertices)
return { deg: p, knots: U, poles: P };
}
Area term of a spline
64-bit
The area term of a spline by Green's formula: the enclosed area when the spline is closed.
½ ∫ (x·y′ − y·x′) du, 3-point Gauss–Legendre per knot span (nodes ±√(3 ÷ 5), 0; weights 5/9, 8/9, 5/9)
The integrand of a cubic is a polynomial of degree 5, which the 3-point rule integrates without method error; what remains is rounding. The formula is taken about the sketch origin, so for a curve far from the origin the terms are large and mostly cancel.
sketch-spline.js
lines 275 to 283
PartCraft 0.8.1
/** ½∫(x y′ − y x′) du over the curve: exact for a cubic (degree-5 integrand, 3-point Gauss–Legendre). */
export function areaTerm(c) {
let s = 0;
for (const [u0, u1] of spans(c)) {
const h = (u1 - u0) / 2, m = (u0 + u1) / 2;
for (const [g, w] of GL3) { const [p, d] = derivsAt(c, m + h * g, 1); s += w * h * (p[0] * d[1] - p[1] * d[0]) / 2; }
}
return s;
}
Curvature of a spline
64-bit
The signed curvature of a spline at a parameter.
κ = (x′·y″ − y′·x″) ÷ |C′|³
One 64-bit expression from the first and second derivatives. A zero speed returns 0 instead of dividing by zero.
sketch-spline.js
lines 318 to 323
PartCraft 0.8.1
/** Signed curvature at u (left turn positive). */
export function curvatureAt(c, u) {
const [, d1, d2] = derivsAt(c, u, 2);
const s = Math.hypot(d1[0], d1[1]);
return s > 0 ? (d1[0] * d2[1] - d1[1] * d2[0]) / (s * s * s) : 0;
}
Knot insertion
64-bit
Adds one knot to a spline without changing its shape.
a = (u − Uᵢ) ÷ (Uᵢ₊ₚ − Uᵢ); Qᵢ = (1 − a)·Pᵢ₋₁ + a·Pᵢ
Each new pole is a weighted mean of two old ones with weights between 0 and 1, so the new poles are exact up to one rounding each.
sketch-spline.js
lines 327 to 340
PartCraft 0.8.1
/** Insert knot u once (Boehm); returns a new open flat curve. */
export function insertKnot(c, u) {
const { deg: p, knots: U, poles: P } = c;
const m = P.length, k = findSpan(U, p, m, u);
const Q = [];
for (let i = 0; i <= k - p; i++) Q.push(P[i].slice());
for (let i = k - p + 1; i <= k; i++) {
const a = (u - U[i]) / (U[i + p] - U[i]);
Q.push([(1 - a) * P[i - 1][0] + a * P[i][0], (1 - a) * P[i - 1][1] + a * P[i][1]]);
}
for (let i = k; i < m; i++) Q.push(P[i].slice());
const K = [...U.slice(0, k + 1), u, ...U.slice(k + 1)];
return { deg: p, knots: K, poles: Q };
}
Exact piece of a spline
64-bit
The part of a spline between two parameters as a curve of its own, used by trim.
Both cut parameters are inserted as knots until the curve can be split there; no fitting or sampling is involved. A cut within 1e-12 × max(1, |lo|, |hi|, hi − lo) of an existing knot is moved onto it.
sketch-spline.js
lines 353 to 371
PartCraft 0.8.1
/**
* The exact piece of c between ua < ub, as an open clamped curve (poles, knots). For a periodic curve
* ub may pass the end of the period (the piece runs round the seam).
*/
export function segment(c, ua, ub) {
const [lo, hi] = domain(c);
let src = c;
if (c.periodic && ub > hi + 1e-12 * Math.max(1, Math.abs(hi))) src = periodicCurve(c.P, c.tau, 2);
const eps = 1e-12 * Math.max(1, Math.abs(lo), Math.abs(hi), hi - lo);
let r = raise(src, ua, eps); const A = r.u;
r = raise(r.c, ub, eps); const Bu = r.u;
const { knots: K, poles: P } = r.c, p = c.deg;
let lastA = -1, firstB = -1;
for (let j = 0; j < K.length; j++) { if (K[j] === A) lastA = j; if (K[j] === Bu && firstB < 0) firstB = j; }
const i0 = lastA - p + 1, i1 = firstB + p - 1;
const knots = [A, ...K.slice(i0, i1 + 1), Bu];
const nP = knots.length - p - 1;
return { deg: p, knots, poles: P.slice(i0 - 1, i0 - 1 + nP).map((q) => q.slice()) };
}
Knots rescaled to start at zero
64-bit
Shifts and scales a spline's knots so that they start at 0 and span the number of its pieces.
u′ = (u − a) × (poles − degree) ÷ L
The shape depends only on ratios of knot differences, which the rescaling keeps up to one rounding per knot. The stored knots of a trimmed piece are these rescaled values.
sketch-spline.js
lines 379 to 384
PartCraft 0.8.1
/** Normalise an open curve's knots to start at 0 (keeps the shape; makes stored knots tidy). */
export function normalizeKnots(c) {
const a = c.knots[0], L = c.knots[c.knots.length - 1] - a;
const s = L > 0 ? (c.poles.length - c.deg) / L : 1;
return { deg: c.deg, knots: c.knots.map((u) => (u - a) * s), poles: c.poles };
}
In detail
Method
The NURBS Book algorithms A2.2 and A2.3 (the file header cites them); chord-length parameters with a knot at every fit point; Gaussian elimination with partial pivoting; Boehm knot insertion; Green's formula with 3-point Gauss–Legendre.
Number types
64-bit numbers in plain arrays.
Tolerances
Pivot refused under 1e-300. A cut parameter snaps to an existing knot within 1e-12 × max(1, |lo|, |hi|, hi − lo). Fit points must not coincide (distance > 0).
Limits
The linear solve has no refinement step and its size is the number of fit points plus 2. Rescaling the knots of a cut piece rounds each knot once.
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.