Slerp¶

An orientation drawn between two steps, alpha from 0 (the step before) to 1 (the last), along the shorter great arc between them: the ends' quaternions as unit 4-vectors, the after end negated where their dot is below 0 (q and -q are the same turn), theta = acos(|d|), and q(alpha) = (sin((1 - alpha) theta) before + sin(alpha theta) after) / sin(theta), a constant angular speed through the arc; ends nearly equal (|d| above 0.9995) are drawn by a normalised lerp, where sin(theta) is too small to divide by.

Drawing only (lusine.render.lerp is the position's), nothing goes back to the simulation. Source: Shoemake, Animating Rotation with Quaternion Curves (SIGGRAPH 1985); Eberly, Quaternion Algebra and Calculus (2002), the nlerp fallback. Checked against a reference in numpy, which has its own tests.

Source

Shoemake, Animating Rotation with Quaternion Curves (SIGGRAPH 1985).

Proof

reference, curve, limit.

slerp.trait

namespace lusine.render;

vector<float, 4> slerp(vector<float, 4> before, vector<float, 4> after, float alpha)
{
  let d0 = dot(before, after);
  let near = d0 < 0.0 ? -after : after;
  let d = abs(d0);
  if (d > 0.9995) return normalize(before + (near - before) * alpha);
  let theta = acos(d);
  let s = sin(theta);
  return before * (sin((1.0 - alpha) * theta) / s) + near * (sin(alpha * theta) / s);
}

An orientation drawn between two steps along the shorter great arc, at a constant angular speed, so a ball turns smoothly between the steps a slow simulation takes. Drawing only.

In [1]:

import os, sys
root = os.getcwd()
while not os.path.isfile(os.path.join(root, 'docsite.json')):
    root = os.path.dirname(root)
sys.path.insert(0, os.path.join(root, 'tools'))
import numpy as np
import matplotlib.pyplot as plt
from lusine import trait_run
F32 = np.float32
U = float(np.finfo(F32).eps) / 2        # binary32's unit roundoff, 2^-24
rng = np.random.default_rng(2026)

import os, sys root = os.getcwd() while not os.path.isfile(os.path.join(root, 'docsite.json')): root = os.path.dirname(root) sys.path.insert(0, os.path.join(root, 'tools')) import numpy as np import matplotlib.pyplot as plt from lusine import trait_run F32 = np.float32 U = float(np.finfo(F32).eps) / 2 # binary32's unit roundoff, 2^-24 rng = np.random.default_rng(2026)

The kernel¶

The law as the library states it.

In [2]:

print(trait_run.source('render/slerp.trait'))

print(trait_run.source('render/slerp.trait'))

// lusine.render.slerp: an orientation drawn between two steps, alpha from 0 (the step before) to 1 (the last), along
// the shorter great arc between them: the ends' quaternions as unit 4-vectors, the after end negated where their dot
// is below 0 (q and -q are the same turn), theta = acos(|d|), and
//   q(alpha) = (sin((1 - alpha) theta) before + sin(alpha theta) after) / sin(theta),
// a constant angular speed through the arc; ends nearly equal (|d| above 0.9995) are drawn by a normalised lerp, where
// sin(theta) is too small to divide by. Drawing only (lusine.render.lerp is the position's), nothing goes back to
// the simulation.
// Source: Shoemake, Animating Rotation with Quaternion Curves (SIGGRAPH 1985); Eberly, Quaternion Algebra and Calculus
// (2002), the nlerp fallback.
// Checked against a reference in numpy, which has its own tests.

namespace lusine.render;

vector<float, 4> slerp(vector<float, 4> before, vector<float, 4> after, float alpha)
{
  let d0 = dot(before, after);
  let near = d0 < 0.0 ? -after : after;
  let d = abs(d0);
  if (d > 0.9995) return normalize(before + (near - before) * alpha);
  let theta = acos(d);
  let s = sin(theta);
  return before * (sin((1.0 - alpha) * theta) / s) + near * (sin(alpha * theta) / s);
}

The program¶

Each section below runs the kernel over a lattice of items, one a case: the host sets every parameter, and main() binds the kernel over the items.

In [3]:

print(trait_run.law_program('render/slerp.trait', 'slerp'))

print(trait_run.law_program('render/slerp.trait', 'slerp'))

namespace probe;

lusine.version version = { { number = 1 } };

struct item {
  host vector<float, 4> before;
  host vector<float, 4> after;
  host float alpha;
  vector<float, 4> result_0;
};
lattice<item, 1> items(spacing = 1 m, boundary = wall);

void main() { items.result_0 = lusine.render.slerp(items.before, items.after, items.alpha); }

In [4]:

from lusine.kernels.render.slerp import slerp as reference
from lusine.kernels.kinematics import _quaternion as Q

from lusine.kernels.render.slerp import slerp as reference from lusine.kernels.kinematics import _quaternion as Q

Reference: arcs at random¶

Ends and alphas at random: the kernel's within binary32's roundings of the reference's (acos and sin are the library's polynomials, each within a few ULP).

