Finds the rigid move (rotation and shift) that lines one scan up with others, from picked point pairs or by iterative closest point, and reports the remaining error.
How the precision is kept
Transforms are 4 × 4 matrices of 64-bit numbers and the work is done in coordinates shifted to the first scan's origin, so large map coordinates keep their small digits. The pair fit is Horn's closed form with a Jacobi eigen solver. Two 32-bit steps reach results: scan points are read from a Float32Array of offsets, and the surface normals used by the point-to-plane step and its reported plane error are stored in a Float32Array.
Open item. Keep the surface normals in a Float64Array so the point-to-plane step and the reported plane error are computed from 64-bit normals.
Calculations under this policy
15 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.
register.jsCadCraft 0.8.2 32-bit in part15 calculationsReviewed
Tests bin/test-register.mjs
Rigid fit from point pairs
64-bit
The rotation and shift that best take the picked points of one scan onto their targets, with the error left at each pair.
R from the top eigenvector (a quaternion) of Horn's 4 × 4 matrix; t = c_dst − R·c_src; rms = √(Σ rᵢ² ÷ n)
Centroids are summed as xᵢ ÷ n term by term, and the covariances are sums of products of coordinates with the centroid taken off. Points on a line are refused when the second eigenvalue of the source covariance is below 1e-9 times the first. Residuals and rms are measured with the fitted matrix applied to the source points.
register.js
lines 108 to 138
CadCraft 0.8.2
/** Horn's closed form: the rotation is the top eigenvector (a quaternion) of a 4 × 4 built from the pairs' cross-covariance. */
export function fitRigid(src, dst) {
const n = Math.min(src.length, dst.length);
if (n < 3) throw fail('few', { n });
const cs = [0, 0, 0], cd = [0, 0, 0];
for (let i = 0; i < n; i++) for (let a = 0; a < 3; a++) { cs[a] += src[i][a] / n; cd[a] += dst[i][a] / n; }
const S = [[0, 0, 0], [0, 0, 0], [0, 0, 0]], C = [[0, 0, 0], [0, 0, 0], [0, 0, 0]];
for (let i = 0; i < n; i++) {
const s = [src[i][0] - cs[0], src[i][1] - cs[1], src[i][2] - cs[2]], d = [dst[i][0] - cd[0], dst[i][1] - cd[1], dst[i][2] - cd[2]];
for (let a = 0; a < 3; a++) for (let b = 0; b < 3; b++) { S[a][b] += s[a] * d[b]; C[a][b] += s[a] * s[b]; }
}
const ev = eigenSym(C).values;
if (!(ev[0] > 0) || ev[1] < 1e-9 * ev[0]) throw fail('collinear'); // a line fixes no turn about itself
const [Sxx, Sxy, Sxz] = S[0], [Syx, Syy, Syz] = S[1], [Szx, Szy, Szz] = S[2];
const N = [
[Sxx + Syy + Szz, Syz - Szy, Szx - Sxz, Sxy - Syx],
[Syz - Szy, Sxx - Syy - Szz, Sxy + Syx, Szx + Sxz],
[Szx - Sxz, Sxy + Syx, -Sxx + Syy - Szz, Syz + Szy],
[Sxy - Syx, Szx + Sxz, Syz + Szy, -Sxx - Syy + Szz]];
const [w, x, y, z] = eigenSym(N).vectors[0];
const R = [
1 - 2 * (y * y + z * z), 2 * (x * y - w * z), 2 * (x * z + w * y),
2 * (x * y + w * z), 1 - 2 * (x * x + z * z), 2 * (y * z - w * x),
2 * (x * z - w * y), 2 * (y * z + w * x), 1 - 2 * (x * x + y * y)];
const t = [0, 1, 2].map((r) => cd[r] - (R[r * 3] * cs[0] + R[r * 3 + 1] * cs[1] + R[r * 3 + 2] * cs[2]));
const M = [R[0], R[1], R[2], t[0], R[3], R[4], R[5], t[1], R[6], R[7], R[8], t[2], 0, 0, 0, 1];
const residuals = [];
let ss = 0;
for (let i = 0; i < n; i++) { const p = apply4(M, src[i]); const r = Math.hypot(p[0] - dst[i][0], p[1] - dst[i][1], p[2] - dst[i][2]); residuals.push(r); ss += r * r; }
return { M, rms: Math.sqrt(ss / n), residuals };
}
Eigenvalues of a symmetric matrix
64-bit
Eigenvalues and eigenvectors of the 3 × 3 and 4 × 4 symmetric matrices the fit and the normals use.
Cyclic Jacobi sweeps, at most 60; the loop ends when the sum of squared off-diagonal entries is below 1e-300, an absolute value far below rounding level for matrices of ordinary size. The tangent is computed in the form that adds, not subtracts, the square root.
register.js
lines 87 to 105
CadCraft 0.8.2
/** Eigenvalues (descending) and eigenvectors (columns of V as arrays, same order) of a symmetric matrix. */
export function eigenSym(A0) {
const n = A0.length, A = A0.map((r) => r.slice()), V = Array.from({ length: n }, (_, i) => Array.from({ length: n }, (_, j) => (i === j ? 1 : 0)));
for (let sweep = 0; sweep < 60; sweep++) {
let off = 0;
for (let i = 0; i < n; i++) for (let j = i + 1; j < n; j++) off += A[i][j] * A[i][j];
if (off < 1e-300) break;
for (let p = 0; p < n; p++) for (let q = p + 1; q < n; q++) {
if (Math.abs(A[p][q]) < 1e-300) continue;
const th = (A[q][q] - A[p][p]) / (2 * A[p][q]);
const t = Math.sign(th || 1) / (Math.abs(th) + Math.sqrt(th * th + 1)), c = 1 / Math.sqrt(t * t + 1), s = t * c;
for (let k = 0; k < n; k++) { const akp = A[k][p], akq = A[k][q]; A[k][p] = c * akp - s * akq; A[k][q] = s * akp + c * akq; }
for (let k = 0; k < n; k++) { const apk = A[p][k], aqk = A[q][k]; A[p][k] = c * apk - s * aqk; A[q][k] = s * apk + c * aqk; }
for (let k = 0; k < n; k++) { const vkp = V[k][p], vkq = V[k][q]; V[k][p] = c * vkp - s * vkq; V[k][q] = s * vkp + c * vkq; }
}
}
const order = [...Array(n).keys()].sort((a, b) => A[b][b] - A[a][a]);
return { values: order.map((i) => A[i][i]), vectors: order.map((i) => V.map((r) => r[i])) };
}
Combining transforms
64-bit
The product of two 4 × 4 transforms, used to add each iteration's step to the running transform.
o = a · b, each entry a sum of four products
Each entry is four products added in order. The running transform is a chain of such products, one per iteration, and is not re-orthonormalized; after a run the result is checked by isRigid with a tolerance of 1e-6.
register.js
lines 39 to 43
CadCraft 0.8.2
export function mul4(a, b) {
const o = new Array(16).fill(0);
for (let r = 0; r < 4; r++) for (let c = 0; c < 4; c++) { let s = 0; for (let k = 0; k < 4; k++) s += a[r * 4 + k] * b[k * 4 + c]; o[r * 4 + c] = s; }
return o;
}
Inverse of a rigid transform
64-bit
The transform that undoes a rotation and shift.
R⁻¹ = Rᵀ; t′ = −Rᵀ·t
Uses the transpose, which is the inverse only when the rotation part is orthonormal; no determinant is divided by.
Rewrites a transform for coordinates measured from the working origin, and back.
t_local = t + R·o − o; back: t = t_local − R·o + o
Only the shift changes. With a large origin, R·o and o are large and nearly equal, so the absolute transform's shift carries the rounding of numbers the size of the origin.
register.js
lines 54 to 60
CadCraft 0.8.2
/** M expressed in coordinates shifted by `o` (local = absolute − o): out · (p − o) = M · p − o. `back` undoes it. */
export function shiftM(M, o, back = false) {
const R = [[M[0], M[1], M[2]], [M[4], M[5], M[6]], [M[8], M[9], M[10]]];
const Ro = R.map((r) => r[0] * o[0] + r[1] * o[1] + r[2] * o[2]);
const t = [M[3], M[7], M[11]].map((v, k) => (back ? v - Ro[k] + o[k] : v + Ro[k] - o[k]));
return [M[0], M[1], M[2], t[0], M[4], M[5], M[6], t[1], M[8], M[9], M[10], t[2], 0, 0, 0, 1];
}
Is a matrix a rigid transform
64-bit
Checks that a matrix is a pure rotation and shift before it is accepted.
|rᵢ·rⱼ − δᵢⱼ| ≤ tol; |det − 1| ≤ tol; last row 0 0 0 1 within tol
The tolerance is absolute on dimensionless entries: 1e-6 for transforms the app makes and 1e-5 for a project file read from outside.
register.js
lines 61 to 72
CadCraft 0.8.2
/** Is M a rigid transform (orthonormal rotation, determinant 1, last row 0 0 0 1)? */
export function isRigid(M, tol = 1e-6) {
if (!Array.isArray(M) || M.length !== 16 || !M.every(Number.isFinite)) return false;
if (Math.abs(M[12]) > tol || Math.abs(M[13]) > tol || Math.abs(M[14]) > tol || Math.abs(M[15] - 1) > tol) return false;
const r = [[M[0], M[1], M[2]], [M[4], M[5], M[6]], [M[8], M[9], M[10]]];
for (let a = 0; a < 3; a++) for (let b = a; b < 3; b++) {
const d = r[a][0] * r[b][0] + r[a][1] * r[b][1] + r[a][2] * r[b][2];
if (Math.abs(d - (a === b ? 1 : 0)) > tol) return false;
}
const det = r[0][0] * (r[1][1] * r[2][2] - r[1][2] * r[2][1]) - r[0][1] * (r[1][0] * r[2][2] - r[1][2] * r[2][0]) + r[0][2] * (r[1][0] * r[2][1] - r[1][1] * r[2][0]);
return Math.abs(det - 1) <= tol;
}
Rotation about an axis
64-bit
The 4 × 4 matrix for a turn by an angle about an axis through a center.
R = c·I + s·[axis]ₓ + (1 − c)·axis·axisᵀ; t = center − R·center
Rodrigues' formula with one cosine and one sine. The axis is used as given and is not normalized here; for a very small angle 1 − cos loses digits, and that term is then itself very small.
register.js
lines 73 to 79
CadCraft 0.8.2
/** A rotation of `angle` radians about the (unit) axis, as a 4 × 4. */
export function rotation4(axis, angle, centre = [0, 0, 0]) {
const [x, y, z] = axis, c = Math.cos(angle), s = Math.sin(angle), k = 1 - c;
const R = [c + x * x * k, x * y * k - z * s, x * z * k + y * s, y * x * k + z * s, c + y * y * k, y * z * k - x * s, z * x * k - y * s, z * y * k + x * s, c + z * z * k];
const t = [0, 1, 2].map((r) => centre[r] - (R[r * 3] * centre[0] + R[r * 3 + 1] * centre[1] + R[r * 3 + 2] * centre[2]));
return [R[0], R[1], R[2], t[0], R[3], R[4], R[5], t[1], R[6], R[7], R[8], t[2], 0, 0, 0, 1];
}
Angle and shift of a transform in words
64-bit
The rotation angle in degrees and the length of the shift that a result is reported with.
angle = acos((trace − 1) ÷ 2); shift = |t|
The cosine is held to [−1, 1] before acos. Near zero rotation the cosine is within rounding of 1, so a turn below about 1.5e-8 rad reads as zero or as a coarse step.
register.js
lines 80 to 84
CadCraft 0.8.2
/** Rotation angle (degrees) and translation length of a rigid transform: for reading a result in words. */
export function describeM(M) {
const tr = M[0] + M[5] + M[10];
return { angle: Math.acos(Math.min(1, Math.max(-1, (tr - 1) / 2))) * 180 / Math.PI, shift: Math.hypot(M[3], M[7], M[11]) };
}
Working points of a scan
32-bit in part
The thinned copy of a scan's points, measured from the session origin, that registration works on.
work = pos + cloud origin − session origin, every stride-th point
cloud.pos is a Float32Array of offsets from the cloud's own origin (pointcloud.js), so each point enters with about 7 digits relative to that origin; the sums are then 64-bit in a Float64Array. The session origin is the first cloud's origin rounded to whole units.
register.js
lines 377 to 386
CadCraft 0.8.2
addScan({ name, size = 0, head = '', cloud }) {
if (this.scans.length >= LIMITS.scans) throw fail('too_many_scans');
if (!this.origin) this.origin = cloud.origin.map((v) => Math.round(v));
const n = cloud.n, stride = Math.max(1, Math.ceil(n / this.workPoints)), m = Math.floor((n + stride - 1) / stride);
const work = new Float64Array(m * 3), o = this.origin;
for (let k = 0; k < m; k++) {
const i = k * stride;
work[k * 3] = cloud.pos[i * 3] + cloud.origin[0] - o[0]; work[k * 3 + 1] = cloud.pos[i * 3 + 1] + cloud.origin[1] - o[1]; work[k * 3 + 2] = cloud.pos[i * 3 + 2] + cloud.origin[2] - o[2];
}
this.scans.push({ name, size, head, cloud, M: I4(), hist: [], pairs: [], work });
Nearest point in the fixed scans
64-bit
For a moved point, the closest point of the fixed scans within the match distance.
d² = Δx² + Δy² + Δz², smallest over the grid cells searched ring by ring
Squared distances are compared in 64-bit and the square root is taken once at the end. The search widens one ring of cells at a time and stops when the best distance found is within the rings already covered.
register.js
lines 189 to 200
CadCraft 0.8.2
nearest(x, y, z, maxR, out = {}) {
const cx = this.cellOf(x, 0), cy = this.cellOf(y, 1), cz = this.cellOf(z, 2);
const R = Math.min(Math.ceil(maxR / this.cell), Math.max(...this.dims));
let best = -1, bd = Infinity;
const p = this.pts;
for (let r = 0; r <= R; r++) {
this.ring(cx, cy, cz, r, (i) => { const dx = p[i * 3] - x, dy = p[i * 3 + 1] - y, dz = p[i * 3 + 2] - z, d = dx * dx + dy * dy + dz * dz; if (d < bd) { bd = d; best = i; } });
if (best >= 0 && Math.sqrt(bd) <= r * this.cell) break;
}
out.d = Math.sqrt(bd);
return best >= 0 && out.d <= maxR ? best : -1;
}
Surface normals
32-bit in part
The direction square to the surface at each fixed point, from its 10 nearest neighbors.
normal = eigenvector of the smallest eigenvalue of the neighbors' covariance
The covariance and its eigenvectors are 64-bit, but each normal is written into a Float32Array, about 7 digits. A point with fewer than 5 neighbors, or whose neighbors lie on a line (second eigenvalue below 1e-6 times the first), gets no normal.
register.js
lines 220 to 236
CadCraft 0.8.2
/** Surface normals of a point set by the smallest axis of its k nearest neighbours; null where the points are too few or lie on a line. */
export function estimateNormals(grid, k = 10, into = null, from = 0, to = grid.n) {
const n = grid.n, out = into?.normals || new Float32Array(n * 3), ok = into?.ok || new Uint8Array(n), p = grid.pts;
for (let i = from; i < to; i++) {
const nb = grid.knn(p[i * 3], p[i * 3 + 1], p[i * 3 + 2], k);
if (nb.length < 5) continue;
const c = [0, 0, 0];
for (const j of nb) for (let a = 0; a < 3; a++) c[a] += p[j * 3 + a] / nb.length;
const C = [[0, 0, 0], [0, 0, 0], [0, 0, 0]];
for (const j of nb) { const d = [p[j * 3] - c[0], p[j * 3 + 1] - c[1], p[j * 3 + 2] - c[2]]; for (let a = 0; a < 3; a++) for (let b = 0; b < 3; b++) C[a][b] += d[a] * d[b]; }
const e = eigenSym(C);
if (!(e.values[0] > 0) || e.values[1] < 1e-6 * e.values[0]) continue; // a line: no normal
const v = e.vectors[2];
out[i * 3] = v[0]; out[i * 3 + 1] = v[1]; out[i * 3 + 2] = v[2]; ok[i] = 1;
}
return { normals: out, ok };
}
Matching and the error of a step
32-bit in part
Pairs each moving point with its nearest fixed point and reports the root-mean-square distance, and the distance along the normals.
rms is 64-bit throughout. The plane error multiplies 64-bit offsets by normals read back from the Float32Array, so the reported plane error carries about 7 digits. Fewer than 6 matches stops the run.
register.js
lines 305 to 316
CadCraft 0.8.2
for (let i = 0; i < n_src; i++) {
const x = src[i * 3], y = src[i * 3 + 1], z = src[i * 3 + 2];
const px = T[0] * x + T[1] * y + T[2] * z + T[3], py = T[4] * x + T[5] * y + T[6] * z + T[7], pz = T[8] * x + T[9] * y + T[10] * z + T[11];
const j = grid.nearest(px, py, pz, cut, hit);
if (j < 0) continue;
if (phase === 'plane' && !nrm.ok[j]) continue;
matched++; ss += hit.d * hit.d;
P.push(px, py, pz); Q.push(dst[j * 3], dst[j * 3 + 1], dst[j * 3 + 2]);
if (nrm) { const nx = nrm.normals[j * 3], ny = nrm.normals[j * 3 + 1], nz = nrm.normals[j * 3 + 2]; N.push(nx, ny, nz); const d = (px - dst[j * 3]) * nx + (py - dst[j * 3 + 1]) * ny + (pz - dst[j * 3 + 2]) * nz; sp += d * d; }
}
if (matched < 6) throw fail('no_overlap', { matched });
const rms = Math.sqrt(ss / matched), planeRms = nrm ? Math.sqrt(sp / matched) : null;
Point-to-plane step
32-bit in part
The small turn and shift that reduce the distances along the normals in one iteration.
minimize Σ ((ω × p + t + p − q)·n)²: 6 × 6 normal equations for (ω, t); R = rotation by |ω| about ω
The 6 × 6 system is built from the 32-bit normals, so the transform it yields is computed from 7-digit directions. Before solving, 1e-9 × trace ÷ 6 + 1e-300 is added to the diagonal. The small-angle answer is turned into an exact rotation by rotation4 about the centroid of the matched points; a turn under 1e-15 rad is taken as none.
register.js
lines 325 to 341
CadCraft 0.8.2
const A = Array.from({ length: 6 }, () => new Array(6).fill(0)), bb = new Array(6).fill(0);
for (let k = 0; k < matched; k++) {
const p = [P[k * 3] - c[0], P[k * 3 + 1] - c[1], P[k * 3 + 2] - c[2]], nn = [N[k * 3], N[k * 3 + 1], N[k * 3 + 2]];
const row = [p[1] * nn[2] - p[2] * nn[1], p[2] * nn[0] - p[0] * nn[2], p[0] * nn[1] - p[1] * nn[0], nn[0], nn[1], nn[2]];
const r = (P[k * 3] - Q[k * 3]) * nn[0] + (P[k * 3 + 1] - Q[k * 3 + 1]) * nn[1] + (P[k * 3 + 2] - Q[k * 3 + 2]) * nn[2];
for (let u = 0; u < 6; u++) { for (let v = 0; v < 6; v++) A[u][v] += row[u] * row[v]; bb[u] -= row[u] * r; }
}
const tr = A.reduce((s, r, i) => s + r[i], 0);
for (let u = 0; u < 6; u++) A[u][u] += 1e-9 * tr / 6 + 1e-300;
const x = solveLinear(A, bb);
if (!x || !x.every(Number.isFinite)) throw fail('degenerate');
const ang = Math.hypot(x[0], x[1], x[2]);
const Rm = ang < 1e-15 ? I4() : rotation4([x[0] / ang, x[1] / ang, x[2] / ang], ang);
// the turn is about the matched centroid: p' = Rδ (p − c) + c + t
D = Rm.slice();
const Rc = [Rm[0] * c[0] + Rm[1] * c[1] + Rm[2] * c[2], Rm[4] * c[0] + Rm[5] * c[1] + Rm[6] * c[2], Rm[8] * c[0] + Rm[9] * c[1] + Rm[10] * c[2]];
D[3] = c[0] - Rc[0] + x[3]; D[7] = c[1] - Rc[1] + x[4]; D[11] = c[2] - Rc[2] + x[5];
Solving the 6 × 6 system
64-bit
Solves the linear equations of the point-to-plane step.
Gaussian elimination with partial pivoting, then back substitution
The largest entry of each column is used as the pivot. The system is called singular only when a pivot is below 1e-300, an absolute value that an ill-conditioned system will not reach.
register.js
lines 238 to 250
CadCraft 0.8.2
/** Solve A x = b (n × n, Gaussian elimination with partial pivoting); null when singular. */
function solveLinear(A, b) {
const n = b.length, M = A.map((r, i) => [...r, b[i]]);
for (let c = 0; c < n; c++) {
let p = c; for (let r = c + 1; r < n; r++) if (Math.abs(M[r][c]) > Math.abs(M[p][c])) p = r;
if (Math.abs(M[p][c]) < 1e-300) return null;
[M[c], M[p]] = [M[p], M[c]];
for (let r = c + 1; r < n; r++) { const f = M[r][c] / M[c][c]; for (let k = c; k <= n; k++) M[r][k] -= f * M[c][k]; }
}
const x = new Array(n).fill(0);
for (let r = n - 1; r >= 0; r--) { let s = M[r][n]; for (let k = r + 1; k < n; k++) s -= M[r][k] * x[k]; x[r] = s / M[r][r]; }
return x;
}
When the iteration stops
64-bit
Decides when a phase has settled and narrows the match distance for the next step.
stop when |previous − error| ≤ tol × max(error, 1e-12) or error < 1e-12; cut = min(start, max(2.5 × rms, 1e-6 × diagonal))
tol is 1e-5 unless set, relative to the current error, and each phase has an iteration cap (20 unless set). The result says whether the last phase settled before its cap. The match distance never exceeds its starting value, 5 % of the fixed cloud's diagonal unless given.
Horn's quaternion method for the rigid fit (named in the code's comment); cyclic Jacobi rotations for symmetric eigenproblems; point-to-point and point-to-plane iterative closest point; Gaussian elimination with partial pivoting.
Number types
64-bit numbers and Float64Array for points and matrices; Float32Array for point offsets as read and for surface normals.
Tolerances
Collinear pairs: second eigenvalue < 1e-9 × first. Normals: none when second eigenvalue < 1e-6 × first. Rigid check 1e-6 (1e-5 for a project file), absolute on matrix entries. Jacobi stops at off-diagonal sum of squares < 1e-300 or 60 sweeps. Each phase stops when the error changes by ≤ 1e-5 × the error, or after 20 iterations by default. Match distance starts at 5 % of the fixed cloud's diagonal.
Limits
Iterative closest point uses at most 20 000 moving and 60 000 fixed points by default. The point-to-plane step is a small-angle solution per iteration.
Approach: 32-bit step, flagged
One step stores or exchanges 32-bit values that reach a number a person reads. It is listed here with its open item, and the step is either moved to 64-bit or the app will warn where the number is read.