Resources Numerical accuracy Calculation policy
The vector, plane, ray and 4×4 matrix helpers that the 3D modeling kernel uses for every point it stores.
Points are plain arrays of three 64-bit numbers in meters and matrices are 16 64-bit numbers. Every helper is a closed-form expression. Lengths use Math.hypot. The tolerances are fixed numbers: 1e-6 m for lengths, and several absolute thresholds for parallel or degenerate cases that do not scale with the size of the input.
Open item. Scale the parallel test in rayToSegment by the two squared lengths, as lineLine already does.
12 calculations in 1 source file. The code is read from the released app when this page loads. It is the calculation people run, not a copy kept for this page.
Tests bin/test-modelcraft.mjs, bin/test-cadcraft-3d.mjs
A vector scaled to length 1.
n = a ÷ |a|; [0, 0, 0] when |a| ≤ 1e-15
The length comes from Math.hypot and each component is divided by it. The cut-off of 1e-15 is an absolute length, not scaled to the input.
export function norm(a) {
const l = len(a);
return l > 1e-15 ? [a[0] / l, a[1] / l, a[2] / l] : [0, 0, 0];
}The unit normal of a polygon given by its corner points.
N = Σ (pᵢ − pᵢ₊₁) × (pᵢ + pᵢ₊₁) taken per axis, then normalized
Three sums of products, added in corner order. Each term multiplies a coordinate difference by a coordinate sum, so for a small polygon far from the origin the terms are large and mostly cancel.
/** Newell normal of a polygon (array of points): unit, CCW = right-hand rule; [0,0,0] if degenerate. */
export function newell(pts) {
let x = 0, y = 0, z = 0;
for (let i = 0; i < pts.length; i++) {
const p = pts[i], q = pts[(i + 1) % pts.length];
x += (p[1] - q[1]) * (p[2] + q[2]);
y += (p[2] - q[2]) * (p[0] + q[0]);
z += (p[0] - q[0]) * (p[1] + q[1]);
}
return norm([x, y, z]);
}The area of a flat polygon in space.
area = |Σ (pᵢ − pᵢ₊₁) × (pᵢ + pᵢ₊₁)| ÷ 2
The same three sums as the polygon normal, then Math.hypot and a division by 2. The sums use absolute coordinates, so the same cancellation applies far from the origin.
/** Twice-area vector length (Newell), unnormalised: the polygon's area = len / 2. */
export function polyArea3(pts) {
let x = 0, y = 0, z = 0;
for (let i = 0; i < pts.length; i++) {
const p = pts[i], q = pts[(i + 1) % pts.length];
x += (p[1] - q[1]) * (p[2] + q[2]);
y += (p[2] - q[2]) * (p[0] + q[0]);
z += (p[0] - q[0]) * (p[1] + q[1]);
}
return Math.hypot(x, y, z) / 2;
}A plane as a unit normal and an offset, with the normal's sign made canonical.
n·p = d; first component of n with |n| > 1e-9 is made positive
The normal is normalized first and the offset is one dot product with the point. The sign test uses an absolute 1e-9 on a component of the unit normal.
/** A plane {n, d} with n·p = d and n made canonical (first non-zero component positive). */
export function planeFrom(n, p) {
let m = norm(n);
const k = Math.abs(m[0]) > 1e-9 ? 0 : Math.abs(m[1]) > 1e-9 ? 1 : 2;
if (m[k] < 0) m = neg(m);
return { n: m, d: dot(m, p) };
}Whether two planes are the same plane.
|n₁·n₂| > 1 − 1e-9 and |d₁ × sign(n₁·n₂) − d₂| ≤ 1e-6
The cosine test passes normals up to about 4.5e-5 rad apart (1 − cos θ ≈ θ² ÷ 2). The offset test is an absolute 1e-6 m, whatever the distance of the plane from the origin.
export const samePlane = (a, b) => Math.abs(dot(a.n, b.n)) > 1 - 1e-9 && Math.abs(a.d * Math.sign(dot(a.n, b.n)) - b.d) <= EPS;The nearest points between a ray's line and an edge.
den = (d·d)(u·u) − (d·u)²; t = (B·E − C·D) ÷ den; s = (A·E − B·D) ÷ den, clamped to [0, 1]
The determinant subtracts two products that are nearly equal when the lines are close to parallel. It is compared with an absolute 1e-14, which is not scaled by the lengths, so a very short edge is treated as parallel at a larger angle than a long one.
/** Closest points between the line o + t·d (a ray, t ≥ 0 not enforced) and the segment a–b. */
export function rayToSegment(o, d, a, b) {
const u = sub(b, a), w = sub(o, a);
const A = dot(d, d), B = dot(d, u), C = dot(u, u), D = dot(d, w), E = dot(u, w);
const den = A * C - B * B;
let t, s;
if (den < 1e-14) { s = 0; t = -D / A; }
else { t = (B * E - C * D) / den; s = (A * E - B * D) / den; }
if (s < 0) { s = 0; t = -D / A; }
else if (s > 1) { s = 1; t = (B - D) / A; }
return { t, s, onRay: add(o, mul(d, t)), onSeg: add(a, mul(u, s)) };
}The nearest points between two infinite lines, or nothing when they are parallel.
den = A·C − B²; parallel when |den| < 1e-14 × A × C
The same 2×2 solve as the segment case, but the parallel test is scaled by the two squared lengths. The determinant still loses digits as the lines approach parallel.
/** Closest points between two infinite lines; null when parallel. */
export function lineLine(o1, d1, o2, d2) {
const w = sub(o1, o2);
const A = dot(d1, d1), B = dot(d1, d2), C = dot(d2, d2), D = dot(d1, w), E = dot(d2, w);
const den = A * C - B * B;
if (Math.abs(den) < 1e-14 * A * C) return null;
const t = (B * E - C * D) / den, s = (A * E - B * D) / den;
return { t, s, p1: add(o1, mul(d1, t)), p2: add(o2, mul(d2, s)) };
}The point where a ray's line meets a plane.
t = (d_plane − n·o) ÷ (n·dir); p = o + t × dir
One division. A dot product of normal and direction below 1e-12 (absolute) counts as parallel and returns nothing. The numerator subtracts the plane offset from n·o, which loses digits when the origin is far from the plane's foot point relative to the gap.
export function rayPlane(o, d, pl) {
const den = dot(pl.n, d);
if (Math.abs(den) < 1e-12) return null;
const t = (pl.d - dot(pl.n, o)) / den;
return { t, p: add(o, mul(d, t)) };
}The distance along a ray to a triangle, or nothing when it misses.
Möller–Trumbore: det = e₁·(d × e₂); u, v, t from three scalar triple products ÷ det
A determinant below 1e-14 (absolute, so it depends on the triangle's size) counts as parallel. The barycentric coordinates are accepted up to 1e-9 outside the triangle.
/** Möller–Trumbore; returns t or null (both sides). */
export function rayTri(o, d, a, b, c) {
const e1 = sub(b, a), e2 = sub(c, a), p = cross(d, e2), det = dot(e1, p);
if (Math.abs(det) < 1e-14) return null;
const inv = 1 / det, s = sub(o, a), u = dot(s, p) * inv;
if (u < -1e-9 || u > 1 + 1e-9) return null;
const q = cross(s, e1), v = dot(d, q) * inv;
if (v < -1e-9 || u + v > 1 + 1e-9) return null;
const t = dot(e2, q) * inv;
return t > 0 ? t : null;
}The 4×4 matrix that turns geometry by an angle about an axis through a given point.
M = T(about) × R(axis, θ) × T(−about), R from Rodrigues' formula
The axis is normalized and the matrix is built from one cosine and one sine. The translation part comes out of two matrix products as about − R × about, which subtracts nearly equal numbers when the point is far from the origin.
/** Rotation by angle (radians, right-hand) about the axis through `about`. */
export function rotation(axis, angle, about = [0, 0, 0]) {
const [x, y, z] = norm(axis), c = Math.cos(angle), s = Math.sin(angle), t = 1 - c;
const R = [
t * x * x + c, t * x * y + s * z, t * x * z - s * y, 0,
t * x * y - s * z, t * y * y + c, t * y * z + s * x, 0,
t * x * z + s * y, t * y * z - s * x, t * z * z + c, 0,
0, 0, 0, 1,
];
return mmul(translation(about), mmul(R, translation(neg(about))));
}A point moved by a 4×4 matrix.
p′ = (M × [x, y, z, 1]) ÷ w
Three sums of four terms and a division by w. When w works out to exactly 0 it is replaced by 1.
export function apply(m, p) {
const x = p[0], y = p[1], z = p[2];
const w = m[3] * x + m[7] * y + m[11] * z + m[15] || 1;
return [(m[0] * x + m[4] * y + m[8] * z + m[12]) / w, (m[1] * x + m[5] * y + m[9] * z + m[13]) / w, (m[2] * x + m[6] * y + m[10] * z + m[14]) / w];
}The matrix that undoes a transform, used to move points into a group's own coordinates.
M⁻¹ = adj(M) ÷ det(M), with det from twelve 2×2 minors
A cofactor expansion without pivoting. Only a determinant of exactly zero is caught, and then the identity is returned without a flag; a nearly singular matrix is inverted as it stands.
export function invert(m) {
const a = m, r = new Array(16);
const [a00, a01, a02, a03, a10, a11, a12, a13, a20, a21, a22, a23, a30, a31, a32, a33] = a;
const b00 = a00 * a11 - a01 * a10, b01 = a00 * a12 - a02 * a10, b02 = a00 * a13 - a03 * a10;
const b03 = a01 * a12 - a02 * a11, b04 = a01 * a13 - a03 * a11, b05 = a02 * a13 - a03 * a12;
const b06 = a20 * a31 - a21 * a30, b07 = a20 * a32 - a22 * a30, b08 = a20 * a33 - a23 * a30;
const b09 = a21 * a32 - a22 * a31, b10 = a21 * a33 - a23 * a31, b11 = a22 * a33 - a23 * a32;
let det = b00 * b11 - b01 * b10 + b02 * b09 + b03 * b08 - b04 * b07 + b05 * b06;
if (!det) return ident();
det = 1 / det;
r[0] = (a11 * b11 - a12 * b10 + a13 * b09) * det;
r[1] = (a02 * b10 - a01 * b11 - a03 * b09) * det;
r[2] = (a31 * b05 - a32 * b04 + a33 * b03) * det;
r[3] = (a22 * b04 - a21 * b05 - a23 * b03) * det;
r[4] = (a12 * b08 - a10 * b11 - a13 * b07) * det;
r[5] = (a00 * b11 - a02 * b08 + a03 * b07) * det;
r[6] = (a32 * b02 - a30 * b05 - a33 * b01) * det;
r[7] = (a20 * b05 - a22 * b02 + a23 * b01) * det;
r[8] = (a10 * b10 - a11 * b08 + a13 * b06) * det;
r[9] = (a01 * b08 - a00 * b10 - a03 * b06) * det;
r[10] = (a30 * b04 - a31 * b02 + a33 * b00) * det;
r[11] = (a21 * b02 - a20 * b04 - a23 * b00) * det;
r[12] = (a11 * b07 - a10 * b09 - a12 * b06) * det;
r[13] = (a00 * b09 - a01 * b07 + a02 * b06) * det;
r[14] = (a31 * b01 - a30 * b03 - a32 * b00) * det;
r[15] = (a20 * b03 - a21 * b01 + a22 * b00) * det;
return r;
}Every value is a 64-bit number from input to result. Formulas are closed form or a direct solve, and every comparison uses a stated tolerance.