Works out the control points and knots of a cubic spline that a drawing stores only by its fit points and end tangents.
How the precision is kept
The fit points are interpolated by a clamped cubic B-spline: one linear system, solved directly in 64-bit numbers. The matrix rows are Float64Array and the right-hand sides are plain 64-bit numbers. The solve is Gaussian elimination with partial pivoting inside the band; a pivot is refused only when it is exactly zero.
Open item. Refuse a pivot that is tiny compared with the largest entry of the matrix, not only one that is exactly zero.
Calculations under this policy
5 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.
Each chord is one Math.hypot call and the parameters are a running sum in fit point order. For knot parameterization 1 the chord's square root is added, for 2 the value 1. The result is refused when the total L is not greater than 0.
dwg/spline-fit.js
lines 17 to 26
CadCraft 0.8.2
const u = [0];
for (let i = 1; i < n; i++) {
const a = fit[i - 1], b = fit[i];
let d = Math.hypot(b.x - a.x, b.y - a.y, (b.z || 0) - (a.z || 0));
if (knotParam === 1) d = Math.sqrt(d); else if (knotParam === 2) d = 1;
u.push(u[i - 1] + d);
}
const L = u[n - 1];
if (!(L > 0)) return null;
const knots = [0, 0, 0, 0, ...u.slice(1, n - 1), L, L, L, L];
Cubic basis functions at a parameter
64-bit
The four B-spline basis values that weigh the control points at one fit point.
Cox–de Boor recurrence, degree 3
The recurrence divides by the difference of two knots without a guard. The span search moves past knots equal to the parameter, so the first-level difference is positive except when the last chord has zero length.
dwg/spline-fit.js
lines 29 to 44
CadCraft 0.8.2
const basis = (t) => {
let span = 3;
while (span < m - 1 && t >= knots[span + 1]) span++;
const N = [1, 0, 0, 0], left = [], right = [];
for (let j = 1; j <= 3; j++) {
left[j] = t - knots[span + 1 - j]; right[j] = knots[span + j] - t;
let saved = 0;
for (let r = 0; r < j; r++) {
const tmp = N[r] / (right[r + 1] + left[j - r]);
N[r] = saved + right[r + 1] * tmp;
saved = left[j - r] * tmp;
}
N[j] = saved;
}
return { span, N };
};
Start condition
64-bit
The rows that pin the first control point and set the spline's start tangent or a natural start.
A tangent counts as given when any of its components is not zero; it is used as it is, without being normalized here. The natural-end row divides by the first two interior knot values, which are positive when the fit points are distinct.
The rows that pin the last control point and set the spline's end tangent or a natural end.
Pₘ₋₁ = Qₙ₋₁; Pₘ₋₁ − Pₘ₋₂ = T₁ × (L − uₘ₋₁) ÷ 3 or the zero second derivative row
The two knot differences e1 = L − u(m−1) and e2 = L − u(m−2) are formed by subtraction and then divided by. They are positive when the last fit points are distinct; there is no guard when they are not.
Solves the banded system for the x, y and z of every control point.
A · P = B
Gaussian elimination with partial pivoting among the next four rows, then back substitution over the next five columns. A pivot is refused only when it equals 0 exactly, so a nearly singular system is solved and returned without a warning.
dwg/spline-fit.js
lines 67 to 85
CadCraft 0.8.2
// Gaussian elimination with partial pivoting (the matrix is banded; m is small)
for (let c = 0; c < m; c++) {
let piv = c;
for (let r = c + 1; r < Math.min(m, c + 5); r++) if (Math.abs(A[r][c]) > Math.abs(A[piv][c])) piv = r;
if (!A[piv][c]) return null;
if (piv !== c) { [A[piv], A[c]] = [A[c], A[piv]]; [B[piv], B[c]] = [B[c], B[piv]]; }
for (let r = c + 1; r < Math.min(m, c + 5); r++) {
const f = A[r][c] / A[c][c];
if (!f) continue;
for (let k = c; k < Math.min(m, c + 6); k++) A[r][k] -= f * A[c][k];
for (let k = 0; k < 3; k++) B[r][k] -= f * B[c][k];
}
}
const P = new Array(m);
for (let r = m - 1; r >= 0; r--) {
const v = B[r].slice();
for (let k = r + 1; k < Math.min(m, r + 6); k++) if (A[r][k]) for (let q = 0; q < 3; q++) v[q] -= A[r][k] * P[k][q];
P[r] = v.map((x) => x / A[r][r]);
}
In detail
Method
Clamped cubic B-spline interpolation at chord-length parameters (square root of the chord for knot parameterization 1, uniform for 2), with the given unit end tangents or a zero second derivative where a tangent is (0,0,0). The file's comment says the method was found by matching another converter's DXF output to 1e-9.
Number types
64-bit numbers; matrix rows in Float64Array.
Tolerances
None. The total parameter length must be greater than 0 and a pivot must be different from 0, both exact tests.
Limits
Degree 3 only, 2 to 2000 fit points, fit tolerance 0. Returns nothing for other input. There is no test for a nearly singular system, such as two fit points at the same place.
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.