NURBS lengths, areas, closest points and tessellation
Measures a curve's length and a surface's area, finds the point of a curve or surface nearest a given point, and turns curves and surfaces into polylines and triangle grids within a tolerance.
How the precision is kept
The arithmetic is 64-bit, but each result comes from a numerical method with its own accuracy: Gauss–Legendre quadrature with a fixed number of points, Newton iteration from sampled starting points, or subdivision to a chord tolerance. None of the quadratures estimates its own error.
Open item. Report when a tessellation stops at its level, segment or per-span limit before the tolerance is met.
Calculations under this policy
10 calculations in 3 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.
The parameter, point and distance of the point of a curve nearest a given point.
f(u) = (C − P)·C′; u ← u − f ÷ (C″·(C − P) + C′·C′)
The start is the nearest of 4 × degree + 5 samples on each knot span. Newton runs at most 30 steps, clamped to the domain, and stops when the step is 1e-15 × (domain length + 1) or less. If the result is farther away than the best sample, the sample is returned.
nurbs/curve.js
lines 343 to 365
ModelCraft 0.8.1
export function closestParam(c, p, { u0 = null } = {}) {
const [a, b] = domain(c);
let best = u0, bd = u0 == null ? Infinity : V.dist(point(c, u0), p);
if (u0 == null) {
const ks = distinctKnots(c.knots), per = 4 * c.deg + 4;
for (let s = 0; s + 1 < ks.length; s++) for (let i = 0; i <= per; i++) {
const u = ks[s].u + (ks[s + 1].u - ks[s].u) * i / per, d = V.dist(point(c, u), p);
if (d < bd) { bd = d; best = u; }
}
}
let u = best;
for (let it = 0; it < 30; it++) {
const [C, d1, d2] = derivs(c, u, 2), r = V.sub(C, p);
const f = V.dot(d1, r), df = V.dot(d2, r) + V.dot(d1, d1);
if (Math.abs(df) < 1e-300) break;
let un = u - f / df;
un = un < a ? a : un > b ? b : un;
if (Math.abs(un - u) <= 1e-15 * (b - a + 1)) { u = un; break; }
u = un;
}
const q = point(c, u), d = V.dist(q, p);
return d <= bd ? { u, point: q, d } : { u: best, point: point(c, best), d: bd };
}
Length of a curve
64-bit
The length of a curve, or of a part of it between two parameters.
s = Σ pieces Σ₈ wᵍ·h·|C′(x₀ + h·(1 + xᵍ))|, 4 pieces per knot span
Gauss–Legendre with 8 points per piece; the nodes and weights are written to 16 digits. The speed |C′| is a square root and not a polynomial, so the sum is an approximation whose error depends on the curve; there is no error estimate and no refinement. Terms are added in order.
nurbs/curve.js
lines 367 to 382
ModelCraft 0.8.1
const GL_X = [-0.9602898564975363, -0.7966664774136267, -0.5255324099163290, -0.1834346424956498, 0.1834346424956498, 0.5255324099163290, 0.7966664774136267, 0.9602898564975363];
const GL_W = [0.1012285362903763, 0.2223810344533745, 0.3137066458778873, 0.3626837833783620, 0.3626837833783620, 0.3137066458778873, 0.2223810344533745, 0.1012285362903763];
/** Gauss–Legendre (8 points) on 4 pieces of each knot span. */
export function arcLength(c, a = null, b = null) {
const [ca, cb] = domain(c);
a = a == null ? ca : a; b = b == null ? cb : b;
if (b < a) return arcLength(c, b, a);
const cuts = [a, ...distinctKnots(c.knots).map((k) => k.u).filter((u) => u > a && u < b), b];
let s = 0;
for (let k = 0; k + 1 < cuts.length; k++) for (let q = 0; q < 4; q++) {
const x0 = cuts[k] + (cuts[k + 1] - cuts[k]) * q / 4, x1 = cuts[k] + (cuts[k + 1] - cuts[k]) * (q + 1) / 4, h = (x1 - x0) / 2;
for (let g = 0; g < 8; g++) s += GL_W[g] * h * V.len(derivs(c, x0 + h * (1 + GL_X[g]), 1)[1]);
}
return s;
}
Parameter at a length
64-bit
The parameter at which the curve has run a given length from its start.
f(u) = length(a … u) − s; u ← u − f ÷ |C′(u)|, or the middle of the bracket
Newton steps kept inside a bracket, falling back to bisection, for at most 50 steps; it stops when |f| is within 1e-12 of the total length. Each step calls arcLength again, so the result carries that quadrature's error.
nurbs/curve.js
lines 384 to 400
ModelCraft 0.8.1
/** The parameter at arc length s from the start (Newton on the length, bracketed). */
export function paramAtLength(c, s) {
const [a, b] = domain(c), total = arcLength(c);
if (s <= 0) return a;
if (s >= total) return b;
let lo = a, hi = b, u = a + (b - a) * s / total;
for (let it = 0; it < 50; it++) {
const f = arcLength(c, a, u) - s;
if (Math.abs(f) <= 1e-12 * total) break;
if (f > 0) hi = u; else lo = u;
const sp = V.len(derivs(c, u, 1)[1]);
let un = sp > 1e-300 ? u - f / sp : (lo + hi) / 2;
if (!(un > lo && un < hi)) un = (lo + hi) / 2;
u = un;
}
return u;
}
Curve to polyline
64-bit
Points along a curve such that the curve stays within a tolerance of the chords between them.
a segment is kept when the curve at 1/4, 1/2 and 3/4 of its parameter range lies within tol of the chord; else it is halved
The deviation is tested at three points per segment, not along all of it. Halving stops at 20 levels or when the segment budget (4096 by default) runs out, and the segment is then kept whether or not it meets the tolerance (default 0.01).
Decides whether one curve lies on another within a tolerance.
33 points of a, each within tol of its closest point on b
A sampled test in one direction: points of the first curve are checked against the second, evenly in the first curve's parameter (samples = 32, tolerance 1e-6 by default). Parts of the second curve beyond the first are not checked.
nurbs/curve.js
lines 432 to 439
ModelCraft 0.8.1
export function sameCurve(a, b, tol = 1e-6, samples = 32) {
const [a0, a1] = domain(a);
for (let i = 0; i <= samples; i++) {
const p = point(a, a0 + (a1 - a0) * i / samples);
if (closestParam(b, p).d > tol) return false;
}
return true;
}
A = Σ cells Σ₅ₓ₅ wₐ·w_b·h_u·h_v·|S_u × S_v|, 2 × 2 pieces per knot-span cell
Gauss–Legendre with 5 × 5 points per piece, nodes and weights written to 16 digits. The integrand is the length of a cross product and not a polynomial, so the sum is an approximation with no error estimate. Terms are added in order.
nurbs/surface.js
lines 308 to 324
ModelCraft 0.8.1
const GX = [-0.9061798459386640, -0.5384693101056831, 0, 0.5384693101056831, 0.9061798459386640];
const GW = [0.2369268850561891, 0.4786286704993665, 0.5688888888888889, 0.4786286704993665, 0.2369268850561891];
/** Area by 5 × 5 Gauss–Legendre on 2 × 2 pieces of every knot-span cell. */
export function area(s) {
const us = distinctKnots(s.knotsU).map((k) => k.u), vs = distinctKnots(s.knotsV).map((k) => k.u);
let A = 0;
const sub = (lst) => { const o = []; for (let i = 0; i + 1 < lst.length; i++) o.push([lst[i], (lst[i] + lst[i + 1]) / 2], [(lst[i] + lst[i + 1]) / 2, lst[i + 1]]); return o; };
for (const [u0, u1] of sub(us)) for (const [v0, v1] of sub(vs)) {
const hu = (u1 - u0) / 2, hv = (v1 - v0) / 2;
for (let a = 0; a < 5; a++) for (let b = 0; b < 5; b++) {
const D = derivs(s, u0 + hu * (1 + GX[a]), v0 + hv * (1 + GX[b]), 1);
A += GW[a] * GW[b] * hu * hv * V.len(V.cross(D[1][0], D[0][1]));
}
}
return A;
}
Closest point on a surface
64-bit
The (u, v), point and distance of the point of a surface nearest a given point.
f = (S − P)·S_u, g = (S − P)·S_v; Newton on the 2 × 2 system with J = [S_u·S_u + r·S_uu, S_u·S_v + r·S_uv; ·, S_v·S_v + r·S_vv]
The start is the nearest of (degU + 3) × (degV + 3) samples per knot cell, leaving out samples on a pole unless they are the point itself. Newton runs at most 40 steps; a step that makes the distance grow is halved once; it stops when the move is below 1e-15 of the domain. A direction whose two edges coincide within 1e-9 × the surface's size wraps round instead of stopping at the edge.
nurbs/surface.js
lines 326 to 376
ModelCraft 0.8.1
export function closestPoint(s, p, { uv0 = null } = {}) {
const [ua, ub] = domainU(s), [va, vb] = domainV(s);
let bu, bv, bd = Infinity;
if (uv0) { [bu, bv] = uv0; bd = V.dist(point(s, bu, bv), p); }
else {
const us = distinctKnots(s.knotsU).map((k) => k.u), vs = distinctKnots(s.knotsV).map((k) => k.u);
const pu = s.degU + 2, pv = s.degV + 2;
// an iso-curve that is one point (a pole: a revolved profile on its axis) gives the same point for every value along it, so as
// a first guess it says nothing about that parameter and Newton then settles on the wrong side; such samples only count when
// they are the point itself
const size = Math.max(1e-300, V.dist(point(s, ua, va), point(s, ub, vb)), V.dist(point(s, ua, vb), point(s, ub, va)));
const one = (q0, q1, q2) => V.dist(q0, q1) < 1e-12 * size && V.dist(q0, q2) < 1e-12 * size;
const poleU = (u) => one(point(s, u, va), point(s, u, (va + vb) / 2), point(s, u, vb));
const poleV = (v) => one(point(s, ua, v), point(s, (ua + ub) / 2, v), point(s, ub, v));
const atPole = new Map();
const isPole = (u, v) => { const k = u + ',' + v; if (!atPole.has(k)) atPole.set(k, ((u === ua || u === ub) && poleU(u)) || ((v === va || v === vb) && poleV(v))); return atPole.get(k); };
for (let a = 0; a + 1 < us.length; a++) for (let i = 0; i <= pu; i++) {
const u = us[a] + (us[a + 1] - us[a]) * i / pu;
for (let b = 0; b + 1 < vs.length; b++) for (let j = 0; j <= pv; j++) {
const v = vs[b] + (vs[b + 1] - vs[b]) * j / pv, d = V.dist(point(s, u, v), p);
if (d < bd && (d < 1e-12 * size || !isPole(u, v))) { bd = d; bu = u; bv = v; }
}
}
if (bu === undefined) { bu = (ua + ub) / 2; bv = (va + vb) / 2; bd = V.dist(point(s, bu, bv), p); }
}
// a direction whose two ends are the same curve (a full revolution's seam) wraps instead of stopping at the end: a point just
// across the seam is reached by going round, not by sticking at the edge
const scale = Math.max(1e-300, V.dist(point(s, ua, va), point(s, ub, vb)), V.dist(point(s, ua, vb), point(s, ub, va)), V.dist(point(s, (ua + ub) / 2, va), point(s, (ua + ub) / 2, vb)));
const meets = (f) => [0, 0.37, 0.71, 1].every((x) => f(x) < 1e-9 * scale);
const wrapU = meets((x) => V.dist(point(s, ua, va + (vb - va) * x), point(s, ub, va + (vb - va) * x)));
const wrapV = meets((x) => V.dist(point(s, ua + (ub - ua) * x, va), point(s, ua + (ub - ua) * x, vb)));
const stepU = (x) => (wrapU ? ua + ((((x - ua) % (ub - ua)) + (ub - ua)) % (ub - ua)) : clamp(x, [ua, ub]));
const stepV = (x) => (wrapV ? va + ((((x - va) % (vb - va)) + (vb - va)) % (vb - va)) : clamp(x, [va, vb]));
let u = bu, v = bv;
for (let it = 0; it < 40; it++) {
const D = derivs(s, u, v, 2), r = V.sub(D[0][0], p), Su = D[1][0], Sv = D[0][1];
const f = V.dot(r, Su), g = V.dot(r, Sv);
const J00 = V.dot(Su, Su) + V.dot(r, D[2][0]), J01 = V.dot(Su, Sv) + V.dot(r, D[1][1]), J11 = V.dot(Sv, Sv) + V.dot(r, D[0][2]);
const det = J00 * J11 - J01 * J01;
if (Math.abs(det) < 1e-300) break;
let du = -(J11 * f - J01 * g) / det, dv = -(J00 * g - J01 * f) / det;
let un = stepU(u + du), vn = stepV(v + dv);
// keep the step only when it does not make things worse (Newton can wander near the edges)
if (V.dist(point(s, un, vn), p) > V.len(r) * (1 + 1e-12)) { un = stepU(u + du / 2); vn = stepV(v + dv / 2); }
const moved = Math.min(Math.abs(un - u), wrapU ? (ub - ua) - Math.abs(un - u) : Infinity) / (ub - ua) + Math.min(Math.abs(vn - v), wrapV ? (vb - va) - Math.abs(vn - v) : Infinity) / (vb - va);
u = un; v = vn;
if (moved < 1e-15) break;
}
const q = point(s, u, v), d = V.dist(q, p);
return d <= bd ? { u, v, point: q, d } : { u: bu, v: bv, point: point(s, bu, bv), d: bd };
}
Tessellation steps per knot span
64-bit
How many parameter steps each knot span gets so that the grid stays within the tolerance.
n = ceil(h·√(M ÷ (4·tol))), 2 ≤ n ≤ maxPerSpan, M = largest |second derivative| at 5 × 7 sample points
A chord over a step h deviates by at most h²·M ÷ 8 when M bounds the second derivative; here M is the largest of 35 sampled values, so it can be below the true maximum. The cap of 64 steps per span can also leave a span coarser than the tolerance asks.
nurbs/surface.js
lines 380 to 399
ModelCraft 0.8.1
function spanParams(s, dir, tol, maxPerSpan) {
const knots = dir === 'u' ? s.knotsU : s.knotsV, ks = distinctKnots(knots).map((k) => k.u);
const [oa, ob] = dir === 'u' ? domainV(s) : domainU(s);
const deg = dir === 'u' ? s.degU : s.degV;
const out = [ks[0]];
for (let k = 0; k + 1 < ks.length; k++) {
const a = ks[k], b = ks[k + 1], h = b - a;
let M = 0;
if (deg > 1) for (let i = 0; i <= 4; i++) for (let j = 0; j <= 6; j++) {
const t = a + h * i / 4, o = oa + (ob - oa) * j / 6;
const D = dir === 'u' ? derivs(s, t, o, 2)[2][0] : derivs(s, o, t, 2)[0][2];
M = Math.max(M, V.len(D));
}
// a chord over parameter step h deviates from the curve by at most h²·M/8; half the tolerance each way, so a cell's
// middle (both directions at once) stays within it
const n = deg > 1 ? Math.min(maxPerSpan, Math.max(2, Math.ceil(h * Math.sqrt(M / (4 * tol))))) : 1;
for (let i = 1; i <= n; i++) out.push(a + h * i / n);
}
return out;
}
Surface to triangle grid
64-bit
A grid of points, normals and parameters on the surface with two triangles per cell.
Every vertex is evaluated on the exact surface and its normal is the surface's own; positions, normals and parameters are Float64Arrays. The triangles between vertices are flat, so they follow the surface only within the tessellation tolerance.
nurbs/surface.js
lines 401 to 421
ModelCraft 0.8.1
export function tessellate(s, { tol = 0.01, maxPerSpan = 64, uParams = null, vParams = null, curvature: withK = false } = {}) {
const up = uParams || spanParams(s, 'u', tol, maxPerSpan), vp = vParams || spanParams(s, 'v', tol, maxPerSpan);
const nu = up.length, nv = vp.length;
const positions = new Float64Array(nu * nv * 3), normals = new Float64Array(nu * nv * 3), uvs = new Float64Array(nu * nv * 2);
const K = withK ? new Float64Array(nu * nv) : null, H = withK ? new Float64Array(nu * nv) : null;
for (let i = 0; i < nu; i++) for (let j = 0; j < nv; j++) {
const k = i * nv + j, q = point(s, up[i], vp[j]), n = normal(s, up[i], vp[j]);
positions.set(q, k * 3); normals.set(n, k * 3); uvs[k * 2] = up[i]; uvs[k * 2 + 1] = vp[j];
if (withK) { const c = curvature(s, up[i], vp[j]); K[k] = c.K; H[k] = c.H; }
}
const indices = new Uint32Array((nu - 1) * (nv - 1) * 6);
let t = 0;
for (let i = 0; i + 1 < nu; i++) for (let j = 0; j + 1 < nv; j++) {
const a = i * nv + j, b = a + nv, c = b + 1, d = a + 1;
indices[t++] = a; indices[t++] = b; indices[t++] = c;
indices[t++] = a; indices[t++] = c; indices[t++] = d;
}
const out = { positions, normals, uvs, indices, nu, nv, uParams: up, vParams: vp };
if (withK) { out.K = K; out.H = H; }
return out;
}
Each distance comes from surface.closestPoint, an iteration from sampled starts, and is kept in a Float64Array. The sums are added in point order.
nurbs/analysis.js
lines 61 to 67
ModelCraft 0.8.1
export function deviation(s, points) {
const d = new Float64Array(points.length);
let max = 0, sum = 0, sq = 0, at = 0;
points.forEach((p, i) => { const x = S.closestPoint(s, p).d; d[i] = x; sum += x; sq += x * x; if (x > max) { max = x; at = i; } });
const n = Math.max(1, points.length);
return { max, mean: sum / n, rms: Math.sqrt(sq / n), at, d };
}
In detail
Method
Length: 8-point Gauss–Legendre on 4 pieces of each knot span. Area: 5 × 5 Gauss–Legendre on 2 × 2 pieces of each knot cell. Closest point: Newton on (C − P)·C′ = 0, or on the 2 × 2 system for a surface (The NURBS Book § 6.1, as cited). Tessellation: halving to a chord tolerance for curves; for surfaces a per-span count from a second-derivative estimate.
Number types
64-bit numbers; quadrature nodes and weights written to 16 digits; Float64Array positions, normals and parameters; Uint32Array indices.
Tolerances
Newton stops at a step of 1e-15 of the domain (curves: 1e-15 × (domain + 1)), after 30 steps for curves and 40 for surfaces, or when the derivative or determinant is below 1e-300. Length to parameter: 1e-12 of the total length, 50 steps. Tessellation tolerance 0.01 by default; curves stop at 20 levels or 4096 segments, surfaces at 64 steps per knot span.
Limits
Closest point finds the minimum near the best sample; a nearer point elsewhere between samples can be missed. When a tessellation limit is reached the tolerance is not met and nothing says so. Surface step counts rest on a sampled, not a proven, bound of the second derivative.
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.