Resources Numerical accuracy Calculation policy
Turns the 3D solids stored in DWG and DXF files (ACIS data) into triangle meshes: shown on CadCraft's 3D tab, and brought into ModelCraft as groups when such a file is imported.
Surfaces and curves are evaluated in 64-bit numbers and the mesh positions are kept in a Float64Array. The accuracy of the result is set by the sampling: curved edges and faces are cut at steps of 7.5 degrees (or a chord height of 0.0002 of the model size, whichever is finer), and spline surfaces at steps of one sixteenth of their parameter range. Edge points must lie on their curve and surface within a tolerance of at least 1e-5 of the model size, otherwise the face is counted as not drawn.
Open item. Sum the volume about the mesh's own center instead of the coordinate origin, so a small solid far from the origin keeps its digits.
26 calculations in 2 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.
Tests tools/acis/test.mjs, tools/acis/test-apps.mjs
A direction scaled to length 1, or nothing when it has no length.
l = √(x² + y² + z²); l > 1e-300 → (x, y, z) ÷ l
Math.hypot for the length, then three divisions. The guard is 1e-300, so it rejects only a vector that is zero or nearly underflowed, not a short one.
export const unit = (a) => { const l = len(a); return l > 1e-300 ? [a[0] / l, a[1] / l, a[2] / l] : null; };A point of a plane from its (u, v), and back.
P = O + U·u + V·v; u = (p − O)·U; v = (p − O)·V
U is the file's u direction made perpendicular to the unit normal, and V = N × U. Evaluation and inverse are sums of three products each, with no tolerance.
case 'plane': {
const N = unit(s.normal);
if (!N) return null;
const U = ortho(N, s.udir && unit(s.udir)), V = cross(N, U), O = s.root;
return {
type: 'plane', curved: false, periodicU: 0, periodicV: 0, poles: null, scaleU: 1, scaleV: 1,
eval: (u, v) => [O[0] + U[0] * u + V[0] * v, O[1] + U[1] * u + V[1] * v, O[2] + U[2] * u + V[2] * v],
inv: (p) => { const d = sub(p, O); return [dot(d, U), dot(d, V)]; },
normal: () => N,
};
}A point of a cone or cylinder (circular or elliptical base) from its angle and its height along the axis, and back.
t = ± sin ÷ cos; q = 1 + h·t ÷ R; P = C + X·R·q·cos u + Y·R_m·q·sin u + N·h; u = atan2((d·Y) ÷ R_m, (d·X) ÷ R); apex at h = −R ÷ t
The surface is refused when the radius is not above 1e-300, or the ratio or the cosine is not above 1e-12 in size. It is treated as a cylinder when |t| < 1e-12. The sign of the half angle is tried both ways by the caller, because the file format's convention is not written down (the comment says so).
case 'cone': {
// base ellipse: centre C, axis N, major axis M (|M| = R), minor = ratio·R along N × M; the radius grows by
// tan(angle) per unit along the axis — with `sign` = ±1 because the angle's sign convention is a spec gap (§ 6.3)
const N = unit(s.normal), R = len(s.major);
if (!N || !(R > 1e-300) || !(Math.abs(s.ratio) > 1e-12)) return null;
const X = ortho(N, unit(s.major)), Y = cross(N, X), C = s.centre, Rm = R * s.ratio;
if (!(Math.abs(s.cos) > 1e-12)) return null;
const t = sign * s.sin / s.cos;
const k = (h) => 1 + (h * t) / R;
const cyl = Math.abs(t) < 1e-12;
const apex = cyl ? null : -R / t;
return {
type: 'cone', curved: true, periodicU: TAU, periodicV: 0, scaleU: R, scaleV: 1,
poles: cyl ? null : (apex < 0 ? [apex, null] : [null, apex]),
eval: (u, h) => { const q = k(h), c = Math.cos(u) * R * q, d = Math.sin(u) * Rm * q; return [C[0] + X[0] * c + Y[0] * d + N[0] * h, C[1] + X[1] * c + Y[1] * d + N[1] * h, C[2] + X[2] * c + Y[2] * d + N[2] * h]; },
inv: (p) => { const d = sub(p, C), h = dot(d, N); return [Math.atan2(dot(d, Y) / Rm, dot(d, X) / R), h]; },
// outward: away from the axis
normal: (u, h) => {
const q = k(h);
const Pu = add(mul(X, -Math.sin(u) * R * q), mul(Y, Math.cos(u) * Rm * q));
const Ph = add(N, mul(add(mul(X, Math.cos(u) * R), mul(Y, Math.sin(u) * Rm)), t / R));
let n = unit(cross(Pu, Ph));
if (!n) return N;
const radial = add(mul(X, Math.cos(u)), mul(Y, Math.sin(u)));
if (dot(n, radial) * Math.sign(q || 1) < 0) n = mul(n, -1);
return n;
},
radiusAt: (h) => R * Math.abs(k(h)),
};
}A point of a sphere from longitude and latitude, and back.
P = C + r·cos v·(X·cos u + Y·sin u) + Z·r·sin v; u = atan2(y, x); v = atan2(z, √(x² + y²))
The latitude is taken with atan2 of z and the horizontal distance, which keeps its digits near the poles where an arcsine would not. A point counts as a pole when its horizontal distance is below 1e-9 × r.
case 'sphere': {
const r = Math.abs(s.radius);
if (!(r > 1e-300)) return null;
const Z = unit(s.pole || [0, 0, 1]) || [0, 0, 1];
const X = ortho(Z, s.origin && unit(s.origin)), Y = cross(Z, X), C = s.centre;
const inside = s.radius < 0 ? -1 : 1;
return {
type: 'sphere', curved: true, periodicU: TAU, periodicV: 0, scaleU: r, scaleV: r, poles: [-Math.PI / 2, Math.PI / 2],
eval: (u, v) => { const cv = Math.cos(v) * r, a = Math.cos(u) * cv, b = Math.sin(u) * cv, z = Math.sin(v) * r; return [C[0] + X[0] * a + Y[0] * b + Z[0] * z, C[1] + X[1] * a + Y[1] * b + Z[1] * z, C[2] + X[2] * a + Y[2] * b + Z[2] * z]; },
inv: (p) => { const d = sub(p, C), x = dot(d, X), y = dot(d, Y), z = dot(d, Z); return [Math.atan2(y, x), Math.atan2(z, Math.hypot(x, y))]; },
normal: (u, v) => mul(add(add(mul(X, Math.cos(u) * Math.cos(v)), mul(Y, Math.sin(u) * Math.cos(v))), mul(Z, Math.sin(v))), inside),
atPole: (p) => { const d = sub(p, C); return Math.hypot(dot(d, X), dot(d, Y)) < r * 1e-9; },
};
}A point of a torus from its two angles, and back.
w = A + B·cos v; P = C + X·w·cos u + Y·w·sin u + N·B·sin v; u = atan2(y, x); v = atan2(z, √(x² + y²) − A)
Both angles come back through atan2. The surface is refused when the minor radius is not above 1e-300.
case 'torus': {
const N = unit(s.normal), A = Math.abs(s.major), B = Math.abs(s.minor);
if (!N || !(B > 1e-300)) return null;
const X = ortho(N, s.origin && unit(s.origin)), Y = cross(N, X), C = s.centre;
const inside = s.minor < 0 ? -1 : 1;
return {
type: 'torus', curved: true, periodicU: TAU, periodicV: TAU, scaleU: A + B, scaleV: B, poles: null,
eval: (u, v) => { const w = A + B * Math.cos(v), a = Math.cos(u) * w, b = Math.sin(u) * w, z = B * Math.sin(v); return [C[0] + X[0] * a + Y[0] * b + N[0] * z, C[1] + X[1] * a + Y[1] * b + N[1] * z, C[2] + X[2] * a + Y[2] * b + N[2] * z]; },
inv: (p) => { const d = sub(p, C), x = dot(d, X), y = dot(d, Y), z = dot(d, N), q = Math.hypot(x, y); return [Math.atan2(y, x), Math.atan2(z, q - A)]; },
normal: (u, v) => mul(add(add(mul(X, Math.cos(u) * Math.cos(v)), mul(Y, Math.sin(u) * Math.cos(v))), mul(N, Math.sin(v))), inside),
};
}A point of a B-spline surface from (u, v), its normal, and the (u, v) nearest to a given point.
first guess: nearest of 25 × 25 samples; then Gauss–Newton on |p − S(u, v)|² with partial derivatives by central differences of 1e-6 of the range
The inverse stops after 30 steps or when both steps are below 1e-12 of the parameter range, and stops early if the 2 × 2 determinant is not above 1e-300. The partial derivatives and the normal are finite differences, so they carry the truncation of a 1e-6 step and not full 64-bit digits. Parameters are clamped to the knot range.
case 'spline': {
const bs = s.bs;
if (!bs || bs.degV === undefined) return null;
const u0 = bs.knotsU[bs.degU], u1 = bs.knotsU[bs.knotsU.length - 1 - bs.degU];
const v0 = bs.knotsV[bs.degV], v1 = bs.knotsV[bs.knotsV.length - 1 - bs.degV];
if (!(u1 > u0) || !(v1 > v0)) return null;
const ev = (u, v) => bsSurfEval(bs, Math.min(u1, Math.max(u0, u)), Math.min(v1, Math.max(v0, v)));
// a sample grid for the first guess of the inverse and for the size
const G = 24, grid = [];
let size = 0;
for (let i = 0; i <= G; i++) for (let j = 0; j <= G; j++) {
const u = u0 + ((u1 - u0) * i) / G, v = v0 + ((v1 - v0) * j) / G;
grid.push([u, v, ev(u, v)]);
}
for (const g of grid) size = Math.max(size, dist(g[2], grid[0][2]));
const du = (u1 - u0) * 1e-6, dv = (v1 - v0) * 1e-6;
const partials = (u, v) => [mul(sub(ev(u + du, v), ev(u - du, v)), 1 / (2 * du)), mul(sub(ev(u, v + dv), ev(u, v - dv)), 1 / (2 * dv))];
const nrm = (u, v) => { const [a, b] = partials(u, v); return unit(cross(a, b)) || [0, 0, 1]; };
const lenU = Math.max(1e-12, dist(ev(u0, (v0 + v1) / 2), ev(u1, (v0 + v1) / 2))) / (u1 - u0);
const lenV = Math.max(1e-12, dist(ev((u0 + u1) / 2, v0), ev((u0 + u1) / 2, v1))) / (v1 - v0);
return {
type: 'spline', curved: true, periodicU: 0, periodicV: 0, poles: null, scaleU: lenU, scaleV: lenV, domain: [u0, u1, v0, v1], size,
eval: ev,
inv: (p) => {
let best = grid[0], bd = Infinity;
for (const g of grid) { const d = dist(g[2], p); if (d < bd) { bd = d; best = g; } }
let u = best[0], v = best[1];
// Newton on the distance (Gauss–Newton with the partials)
for (let it = 0; it < 30; it++) {
const P = ev(u, v), [a, b] = partials(u, v), r = sub(p, P);
const aa = dot(a, a), ab = dot(a, b), bb = dot(b, b), ra = dot(r, a), rb = dot(r, b);
const det = aa * bb - ab * ab;
if (!(Math.abs(det) > 1e-300)) break;
const su = (ra * bb - rb * ab) / det, sv = (rb * aa - ra * ab) / det;
u = Math.min(u1, Math.max(u0, u + su)); v = Math.min(v1, Math.max(v0, v + sv));
if (Math.abs(su) < (u1 - u0) * 1e-12 && Math.abs(sv) < (v1 - v0) * 1e-12) break;
}
return [u, v];
},
normal: (u, v) => { const n = nrm(u, v); return s.reversed ? mul(n, -1) : n; },
};
}Which knot interval a parameter falls in.
A binary search on the knot vector with plain comparisons; the shift operator is applied to the index only. Parameters at or beyond the ends are placed in the first or last span.
/** knot span index for t (clamped vector, degree p, n control points) */
function span(U, p, n, t) {
if (t >= U[n]) return n - 1;
if (t <= U[p]) return p;
let lo = p, hi = n;
while (hi - lo > 1) { const m = (lo + hi) >> 1; if (t < U[m]) hi = m; else lo = m; }
return lo;
}A point of a B-spline at a parameter, in homogeneous coordinates.
a = (t − Uᵢ) ÷ (Uᵢ₊ₚ₋ᵣ₊₁ − Uᵢ); d_j = (1 − a)·d_j₋₁ + a·d_j
Control points are multiplied by their weight first and the blend is the convex form (1 − a)·x + a·y in 64-bit numbers. A knot interval of zero length gives a = 0 instead of a division by zero.
/** de Boor (Piegl–Tiller A2.1 / A4.1) on homogeneous points [x w, y w, z w, w] */
function deBoor(U, p, P, k, t) {
const d = [];
for (let j = 0; j <= p; j++) { const q = P[k - p + j]; d.push([q[0] * q[3], q[1] * q[3], q[2] * q[3], q[3]]); }
for (let r = 1; r <= p; r++) {
for (let j = p; j >= r; j--) {
const i = k - p + j, den = U[i + p - r + 1] - U[i];
const a = den > 0 ? (t - U[i]) / den : 0;
const x = d[j - 1], y = d[j];
d[j] = [(1 - a) * x[0] + a * y[0], (1 - a) * x[1] + a * y[1], (1 - a) * x[2] + a * y[2], (1 - a) * x[3] + a * y[3]];
}
}
return d[p];
}A point of a rational B-spline curve.
P = (X, Y, Z) ÷ W from de Boor
One de Boor evaluation and three divisions by the weight. A weight of exactly 0 is replaced by 1.
export function bsEval(bs, t) {
const n = bs.nu, p = bs.degU, U = bs.knotsU;
const h = deBoor(U, p, bs.ctrl, span(U, p, n, t), t);
const w = h[3] || 1;
return [h[0] / w, h[1] / w, h[2] / w];
}A point of a rational B-spline surface.
de Boor along v on each of the p + 1 rows, then de Boor along u on the results
Each row result is divided by its weight and handed on with that weight, so the second pass multiplies by it again. A weight of exactly 0 is replaced by 1.
export function bsSurfEval(bs, u, v) {
const { nu, nv, degU: p, degV: q, knotsU: U, knotsV: V } = bs;
const ku = span(U, p, nu, u), kv = span(V, q, nv, v);
// first along v for each of the p + 1 rows that matter, then along u
const rows = [];
for (let i = ku - p; i <= ku; i++) {
const row = bs.ctrl.slice(i * nv, i * nv + nv);
const h = deBoor(V, q, row, kv, v);
rows[i] = [h[0] / (h[3] || 1), h[1] / (h[3] || 1), h[2] / (h[3] || 1), h[3]];
}
const h = deBoor(U, p, rows, ku, u);
const w = h[3] || 1;
return [h[0] / w, h[1] / w, h[2] / w];
}Tests tools/acis/test.mjs, tools/acis/test-apps.mjs
The size of the whole solid and the tolerances every later test uses.
size = max(|hi − lo|, 2·rad, 1e-9); tol = max(10 × resabs, 1e-5 × size); weld = 1e-9 × size
The size is the diagonal of the box round every point, center and root record, or twice the largest radius. Both tolerances scale with it; resabs is the file's own value, taken as 1e-6 when the file has none.
let size = Number.isFinite(lo[0]) ? dist(lo, hi) : 0;
size = Math.max(size, 2 * rad, 1e-9);
const tol = Math.max((dec.resabs || 1e-6) * 10, size * 1e-5);
const qWeld = size * 1e-9;The points along an edge that is part of a circle or an ellipse.
step = min(7.5°, 2·acos(1 − min(1, 0.0002·size ÷ R))); n = min(4096, ⌈|t₁ − t₀| ÷ max(step, 0.001)⌉)
The first and last points are the file's own vertices, not recomputed ones; the points between are evaluated on the ellipse at equal parameter steps. The edge is refused when either vertex is further than tol from the curve. The acos argument is kept within −1 to 1.
if (c.type === 'ellipse') {
const N = unit(c.normal), R = len(c.major);
if (!N || !(R > 0) || !(Math.abs(c.ratio) > 0)) return { why: 'curve', what: 'ellipse-curve' };
const X = unit(sub(c.major, mul(N, dot(c.major, N)))), Y = cross(N, X), Rm = R * c.ratio, C = c.centre;
const at = (t) => add(C, add(mul(X, R * Math.cos(t)), mul(Y, Rm * Math.sin(t))));
const par = (p) => { const d = sub(p, C); return Math.atan2(dot(d, Y) / Rm, dot(d, X) / R); };
const t0 = par(p0);
let t1 = par(p1);
const dir = e.reversed ? -1 : 1;
if (same) t1 = t0 + dir * TAU;
else if (dir > 0) { while (t1 <= t0) t1 += TAU; while (t1 - t0 > TAU) t1 -= TAU; }
else { while (t1 >= t0) t1 -= TAU; while (t0 - t1 > TAU) t1 += TAU; }
if (dist(at(t0), p0) > tol || dist(at(t1), p1) > tol) return { why: 'off_curve', what: 'ellipse-curve' };
const st = Math.min(step, R > 0 ? 2 * Math.acos(Math.max(-1, 1 - Math.min(1, size * 2e-4 / Math.max(R, Math.abs(Rm))))) : step);
const n = Math.min(4096, Math.max(same ? 3 : 1, Math.ceil(Math.abs(t1 - t0) / Math.max(st, 1e-3))));
const pts = [p0];
for (let k = 1; k < n; k++) pts.push(at(t0 + ((t1 - t0) * k) / n));
if (!same) pts.push(p1);
return { pts };
}The points along an edge that follows a B-spline curve.
N = min(4000, max(64, 24 × spans)) samples; end parameters by nearest sample, then 60 steps of a three-way search; n = ⌈|t₁ − t₀| ÷ (b − a) × N ÷ 4⌉ pieces, 2 to 4096
The parameter of each end vertex is found by a search that keeps two thirds of the interval 60 times, which leaves about 3e-11 of the starting interval (two sample spacings). The edge is refused when a vertex is further than tol from the curve. The end points written are the file's vertices.
if (c.type === 'intcurve') {
const bs = c.bs;
if (!bs) return { why: 'curve', what: 'intcurve-curve (' + (c.subtype || '?') + ')' };
const a = bs.knotsU[bs.degU], b = bs.knotsU[bs.knotsU.length - 1 - bs.degU];
if (!(b > a)) return { why: 'curve', what: 'intcurve-curve' };
const N = Math.min(4000, Math.max(64, (bs.nu - bs.degU) * 24));
const ts = [], ps = [];
for (let k = 0; k <= N; k++) { const t = a + ((b - a) * k) / N; ts.push(t); ps.push(bsEval(bs, t)); }
const near = (p) => {
let bi = 0, bd = Infinity;
for (let k = 0; k <= N; k++) { const d = dist(ps[k], p); if (d < bd) { bd = d; bi = k; } }
// refine by halving around the nearest sample
let lo2 = ts[Math.max(0, bi - 1)], hi2 = ts[Math.min(N, bi + 1)];
for (let it = 0; it < 60; it++) {
const m1 = lo2 + (hi2 - lo2) / 3, m2 = hi2 - (hi2 - lo2) / 3;
if (dist(bsEval(bs, m1), p) < dist(bsEval(bs, m2), p)) hi2 = m2; else lo2 = m1;
}
const t = (lo2 + hi2) / 2;
return { t, d: dist(bsEval(bs, t), p) };
};
const A = near(p0), B = same ? A : near(p1);
if (A.d > tol || B.d > tol) return { why: 'off_curve', what: 'intcurve-curve' };
let t0 = A.t, t1 = B.t;
if (same) { t0 = a; t1 = b; }
const n = Math.max(2, Math.min(4096, Math.ceil((Math.abs(t1 - t0) / (b - a)) * N / 4)));
const pts = [p0];
for (let k = 1; k < n; k++) pts.push(bsEval(bs, t0 + ((t1 - t0) * k) / n));
if (!same) pts.push(p1);
// the edge runs from p0 to p1 whichever way the curve's parameter goes: the points above already do
return { pts };
}Gives points that coincide one shared vertex, so a closed solid stays closed.
key = round(x ÷ weld), round(y ÷ weld), round(z ÷ weld)
Points merge when they round to the same cell of a grid of 1e-9 × size. Two points closer than that but on either side of a cell boundary do not merge. The first point's coordinates are kept unrounded.
const vid = (p) => {
const k = Math.round(p[0] / qWeld) + ',' + Math.round(p[1] / qWeld) + ',' + Math.round(p[2] / qWeld);
let id = mesh.weld.get(k);
if (id === undefined) { id = mesh.verts.length; mesh.verts.push(p); mesh.weld.set(k, id); }
return id;
};Checks that every edge point of a face lies on the face's surface before the face is used.
|S(inv(p)) − p| ≤ tol for every loop point
Each point is mapped to (u, v) and back and the distance is compared with tol. A cone is tried with both signs of its angle and a spline with its control points read both ways; if none fits, the face is counted as off the surface and left out.
const tries = s.type === 'cone' ? [surfaceOf(s, 1), surfaceOf(s, -1)]
: s.type === 'spline' ? [surfaceOf(s), surfaceOf({ ...s, bs: bsTransposed(s.bs) })] : [surfaceOf(s)];
let surf = null;
for (const t of tries) {
if (!t) continue;
if (loops.every((pts) => pts.every((p) => dist(t.eval(...t.inv(p)), p) <= tol))) { surf = t; break; }
}
if (!surf) { miss(tries.some(Boolean) ? 'off_surface' : 'surface', sname); return null; }Turns a face's triangles so they point the way the face's outward normal does.
n = (t₁ − t₀) × (t₂ − t₀); vote by the sign of n · normal at the triangle's center
A majority vote over the face: if more triangles disagree with the surface normal than agree, all are flipped. Triangles with a dot product of exactly 0 do not vote.
let agree = 0, against = 0;
for (const t of r.tris) {
const n = cross(sub(t[1], t[0]), sub(t[2], t[0]));
const c = mul(add(add(t[0], t[1]), t[2]), 1 / 3);
const want = mul(surf.normal(...surf.inv(c)), fsign);
const d = dot(n, want);
if (d > 0) agree++; else if (d < 0) against++;
}
if (against > agree) for (const t of r.tris) { const x = t[1]; t[1] = t[2]; t[2] = x; }Triangles of a flat face with its holes.
The loop points are projected to the plane's (u, v) for the triangulation only; the triangles use the original 3D points, so no point is moved.
function planarFace(surf, loops, ctx) {
const rings = loops.filter((l) => l.length >= 3).map((pts) => ({ pts, uv: pts.map((p) => surf.inv(p)) }));
if (!rings.length) return null;
const tris = [], outLoops = [];
for (const g of nest(rings)) {
const all = [g.outer, ...g.holes];
const flat = [], spans = [], P = [];
for (const r of all) {
spans.push([P.length, r.uv.length]);
for (let i = 0; i < r.uv.length; i++) { flat.push(r.uv[i][0], r.uv[i][1]); P.push(r.pts[i]); }
}
const idx = cdt(flat, spans);
for (let i = 0; i < idx.length; i += 3) tris.push([P[idx[i]], P[idx[i + 1]], P[idx[i + 2]]]);
outLoops.push(all.map((r) => r.pts));
}
return { tris, loops: outLoops };
}Sorts the loops of a face into outlines and the holes inside them.
A = |Σ (uᵢ·vᵢ₊₁ − uᵢ₊₁·vᵢ)| ÷ 2; depth = number of larger loops that contain the middle of the loop's first edge
Even depths are outlines and odd depths are holes of the smallest loop one level out. Containment is the even-odd test in plain 64-bit arithmetic with no tolerance.
/** rings ({uv}) → [{ outer, holes }] by containment depth */
function nest(rings) {
const area = (uv) => { let a = 0; for (let i = 0; i < uv.length; i++) { const p = uv[i], q = uv[(i + 1) % uv.length]; a += p[0] * q[1] - q[0] * p[1]; } return a / 2; };
const info = rings.map((r) => ({ r, A: Math.abs(area(r.uv)) }));
for (const a of info) {
a.in = info.filter((b) => b !== a && b.A > a.A && inside(probe(a.r.uv), b.r.uv));
a.depth = a.in.length;
}
const groups = [];
for (const a of info) if (a.depth % 2 === 0) groups.push({ outer: a.r, a, holes: [] });
for (const a of info) if (a.depth % 2 === 1) {
const host = a.in.filter((b) => b.depth === a.depth - 1).sort((x, y) => x.A - y.A)[0];
const g = host && groups.find((q) => q.a === host);
if (g) g.holes.push(a.r);
}
return groups;
}The (u, v) of a loop's points with the angles unwound along the loop, and how many times it goes round the surface.
du wrapped to ±half a turn and added up; windings = round(total ÷ period); pole when cos v · r ≤ tol (sphere) or the cone radius ≤ tol
A jump within 1e-6 rad of exactly half a turn is given the sign that keeps the running total nearest zero. Pole points are split into two (u, v) points that map to the same 3D point.
/**
* A loop's points → [{u, v, p}] with angles unwound along it, poles / the apex split into two points at the angles on
* either side, and how many times it winds round in u and in v.
*/
function loopUV(surf, pts, tol) {
const n = pts.length;
const isPole = (p, uv) => {
if (surf.type === 'sphere') return Math.cos(uv[1]) * surf.scaleU <= tol;
if (surf.type === 'cone' && surf.poles) return surf.radiusAt(uv[1]) <= tol;
return false;
};
const raw = pts.map((p) => { const uv = surf.inv(p); return { u: uv[0], v: uv[1], p, pole: isPole(p, uv) }; });
const Pu = surf.periodicU, Pv = surf.periodicV;
const wrap = (d, P) => { if (!P) return d; while (d > P / 2) d -= P; while (d < -P / 2) d += P; return d; };
const good = raw.filter((x) => !x.pole);
if (!good.length) return null;
// unwind u (and v on a torus) along the non-pole points; a jump of exactly half a turn (through a pole) is given the
// sign that keeps the running total nearest nought
let total = 0, totalV = 0;
for (let i = 0; i < good.length; i++) {
if (i === 0) continue;
let du = wrap(good[i].u - good[i - 1].u, Pu);
if (Pu && Math.abs(Math.abs(du) - Pu / 2) < 1e-6) du = total > 0 ? -Pu / 2 : Pu / 2;
good[i].u = good[i - 1].u + du; total += du;
const dv = wrap(good[i].v - good[i - 1].v, Pv);
good[i].v = good[i - 1].v + dv; totalV += dv;
}
let du = wrap(good[0].u - good[good.length - 1].u, Pu);
if (Pu && Math.abs(Math.abs(du) - Pu / 2) < 1e-6) du = total > 0 ? -Pu / 2 : Pu / 2;
total += du;
totalV += wrap(good[0].v - good[good.length - 1].v, Pv);
const wu = Pu ? Math.round(total / Pu) : 0, wv = Pv ? Math.round(totalV / Pv) : 0;
// poles: (u before, v pole), (u after, v pole); both map to the pole itself
const out = [];
for (let i = 0; i < n; i++) {
const x = raw[i];
if (!x.pole) { out.push({ u: x.u, v: x.v, p: x.p }); continue; }
let a = i - 1, b = i + 1;
while (raw[(a + n) % n].pole && a > i - n) a--;
while (raw[b % n].pole && b < i + n) b++;
const A = raw[(a + n) % n], B = raw[b % n];
let ua = A.u, ub = B.u;
// B comes after A along the loop; at the end of the list it is the first point again, one winding further on
if (b >= n) ub += wu * Pu;
if (a < 0) ua -= wu * Pu;
const vp = surf.type === 'sphere' ? (x.v > 0 ? Math.PI / 2 : -Math.PI / 2) : surf.inv(x.p)[1];
out.push({ u: ua, v: vp, p: x.p }, { u: ub, v: vp, p: x.p });
}
return { pts: out, wu, wv };
}How finely a curved face is divided.
angles: 7.5°; spline: range ÷ 16 each way; along a cone's axis: max(3 × 7.5° in radians × R, 1e-12)
These steps set how far the facets stand off the true surface: a chord of 7.5 degrees is about 0.2 percent of the radius away from its arc at the middle (1 − cos 3.75°).
const su = surf.scaleU || 1, sv = surf.scaleV || 1;
const stepU = surf.type === 'spline' ? (surf.domain[1] - surf.domain[0]) / 16 : step;
const stepV = surf.type === 'spline' ? (surf.domain[3] - surf.domain[2]) / 16 : surf.type === 'cone' ? Math.max(step * su * 3, 1e-12) : step;
const segV = (a, b) => Math.max(1, Math.ceil(Math.abs(b - a) / stepV));
const segU = (a, b) => Math.max(1, Math.ceil(Math.abs(b - a) / stepU));The points added inside a curved face so that its triangles follow the surface.
one point per grid cell at (i + 0.5 + 0.3·jitter)·cell; kept when inside the face and at least 0.35 × the smaller cell side from every boundary edge
The 3D position of each inside point is evaluated on the surface in 64-bit numbers. The boundary keeps only its edge points, so neighboring faces meet on the same vertices. A face with more than 200,000 grid cells gets no inside points.
const gu = g.axisV ? stepV : stepU, gv = g.axisV ? stepU : stepV;
const nu = Math.max(1, Math.round((u1 - u0) / gu)), nv = Math.max(1, Math.round((v1 - v0) / gv));
if (nu * nv > 1 && nu * nv <= MAX_GRID) {
const cu = (u1 - u0) / nu, cv = (v1 - v0) / nv;
const segs = [];
for (const r of rings) for (let i = 0; i < r.length; i++) { const a = r[i], b = r[(i + 1) % r.length]; segs.push([a.u * sx, a.v * sy, b.u * sx, b.v * sy]); }
const minD = 0.35 * Math.min(cu * sx, cv * sy);
const polys = rings.map((r) => r.map((x) => [x.u * sx, x.v * sy]));
// only the segments near a row can be near a point on it
const rowSegs = [];
for (let j = 0; j < nv; j++) {
const y0 = (v0 + cv * (j - 0.5)) * sy - minD, y1 = (v0 + cv * (j + 1.5)) * sy + minD;
rowSegs.push(segs.filter((s) => Math.max(s[1], s[3]) >= Math.min(y0, y1) && Math.min(s[1], s[3]) <= Math.max(y0, y1)));
}
for (let i = 0; i < nu; i++) for (let j = 0; j < nv; j++) {
const u = u0 + cu * (i + 0.5 + 0.3 * jitter(i, j)), v = v0 + cv * (j + 0.5 + 0.3 * jitter(j + 7, i));
const q = [u * sx, v * sy];
let ins = false;
for (const pl of polys) if (inside(q, pl)) ins = !ins;
if (!ins) continue;
let near = false;
for (const s of rowSegs[j]) if (segDist(q, s) < minD) { near = true; break; }
if (near) continue;
flat.push(q[0], q[1]);
P3.push(back(u, v));
}
}
const idx = cdt(flat, spans);A repeatable offset between −0.5 and 0.5 for each grid cell.
s = sin(12.9898·i + 78.233·j) × 43758.5453; s − ⌊s⌋ − 0.5
Not a measured value: it only moves inside points off the grid lines so that no four lie on one circle. The same cell gives the same offset every time in one JavaScript engine; Math.sin is not required to be correctly rounded, so the offset can differ between engines.
/** a fixed pseudo-random number in [-0.5, 0.5) for grid cell (i, j): the same every time, so a face always comes out the same */
function jitter(i, j) { const s = Math.sin(i * 12.9898 + j * 78.233) * 43758.5453; return s - Math.floor(s) - 0.5; }How far an inside point is from a boundary edge, in scaled (u, v).
t = clamp(((q − s₀)·d) ÷ |d|², 0, 1); distance = |q − s₀ − t·d|
The division is skipped for a segment of zero length. Math.hypot gives the distance.
function segDist(q, s) {
const dx = s[2] - s[0], dy = s[3] - s[1];
const l2 = dx * dx + dy * dy;
let t = l2 > 0 ? ((q[0] - s[0]) * dx + (q[1] - s[1]) * dy) / l2 : 0;
t = Math.max(0, Math.min(1, t));
return Math.hypot(q[0] - s[0] - t * dx, q[1] - s[1] - t * dy);
}Moves the mesh by the transform stored with the body.
p′ = s·(p · M) + t; mirrored when det(M)·s³ < 0
Row vector times a 3 × 3 matrix, a uniform scale and a translation, each coordinate in one 64-bit expression. A mirroring transform swaps two indexes of every triangle.
const [a, b, c] = tr.rows, s = tr.s, t = tr.t;
// p' = s · (p · M) + t (row vectors; docs § 4)
V = V.map((p) => [s * (p[0] * a[0] + p[1] * b[0] + p[2] * c[0]) + t[0], s * (p[0] * a[1] + p[1] * b[1] + p[2] * c[1]) + t[1], s * (p[0] * a[2] + p[1] * b[2] + p[2] * c[2]) + t[2]]);
const det = a[0] * (b[1] * c[2] - b[2] * c[1]) - a[1] * (b[0] * c[2] - b[2] * c[0]) + a[2] * (b[0] * c[1] - b[1] * c[0]);
flip = det * s * s * s < 0;Decides whether the mesh is a closed solid: every edge used once in each direction.
key(a, b) = a × 2³² + b
An exact count on whole numbers. The key is exact in a 64-bit number while the first vertex index is below 2²¹ = 2,097,152; above that two different edges could share a key.
const edges = new Map();
const key = (a, b) => a * 4294967296 + b;
for (let i = 0; i < T.length; i += 3) for (let k = 0; k < 3; k++) {
const a = T[i + k], b = T[i + (k + 1) % 3];
edges.set(key(a, b), (edges.get(key(a, b)) || 0) + 1);
}
let closed = T.length > 0;
for (const [k, n] of edges) {
const a = Math.floor(k / 4294967296), b = k % 4294967296;
if (n !== 1 || edges.get(key(b, a)) !== 1) { closed = false; break; }
}The volume enclosed by the triangles and their total area.
V = Σ p · (q × r) ÷ 6; A = Σ |(q − p) × (r − p)| ÷ 2
The volume adds signed tetrahedra from the coordinate origin, triangle by triangle, so for a small solid far from the origin large terms cancel and digits are lost. The area uses edge differences and does not have that problem. Both are those of the facets, not of the curved surfaces, and the volume is reported only when the mesh is closed.
let vol = 0, area = 0;
for (let i = 0; i < T.length; i += 3) {
const p = V[T[i]], q = V[T[i + 1]], r = V[T[i + 2]];
vol += dot(p, cross(q, r)) / 6;
area += len(cross(sub(q, p), sub(r, p))) / 2;
}
if (closed && vol < 0) { for (let i = 0; i < T.length; i += 3) { const x = T[i + 1]; T[i + 1] = T[i + 2]; T[i + 2] = x; } vol = -vol; }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.