Turning a vector¶
A vector turned by a body's orientation, R(q) v, from the body's frame into the world's: q a unit quaternion (x, y, z, w), its vector part u = (x, y, z) and scalar part w, and with t = 2 (u x v), R(q) v = v + w t + u x t.
Its inverse is lusine.kinematics.unrotate. Any quantity's vector turns so (a position in metres, a gradient per metre). q is a unit quaternion, the caller's to keep so (lusine.integrator.orientation does). Source: Shoemake, Animating Rotation with Quaternion Curves (SIGGRAPH 1985); Kuipers, Quaternions and Rotation Sequences (1999), 5.13, the rotation operator q v q*. Checked against a reference in numpy, which has its own tests.
- Source
- Kuipers, Quaternions and Rotation Sequences (1999), 5.13.
- Proof
- reference, invariant, limit.
rotate.trait
namespace lusine.kinematics;
vector<float<_>, 3> rotate(vector<float, 4> orientation, vector<float<_>, 3> value)
{
let u = vector<float, 3>(orientation.x, orientation.y, orientation.z);
let t = cross(u, value) * 2.0;
return (value + t * orientation.w) + cross(u, t);
}
A vector turned by a body's orientation from the body's frame into the world's, R(q) v, q a unit quaternion (x, y, z, w). The render turns a shape's layer gradients back to the world by it.
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.
print(trait_run.source('kinematics/rotate.trait'))
// lusine.kinematics.rotate: a vector turned by a body's orientation, R(q) v, from the body's frame into the world's: q
// a unit quaternion (x, y, z, w), its vector part u = (x, y, z) and scalar part w, and with t = 2 (u x v),
// R(q) v = v + w t + u x t. Its inverse is lusine.kinematics.unrotate. Any quantity's vector turns so (a position in
// metres, a gradient per metre). q is a unit quaternion, the caller's to keep so (lusine.integrator.orientation does).
// Source: Shoemake, Animating Rotation with Quaternion Curves (SIGGRAPH 1985); Kuipers, Quaternions and Rotation
// Sequences (1999), 5.13, the rotation operator q v q*.
// Checked against a reference in numpy, which has its own tests.
namespace lusine.kinematics;
vector<float<_>, 3> rotate(vector<float, 4> orientation, vector<float<_>, 3> value)
{
let u = vector<float, 3>(orientation.x, orientation.y, orientation.z);
let t = cross(u, value) * 2.0;
return (value + t * orientation.w) + cross(u, t);
}
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.
print(trait_run.law_program('kinematics/rotate.trait', 'rotate'))
namespace probe;
lusine.version version = { { number = 1 } };
struct item {
host vector<float, 4> orientation;
host vector<float<_>, 3> value;
vector<float<_>, 3> result_0;
};
lattice<item, 1> items(spacing = 1 m, boundary = wall);
void main() { items.result_0 = lusine.kinematics.rotate(items.orientation, items.value); }
from lusine.kernels.kinematics.rotate import rotate as reference
from lusine.kernels.kinematics import _quaternion as Q
qs = F32(Q.random(rng, 400)).astype(float)
vs = F32(rng.normal(0.0, 2.0, (400, 3))).astype(float)
Reference: turns at random¶
The kernel's turned vector within binary32's roundings of the reference's, which goes through the rotation matrix, not the kernel's cross products.
got = trait_run.law('kinematics/rotate.trait', 'rotate', orientation=qs, value=vs)
err = np.max(np.linalg.norm(got - reference(qs, vs), axis=1) / np.linalg.norm(vs, axis=1))
assert err <= 32 * U, err
print(f'{len(qs)} turns within {err / U:.1f} unit roundoffs')
400 turns within 3.7 unit roundoffs
Invariant: lengths kept, and q and -q one turn¶
A turn keeps a vector's length, and a quaternion and its negative turn alike.
assert np.max(np.abs(np.linalg.norm(got, axis=1) / np.linalg.norm(vs, axis=1) - 1)) <= 32 * U
neg = trait_run.law('kinematics/rotate.trait', 'rotate', orientation=-qs, value=vs)
assert np.max(np.abs(neg - got)) <= 1e-6
print('lengths kept; q and -q the same turn')
lengths kept; q and -q the same turn
Limit: a quarter turn, and none¶
A quarter turn about +z takes +x to +y; the identity leaves a vector as it is.
quarter = Q.about([0, 0, 1], np.pi / 2)
r = trait_run.law('kinematics/rotate.trait', 'rotate', orientation=np.array([quarter, [0, 0, 0, 1.0]]),
value=np.array([[1.0, 0, 0], [0.3, -2.0, 5.0]]))
assert np.max(np.abs(r[0] - [0, 1, 0])) <= 4 * U and np.array_equal(r[1], F32([0.3, -2.0, 5.0]).astype(float))
print('+x to +y; the identity, nothing')
+x to +y; the identity, nothing