Resources Numerical accuracy Calculation policy

Inverse kinematics and constraints

Moves a bone or a placement so it reaches a target, looks at it, copies it, stays within limits, hangs from it or sits on a path.

How the precision is kept

Matrices are plain arrays of 16 64-bit numbers and turns are quaternions. A two-bone chain is solved in closed form by the law of cosines; longer chains are solved by FABRIK, an iteration that stops at a fixed distance from the target or after a set number of sweeps. A constraint's result is blended into the pose by its influence.

Calculations under this policy

11 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.

anim/constraints.js ModelCraft 0.8.1 64-bit 9 calculations Reviewed

Tests bin/test-animation.mjs

Shortest turn between two directions

64-bit

The quaternion that turns one unit vector onto another.

q = (u × v, 1 + u · v), made unit

Built from the cross product and 1 + the dot product, so no arc cosine is taken and no clamp is needed. Directions with a dot product above 1 − 1e-12 return no turn; below −1 + 1e-12 a half turn about a perpendicular axis is returned.

  • anim/constraints.js
  • lines 44 to 55
  • ModelCraft 0.8.1
/** The shortest turn taking unit vector u onto unit vector v (a quaternion). */
export function swing(u, v) {
  const d = dot(u, v);
  if (d > 1 - 1e-12) return [0, 0, 0, 1];
  if (d < -1 + 1e-12) {
    // half a turn about any axis square to u
    const ax = norm(Math.abs(u[0]) < 0.9 ? cross(u, [1, 0, 0]) : cross(u, [0, 1, 0]));
    return [ax[0], ax[1], ax[2], 0];
  }
  const c = cross(u, v), w = 1 + d, l = Math.hypot(c[0], c[1], c[2], w);
  return [c[0] / l, c[1] / l, c[2] / l, w / l];
}

Two-bone reach

64-bit

Where the elbow and the end of a two-bone chain go so the end reaches a target.

cos A = (L₁² + d² − L₂²) ÷ (2 L₁ d); sin A = √(1 − cos² A); elbow = a + dir × L₁ cos A + n × L₁ sin A

The distance is clamped to [|L₁ − L₂| + 1e-9, L₁ + L₂ − 1e-9] and cos A to [−1, 1]. √(1 − cos² A) subtracts nearly equal numbers when the chain is almost straight or folded. The division by 2 L₁ d is not guarded against a first bone of zero length.

  • anim/constraints.js
  • lines 63 to 75
  • ModelCraft 0.8.1
/** Two bones from a with lengths L1, L2 reaching for target: the elbow (bent towards poleDir) and the end. */
export function twoBoneIK(a, L1, L2, target, poleDir) {
  const toT = sub(target, a), dRaw = len(toT);
  const dir = dRaw > 1e-12 ? scl(toT, 1 / dRaw) : [0, 0, 1];
  const d = Math.min(Math.max(dRaw, Math.abs(L1 - L2) + 1e-9), L1 + L2 - 1e-9);
  // the bend: poleDir made square to the line, else any square direction
  let n = sub(poleDir || [0, 0, 0], scl(dir, dot(poleDir || [0, 0, 0], dir)));
  if (len(n) < 1e-9) n = Math.abs(dir[2]) < 0.9 ? cross(dir, [0, 0, 1]) : cross(dir, [1, 0, 0]);
  n = norm(n);
  const cosA = Math.min(1, Math.max(-1, (L1 * L1 + d * d - L2 * L2) / (2 * L1 * d))), sinA = Math.sqrt(1 - cosA * cosA);
  const elbow = add(a, add(scl(dir, L1 * cosA), scl(n, L1 * sinA)));
  return { elbow, end: add(a, scl(dir, d)) };
}

Chain of three or more bones (FABRIK)

64-bit

Joint positions of a longer chain whose tip reaches for a target.

backward sweep from the target, forward sweep from the root, each joint kept its segment length from its neighbor

Stops when the tip is within 1e-7 of the target (absolute) or after the given number of sweeps, 12 by default. A target at or beyond the summed bone lengths gives a straight chain pointing at it.

  • anim/constraints.js
  • lines 76 to 94
  • ModelCraft 0.8.1
