Turning a vector back¶

A vector turned back by a body's orientation, R(q)^T v, from the world's frame into the body's: lusine.kinematics.rotate by q's conjugate (-x, -y, -z, w), with t = 2 (-u x v), R^T v = v + w t - u x t.

A point p of the world is at R^T (p - c) in the frame of a body centred at c. Source: Kuipers, Quaternions and Rotation Sequences (1999), 5.13. Checked against a reference in numpy, which has its own tests.

Source

Kuipers, Quaternions and Rotation Sequences (1999), 5.13.

Proof

reference, invariant.

unrotate.trait

namespace lusine.kinematics;

vector<float<_>, 3> unrotate(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 back by a body's orientation, from the world's frame into the body's, R(q)^T v: the render evaluates a shape's layers at the hit point less its centre turned back so, and they turn with the ball.

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('kinematics/unrotate.trait'))

print(trait_run.source('kinematics/unrotate.trait'))

// lusine.kinematics.unrotate: a vector turned back by a body's orientation, R(q)^T v, from the world's frame into the
// body's: lusine.kinematics.rotate by q's conjugate (-x, -y, -z, w), with t = 2 (-u x v), R^T v = v + w t - u x t. A
// point p of the world is at R^T (p - c) in the frame of a body centred at c.
// Source: Kuipers, Quaternions and Rotation Sequences (1999), 5.13.
// Checked against a reference in numpy, which has its own tests.

namespace lusine.kinematics;

vector<float<_>, 3> unrotate(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.

In [3]:

print(trait_run.law_program('kinematics/unrotate.trait', 'unrotate'))

print(trait_run.law_program('kinematics/unrotate.trait', 'unrotate'))

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.unrotate(items.orientation, items.value); }

In [4]:

from lusine.kernels.kinematics.unrotate import unrotate 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)

from lusine.kernels.kinematics.unrotate import unrotate 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¶

Within binary32's roundings of the reference's transpose of the rotation matrix.

In [5]:

got = trait_run.law('kinematics/unrotate.trait', 'unrotate', 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')

got = trait_run.law('kinematics/unrotate.trait', 'unrotate', 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.4 unit roundoffs

Invariant: it undoes rotate¶

Turning a vector into the world and back gives it again, within the roundings of the two turns.

In [6]:

there = trait_run.law('kinematics/rotate.trait', 'rotate', orientation=qs, value=vs)
back = trait_run.law('kinematics/unrotate.trait', 'unrotate', orientation=qs, value=there)
err = np.max(np.linalg.norm(back - vs, axis=1) / np.linalg.norm(vs, axis=1))
assert err <= 64 * U, err
print(f'there and back within {err / U:.1f} unit roundoffs')

there = trait_run.law('kinematics/rotate.trait', 'rotate', orientation=qs, value=vs) back = trait_run.law('kinematics/unrotate.trait', 'unrotate', orientation=qs, value=there) err = np.max(np.linalg.norm(back - vs, axis=1) / np.linalg.norm(vs, axis=1)) assert err <= 64 * U, err print(f'there and back within {err / U:.1f} unit roundoffs')

there and back within 7.6 unit roundoffs