In [5]:

a = F32(Q.random(rng, 400)).astype(float)
b = F32(Q.random(rng, 400)).astype(float)
t = F32(rng.random(400)).astype(float)
got = trait_run.law('render/slerp.trait', 'slerp', before=a, after=b, alpha=t)
want = reference(a, b, t[:, None])
err = np.max(np.linalg.norm(got - want, axis=1))
assert err <= 2e-5, err
print(f'{len(a)} arcs within {err:.2g}')

a = F32(Q.random(rng, 400)).astype(float) b = F32(Q.random(rng, 400)).astype(float) t = F32(rng.random(400)).astype(float) got = trait_run.law('render/slerp.trait', 'slerp', before=a, after=b, alpha=t) want = reference(a, b, t[:, None]) err = np.max(np.linalg.norm(got - want, axis=1)) assert err <= 2e-5, err print(f'{len(a)} arcs within {err:.2g}')

400 arcs within 1.4e-07

Curve: a constant angular speed¶

Along a turn of 2 rad about +y, the angle drawn at alpha is 2 alpha.

In [6]:

alphas = np.linspace(0.0, 1.0, 51)
end = Q.about([0, 1, 0], 2.0)
q = trait_run.law('render/slerp.trait', 'slerp', before=np.array([0.0, 0, 0, 1]), after=end, alpha=alphas)
angle = 2 * np.arccos(np.clip(q[:, 3], -1, 1))
assert np.max(np.abs(angle - 2 * alphas)) <= 1e-5
fig, ax = plt.subplots(figsize=(6, 3)); ax.plot(alphas, angle)
ax.set_xlabel('alpha'); ax.set_ylabel('angle drawn (rad)'); plt.show()

alphas = np.linspace(0.0, 1.0, 51) end = Q.about([0, 1, 0], 2.0) q = trait_run.law('render/slerp.trait', 'slerp', before=np.array([0.0, 0, 0, 1]), after=end, alpha=alphas) angle = 2 * np.arccos(np.clip(q[:, 3], -1, 1)) assert np.max(np.abs(angle - 2 * alphas)) <= 1e-5 fig, ax = plt.subplots(figsize=(6, 3)); ax.plot(alphas, angle) ax.set_xlabel('alpha'); ax.set_ylabel('angle drawn (rad)'); plt.show()

Limit: the ends, and the shorter arc¶

At 0 the step before's, at 1 the last's; an end given as -q draws the same turns as q.

In [7]:

e0 = trait_run.law('render/slerp.trait', 'slerp', before=a, after=b, alpha=0.0)
e1 = trait_run.law('render/slerp.trait', 'slerp', before=a, after=b, alpha=1.0)
sign = np.sign(np.sum(a * b, axis=1, keepdims=True))
assert np.max(np.abs(e0 - a)) <= 1e-6 and np.max(np.abs(e1 - b * sign)) <= 1e-6
flip = trait_run.law('render/slerp.trait', 'slerp', before=a, after=-b, alpha=t)
assert np.max(np.abs(flip - got)) <= 1e-6
print('the ends; the shorter arc either way')

e0 = trait_run.law('render/slerp.trait', 'slerp', before=a, after=b, alpha=0.0) e1 = trait_run.law('render/slerp.trait', 'slerp', before=a, after=b, alpha=1.0) sign = np.sign(np.sum(a * b, axis=1, keepdims=True)) assert np.max(np.abs(e0 - a)) <= 1e-6 and np.max(np.abs(e1 - b * sign)) <= 1e-6 flip = trait_run.law('render/slerp.trait', 'slerp', before=a, after=-b, alpha=t) assert np.max(np.abs(flip - got)) <= 1e-6 print('the ends; the shorter arc either way')

the ends; the shorter arc either way