/** FABRIK: joints [p0 … pn] (p0 stays), segment lengths, a target for pn. → new joints. */
export function fabrik(joints, lengths, target, iterations = 12) {
  const p = joints.map((x) => x.slice()), n = p.length - 1, root = p[0].slice();
  const total = lengths.reduce((s, x) => s + x, 0);
  if (len(sub(target, root)) >= total) {
    // out of reach: straight at it
    const d = norm(sub(target, root));
    for (let i = 1; i <= n; i++) p[i] = add(p[i - 1], scl(d, lengths[i - 1]));
    return p;
  }
  for (let it = 0; it < iterations; it++) {
    p[n] = target.slice();
    for (let i = n - 1; i >= 0; i--) p[i] = add(p[i + 1], scl(norm(sub(p[i], p[i + 1])), lengths[i]));
    p[0] = root.slice();
    for (let i = 1; i <= n; i++) p[i] = add(p[i - 1], scl(norm(sub(p[i], p[i - 1])), lengths[i - 1]));
    if (len(sub(p[n], target)) < 1e-7) break;
  }
  return p;
}

Point along a path by length

64-bit

The point and direction a given fraction of the way along a polyline.

s = k × Σ segment lengths; walk the segments subtracting each length until s falls inside one

Segment lengths are added in order with Math.hypot, and the remaining distance is found by subtracting them one by one. The fraction is clamped to [0, 1] and the last segment takes any remainder left by rounding.

  • anim/constraints.js
  • lines 95 to 113
  • ModelCraft 0.8.1
/** A point k (0…1) of the way along a polyline by length, and the direction there. */
export function alongPath(pts, closed, k) {
  const n = pts.length / 3;
  if (n < 2) return null;
  const P = (i) => [pts[(i % n) * 3], pts[(i % n) * 3 + 1], pts[(i % n) * 3 + 2]];
  const segs = closed ? n : n - 1, L = [];
  let total = 0;
  for (let i = 0; i < segs; i++) { const l = len(sub(P(i + 1), P(i))); L.push(l); total += l; }
  if (total <= 0) return { p: P(0), d: [0, 1, 0] };
  let s = Math.min(1, Math.max(0, k)) * total;
  for (let i = 0; i < segs; i++) {
    if (s <= L[i] || i === segs - 1) {
      const f = L[i] > 0 ? Math.min(1, s / L[i]) : 0, a = P(i), b = P(i + 1);
      return { p: add(a, scl(sub(b, a), f)), d: norm(sub(b, a)) };
    }
    s -= L[i];
  }
  return null;
}

Point an axis at a direction

64-bit

A matrix whose chosen axis points along a direction, with its up axis leaning to world +Z.

Two shortest turns as quaternions, applied to the decomposed matrix and composed again. The roll correction is skipped when either projected up vector is shorter than 1e-6.

  • anim/constraints.js
  • lines 114 to 124
  • ModelCraft 0.8.1
/** A matrix whose `axis` points along d (up kept as near `upAxis` → world +Z as it can), its place and scale from m. */
function aim(m, d, axis = '+y', upAxis = '+z') {
  const D = decompose(m), cur = norm(point(compose({ t: [0, 0, 0], q: D.q, s: [1, 1, 1] }), AXIS[axis]));
  const q1 = swing(cur, norm(d));
  let q = qmul(q1, D.q);
  // then turn about d so the up axis leans to world +Z
  const R = compose({ t: [0, 0, 0], q, s: [1, 1, 1] }), up = point(R, AXIS[upAxis]), dn = norm(d);
  const want = sub([0, 0, 1], scl(dn, dot([0, 0, 1], dn))), have = sub(up, scl(dn, dot(up, dn)));
  if (len(want) > 1e-6 && len(have) > 1e-6 && axis[1] !== upAxis[1]) q = qmul(swing(norm(have), norm(want)), q);
  return compose({ t: D.t, q, s: D.s });
}

Blend by influence

64-bit

Mixes two matrices by a factor from 0 to 1: place and size linearly, the turn along the shortest arc.

t = t_a + (t_b − t_a) k; s = s_a + (s_b − s_a) k; q = slerp(q_a, q_b, k)

A factor of 1 or more returns the second matrix unchanged and 0 or less the first, so full influence adds no decompose-and-compose rounding. Other factors go through decompose, slerp and compose from render/anim.js.

  • anim/constraints.js
  • lines 125 to 131
  • ModelCraft 0.8.1
/** Blend two matrices (place, turn, size) by k. */
function blend(a, b, k) {
  if (k >= 1) return b.slice();
  if (k <= 0) return a.slice();
  const A = decompose(a), B = decompose(b);
  return compose({ t: A.t.map((v, i) => v + (B.t[i] - v) * k), q: slerp(A.q, B.q, k), s: A.s.map((v, i) => v + (B.s[i] - v) * k) });
}

IK on a rig

64-bit

Builds the chain for an IK constraint, solves it and swings each bone so it runs to the next joint.

Joints are the bones' heads in the rig's mesh coordinates. Each bone is turned about its head by the shortest turn from its current Y axis to the solved direction, root first, and the world matrices are rebuilt after every bone.

  • anim/constraints.js
  • lines 173 to 192
  • ModelCraft 0.8.1
