Resources Numerical accuracy Calculation policy
Moves the end of a curve or the edge of a surface so it meets another with position, tangent or curvature continuity, and builds blends between them.
A curve's end is set in closed form from the target's point, direction and curvature vector, in 64-bit. A surface edge is fitted row by row by least squares at sampled stations and refined until the measured gap, angle and curvature difference are within tolerance or the refinement limit is reached. Each result carries the continuity that was measured afterward, so a result outside tolerance is reported as such.
7 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.
Tests tools/nurbs/test.mjs
The factors that link the first control points to the first and second derivative at the start.
a₁ = p ÷ (U[p+1] − U[1]); b₁ = p ÷ (U[p+2] − U[2]); c₂ = (p − 1) ÷ (U[p+1] − U[2])
Each is one division by a knot difference. The divisions are not guarded here; the callers reparameterize to 0…1 and insert knots first.
const startCoefs = (U, p) => ({
a1: p / (U[p + 1] - U[1]),
b1: p / (U[p + 2] - U[2]),
c2: p > 1 ? (p - 1) / (U[p + 1] - U[2]) : 0,
});What a curve's end looks like to the curve that must meet it.
t = C′ ÷ |C′|; out = ± t; speed = |C′|
The first derivative is normalized by vec.js; the curvature vector comes from curve.js. Nothing is rounded.
/** Position, unit direction of travel (`out`: away from the curve past its end) and curvature vector at an end. */
function endFrame(c, end) {
const [a, b] = C.domain(c), u = end === 'start' ? a : b, d = C.derivs(c, u, 1);
const t = V.norm(d[1]);
return { p: d[0], out: end === 'end' ? t : V.mul(t, -1), kv: C.curvatureVector(c, u), speed: V.len(d[1]) };
}Moves the first one to three control points so the curve meets the target with G0, G1 or G2.
P₀ = T; P₁ = P₀ + out × speed ÷ a₁; D₂ = out (C″ · out) + κ speed²; P₂ = P₁ + (a₁ (P₁ − P₀) + D₂ ÷ c₂) ÷ b₁
Closed form, with the curve's own speed kept (1 when it is zero). A rational curve is first replaced by a polynomial fit within 1e-7 and the deviation is added to the report. The result's continuity is measured afterward by continuity.js and returned with it.
export function matchCurve(c, end, target, targetEnd, { continuity = 'G2' } = {}) {
const k = level(continuity);
let w = end === 'start' ? C.reparam(c, 0, 1) : C.reverse(C.reparam(c, 0, 1));
let note = null;
if (w.w) { const q = K.polynomial(w, 1e-7); w = q.curve; note = { refit: q.deviation }; }
if (w.deg < k) w = C.elevate(w, k - w.deg);
// enough control points that the other end keeps its own (k + 1 at each end, none shared)
while (w.cvs.length < 2 * (k + 1)) {
const ks = B.distinctKnots(w.knots), widest = ks.slice(1).reduce((m, x, i) => (x.u - ks[i].u > m.w ? { w: x.u - ks[i].u, u: (x.u + ks[i].u) / 2 } : m), { w: 0, u: 0.5 });
w = C.insertKnot(w, widest.u);
}
const T = endFrame(target, targetEnd);
// the joined path runs from the target into this curve: at our start the direction must be T.out
const P = w.cvs.map((q) => q.slice()), U = w.knots, p = w.deg, { a1, b1, c2 } = startCoefs(U, p);
const d0 = C.derivs(w, 0, 2);
P[0] = T.p.slice();
if (k >= 1) {
const speed = V.len(d0[1]) || 1;
P[1] = V.add(P[0], V.mul(T.out, speed / a1));
if (k >= 2 && p >= 2) {
const d1 = V.mul(T.out, speed);
const tang = V.mul(T.out, V.dot(d0[2], T.out));
const D2 = V.add(tang, V.mul(T.kv, speed * speed));
const Q0 = V.mul(V.sub(P[1], P[0]), a1), Q1 = V.add(Q0, V.mul(D2, 1 / c2));
P[2] = V.add(P[1], V.mul(Q1, 1 / b1));
void d1;
}
}
let out = { deg: p, knots: U.slice(), cvs: P, w: null };
if (end === 'end') out = C.reverse(out);
const rep = G.curveJoin(target, targetEnd, out, end, G.TOL);
return { curve: out, report: { ...rep, want: continuity, ...(note || {}) } };
}A Bézier curve joining two curve ends with the asked continuity.
G1: P₁ = A + out × α ÷ 3; G2: P₁ = A + out × α ÷ 5, P₂ = 2 P₁ − A + κ α² ÷ 20; α = chord × bulge
The control points are set directly from the end frames; the divisors 3, 5 and 20 are the Bézier derivative factors for degrees 3 and 5. Ends closer than 1e-12 are refused, an absolute distance.
export function blendCurve(a, endA, b, endB, { continuity = 'G2', bulge = [1, 1] } = {}) {
const k = level(continuity);
const A = endFrame(a, endA), Bf = endFrame(b, endB);
const chord = V.dist(A.p, Bf.p);
if (chord < 1e-12) fail('nurbs.err.blend_ends');
const al = chord * bulge[0], be = chord * bulge[1];
const arrive = V.mul(Bf.out, -1); // the blend arrives at b running into it
let cvs;
if (k === 0) cvs = [A.p, Bf.p];
else if (k === 1) cvs = [A.p, V.add(A.p, V.mul(A.out, al / 3)), V.sub(Bf.p, V.mul(arrive, be / 3)), Bf.p];
else {
const P1 = V.add(A.p, V.mul(A.out, al / 5)), P2 = V.add(V.sub(V.mul(P1, 2), A.p), V.mul(A.kv, al * al / 20));
const P4 = V.sub(Bf.p, V.mul(arrive, be / 5)), P3 = V.add(V.sub(V.mul(P4, 2), Bf.p), V.mul(Bf.kv, be * be / 20));
cvs = [A.p, P1, P2, P3, P4, Bf.p];
}
const deg = cvs.length - 1;
const curve = { deg, knots: [...new Array(deg + 1).fill(0), ...new Array(deg + 1).fill(1)], cvs: cvs.map((q) => q.slice()), w: null };
return { curve, report: { start: G.curveJoin(a, endA, curve, 'start'), end: G.curveJoin(curve, 'end', b, endB), want: continuity } };
}For a point near the target's edge: the nearest edge point, the outward direction in the tangent plane, the normal and the normal curvature.
d = (n × tₑ) ÷ |n × tₑ|, turned to point away from the surface
The nearest point comes from curve.js closestParam; its method and stopping rule are in that file. An edge tangent shorter than 1e-300 is replaced by any vector perpendicular to the normal.
/** At target's edge, for a point near it: the point on the edge, the direction away from the target (in its tangent plane,
* across the edge), the target normal and its normal curvature in that direction. */
function targetFrame(target, tEdge, tCurve, p) {
const cp = C.closestParam(tCurve, p), uv = S.edgeUV(target, tEdge, cp.u);
const n = S.normal(target, ...uv), te = C.derivs(tCurve, cp.u, 1)[1];
let d = V.norm(V.cross(n, V.len(te) > 1e-300 ? te : V.perp(n)));
// pointing away from the target: against the inward cross derivative
const e = S.edgeParams(target, tEdge), D = S.derivs(target, ...uv, 1);
const inward = V.mul(e.dir === 'u' ? D[1][0] : D[0][1], e.inward);
if (V.dot(d, inward) > 0) d = V.mul(d, -1);
return { p: cp.point, d, n, k: S.normalCurvature(target, ...uv, d), u: cp.u };
}Refits the first one to three rows of control points so the surface meets the target along an edge.
stations vᵢ = i ÷ M, M = max(samples, 4 × columns); rows fitted by least squares; accepted when gap ≤ g0, angle ≤ g1 and Δκ ≤ g2 (or relative ≤ g2rel)
The wanted rows are computed at M + 1 stations and fitted with the surface's own knots, so between stations the match is as good as the fit. Every span along the edge is halved and the fit repeated, up to 6 times; the last candidate is returned with its measured gap, angle and curvature difference even when it is outside tolerance.
export function matchSurface(s, edge, target, targetEdge, { continuity = 'G1', tol = G.TOL, maxRefine = 6, samples = 33 } = {}) {
const k = level(continuity);
let w = toU0(s, edge), refit = 0;
if (w.w) {
const r = F.fitSurfaceTo((u, v) => S.point(w, u, v), Math.max(3, w.degU), Math.max(3, w.degV), tol.g0 / 10, { start: [w.cvs.length, w.cvs[0].length] });
w = r.surface; refit = r.deviation;
}
if (w.degU < Math.max(k, 1)) w = S.elevate(w, 'u', Math.max(k, 1) - w.degU);
// rows near the edge free to move without touching the far edge: k + 1 rows here, k + 1 there, none shared
while (w.cvs.length < 2 * (k + 1) + 1) {
const ks = B.distinctKnots(w.knotsU);
w = S.insertKnot(w, 'u', (ks[0].u + ks[1].u) / 2);
}
// the target's edge curve in our v space (degree and knots), so G0 can be exact when the edges are the same curve
const tc = S.edgeCurve(target, targetEdge);
const tcPoly = tc.w ? null : C.reparam(tc, 0, 1);
const ourEdge = S.edgeCurve(w, 'u0');
const forward = V.dist(C.endsOf(ourEdge).start, C.endsOf(tc).start) + V.dist(C.endsOf(ourEdge).end, C.endsOf(tc).end)
<= V.dist(C.endsOf(ourEdge).start, C.endsOf(tc).end) + V.dist(C.endsOf(ourEdge).end, C.endsOf(tc).start);
if (tcPoly) {
const tk = forward ? tcPoly : C.reverse(tcPoly);
if (w.degV < tk.deg) w = S.elevate(w, 'v', tk.deg - w.degV);
const tk2 = tk.deg < w.degV ? C.elevate(tk, w.degV - tk.deg) : tk;
for (const x of B.distinctKnots(tk2.knots)) if (x.u > 0 && x.u < 1) { const have = B.multiplicity(w.knotsV, x.u, 1e-12); if (x.m > have) w = S.insertKnot(w, 'v', x.u, x.m - have); }
}
while (w.cvs[0].length < 8) {
const ks = B.distinctKnots(w.knotsV), widest = ks.slice(1).reduce((m, x, i) => (x.u - ks[i].u > m.w ? { w: x.u - ks[i].u, u: (x.u + ks[i].u) / 2 } : m), { w: 0, u: 0.5 });
w = S.insertKnot(w, 'v', widest.u);
}
let best = null;
for (let round = 0; round <= maxRefine; round++) {
const U = w.knotsU, p = w.degU, { a1, b1, c2 } = startCoefs(U, p);
const M = Math.max(samples, w.cvs[0].length * 4);
const tv = Array.from({ length: M + 1 }, (_, i) => i / M);
const tparam = (v) => (forward ? v : 1 - v);
const r0 = [], r1 = [], r2 = [];
tv.forEach((v) => {
const D = S.derivs(w, 0, v, 2);
const fr = targetFrame(target, targetEdge, tc, tcPoly ? C.point(tcPoly, tparam(v)) : C.point(tc, C.domain(tc)[0] + (C.domain(tc)[1] - C.domain(tc)[0]) * tparam(v)));
r0.push(fr.p);
if (k >= 1) {
const speed = V.len(D[1][0]) || 1;
const D1 = V.mul(fr.d, speed);
r1.push(V.add(fr.p, V.mul(D1, 1 / a1)));
if (k >= 2) {
const tang = V.sub(D[2][0], V.mul(fr.n, V.dot(D[2][0], fr.n)));
const D2 = V.add(tang, V.mul(fr.n, fr.k * speed * speed));
const Q0 = D1, Q1 = V.add(Q0, V.mul(D2, 1 / c2));
r2.push(V.add(V.add(fr.p, V.mul(D1, 1 / a1)), V.mul(Q1, 1 / b1)));
}
}
});
const fitRow = (Q) => F.fitFixedKnots(Q, tv, w.degV, w.knotsV).cvs;
const cvs = w.cvs.map((r) => r.map((q) => q.slice()));
cvs[0] = fitRow(r0);
if (k >= 1) cvs[1] = fitRow(r1);
if (k >= 2) cvs[2] = fitRow(r2);
const cand = { ...w, cvs, w: null };
const back = fromU0(cand, edge, s);
const m = G.surfaceEdge(back, edge, target, targetEdge, { samples, tol });
const good = [m.gap <= tol.g0, m.angle <= tol.g1, m.dk.max <= tol.g2 || m.rel <= tol.g2rel].slice(0, k + 1).every(Boolean);
best = { surface: back, report: { want: continuity, measured: m.cls, gap: m.gap, angle: m.angle, dk: m.dk.max, rel: m.rel, refinements: round, cvs: [cand.cvs.length, cand.cvs[0].length], refit, tol } };
if (good) return best;
// refine along the edge where it is worst: halve every span
const ks = B.distinctKnots(w.knotsV);
for (let i = 0; i + 1 < ks.length; i++) w = S.insertKnot(w, 'v', (ks[i].u + ks[i + 1].u) / 2);
void rowCurve;
}
return best;
}A surface joining two surface edges, built from blend sections along them.
sections at t = i ÷ n, n = 8, 16, … ≤ 128; each a Bézier of degree 1, 3 or 5; interpolated along the edge
Each section is set in closed form from the two edge frames. The sections are interpolated along the edge, and the count doubles until both edges measure within tolerance or 129 sections are reached; the measured values are returned either way.
export function blendSurface(sA, edgeA, sB, edgeB, { continuity = 'G1', bulge = [1, 1], tol = G.TOL, maxSections = 129 } = {}) {
const k = level(continuity);
const cA = S.edgeCurve(sA, edgeA), cB = S.edgeCurve(sB, edgeB);
const [a0, a1] = C.domain(cA), [b0, b1] = C.domain(cB);
const eA = C.endsOf(cA), eB = C.endsOf(cB);
const forward = V.dist(eA.start, eB.start) + V.dist(eA.end, eB.end) <= V.dist(eA.start, eB.end) + V.dist(eA.end, eB.start);
const section = (t) => {
const pa = C.point(cA, a0 + (a1 - a0) * t), pb = C.point(cB, b0 + (b1 - b0) * (forward ? t : 1 - t));
const fa = targetFrame(sA, edgeA, cA, pa), fb = targetFrame(sB, edgeB, cB, pb);
const chord = V.dist(fa.p, fb.p);
if (chord < 1e-12) fail('nurbs.err.blend_ends');
const al = chord * bulge[0], be = chord * bulge[1];
if (k === 0) return [fa.p, fb.p];
if (k === 1) return [fa.p, V.add(fa.p, V.mul(fa.d, al / 3)), V.add(fb.p, V.mul(fb.d, be / 3)), fb.p];
const P1 = V.add(fa.p, V.mul(fa.d, al / 5)), P2 = V.add(V.sub(V.mul(P1, 2), fa.p), V.mul(fa.n, fa.k * al * al / 20));
const P4 = V.add(fb.p, V.mul(fb.d, be / 5)), P3 = V.add(V.sub(V.mul(P4, 2), fb.p), V.mul(fb.n, fb.k * be * be / 20));
return [fa.p, P1, P2, P3, P4, fb.p];
};
const deg = k === 0 ? 1 : k === 1 ? 3 : 5;
let n = 8, best = null;
for (;;) {
const ts = Array.from({ length: n + 1 }, (_, i) => i / n), secs = ts.map(section);
// rows across: section point index i (u); along: the sections (v)
const lists = secs[0].map((_, i) => secs.map((sc) => sc[i]));
const r = F.interpolateMany(lists, ts, 3);
const surface = { degU: deg, degV: r.deg, knotsU: [...new Array(deg + 1).fill(0), ...new Array(deg + 1).fill(1)], knotsV: r.knots, cvs: r.rows.map((row) => row.map((q) => q.slice())), w: null };
const mA = G.surfaceEdge(surface, 'u0', sA, edgeA, { tol }), mB = G.surfaceEdge(surface, 'u1', sB, edgeB, { tol });
const okk = (m) => [m.gap <= tol.g0, m.angle <= tol.g1, m.dk.max <= tol.g2 || m.rel <= tol.g2rel].slice(0, k + 1).every(Boolean);
best = { surface, report: { want: continuity, sections: n + 1, a: { measured: mA.cls, gap: mA.gap, angle: mA.angle, dk: mA.dk.max }, b: { measured: mB.cls, gap: mB.gap, angle: mB.angle, dk: mB.dk.max }, tol } };
if ((okk(mA) && okk(mB)) || n + 1 >= maxSections) return best;
n = Math.min(maxSections - 1, n * 2);
}
}The arithmetic is 64-bit, but the method or the data sets the accuracy: a solver tolerance, a mesh, a step size, a datum shift, a table. The limit is stated where the result is shown.