Tessellates faces that lie on cylinders, cones, spheres, tori, swept surfaces and B-spline surfaces into triangles.
How the precision is kept
Boundary points keep the 3D positions computed from the file's edges, and interior points are evaluated on the surface in 64-bit. The face is triangulated in the surface's parameter plane, with a spacing chosen so the surface normal turns by about 0.2 radian between points. The triangle size, and for B-splines a Newton iteration that finds parameters, set the accuracy of the result. The triangulation uses ordinary 64-bit tests, not exact predicates.
Calculations under this policy
11 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.
A point on a cylindrical or conical surface from (u, v), and (u, v) from a point.
p = o + (r + v tan α)(x cos u + y sin u) + z v; u = atan2(d · y, d · x), v = d · z
Both directions are closed forms in 64-bit. The semi-angle is multiplied by the file's angle unit. The inverse uses atan2, which keeps its digits at every angle.
ifc/geo-nurbs.js
lines 149 to 154
BimCraft 0.8.1
if (t === 'IFCCYLINDRICALSURFACE' || t === 'IFCCONICALSURFACE') {
const F = frameOf(P()), r = n('Radius'), k = t === 'IFCCONICALSURFACE' ? Math.tan(n('SemiAngle') * g.angleUnit) : 0;
const at = (u, v) => { const rr = r + v * k; return add(F.o, add(add(scale(F.x, rr * Math.cos(u)), scale(F.y, rr * Math.sin(u))), scale(F.z, v))); };
const inv = (p) => { const d = sub(p, F.o); return [Math.atan2(dot(d, F.y), dot(d, F.x)), dot(d, F.z)]; };
return { kind: 'cylinder', at, inv, period: [TAU, 0], curv: [1 / Math.max(r, 1e-9), 0], scaleUV: [r, 1] };
}
Swept curve by length
64-bit
The generating curve of an extruded or revolved surface as a polyline with a length parameter.
s_i = s_{i−1} + |p_i − p_{i−1}|; point at s by linear interpolation
The curve is the sampled polyline of the profile, so the surface between samples is ruled, not curved. The nearest-parameter search projects on every chord with a clamped parameter and a divisor floored at 1e-300.
ifc/geo-nurbs.js
lines 200 to 232
BimCraft 0.8.1
/** A swept curve (an open profile's curve) sampled and parameterised by length: at(s, M), near(p, M) → s. */
function sweptCurve(g, prof) {
let pts = null;
const t = g.t(prof);
if (t === 'IFCARBITRARYOPENPROFILEDEF' || t === 'IFCARBITRARYCLOSEDPROFILEDEF') pts = curvePoints(g, g.ref(g.a(prof, t === 'IFCARBITRARYOPENPROFILEDEF' ? 'Curve' : 'OuterCurve')));
else if (t && g.s.isA(t, 'IFCCURVE')) pts = curvePoints(g, prof);
else if (t && g.s.isA(t, 'IFCPROFILEDEF')) { const sh = profileShapes(g, prof); if (sh.length) pts = sh[0].outer.map((p) => [p[0], p[1], 0]); }
if (!pts || pts.length < 2) return null;
const s = [0];
for (let i = 1; i < pts.length; i++) s.push(s[i - 1] + len(sub(pts[i], pts[i - 1])));
const total = s[s.length - 1] || 1;
const at = (x, M) => {
const v = Math.max(0, Math.min(total, x));
let i = 1;
while (i < s.length - 1 && s[i] < v) i++;
const f = (v - s[i - 1]) / ((s[i] - s[i - 1]) || 1);
const p = add(pts[i - 1], scale(sub(pts[i], pts[i - 1]), f));
return apply(M, p[0], p[1], p[2]);
};
const W = pts.map((p) => p);
const near = (q, M) => {
let best = 0, bd = Infinity;
for (let i = 1; i < W.length; i++) {
const a = apply(M, ...W[i - 1]), b = apply(M, ...W[i]);
const ab = sub(b, a), l2 = dot(ab, ab) || 1e-300;
const f = Math.max(0, Math.min(1, dot(sub(q, a), ab) / l2));
const d = len(sub(q, add(a, scale(ab, f))));
if (d < bd) { bd = d; best = s[i - 1] + f * (s[i] - s[i - 1]); }
}
return best;
};
return { at, near, total, mid: total / 2, curv: pts.length > 2 ? 4 / total : 0 };
}
B-spline basis functions
64-bit
The non-zero basis functions of a B-spline at a parameter.
The standard triangular recurrence in 64-bit; a zero denominator contributes 0. The functions are sums of positive terms, so no cancellation occurs inside a span.
ifc/geo-nurbs.js
lines 248 to 263
BimCraft 0.8.1
const basis = (K, p, k, x) => {
// the p + 1 non-zero B-spline basis functions at x in span k (Cox–de Boor)
const N = [1], left = [], right = [];
for (let j = 1; j <= p; j++) {
left[j] = x - K[k + 1 - j]; right[j] = K[k + j] - x;
let saved = 0;
for (let r = 0; r < j; r++) {
const den = right[r + 1] + left[j - r];
const tmp = den ? N[r] / den : 0;
N[r] = saved + right[r + 1] * tmp;
saved = left[j - r] * tmp;
}
N[j] = saved;
}
return N;
};
Point on a B-spline surface
64-bit
The point of a B-spline surface, rational or not, at (u, v).
p = Σ N_i(u) N_j(v) (w P)_ij ÷ Σ N_i(u) N_j(v) w_ij
Parameters are clamped to the domain, or wrapped in a closed direction. The weighted sum and the weight sum are accumulated in 64-bit and divided once; the division is not guarded against a zero weight sum.
ifc/geo-nurbs.js
lines 266 to 278
BimCraft 0.8.1
const wrapU = (u) => (closedU ? u0 + ((((u - u0) % (u1 - u0)) + (u1 - u0)) % (u1 - u0)) : Math.max(u0, Math.min(u1, u)));
const wrapV = (v) => (closedV ? v0 + ((((v - v0) % (v1 - v0)) + (v1 - v0)) % (v1 - v0)) : Math.max(v0, Math.min(v1, v)));
const at = (u, v) => {
u = wrapU(u); v = wrapV(v);
const ku = span(U, pu, nu, u), kv = span(V, pv, nv, v);
const Nu = basis(U, pu, ku, u), Nv = basis(V, pv, kv, v);
let x = 0, y = 0, z = 0, w = 0;
for (let i = 0; i <= pu; i++) for (let j = 0; j <= pv; j++) {
const h = H[ku - pu + i][kv - pv + j], b = Nu[i] * Nv[j];
x += b * h[0]; y += b * h[1]; z += b * h[2]; w += b * h[3];
}
return [x / w, y / w, z / w];
};
Parameters of a point on a B-spline surface
64-bit
Finds (u, v) of a 3D point on a B-spline surface.
start at the nearest of 25 × 25 samples; Newton: solve (Sᵤ·Sᵤ, Sᵤ·Sᵥ; Sᵤ·Sᵥ, Sᵥ·Sᵥ) Δ = (r·Sᵤ, r·Sᵥ)
Derivatives are one-sided differences with a step of 1e-6 of the parameter range, so they are approximate to about that step and not to full precision. The iteration stops when both steps are below that same size or after 12 steps, and it is not reported when it stops at 12. A determinant under 1e-300 ends it.
ifc/geo-nurbs.js
lines 280 to 308
BimCraft 0.8.1
const G = 24, grid = [];
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, at(u, v)]); }
const du = (u1 - u0) * 1e-6 || 1e-9, dv = (v1 - v0) * 1e-6 || 1e-9;
const inv = (p, hint) => {
let u, v;
if (hint) { [u, v] = hint; }
else { let bd = Infinity; for (const q of grid) { const d = len(sub(q[2], p)); if (d < bd) { bd = d; u = q[0]; v = q[1]; } } }
for (let it = 0; it < 12; it++) {
// one-sided differences that stay inside the domain (a step past its edge would be clamped to a zero derivative)
const s = at(u, v);
const hu = closedU || u + du <= u1 ? du : -du, hv = closedV || v + dv <= v1 ? dv : -dv;
const su = scale(sub(at(u + hu, v), s), 1 / hu), sv = scale(sub(at(u, v + hv), s), 1 / hv);
const r = sub(p, s);
const a = dot(su, su), b = dot(su, sv), c = dot(sv, sv), e = dot(r, su), f = dot(r, sv);
const det = a * c - b * b;
if (Math.abs(det) < 1e-300) break;
const xu = (e * c - f * b) / det, xv = (a * f - b * e) / det;
u = closedU ? u + xu : Math.max(u0, Math.min(u1, u + xu)); v = closedV ? v + xv : Math.max(v0, Math.min(v1, v + xv));
if (Math.abs(xu) < du && Math.abs(xv) < dv) break;
}
if (closedU) u = wrapU(u);
if (closedV) v = wrapV(v);
if (hint) {
// a hint can trap Newton on the far side of a fold: keep the better of it and a fresh grid start
const d1 = len(sub(at(u, v), p));
if (d1 > 1e-6 * (1 + len(p))) { const fresh = inv(p, null); if (len(sub(at(fresh[0], fresh[1]), p)) < d1) return fresh; }
}
return [u, v];
};
Boundary points put back on the surface
64-bit
Moves a boundary point onto the surface when the file's edge curve leaves it.
tolerance = max(10⁻⁹, 0.001 × diagonal of the loops' box)
A point within the tolerance keeps its own position, so neighboring faces still meet exactly; a point farther away is replaced by the surface point at its (u, v). The tolerance scales with the face.
ifc/geo-nurbs.js
lines 345 to 346
BimCraft 0.8.1
const tol = Math.max(1e-9, len(sub(hi3, lo3)) * 1e-3);
for (const L of closedL) L.p = L.p.map((q, i) => { const on = S.at(L.uv[i][0], L.uv[i][1]); return len(sub(on, q)) > tol ? on : q; });
Spacing of interior points
64-bit
The spacing of the points added inside a face, from how fast the surface normal turns.
h = clamp(0.2 ÷ turn, extent ÷ 96, extent ÷ 4), turn = largest change of the normal per unit over a 9 × 9 grid
Normals are cross products of two finite differences with a step of 1e-4 of the face's extent, and the angle between neighbors is an acos of a clamped dot product. The spacing aims at a turn of 0.2 rad between points, never finer than 1/96 of the extent.
ifc/geo-nurbs.js
lines 393 to 403
BimCraft 0.8.1
const K = 8;
let turn = 0;
for (let i = 0; i <= K; i++) for (let j = 0; j <= K; j++) {
const x = lo[0] + (hi[0] - lo[0]) * i / K, y = lo[1] + (hi[1] - lo[1]) * j / K;
const n0 = nrmAt(x, y);
if (!n0) continue;
const ang = (n1) => (n1 ? Math.acos(Math.max(-1, Math.min(1, dot(n0, n1)))) : 0);
if (i < K) turn = Math.max(turn, ang(nrmAt(lo[0] + (hi[0] - lo[0]) * (i + 1) / K, y)) / ((hi[0] - lo[0]) / K || 1));
if (j < K) turn = Math.max(turn, ang(nrmAt(x, lo[1] + (hi[1] - lo[1]) * (j + 1) / K)) / ((hi[1] - lo[1]) / K || 1));
}
const h = Math.max(ext / 96, Math.min(ext / 4, turn > 0 ? 0.2 / turn : ext / 4));
Scale of the parameter plane
64-bit
Scale factors that make distances in (u, v) roughly equal to distances on the surface.
s_u = mean |S(u + e, v) − S(u, v)| ÷ e, e = 10⁻⁴ (1 + |u|)
Analytic surfaces use fixed factors (a radius or 1). Others average finite differences at up to 24 loop points. The factors only shape the triangulation; the 3D positions do not depend on them.
ifc/geo-nurbs.js
lines 643 to 655
BimCraft 0.8.1
/** Scale factors for (u, v) so the parameter plane is roughly isometric to the surface near the loops. */
function metric(S, uv) {
if (S.scaleUV) return S.scaleUV;
let su = 0, sv = 0, n = 0;
const step = Math.max(1, Math.floor(uv.length / 24));
for (let i = 0; i < uv.length; i += step) {
const [u, v] = uv[i];
const e = 1e-4 * (1 + Math.abs(u)), e2 = 1e-4 * (1 + Math.abs(v));
const p = S.at(u, v);
su += len(sub(S.at(u + e, v), p)) / e; sv += len(sub(S.at(u, v + e2), p)) / e2; n++;
}
return [su / n || 1, sv / n || 1];
}
Circumcircle of a triangle
64-bit
The center and squared radius of a triangle's circumcircle, for the Delaunay test.
D = 2 (a_x (b_y − c_y) + b_x (c_y − a_y) + c_x (a_y − b_y)); center = Σ |p|² × (difference) ÷ D
The textbook formula with squared absolute coordinates, so for a thin triangle or points far from the origin the sums nearly cancel. A zero D is replaced by 1e-300, which gives a very large circle instead of a division by zero.
Triangulates the face's points in the scaled parameter plane.
a point is inside a circumcircle when (x − c_x)² + (y − c_y)² < r² (1 + 10⁻¹²)
Bowyer-Watson insertion with the in-circle test in ordinary 64-bit arithmetic and a relative slack of 1e-12, not an exact predicate. The enclosing triangle is 20 times the points' extent. Points are inserted in grid-column order.
ifc/geo-nurbs.js
lines 719 to 773
BimCraft 0.8.1
/**
* Bowyer–Watson Delaunay triangulation of 2D points → [i, j, k, …] counter-clockwise. Triangles are indexed in a grid by
* their circumcircle's box, so an insertion only tests the triangles whose circle can hold the new point.
*/
export function delaunay(P) {
const n = P.length;
if (n < 3) return [];
let lo = [Infinity, Infinity], hi = [-Infinity, -Infinity];
for (const p of P) { if (p[0] < lo[0]) lo[0] = p[0]; if (p[1] < lo[1]) lo[1] = p[1]; if (p[0] > hi[0]) hi[0] = p[0]; if (p[1] > hi[1]) hi[1] = p[1]; }
const d = Math.max(hi[0] - lo[0], hi[1] - lo[1]) || 1, cx = (lo[0] + hi[0]) / 2, cy = (lo[1] + hi[1]) / 2;
const pts = P.concat([[cx - 20 * d, cy - d], [cx + 20 * d, cy - d], [cx, cy + 20 * d]]);
const G = Math.max(1, Math.min(256, Math.ceil(Math.sqrt(n / 2))));
const cw = (hi[0] - lo[0]) / G || 1, ch = (hi[1] - lo[1]) / G || 1;
const cells = Array.from({ length: G * G }, () => []);
const big = [];
const cellOf = (x, y) => [Math.min(G - 1, Math.max(0, Math.floor((x - lo[0]) / cw))), Math.min(G - 1, Math.max(0, Math.floor((y - lo[1]) / ch)))];
const add = (t) => {
const r = Math.sqrt(t.r2);
const [x0, y0] = cellOf(t.cx - r, t.cy - r), [x1, y1] = cellOf(t.cx + r, t.cy + r);
if ((x1 - x0 + 1) * (y1 - y0 + 1) > 64) { big.push(t); return; }
for (let x = x0; x <= x1; x++) for (let y = y0; y <= y1; y++) cells[y * G + x].push(t);
};
add(mk(pts, n, n + 1, n + 2));
const order = [...Array(n).keys()].sort((a, b) => (Math.floor((P[a][0] - lo[0]) / cw) - Math.floor((P[b][0] - lo[0]) / cw)) || P[a][1] - P[b][1]);
for (const i of order) {
const [x, y] = pts[i];
const [gx, gy] = cellOf(x, y);
const cand = cells[gy * G + gx];
const bad = [];
const test = (list) => {
for (let k = 0; k < list.length; k++) {
const t = list[k];
if (t.dead) { list[k] = list[list.length - 1]; list.pop(); k--; continue; }
if (!t.mark && (x - t.cx) ** 2 + (y - t.cy) ** 2 < t.r2 * (1 + 1e-12)) { t.mark = true; bad.push(t); }
}
};
test(cand); test(big);
const edges = new Map();
for (const t of bad) {
t.dead = true;
for (const [a, b] of [[t.a, t.b], [t.b, t.c], [t.c, t.a]]) {
const k = a < b ? a * 1048576 + b : b * 1048576 + a;
if (edges.has(k)) edges.delete(k); else edges.set(k, [a, b]);
}
}
for (const [a, b] of edges.values()) add(mk(pts, a, b, i));
}
const out = [], seen = new Set();
for (const list of [...cells, big]) for (const t of list) {
if (t.dead || seen.has(t)) continue;
seen.add(t);
if (t.a < n && t.b < n && t.c < n) out.push(t.a, t.b, t.c);
}
return out;
}
Orientation and crossing tests
64-bit
Whether three points turn left or right, and whether two segments cross, for recovering the boundary edges.
A plain 64-bit cross product of coordinate differences with no error bound; for nearly collinear points its sign can be wrong. Segments are taken to cross only when all four orientation values are strictly non-zero and of opposite signs.
Closed-form mapping to (u, v) for analytic surfaces; Cox-de Boor basis functions for B-spline surfaces; nearest of a 25 × 25 grid then Newton steps for the inverse; Bowyer-Watson Delaunay triangulation with edge recovery by diagonal flips (Sloan 1993, as the code's comment cites).
Number types
64-bit numbers; Float64Array for accumulated normals.
Tolerances
Newton stops when the step is under 1e-6 of the parameter range, after at most 12 steps. A boundary point farther than 0.001 × the loops' box diagonal from the surface is moved onto it. The in-circle test allows 1e-12 of the squared radius. Repeated parameter points are merged at 1e-9.
Limits
Interior spacing is between 1/96 and 1/4 of the face's extent. Orientation and crossing tests are plain 64-bit cross products, so nearly degenerate point sets can be decided either way. Edge keys allow at most 1 048 576 points in a face.
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.