if (c.type === 'ik') {
  const T = tgtPoint();
  if (!T) continue;
  const chain = [b];
  while (chain.length < c.chain && chain[0].parent && byId.has(chain[0].parent)) chain.unshift(byId.get(chain[0].parent));
  const joints = [...chain.map((x) => head(x.id)), tip(b.id)], lengths = chain.map(lengthOf);
  let next;
  if (chain.length === 2) {
    const pole = ptOf(c.pole);
    const bendNow = sub(joints[1], add(joints[0], scl(norm(sub(joints[2], joints[0])), dot(sub(joints[1], joints[0]), norm(sub(joints[2], joints[0]))))));
    const r = twoBoneIK(joints[0], lengths[0], lengths[1], T, pole ? sub(pole, joints[0]) : bendNow);
    next = [joints[0], r.elbow, r.end];
  } else if (chain.length === 1) next = [joints[0], add(joints[0], scl(norm(sub(T, joints[0])), lengths[0]))];
  else next = fabrik(joints, lengths, T, c.iterations);
  // swing each bone of the chain, root first, so its Y runs to the next joint
  chain.forEach((x, i) => {
    const W = world.get(x.id), from = norm(col(W, 1)), to = norm(sub(next[i + 1], head(x.id)));
    setWorld(x, turnAbout(W, swing(from, to), head(x.id)), k);
  });
} else if (c.type === 'lookat') {

World matrix back to a local pose

64-bit

Works a bone's local matrix back from a wanted world matrix through its parent.

local = (parent world × parent rest⁻¹ × rest)⁻¹ × wanted world

Uses the general 4×4 inverse from anim/rig.js. Every constraint adds one inverse and several products to the chain of matrices, all in 64-bit.

  • anim/constraints.js
  • lines 147 to 154
  • ModelCraft 0.8.1
const parentRel = (b, world) => (b.parent && world.has(b.parent) ? mul(world.get(b.parent), mul(inv(rest.get(b.parent)), rest.get(b.id))) : rest.get(b.id));
const worlds = () => {
  const w = new Map();
  for (const b of ord) w.set(b.id, mul(parentRel(b, w), locals.get(b.id)));
  return w;
};
let world = worlds();
/** Give bone b a new world matrix (by influence k): its local is worked back from its parent's. */

Constraints on placements

64-bit

Applies look-at, copy, limit, child-of and follow-path constraints to placements in world coordinates and returns their new local matrices.

new local = blend(local, parent world⁻¹ × constrained world, influence)

Owners whose targets are other constrained placements are solved after them. The constrained world matrix is carried back through the inverse of the parent's world matrix in 64-bit.

  • anim/constraints.js
  • lines 215 to 260
  • ModelCraft 0.8.1
// ---- placements ----
export function solveObjects(anim, poses, vals, ctx) {
  const out = new Map();
  const all = (anim.constraints || []).filter((c) => c.on && c.owner.startsWith('obj:'));
  if (!all.length) return out;
  // owners whose targets are other owners come after them (a placement looking at one that follows a path sees it on the
  // path); each owner's own constraints keep their order; a loop keeps the list's order
  const owners = [...new Set(all.map((c) => c.owner.slice(4)))], deps = new Map(owners.map((o) => [o, new Set()]));
  for (const c of all) for (const x of [c.target, c.pole]) if (x && x.kind === 'obj') { const k = pathKey(x.path); if (deps.has(k) && k !== c.owner.slice(4)) deps.get(c.owner.slice(4)).add(k); }
  const sorted = [], state = new Map();
  const visit = (o) => { if (state.get(o) === 2 || state.get(o) === 1) return; state.set(o, 1); for (const d of deps.get(o)) visit(d); state.set(o, 2); sorted.push(o); };
  for (const o of owners) visit(o);
  const cons = sorted.flatMap((o) => all.filter((c) => c.owner.slice(4) === o));
  const cur = new Map(poses);
  const localOf = (key, path) => cur.get(key) || ctx.local(path);
  for (const c of cons) {
    const key = c.owner.slice(4), path = key.split('/').map(Number);
    const k = vals && vals.has('cns:' + c.id + '|influence') ? vals.get('cns:' + c.id + '|influence')[0] : c.influence;
    if (!(k > 0)) continue;
    const L = localOf(key, path);
    if (!L) continue;
    const PW = ctx.parentWorld(path, cur), W = mul(PW, L);
    let N = null;
    if (c.type === 'lookat') { const T = ctx.point(c.target, cur); if (T) N = aim(W, sub(T, col(W, 3)), c.axis, c.up); }
    else if (c.type === 'copy') {
      const M = ctx.matrix(c.target, cur);
      if (M) { const A = decompose(W), B = decompose(M); N = compose({ t: c.mix.location ? B.t : A.t, q: c.mix.rotation ? B.q : A.q, s: c.mix.scale ? B.s : A.s }); }
    } else if (c.type === 'limit') {
      const x = clampTRS(L, c.limits, null);
      cur.set(key, blend(L, compose({ t: x.location, q: eulerToQuat(x.rotation), s: x.scale }), k));
      out.set(key, cur.get(key));
      continue;
    } else if (c.type === 'childof') { const M = ctx.matrix(c.target, cur); if (M) N = mul(M, mul(c.inverse || inv(M), W)); }
    else if (c.type === 'path') {
      const cv = c.target && c.target.kind === 'obj' ? ctx.curve(c.target.path, cur) : null;
      const off = vals && vals.has('cns:' + c.id + '|offset') ? vals.get('cns:' + c.id + '|offset')[0] : c.offset;
      const at = cv && alongPath(cv.pts, cv.closed, off);
      if (at) { const D = decompose(W); N = compose({ t: at.p, q: D.q, s: D.s }); if (c.follow) N = aim(N, at.d, c.axis, c.up); }
    }
    if (!N) continue;
    const next = blend(L, mul(inv(PW), N), k);
    cur.set(key, next);
    out.set(key, next);
  }
  return out;
}
anim/cns-model.js ModelCraft 0.8.1 64-bit 2 calculations Reviewed

Tests bin/test-animation.mjs

World matrix of a placement

64-bit

The product of the matrices of a placement and everything it sits in.

W = M₁ × M₂ × … × M_own

One 4×4 product per level, in 64-bit, starting from the identity. A constrained placement starts from its stored rest matrix, not from the matrix the last frame left on the model.

  • anim/cns-model.js
  • lines 22 to 23
  • ModelCraft 0.8.1
const parentWorld = (path, cur) => { let W = IDENT; for (let i = 1; i < path.length; i++) { const m = own(path.slice(0, i), cur); if (m) W = mul(W, m); } return W; };
const worldOf = (path, cur) => { const m = own(path, cur); return m ? mul(parentWorld(path, cur), m) : null; };

A curve object's line in world coordinates

64-bit

The polyline a follow-path constraint runs along.

Each point of the curve's line segments is multiplied by the world matrix and kept in a Float64Array. The curve counts as closed when its first and last points are within 1e-9 (absolute).

  • anim/cns-model.js
  • lines 40 to 50
  • ModelCraft 0.8.1
curve(path, cur) {
  const p = model.placementAt(path), def = p && model.defs.get(p.inst.def);
  const L = def && def.curve && MESH_HOOKS.lines ? MESH_HOOKS.lines(def) : null;
  if (!L || L.length < 6) return null;
  const W = worldOf(path, cur), tr = (i) => [W[0] * L[i] + W[4] * L[i + 1] + W[8] * L[i + 2] + W[12], W[1] * L[i] + W[5] * L[i + 1] + W[9] * L[i + 2] + W[13], W[2] * L[i] + W[6] * L[i + 1] + W[10] * L[i + 2] + W[14]];
  // segments a-b, a-b …: chain them (the next segment starts where one ends)
  const pts = [...tr(0)];
  for (let i = 0; i < L.length; i += 6) pts.push(...tr(i + 3));
  const n = pts.length / 3, closed = n > 2 && Math.hypot(pts[0] - pts[(n - 1) * 3], pts[1] - pts[(n - 1) * 3 + 1], pts[2] - pts[(n - 1) * 3 + 2]) < 1e-9;
  return { pts: Float64Array.from(closed ? pts.slice(0, -3) : pts), closed };
},

In detail

Method
Law of cosines for two bones; FABRIK for longer chains; the shortest turn between two unit vectors as a quaternion without an arc cosine; position along a polyline by summed segment lengths.
Number types
64-bit numbers; a curve's polyline is a Float64Array.
Tolerances
FABRIK stops when the tip is within 1e-7 of the target (absolute, model units) or after the constraint's number of sweeps (12 by default). Two-bone reach is kept 1e-9 inside its limits (absolute). Parallel and opposite directions are detected at 1 − 1e-12 of the dot product. A curve counts as closed when its ends are within 1e-9.
Checked by
bin/test-animation.mjs loads the development build's copy of anim/constraints.js and calls twoBoneIK, fabrik, swing, alongPath, solveRig and solveObjects.
Limits
FABRIK is an iteration: the tip is as close to the target as its stop rule, not exact. The absolute tolerances do not scale with the size of the rig.

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.

Every policy that takes this approach