Builds one B-spline surface that passes through a grid of profile curves and guide curves, and reports how far it strays from them.
How the precision is kept
All knots, parameters and poles are 64-bit numbers; the band matrices are Float64Array. Degree raising, knot merging and the final sum are algebraically exact steps carried out as 64-bit band solves. The accuracy of the result is set by the method: curve crossings are refined with numeric derivatives and accepted within a gap tolerance, and curves that are not polynomial are fitted until a sampled error is under the fit tolerance (by default 1e-7 of the network's size).
Calculations under this policy
19 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.
The point of a curve at a parameter, from its poles and knots.
C(t) = Σ Nⱼ(t) × Pⱼ over the p + 1 poles of the span
A sum of p + 1 products in 64-bit numbers, added in pole order. The span and the basis values come from findSpan and basisFuns in sketch-spline.js, which this review did not cover.
gordon.js
lines 53 to 58
PartCraft 0.8.1
export function evalCurve(c, t) {
const n = c.P.length, k = SP.findSpan(c.U, c.p, n, t), N = SP.basisFuns(k, t, c.p, c.U);
const d = c.P[0].length, out = new Array(d).fill(0);
for (let j = 0; j <= c.p; j++) { const q = c.P[k - c.p + j]; for (let a = 0; a < d; a++) out[a] += N[j] * q[a]; }
return out;
}
Point on a B-spline surface
64-bit
The point of the surface at (u, v), from its pole grid and two knot vectors.
S(u, v) = Σᵢ Σⱼ Nᵢ(u) × Nⱼ(v) × Pᵢⱼ
The two basis values are multiplied first, then applied to each coordinate; (pu + 1) × (pv + 1) terms added in order in 64-bit numbers.
gordon.js
lines 60 to 70
PartCraft 0.8.1
export function evalSurface(s, u, v) {
const nu = s.P.length, nv = s.P[0].length;
const ku = SP.findSpan(s.U, s.pu, nu, u), kv = SP.findSpan(s.V, s.pv, nv, v);
const Nu = SP.basisFuns(ku, u, s.pu, s.U), Nv = SP.basisFuns(kv, v, s.pv, s.V);
const out = [0, 0, 0];
for (let i = 0; i <= s.pu; i++) for (let j = 0; j <= s.pv; j++) {
const w = Nu[i] * Nv[j], q = s.P[ku - s.pu + i][kv - s.pv + j];
out[0] += w * q[0]; out[1] += w * q[1]; out[2] += w * q[2];
}
return out;
}
Interpolation knots by averaging
64-bit
The knot vector used to interpolate points at given parameters.
inner knot j = (t_j + … + t_(j+p−1)) ÷ p
Each inner knot is a sum of p parameters divided once by p; the end knots are copies of the first and last parameter, so the ends are clamped without rounding.
gordon.js
lines 72 to 79
PartCraft 0.8.1
/** Interpolation knots by averaging (NURBS Book eq. 9.8) for n params and degree p. */
export function averagedKnots(t, p) {
const n = t.length, U = [];
for (let i = 0; i <= p; i++) U.push(t[0]);
for (let j = 1; j <= n - 1 - p; j++) { let s = 0; for (let i = j; i < j + p; i++) s += t[i]; U.push(s / p); }
for (let i = 0; i <= p; i++) U.push(t[n - 1]);
return U;
}
Band solve for interpolation
64-bit
Factors the collocation matrix once and solves for the poles of every curve that shares its parameters.
N × P = Q, with N banded; LU without pivoting, then forward and back substitution
The band is a Float64Array and the right-hand sides are 64-bit arrays. There is no pivoting; the code's comment gives total positivity of B-spline collocation matrices as the reason. A pivot whose magnitude is not above 1e-13 (absolute; basis values lie in 0 to 1) stops the build with an error instead of dividing.
gordon.js
lines 83 to 114
PartCraft 0.8.1
/**
* The collocation matrix of degree p / knots U (n functions) at n params, factored as a band LU without pivoting
* (B-spline collocation matrices are totally positive: de Boor). → solve(points) → poles.
*/
function bandLU(p, U, n, params) {
if (params.length !== n) fail('Internal: the interpolation is not square.'); // i18n-ignore: internal check, never reached by a person
const rows = params.map((t) => { const k = SP.findSpan(U, p, n, t); return { lo: k - p, N: SP.basisFuns(k, t, p, U) }; });
let kl = 0, ku = 0;
rows.forEach((r, i) => { kl = Math.max(kl, i - r.lo); ku = Math.max(ku, r.lo + p - i); });
const w = kl + ku + 1, A = new Float64Array(n * w);
rows.forEach((r, i) => r.N.forEach((v, j) => { A[i * w + (r.lo + j - i + kl)] = v; }));
for (let c = 0; c < n; c++) {
const piv = A[c * w + kl];
if (!(Math.abs(piv) > 1e-13)) fail(tr('model.gordon.err.too_close'));
for (let r = c + 1; r <= Math.min(n - 1, c + kl); r++) {
const ir = r * w + (c - r + kl), f = A[ir] / piv;
if (f === 0) continue;
A[ir] = f;
for (let j = c + 1; j <= Math.min(n - 1, c + ku); j++) A[r * w + (j - r + kl)] -= f * A[c * w + (j - c + kl)];
}
}
return (Q) => {
const d = Q[0].length, X = Q.map((q) => q.slice());
for (let r = 0; r < n; r++) for (let c = Math.max(0, r - kl); c < r; c++) { const f = A[r * w + (c - r + kl)]; if (f) for (let a = 0; a < d; a++) X[r][a] -= f * X[c][a]; }
for (let r = n - 1; r >= 0; r--) {
for (let j = r + 1; j <= Math.min(n - 1, r + ku); j++) { const v = A[r * w + (j - r + kl)]; if (v) for (let a = 0; a < d; a++) X[r][a] -= v * X[j][a]; }
const piv = A[r * w + kl];
for (let a = 0; a < d; a++) X[r][a] /= piv;
}
return X;
};
}
Degree raising and knot refinement
64-bit
Re-expresses a curve with a higher degree or more knots without changing its shape.
new poles = solve(N_fine(g), C(g)), g = Greville points of the finer knot vector
The curve is evaluated at the Greville points of the finer space and interpolated there. The algebra reproduces the curve exactly; in 64-bit numbers the new poles carry the rounding of one evaluation and one band solve. When degree and knots already match within 1e-12 the poles are copied unchanged.
gordon.js
lines 149 to 158
PartCraft 0.8.1
/** Poles in (p, U) → poles of the same curve in the finer space (d, W): interpolation at W's Greville points (exact). */
function projector(p, U, n, d, W) {
if (sameSpace(p, U, d, W)) return (P) => P.map((q) => q.slice());
const m = W.length - d - 1, g = greville(d, W, m), solve = bandLU(d, W, m, g);
const rows = g.map((t) => { const k = SP.findSpan(U, p, n, t); return { lo: k - p, N: SP.basisFuns(k, t, p, U) }; });
return (P) => {
const dim = P[0].length;
return solve(rows.map((r) => { const o = new Array(dim).fill(0); r.N.forEach((v, j) => { const q = P[r.lo + j]; for (let a = 0; a < dim; a++) o[a] += v * q[a]; }); return o; }));
};
}
Merged knot vector
64-bit
The smallest knot vector that contains the knot vectors of several curves.
Knots closer than 1e-9 are treated as one knot and keep the first value met; within one curve the limit is 1e-10 (knotsToMults). Both are absolute on the 0 to 1 range. The first and last knots are set to exactly 0 and 1.
gordon.js
lines 136 to 145
PartCraft 0.8.1
/** The smallest degree-d knot vector containing each of the (degree-d) knot vectors. */
function unionKnots(list, d) {
const all = [];
for (const U of list) { const { k, m } = knotsToMults(U); k.forEach((u, i) => all.push([u, m[i]])); }
all.sort((a, b) => a[0] - b[0]);
const k = [], m = [];
for (const [u, x] of all) { if (k.length && u - k[k.length - 1] <= 1e-9) m[m.length - 1] = Math.max(m[m.length - 1], x); else { k.push(u); m.push(x); } }
k[0] = 0; k[k.length - 1] = 1;
return multsToKnots(k, m.map((x, i) => (i === 0 || i === k.length - 1 ? d + 1 : Math.min(x, d))));
}
Joining exact curves end to end
64-bit
Makes one B-spline from a chain of lines and splines, with a corner knot at each joint.
knot = a + (b − a) × u for each piece's inner knots; joint knots repeated d times
Pieces are raised to the highest degree by the band solve, then their inner knots are mapped linearly into the piece's share of 0 to 1. The poles are copied, dropping the repeated joint pole.
gordon.js
lines 302 to 316
PartCraft 0.8.1
/** Exact curves on [0, 1] placed end to end at the breaks ts (C0 joins). */
function concatenate(parts, ts) {
const d = Math.max(...parts.map((c) => c.p));
const el = parts.map((c) => elevate(c, d));
const U = new Array(d + 1).fill(0), P = [];
el.forEach((c, k) => {
const a = ts[k], b = ts[k + 1];
const inner = c.U.slice(d + 1, c.U.length - d - 1).map((u) => a + (b - a) * u);
U.push(...inner);
if (k < el.length - 1) U.push(...new Array(d).fill(b));
P.push(...(k ? c.P.slice(1) : c.P));
});
U.push(...new Array(d + 1).fill(1));
return { p: d, U, P };
}
Sketch chains as world curves
64-bit
Turns the chosen lines, arcs, circles and splines of a sketch into curves over 0 to 1, one per chain.
share of a piece = its length ÷ total length; arc length = r × sweep; circle = 2π × r
Ends are joined when they are within 1e-6 times the largest sketch coordinate (at least 1e-6). The chain parameter is split by arc length, with the last break set to exactly 1. Spline knots are rescaled with (u − lo) ÷ (hi − lo); Math.PI is used for circles.
gordon.js
lines 223 to 298
PartCraft 0.8.1
export function sketchCurves(sk, frame, { entities = null, label = tr('model.gordon.label.sketch') } = {}) {
const want = entities ? new Set(entities) : null;
const W = (p) => V.toWorld(frame, p[0], p[1]);
const pieces = [];
let size = 1;
for (const e of sk.entities) {
if (want ? !want.has(e.id) : e.construction) continue;
if (e.type === 'line') { const a = pointXY(sk, e.p1), b = pointXY(sk, e.p2); pieces.push({ kind: 'line', a, b, len: Math.hypot(b[0] - a[0], b[1] - a[1]) }); }
else if (e.type === 'arc') {
const { c, r, a0, sweep } = arcAngles(sk, e);
pieces.push({ kind: 'arc', c, r, a0, sweep, a: [c[0] + r * Math.cos(a0), c[1] + r * Math.sin(a0)], b: [c[0] + r * Math.cos(a0 + sweep), c[1] + r * Math.sin(a0 + sweep)], len: r * sweep });
} else if (e.type === 'circle') { const c = pointXY(sk, e.c); pieces.push({ kind: 'arc', c, r: e.r, a0: 0, sweep: 2 * Math.PI, closed: true, len: 2 * Math.PI * e.r }); }
else if (e.type === 'spline') {
const bs = splineCurve(sk, e), [lo, hi] = SP.domain(bs);
const at = (u) => SP.pointAt(bs, lo + (hi - lo) * u);
pieces.push({ kind: 'spline', bs, lo, hi, closed: !!e.closed, a: at(0), b: at(1), len: SP.arcLength(bs) });
} else continue;
for (const p of [pieces[pieces.length - 1].a, pieces[pieces.length - 1].b, pieces[pieces.length - 1].c]) if (p) size = Math.max(size, Math.abs(p[0]), Math.abs(p[1]));
}
if (!pieces.length) fail(tr('model.gordon.err.no_curves', { label }));
const tol = 1e-6 * size;
const near = (p, q) => Math.hypot(p[0] - q[0], p[1] - q[1]) <= tol;
const at2 = (pc, u) => { // a piece at its own parameter u ∈ [0, 1], sketch coordinates
if (pc.kind === 'line') return [pc.a[0] + (pc.b[0] - pc.a[0]) * u, pc.a[1] + (pc.b[1] - pc.a[1]) * u];
if (pc.kind === 'arc') { const t = pc.a0 + pc.sweep * u; return [pc.c[0] + pc.r * Math.cos(t), pc.c[1] + pc.r * Math.sin(t)]; }
return SP.pointAt(pc.bs, pc.lo + (pc.hi - pc.lo) * u);
};
const out = [];
const open = [];
for (const pc of pieces) {
if (!pc.closed) { open.push(pc); continue; }
out.push(samplerCurve((t) => W(at2(pc, ((t % 1) + 1) % 1)), { closed: true, label: pc.kind === 'arc' ? tr('model.gordon.label.circle', { label }) : tr('model.gordon.label.closed_spline', { label }) }));
}
const deg = (p) => open.reduce((k, pc) => k + (near(pc.a, p) ? 1 : 0) + (near(pc.b, p) ? 1 : 0), 0);
for (const pc of open) for (const p of [pc.a, pc.b]) if (deg(p) > 2) fail(tr('model.gordon.err.branches', { label, at: p.map((v) => +v.toFixed(3)).join(', ') }));
const used = new Set();
while (used.size < open.length) {
// start at a free end when the chain has one
let first = open.find((pc) => !used.has(pc) && (deg(pc.a) === 1 || deg(pc.b) === 1)) || open.find((pc) => !used.has(pc));
const chain = [{ pc: first, rev: deg(first.a) !== 1 && deg(first.b) === 1 }];
used.add(first);
for (;;) {
const last = chain[chain.length - 1], tip = last.rev ? last.pc.a : last.pc.b;
const nx = open.find((pc) => !used.has(pc) && (near(pc.a, tip) || near(pc.b, tip)));
if (!nx) break;
used.add(nx); chain.push({ pc: nx, rev: !near(nx.a, tip) });
}
const start = chain[0].rev ? chain[0].pc.b : chain[0].pc.a, end = chain[chain.length - 1].rev ? chain[chain.length - 1].pc.a : chain[chain.length - 1].pc.b;
const closed = chain.length > 1 && near(start, end);
const total = chain.reduce((s, x) => s + x.pc.len, 0);
if (!(total > 0)) fail(tr('model.gordon.err.zero_length', { label }));
const ts = [0];
for (const x of chain) ts.push(ts[ts.length - 1] + x.pc.len / total);
ts[ts.length - 1] = 1;
const f = (t0) => {
const t = closed ? ((t0 % 1) + 1) % 1 : Math.min(1, Math.max(0, t0));
let k = 0;
while (k < chain.length - 1 && t > ts[k + 1]) k++;
const u = (t - ts[k]) / (ts[k + 1] - ts[k]);
return W(at2(chain[k].pc, chain[k].rev ? 1 - u : u));
};
// exact form when every piece is polynomial (lines, open clamped splines)
let bs = null;
if (!closed && chain.every((x) => x.pc.kind === 'line' || (x.pc.kind === 'spline' && clamped(x.pc.bs)))) {
const parts = chain.map((x) => {
let c;
if (x.pc.kind === 'line') c = { p: 1, U: [0, 0, 1, 1], P: [W(x.pc.a), W(x.pc.b)] };
else { const b = x.pc.bs, [lo, hi] = SP.domain(b); c = { p: b.deg, U: b.knots.map((u) => (u - lo) / (hi - lo)), P: b.poles.map(W) }; }
return x.rev ? reverseBs(c) : c;
});
bs = concatenate(parts, ts);
}
out.push(samplerCurve(f, { breaks: ts.slice(1, -1), seamBreak: closed, closed, bs, label: chain.length > 1 ? tr('model.gordon.label.chain', { label, n: chain.length }) : label }));
}
return out;
}
Closest points of two segments
64-bit
The nearest pair of points on two straight segments, used to find where two curves may cross.
s = (b·f − c·e) ÷ (a·e − b²), t = (b·s + f) ÷ e, both clamped to 0…1
For nearly parallel segments a·e − b² subtracts two nearly equal numbers; the only guard is that it must exceed 1e-300, an absolute value. The result is only the starting point for the refinement that follows, so this affects where the iteration starts, not the crossing that is kept.
gordon.js
lines 320 to 338
PartCraft 0.8.1
function segSeg(p1, q1, p2, q2) {
// closest points of two segments (Ericson, Real-Time Collision Detection 5.1.9)
const d1 = V.sub(q1, p1), d2 = V.sub(q2, p2), r = V.sub(p1, p2);
const a = V.dot(d1, d1), e = V.dot(d2, d2), f = V.dot(d2, r);
let s, t;
if (a <= 1e-300 && e <= 1e-300) { s = t = 0; }
else if (a <= 1e-300) { s = 0; t = Math.min(1, Math.max(0, f / e)); }
else {
const c = V.dot(d1, r);
if (e <= 1e-300) { t = 0; s = Math.min(1, Math.max(0, -c / a)); }
else {
const b = V.dot(d1, d2), den = a * e - b * b;
s = den > 1e-300 ? Math.min(1, Math.max(0, (b * f - c * e) / den)) : 0;
t = (b * s + f) / e;
if (t < 0) { t = 0; s = Math.min(1, Math.max(0, -c / a)); } else if (t > 1) { t = 1; s = Math.min(1, Math.max(0, (b - c) / a)); }
}
}
return { s, t, d: V.len(V.sub(V.add(p1, V.mul(d1, s)), V.add(p2, V.mul(d2, t)))) };
}
Crossing of two curves
64-bit
The parameters on a profile and a guide where they meet, and the gap between them there.
minimize |A(s) − B(t)|² by Gauss–Newton; point = (A(s) + B(t)) ÷ 2; gap = |A(s) − B(t)|
Up to 40 Gauss–Newton steps with derivatives taken as central differences of step 1e-7 in the parameter (one-sided at the ends of an open curve). It stops when both steps are under 1e-15, when the 2 × 2 determinant is not above 1e-300, or when the gap grows on any step after the first four. The grid point kept is the midpoint of the two curve points.
gordon.js
lines 346 to 363
PartCraft 0.8.1
/** Gauss–Newton on |A(s) − B(t)|² from (s, t). */
function refine(A, B, s, t) {
let a = A.f(s), b = B.f(t);
for (let it = 0; it < 40; it++) {
const r = V.sub(a, b), da = deriv(A, s), db = deriv(B, t);
const m11 = V.dot(da, da), m12 = -V.dot(da, db), m22 = V.dot(db, db);
const g1 = -V.dot(da, r), g2 = V.dot(db, r);
const det = m11 * m22 - m12 * m12;
if (!(Math.abs(det) > 1e-300)) break;
const ds = (g1 * m22 - m12 * g2) / det, dt = (m11 * g2 - m12 * g1) / det;
const s2 = wrap(A, s + ds), t2 = wrap(B, t + dt);
const a2 = A.f(s2), b2 = B.f(t2);
if (V.len(V.sub(a2, b2)) > V.len(r) * (1 + 1e-12) + 1e-300 && it > 3) break;
s = s2; t = t2; a = a2; b = b2;
if (Math.abs(ds) < 1e-15 && Math.abs(dt) < 1e-15) break;
}
return { s, t, d: V.len(V.sub(a, b)), p: V.mul(V.add(a, b), 0.5) };
}
All crossings of a profile and a guide
64-bit
Finds every place two curves come within the gap tolerance, so the network can be checked to be a grid.
Each curve is sampled as a polyline; a segment pair is a candidate when its distance is under tol + 0.25 × the two segment lengths and is a local minimum among its neighbors. A refined crossing is kept when its gap is at most tol, and dropped as a duplicate when within 1e-7 in both parameters or within max(tol, 1e-9) × 0.01 in distance of one already kept.
gordon.js
lines 372 to 400
PartCraft 0.8.1
/** Every point where curves A and B come within tol: [{ s, t, d, p }]. */
function crossings(A, B, tol, KA, KB) {
const pa = polyline(A, KA), pb = polyline(B, KB);
const D = [];
for (let k = 0; k < pa.length - 1; k++) {
const row = [];
for (let l = 0; l < pb.length - 1; l++) row.push(segSeg(pa[k].p, pa[k + 1].p, pb[l].p, pb[l + 1].p));
D.push(row);
}
const out = [];
const lenA = (k) => V.len(V.sub(pa[k + 1].p, pa[k].p)), lenB = (l) => V.len(V.sub(pb[l + 1].p, pb[l].p));
for (let k = 0; k < D.length; k++) for (let l = 0; l < D[k].length; l++) {
const x = D[k][l];
if (x.d > tol + 0.25 * (lenA(k) + lenB(l))) continue;
let min = true;
for (let dk = -1; dk <= 1 && min; dk++) for (let dl = -1; dl <= 1; dl++) {
const y = D[k + dk]?.[l + dl];
if ((dk || dl) && y && (y.d < x.d || (y.d === x.d && (dk < 0 || (dk === 0 && dl < 0))))) { min = false; break; }
}
if (!min) continue;
const s0 = pa[k].t + (pa[k + 1].t - pa[k].t) * x.s, t0 = pb[l].t + (pb[l + 1].t - pb[l].t) * x.t;
const r = refine(A, B, s0, t0);
if (r.d > tol) continue;
const cyc = (c, u, w) => { const d = Math.abs(u - w); return c.closed ? Math.min(d, 1 - d) : d; };
if (out.some((o) => (cyc(A, o.s, r.s) < 1e-7 && cyc(B, o.t, r.t) < 1e-7) || V.len(V.sub(o.p, r.p)) < Math.max(tol, 1e-9) * 1e-2)) continue;
out.push(r);
}
return out;
}
Common parameters of a curve family
64-bit
The shared parameter at which every profile (or guide) meets each curve of the other family.
t_q = mean over the curves of (x_q − x_first) ÷ (x_last − x_first)
Each curve's crossing parameters are normalized by its own span and averaged, dividing by the curve count inside the sum. The ends are set to exactly 0 and 1. A span not above 1e-9 or neighbors not more than 1e-6 apart stop the build with an error.
gordon.js
lines 571 to 581
PartCraft 0.8.1
const common = (rows, what) => {
const m = rows[0].length, t = new Array(m).fill(0);
for (const r of rows) {
const a = r[0], b = r[m - 1];
if (!(b - a > 1e-9)) fail(tr(what === 'profiles' ? 'model.gordon.err.same_point_profiles' : 'model.gordon.err.same_point_guides'));
r.forEach((x, q) => { t[q] += (x - a) / (b - a) / rows.length; });
}
t[0] = 0; t[m - 1] = 1;
for (let q = 1; q < m; q++) if (!(t[q] - t[q - 1] > 1e-6)) fail(tr(what === 'profiles' ? 'model.gordon.err.too_close_profiles' : 'model.gordon.err.too_close_guides'));
return t;
};
Monotone reparametrization
64-bit
A smooth increasing map from the common parameters to one curve's own crossing parameters.
m_k = (w₁ + w₂) ÷ (w₁ ÷ d_(k−1) + w₂ ÷ d_k); cubic Hermite on each interval
Slopes are the weighted harmonic mean of the neighboring secant slopes, zero where they change sign. The cubic is evaluated in the power form of the Hermite basis in 64-bit numbers. Its inverse (invert) is 60 bisection steps on 0 to 1.
gordon.js
lines 404 to 421
PartCraft 0.8.1
/** Monotone piecewise cubic through (x_k, y_k), x and y strictly increasing (Fritsch–Carlson / Fritsch–Butland slopes). */
function pchip(x, y) {
const n = x.length, h = [], d = [];
for (let k = 0; k < n - 1; k++) { h.push(x[k + 1] - x[k]); d.push((y[k + 1] - y[k]) / h[k]); }
const m = new Array(n);
m[0] = d[0]; m[n - 1] = d[n - 2];
for (let k = 1; k < n - 1; k++) {
if (d[k - 1] * d[k] <= 0) { m[k] = 0; continue; }
const w1 = 2 * h[k] + h[k - 1], w2 = h[k] + 2 * h[k - 1];
m[k] = (w1 + w2) / (w1 / d[k - 1] + w2 / d[k]);
}
return (t) => {
let k = 0;
while (k < n - 2 && t > x[k + 1]) k++;
const s = (t - x[k]) / h[k], s2 = s * s, s3 = s2 * s;
return (2 * s3 - 3 * s2 + 1) * y[k] + (s3 - 2 * s2 + s) * h[k] * m[k] + (-2 * s3 + 3 * s2) * y[k + 1] + (s3 - s2) * h[k] * m[k + 1];
};
}
Distance from a point to a curve
64-bit
How far a point of the fitted surface or curve is from the true curve near a given parameter range.
golden-section search for the minimum of |g(t) − q|
The search range is the given range widened by half its width on each side and clamped to 0…1. It runs 60 steps or until the range is under 1e-15, and returns a local minimum of the distance inside that range, not a global one.
gordon.js
lines 424 to 436
PartCraft 0.8.1
/** Distance from q to the curve g between parameters a and b (golden section on |g(t) − q|). */
function distNear(g, q, a, b) {
const R = (Math.sqrt(5) - 1) / 2;
let lo = Math.min(a, b), hi = Math.max(a, b);
const pad = (hi - lo) * 0.5; lo = Math.max(0, lo - pad); hi = Math.min(1, hi + pad);
let x1 = hi - R * (hi - lo), x2 = lo + R * (hi - lo);
let f1 = V.len(V.sub(g(x1), q)), f2 = V.len(V.sub(g(x2), q));
for (let i = 0; i < 60 && hi - lo > 1e-15; i++) {
if (f1 < f2) { hi = x2; x2 = x1; f2 = f1; x1 = hi - R * (hi - lo); f1 = V.len(V.sub(g(x1), q)); }
else { lo = x1; x1 = x2; f1 = f2; x2 = lo + R * (hi - lo); f2 = V.len(V.sub(g(x2), q)); }
}
return Math.min(f1, f2);
}
Making a curve family compatible
64-bit
Gives every profile (or every guide) the same degree and knots, exactly when possible and by fitting otherwise.
exact: raise degree, merge knots; otherwise: cubic interpolation of C(r(s)), points doubled until error ≤ fitTol
The exact path is taken when every curve is an open polynomial B-spline and its crossing parameters equal the common ones within 1e-10; the reported fit error is then 0. Otherwise each curve is interpolated by one cubic and the error is measured at three points per span (0.25, 0.5, 0.75), so it is a sampled maximum. The fit stops with an error message if it would need more than 6000 points.
gordon.js
lines 440 to 493
PartCraft 0.8.1
/**
* One family (profiles or guides) made compatible: { p, U, P: [poles of each curve], approx, err, maps }.
* cs: oriented curves { g(σ) on [0, 1], bs, breaks (in σ), closed }; sig[i][j] = σ where curve i meets column j;
* tbar[j] = the common parameters.
*/
function family(cs, sig, tbar, fitTol, what) {
const aligned = cs.every((c, i) => sig[i].every((x, j) => Math.abs(x - tbar[j]) < 1e-10));
if (aligned && cs.every((c) => c.bs && !c.closed)) {
const d = Math.max(...cs.map((c) => c.bs.p));
const el = cs.map((c) => elevate(c.bs, d));
const U = unionKnots(el.map((c) => c.U), d);
return { p: d, U, P: el.map((c) => projector(d, c.U, c.P.length, d, U)(c.P)), approx: false, err: 0, maps: cs.map(() => (s) => s) };
}
const maps = cs.map((c, i) => (aligned ? (s) => s : pchip(tbar, sig[i])));
const bounds = new Set([0, 1]);
cs.forEach((c, i) => { for (const b of c.breaks) if (b > 1e-9 && b < 1 - 1e-9) bounds.add(aligned ? b : invert(maps[i], b)); });
const bl = [...bounds].sort((a, b) => a - b).filter((b, k, l) => k === 0 || b - l[k - 1] > 1e-9);
const anchors = [...new Set([...bl, ...tbar])].sort((a, b) => a - b).filter((b, k, l) => k === 0 || b - l[k - 1] > 1e-9);
const isBound = (s) => bl.some((b) => Math.abs(b - s) <= 1e-9);
for (let nsub = 4; ; nsub *= 2) {
const params = [];
for (let k = 0; k < anchors.length - 1; k++) for (let q = 0; q < nsub; q++) params.push(anchors[k] + (anchors[k + 1] - anchors[k]) * q / nsub);
params.push(1);
if (params.length > MAX_FIT) fail(tr(what === 'profiles' ? 'model.gordon.err.fit_profiles' : 'model.gordon.err.fit_guides', { max: MAX_FIT, tol: fitTol.toExponential(1) }));
// knots: averaged within each run between corners, full multiplicity at corners
const U = [0, 0, 0, 0];
let run = [params[0]];
for (let k = 1; k < params.length; k++) {
run.push(params[k]);
if (isBound(params[k]) || k === params.length - 1) {
const Kr = averagedKnots(run, 3);
U.push(...Kr.slice(4, Kr.length - 4));
if (k < params.length - 1) U.push(params[k], params[k], params[k]);
run = [params[k]];
}
}
U.push(1, 1, 1, 1);
const solve = bandLU(3, U, params.length, params);
let err = 0;
const P = cs.map((c, i) => {
const r = maps[i];
const poles = solve(params.map((s) => c.g(r(s))));
const bsC = { p: 3, U, P: poles };
for (let k = 0; k < params.length - 1; k++) for (const w of [0.25, 0.5, 0.75]) {
const s = params[k] + (params[k + 1] - params[k]) * w, q = evalCurve(bsC, s);
let e = V.len(V.sub(q, c.g(r(s))));
if (e > fitTol) e = distNear(c.g, q, r(params[k]), r(params[k + 1]));
err = Math.max(err, e);
}
return poles;
});
if (err <= fitTol) return { p: 3, U, P, approx: true, err, maps, points: params.length };
}
}
Boolean sum of the three surfaces
64-bit
Adds the surface through the profiles and the surface through the guides and subtracts the surface through the grid points.
P = Lu + Lv − T, pole by pole
The three surfaces are first brought to one degree and one knot vector by the band solves above; the poles are then combined as (Lu + Lv) − T in 64-bit numbers. Where the three surfaces nearly agree the subtraction cancels leading digits of numbers that each carry the rounding of their own solves.
gordon.js
lines 607 to 608
PartCraft 0.8.1
const [A, B, C] = [Lu, Lv, Tt].map((s) => projectSurface(s, du, Uc, dv, Vc));
let S = { pu: du, pv: dv, U: Uc, V: Vc, P: A.P.map((row, i) => row.map((q, j) => V.sub(V.add(q, B.P[i][j]), C.P[i][j]))) };
Tangency to a plane along a boundary
64-bit
Makes the surface leave a boundary curve tangent to a plane and reports the angle it made with the plane before.
second pole row: q ← a + (d − n × (d·n)), d = q − a; angle = asin(|d·n| ÷ |d|)
The boundary is chosen by its largest distance to the plane at 65 samples, which must be within tol. The reported angle uses a one-sided difference of step 1e-6 in the parameter at 33 samples; the asin argument is capped at 1. The second row of poles is projected into the plane in 64-bit numbers, after an inserted knot placed 0.3 of the way to the next curve.
gordon.js
lines 187 to 212
PartCraft 0.8.1
export function tangentToPlane(S, plane, params, { tol = 0.01, closedU = false, closedV = false, skip = [] } = {}) {
const n = V.norm(plane.n), dist = (q) => Math.abs(V.dot(V.sub(q, plane.o), n));
const sides = [];
if (!closedV) sides.push(['v0', (t) => evalSurface(S, t, 0)], ['v1', (t) => evalSurface(S, t, 1)]);
if (!closedU) sides.push(['u0', (t) => evalSurface(S, 0, t)], ['u1', (t) => evalSurface(S, 1, t)]);
let best = null;
for (const [side, f] of sides) { if (skip.includes(side)) continue; let d = 0; for (let k = 0; k <= 64; k++) d = Math.max(d, dist(f(k / 64))); if (!best || d < best.d) best = { side, d }; }
if (!best || best.d > tol) fail(t('gordon.err.no_boundary', { d: fmt(best ? best.d : Infinity), tol }));
const { side } = best, dir = side[0], end = side[1] === '1';
// the angle the surface makes with the plane along that boundary now
const deriv = (T, t) => { const h = 1e-6, a = dir === 'v' ? [t, end ? 1 - h : h] : [end ? 1 - h : h, t], b = dir === 'v' ? [t, end ? 1 : 0] : [end ? 1 : 0, t]; return V.sub(evalSurface(T, a[0], a[1]), evalSurface(T, b[0], b[1])); };
let angle = 0;
for (let k = 0; k <= 32; k++) { const d = deriv(S, k / 32), l = V.len(d); if (l > 0) angle = Math.max(angle, Math.asin(Math.min(1, Math.abs(V.dot(d, n)) / l)) * 180 / Math.PI); }
// localise: a knot 0.3 of the way to the next curve of that family, unless one is that close already
const par = dir === 'v' ? params.v : params.u, K = dir === 'v' ? S.V : S.U, deg = dir === 'v' ? S.pv : S.pu;
const next = end ? par[par.length - 2] : par[1], tk = end ? 1 - 0.3 * (1 - next) : 0.3 * next;
const firstInner = end ? K[K.length - deg - 2] : K[deg + 1];
let T = S;
if (end ? firstInner < tk - 1e-9 : firstInner > tk + 1e-9) T = withKnot(S, dir, tk);
// the second row (or column) of poles into the plane
const P = T.P.map((row) => row.map((q) => q.slice()));
const into = (a, b) => { const d = V.sub(b, a); return V.add(a, V.sub(d, V.mul(n, V.dot(d, n)))); };
if (dir === 'v') { const j0 = end ? P[0].length - 1 : 0, j1 = end ? j0 - 1 : 1; for (const row of P) row[j1] = into(row[j0], row[j1]); }
else { const i0 = end ? P.length - 1 : 0, i1 = end ? i0 - 1 : 1; P[i1] = P[i1].map((q, j) => into(P[i0][j], q)); }
return { surface: { ...T, P }, side, angle };
}
Deviation from the input curves
64-bit
The largest distance between the finished surface and any profile or guide, shown in the panel.
dev = max over 401 samples per curve of |S(s) − C(r(s))|
Each curve is compared with its iso-curve of the surface at 401 evenly spaced parameters; a sample over the fit tolerance is replaced by the closest-point distance within one sample step either side. The reported value is a sampled maximum. The panel text formats it to 3 significant digits (fmt), which is display only.
gordon.js
lines 620 to 630
PartCraft 0.8.1
let dev = 0;
const check = (c, r, iso, isU) => {
for (let k = 0; k <= 400; k++) {
const s = k / 400, q = isU ? evalSurface(S, s, iso) : evalSurface(S, iso, s);
let e = V.len(V.sub(q, c.g(r(s))));
if (e > ftol) e = distNear(c.g, q, r(Math.max(0, s - 1 / 400)), r(Math.min(1, s + 1 / 400)));
dev = Math.max(dev, e);
}
};
pc.forEach((c, i) => check(c, FP.maps[i], vbar[i], true));
gc.forEach((c, j) => check(c, FG.maps[j], ubar[j], false));
Fold check
64-bit
Refuses a surface whose normal flips between neighboring sample points.
n = (S(u + h, v) − S(u − h, v)) × (S(u, v + h) − S(u, v − h)), h = 1e-7; fold if n₁·n₂ < 0
Normals come from central differences of step 1e-7 in the parameter, up to 4 samples per knot span. The decision is the sign of a 64-bit dot product with no tolerance; a pair where either normal has zero length is not judged.
gordon.js
lines 639 to 661
PartCraft 0.8.1
/**
* Refuse a surface that folds over: its normal S_u × S_v, sampled up to 4 times per knot span each way, must not turn by more
* than 90° between neighbouring samples (pairs across a corner knot — multiplicity = degree — are not compared).
*/
function foldCheck(S) {
const ts = (U, p) => {
const { k, m } = knotsToMults(U), out = [], per = Math.max(1, Math.min(4, Math.floor(160 / (k.length - 1))));
for (let q = 0; q < k.length - 1; q++) for (let w = 0; w < per; w++) out.push({ t: k[q] + (k[q + 1] - k[q]) * (w + 0.5) / per, cut: w === 0 && q > 0 && m[q] >= p });
return out;
};
const us = ts(S.U, S.pu), vs = ts(S.V, S.pv), h = 1e-7;
const N = us.map(({ t: u }) => vs.map(({ t: v }) => {
const su = V.sub(evalSurface(S, Math.min(1, u + h), v), evalSurface(S, Math.max(0, u - h), v));
const sv = V.sub(evalSurface(S, u, Math.min(1, v + h)), evalSurface(S, u, Math.max(0, v - h)));
return V.cross(su, sv);
}));
const bad = (a, b) => { const la = V.len(a), lb = V.len(b); return la > 0 && lb > 0 && V.dot(a, b) < 0; };
for (let i = 0; i < us.length; i++) for (let j = 0; j < vs.length; j++) {
if ((i + 1 < us.length && !us[i + 1].cut && bad(N[i][j], N[i + 1][j])) || (j + 1 < vs.length && !vs[j + 1].cut && bad(N[i][j], N[i][j + 1]))) {
fail(tr('model.gordon.err.folds'));
}
}
}
In detail
Method
Gordon surface as a Boolean sum S = Lu + Lv − T (the file header cites Gordon 1969). Skinning by global interpolation with averaged knots (the header cites the NURBS Book, eq. 9.8); crossings by closest segments (the code cites Ericson, Real-Time Collision Detection 5.1.9) then Gauss–Newton; reparametrization by a monotone cubic (Fritsch–Carlson / Fritsch–Butland slopes).
Number types
64-bit numbers throughout; Float64Array for the band matrix of each interpolation.
Tolerances
Crossing gap tol 0.01 (model units, absolute, an option). Fit tolerance 1e-7 × network size unless given. Knot equality 1e-10 (knotsToMults), 1e-9 (unionKnots), 1e-12 (sameSpace), all absolute on the [0, 1] parameter range. Pivot test |pivot| > 1e-13 absolute. Curve end joining 1e-6 × the largest sketch coordinate. Parameter spacing between neighboring curves must exceed 1e-6. Finite-difference steps 1e-7 (curve derivative, fold check) and 1e-6 (boundary angle) in parameter.
Limits
The fit error and the reported deviation are maxima over sample points, not proven bounds. The band factorization does not pivot. Crossings start from polylines of 96 segments plus 16 per corner, so two crossings closer than the sampling can be merged or missed.
Approach: 64-bit, stated limit
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.