Resources Numerical accuracy Calculation policy
Works out each bone's rest frame from its head, tail and roll, and the matrices that carry the mesh from rest to a pose.
Matrices are arrays of 16 64-bit numbers in column order. A bone's frame comes from its head and tail by cross products; a pose is composed down the bone tree; the skin matrix is the posed frame times the inverse of the rest frame. The inverse is the general 4×4 cofactor form.
6 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 bin/test-animation.mjs
The product of two column-major 4×4 matrices.
Sixteen sums of four 64-bit products each, written to a new plain array.
/** Column-major 4×4 product a × b. */
export function mul(a, b) {
const o = new Array(16);
for (let c = 0; c < 4; c++) for (let r = 0; r < 4; r++) o[c * 4 + r] = a[r] * b[c * 4] + a[4 + r] * b[c * 4 + 1] + a[8 + r] * b[c * 4 + 2] + a[12 + r] * b[c * 4 + 3];
return o;
}The inverse of a 4×4 matrix by cofactors.
M⁻¹ = adjugate(M) ÷ det(M)
The determinant is a sum of products of 2×2 minors and can lose digits when the matrix is nearly singular. Only a determinant of exactly 0 returns the identity; there is no tolerance scaled to the matrix.
/** Inverse of a rigid (rotation + translation) or rotation-scale-translation matrix, the general way. */
export function inv(m) {
const [a00, a01, a02, a03, a10, a11, a12, a13, a20, a21, a22, a23, a30, a31, a32, a33] = m;
const b00 = a00 * a11 - a01 * a10, b01 = a00 * a12 - a02 * a10, b02 = a00 * a13 - a03 * a10, b03 = a01 * a12 - a02 * a11;
const b04 = a01 * a13 - a03 * a11, b05 = a02 * a13 - a03 * a12, b06 = a20 * a31 - a21 * a30, b07 = a20 * a32 - a22 * a30;
const b08 = a20 * a33 - a23 * a30, b09 = a21 * a32 - a22 * a31, b10 = a21 * a33 - a23 * a31, b11 = a22 * a33 - a23 * a32;
let d = b00 * b11 - b01 * b10 + b02 * b09 + b03 * b08 - b04 * b07 + b05 * b06;
if (!d) return [1, 0, 0, 0, 0, 1, 0, 0, 0, 0, 1, 0, 0, 0, 0, 1];
d = 1 / d;
return [
(a11 * b11 - a12 * b10 + a13 * b09) * d, (a02 * b10 - a01 * b11 - a03 * b09) * d, (a31 * b05 - a32 * b04 + a33 * b03) * d, (a22 * b04 - a21 * b05 - a23 * b03) * d,
(a12 * b08 - a10 * b11 - a13 * b07) * d, (a00 * b11 - a02 * b08 + a03 * b07) * d, (a32 * b02 - a30 * b05 - a33 * b01) * d, (a20 * b05 - a22 * b02 + a23 * b01) * d,
(a10 * b10 - a11 * b08 + a13 * b06) * d, (a01 * b08 - a00 * b10 - a03 * b06) * d, (a30 * b04 - a31 * b02 + a33 * b00) * d, (a21 * b02 - a20 * b04 - a23 * b00) * d,
(a11 * b07 - a10 * b09 - a12 * b06) * d, (a00 * b09 - a01 * b07 + a02 * b06) * d, (a31 * b01 - a30 * b03 - a32 * b00) * d, (a20 * b03 - a21 * b01 + a22 * b00) * d,
];
}The matrix from a bone's own frame to the mesh: Y from head to tail, Z as near the mesh's +Z as it can be, turned by the roll.
y = (tail − head) made unit; z = (up − y (up · y)) made unit; x = y × z; roll turns x and z about y
The up reference switches from +Z to +Y when |y_z| is 0.999 or more. A bone whose head and tail coincide gets a zero Y axis, because the length falls back to 1. The roll angle is degrees × π ÷ 180.
/** A bone's rest frame → mesh (see the header). */
export function boneRest(b) {
const y = norm(sub(b.tail, b.head));
const up = Math.abs(y[2]) < 0.999 ? [0, 0, 1] : [0, 1, 0];
let z = norm(sub(up, y.map((v) => v * dot(up, y))));
let x = cross(y, z);
if (b.roll) {
// turn X and Z about Y
const r = (b.roll * Math.PI) / 180, c = Math.cos(r), s = Math.sin(r);
const x2 = [x[0] * c - z[0] * s, x[1] * c - z[1] * s, x[2] * c - z[2] * s], z2 = [x[0] * s + z[0] * c, x[1] * s + z[1] * c, x[2] * s + z[2] * c];
x = x2; z = z2;
}
return [x[0], x[1], x[2], 0, y[0], y[1], y[2], 0, z[0], z[1], z[2], 0, b.head[0], b.head[1], b.head[2], 1];
}Each bone's posed frame and the matrix that moves rest points to the pose.
W = W_parent × R_parent⁻¹ × R × local; skin = W × R⁻¹
Parents are done before children. Each level adds two inverses and three products, so rounding grows with the depth of the bone chain.
export function poseMatrices(rig, poseOf = () => REST_POSE) {
const rest = restMatrices(rig), world = new Map(), skin = new Map();
for (const b of order(rig)) {
const p = poseOf(b.id) || REST_POSE;
const local = compose({ t: p.location || [0, 0, 0], q: eulerToQuat(p.rotation || [0, 0, 0]), s: p.scale || [1, 1, 1] });
const R = rest.get(b.id);
let W;
if (b.parent && world.has(b.parent)) W = mul(mul(world.get(b.parent), mul(inv(rest.get(b.parent)), R)), local);
else W = mul(R, local);
world.set(b.id, W);
skin.set(b.id, mul(W, inv(R)));
}
return { world, skin, rest };
}Each bone's turn and move as a unit dual quaternion, for dual-quaternion skinning.
dual part = ½ × (t, 0) × q
The matrix is decomposed by render/anim.js and its scale is left out. The eight numbers per bone are kept in a Float64Array.
/** Each bone's skin matrix as a unit dual quaternion [qx, qy, qz, qw, dx, dy, dz, dw] (its turn and move; its size is left out). */
export function dualQuats(M, nb) {
const out = new Float64Array(nb * 8);
for (let b = 0; b < nb; b++) {
const m = Array.from(M.subarray(b * 16, b * 16 + 16)), d = decompose(m), q = d.q, t = d.t;
// dual part = ½ · (t, 0) · q
const dx = 0.5 * (t[0] * q[3] + t[1] * q[2] - t[2] * q[1]), dy = 0.5 * (-t[0] * q[2] + t[1] * q[3] + t[2] * q[0]);
const dz = 0.5 * (t[0] * q[1] - t[1] * q[0] + t[2] * q[3]), dw = -0.5 * (t[0] * q[0] + t[1] * q[1] + t[2] * q[2]);
out.set([q[0], q[1], q[2], q[3], dx, dy, dz, dw], b * 8);
}
return out;
}Tests None found by the review
The location or rotation a bone gets when it is dragged in the view.
turn: angle swept round the bone's head on screen, about the view axis carried into the bone's frame; move: pixels × pixel size, along the camera's right and up, carried into the bone's frame
The angle comes from two Math.atan2 calls on pointer positions, so it has the resolution of the pointer. The new pose is always worked out from the pose at the start of the drag, so the steps of one drag do not add up.
/** The pose the drag gives: a turn about the view's axis through the bone's head, or a move in the view's plane. */
posed(d, ev) {
const { rg, s, base, W } = d;
const { world } = poseMatrices(rg, (id) => bonePose(this.tl.values, rg.id, id));
const Wb = world.get(s.bone);
let p;
if (d.move) {
// pixels → mesh units at the bone's head, along the camera's right and up, into the bone's own frame
const head = apply(W, point(Wb, [0, 0, 0])), px = this.h.cam.pixelSize(head);
const cam = this.h.cam, fwd = norm([cam.target[0] - cam.eye[0], cam.target[1] - cam.eye[1], cam.target[2] - cam.eye[2]]);
const upW = Math.abs(fwd[2]) < 0.999 ? [0, 0, 1] : [0, 1, 0];
const right = norm([fwd[1] * upW[2] - fwd[2] * upW[1], fwd[2] * upW[0] - fwd[0] * upW[2], fwd[0] * upW[1] - fwd[1] * upW[0]]);
const up = [right[1] * fwd[2] - right[2] * fwd[1], right[2] * fwd[0] - right[0] * fwd[2], right[0] * fwd[1] - right[1] * fwd[0]];
const dw = [0, 1, 2].map((i) => ((ev.x - d.x) * right[i] - (ev.y - d.y) * up[i]) * px);
const dm = applyDir(invert(W), dw);
// a move in the mesh → the bone's local frame (its pose's location is in its rest frame, parent's pose applied)
const local = invDir(mul(Wb, inv(compose({ t: base.location, q: eulerToQuat(base.rotation), s: base.scale }))), dm);
p = { location: base.location.map((v, i) => v + local[i]), rotation: base.rotation, scale: base.scale };
} else {
// the turn: the angle swept round the bone's head on screen
const h = this.h.cam.project(apply(W, point(Wb, [0, 0, 0])));
if (!h) return null;
const a0 = Math.atan2(d.y - h[1], d.x - h[0]), a1 = Math.atan2(ev.y - h[1], ev.x - h[0]);
const ang = -(a1 - a0); // screen y runs down
const cam = this.h.cam, axisW = norm([cam.eye[0] - cam.target[0], cam.eye[1] - cam.target[1], cam.eye[2] - cam.target[2]]);
const axisM = norm(applyDir(invert(W), axisW));
// the turn in the mesh, carried into the bone's frame: X = Wb⁻¹ · R · Wb (about the bone's own origin, its head)
const R = compose({ t: [0, 0, 0], q: [axisM[0] * Math.sin(ang / 2), axisM[1] * Math.sin(ang / 2), axisM[2] * Math.sin(ang / 2), Math.cos(ang / 2)], s: [1, 1, 1] });
const Wr = Wb.slice(); Wr[12] = Wr[13] = Wr[14] = 0;
const X = mul(inv(Wr), mul(R, Wr));
const L0 = compose({ t: base.location, q: eulerToQuat(base.rotation), s: base.scale });
const L1 = mul(L0, X), dd = decompose(L1);
p = { location: base.location, rotation: quatToEuler(dd.q, base.rotation), scale: base.scale };
}
return keyOrPose(this.tl, d.start, rg.id, s.bone, p);
}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.