Quaternions
A
quaternion is a four dimensional vector that can be used to represent rotations in three dimensions.
Compared to standard vector algebra, quaternions are less intuitive and much less familiar. However,
rotations represented by quaternions are more computationally efficient and stable, and have therefore become
popular in 3D applications such as Houdini.
Quaternions were discovered by William Rowen Hamilton in 1843, an Irish mathematician whose name is
immediately familiar to all physicists. By his own account, Hamilton was out for a walk with his wife in Dublin,
then when crossing a bridge, was struck with an idea that would lead him to develop the theory of quaternions.
He took out his pen-knife to scratch his inspiration into the stone of the bridge:
\( \mathbf{i}^2 = \mathbf{j}^2 = \mathbf{k}^2 = \mathbf{ijk} = -1 \).
His engraving has since disappeared from the Brougham Bridge, but a plaque commemorating the event can now
be found there, and the theory it inspired lives on.
Interesting aside - research into the writings of Lewis Carroll, who was also a mathematician,
has apparently established that
Alice in Wonderland was in large part a
satire of the math being developed in his time. Hamilton was a contemporary of Carroll, and the former's now
fashionable quaternions (one imagines only fashionable amongst mathematicians) were the focus of the story of the Mad Hatter.
The famous tea party was a critique of Hamilton, meant to illustrate the absurdity of his "new" math.
There is a great discussion of quaternions at
euclideanspace.com,
and as always, a detailed reference at
wikipedia.
What exactly is a quaternion
A quaternion is a
hypercomplex number
with four real coefficients. A hypercomplex number is like a complex number, only more so. If
\( a + b\mathbf{i} \) is a normal complex number, then a quaternion can be written
\( \mathbf{q} = a + b\mathbf{i} + c\mathbf{j} + d\mathbf{k} \)
where \( a, b, c, d \) are real numbers, and
\( \mathbf{i}, \mathbf{j}, \mathbf{k} \) are different imaginary number known
as quaternion units.
Axis-Angle
When quaternions are used for rotations, they are defined in an axis-angle
representation. That is, any rotation can be represented by a vector (an axis of rotation) and a rotation
around that vector,
\( \mathbf{q} = \cos(\alpha / 2) + \sin(\alpha / 2) ( \mathbf{i} x + \mathbf{j} y + \mathbf{k} z )\)
where \( \alpha \) is the rotation angle and \( [x,y,z] \) is a unit vector that gives the rotation axis.
Quaternion rotation
Given a vector \( \mathbf{r} \) in three dimensions and a normalised quaternion \( \mathbf{q} \) in the axis-angle representation,
the rotation of \( \mathbf{r} \) by \( \mathbf{q} \) is achieved by the conjugate product,
\(\mathbf{r}^{'} = \mathbf{q}\mathbf{rq}^{-1}\)
where \( \mathbf{r}^{'} \) is the rotated vector. This, the central result of this section, is stated here without
proof, though one could imagine deriving a matrix-vector representation of the above conjugate product, and then
showing the equivalency with the three dimensional rotation matrix.
To evaluate the conjugate product, we must establish the inverse of a quaternion, and define quaternion muliplication.
Quaternion multiplication
The multiplication of two quaternions is called a Hamilton product. To calculate the Hamilton product, we'll need
distributivity, and we'll need to know how the basis elements multiply.
From the Brougham Bridge equations, and the fact that quaternion units are imaginary,
\(\mathbf{i}^{2} = \mathbf{j}^{2} = \mathbf{k}^{2} = -1\).
The "off-diagonal" multiplications can be derived from \( \mathbf{ijk} = -1 \),
\(\mathbf{ij} = \mathbf{k}; \: \mathbf{jk} = \mathbf{i}; \: \mathbf{ki} = \mathbf{j} \)
\( \mathbf{ji} = -\mathbf{k}; \: \mathbf{kj} = -\mathbf{i}; \: \mathbf{ik} = -\mathbf{j} \)
So we see that the rules for multiplying the quaternion units are encoded into the equations Hamilton etched into the bridge.
For two quaternions, \( \mathbf{q}_{1} = a + b\mathbf{i} + c\mathbf{j} + d\mathbf{k} \) and
\( \mathbf{q}_{2} = e + f\mathbf{i} + g\mathbf{j} + h\mathbf{k} \), using distributivity to multiply through the brackets
and the rules for multiplying basis elements, we find
\(\mathbf{q}_{1}\mathbf{q}_{2} = (ae-bf-cg-dh) + (af+be+ch-dg)\mathbf{i} + (ag+ce-bh+df)\mathbf{j} + (ak+de+bg-cf)\mathbf{k}\).
Norm and inverse of a quaternion
If \( \mathbf{q} = a + b\mathbf{i} + c\mathbf{j} + d\mathbf{k} \), the complex conjugate \( \mathbf{q}^{*} \) is
\(\mathbf{q}^{*} = a - b\mathbf{i} - c\mathbf{j} - d\mathbf{k}\)
Using the rules of quaternion multiplication, we find the product \( \mathbf{qq}^{*} \) to be a non-negative real number,
\(\mathbf{qq}^{*} = a^{2}+ b^{2} + c^{2} + d^{2}\)
So we can define the norm of a quaternion \(\mathbf{q}\) to be
\(||\mathbf{q}|| = (\mathbf{qq}^{*})^{1/2}\)
It follows that dividing the quaternion \(\mathbf{q}\) by the norm \(||\mathbf{q}||\) produces a unit quaternion whose norm,
given by the above equation, is one. And since the real number one can be considered a quaternion with real element one
with imaginary elements equal zero, multiplication with any other quaternion shows one to be the identity quaternion
under multiplcation. That is,
\(\mathbf{Iq} = \mathbf{qI} = \mathbf{q}\)
where \(\mathbf{I}\) represents the identity quaternion. It then follows that the multiplicative inverse, \(\mathbf{q}^{-1}\),
\(\mathbf{qq}^{-1} = \mathbf{I} = \mathbf{qq}^{*} / ||\mathbf{q}||^{2}\)
showing that \( \mathbf{q}^{*}/||\mathbf{q}||^{2}\) is the inverse of \(\mathbf{q}\), which is simply \(\mathbf{q}^{*}\) if
\(\mathbf{q}\) is a unit quaternion.
Combining rotations
Suppose we have a vector \( \mathbf{p} \), rotated by the quaternion \( \mathbf{q}_{1} \), followed by a rotation by \( \mathbf{q}_{1} \).
Using our formula for quaternion rotation, we find for the resulting rotated vector \( \mathbf{p}^{'} \)
\( \mathbf{p}^{'} = \mathbf{q}_{2}(\mathbf{q}_{1}\mathbf{p}\mathbf{q}_{1}^{*})\mathbf{q}_{2}^{*} \)
\( = (\mathbf{q}_{2}\mathbf{q}_{1})\mathbf{p}(\mathbf{q}_{1}^{*}\mathbf{q}_{2}^{*}) \)
\( = \mathbf{q}_{3}\mathbf{p}\mathbf{q}_{3}^{* } \)
where \( \mathbf{q}_{3} = \mathbf{q}_{2} \mathbf{q}_{1} \) is the combined rotation, showing that multiplying quaternions combines their rotations.
Simple example, rotation in xz plane
Although we won't prove the rotation formula here, a simple demonstration will clarify how it works, and illustrate how
practical quaternions can be for calculating rotations. As an example, we will rotate the vector i
by θ around the axis j; that is, a rotation in the xz plane.
\( \mathbf{r} = \mathbf{i} \)
\(\mathbf{q} = \cos(\theta/2) + \mathbf{j} \sin(\theta/2) \)
\( \mathbf{q}^{-1} = \mathbf{q}^{*} = \cos(\theta/2) - \mathbf{j} \sin(\theta/2) \)
\( \mathbf{rq}^{-1}= \mathbf{i} \cos(\theta/2) - \mathbf{k} \sin(\theta/2) \)
\( \mathbf{qrq}^{-1} = \mathbf{i} \cos^{2}(\theta/2) - \mathbf{k} \cos(\theta/2) \sin(\theta/2) + \mathbf{ji} \cos(\theta/2) \sin(\theta/2) - \mathbf{jk} \sin^{2}(\theta/2) \)
\(= \mathbf{i}( \cos^{2}(\theta/2) - \sin^{2}(\theta/2) ) - 2 \mathbf{k} \sin(\theta/2) \cos(\theta/2) \)
Recall the half-angle identities:
\( \sin(2\theta) = 2 \sin(\theta) \cos(\theta) \) and \( \cos(2\theta) = \cos^{2}(\theta) - \sin^{2}(\theta) \)
then the conjugate product above becomes
\(\mathbf{r}^{'} = \mathbf{i} \cos(\theta) - \mathbf{k} \sin(\theta) \),
which is evidently a rotation around the j axis, as required.
Houdini and Quaternions
In this section, I've assembled functions and VEX that I use regularly in Houdini.
A great reference on this topic is the always excellent
cgwiki,
where you'll find a lots of examples and tricks for using @orient.
In Houdini, @orient is a vector4 attribute reserved for quaternions. Note that @orient is defined
backward compared to this discussion and what seems to be standard in the math literature. That is,
the scalar part is the last component.
@orient \(= [b,c,d,a] = a + b \mathbf{i} + c \mathbf{j} + d \mathbf{k} \)
Useful VEX for quaternion rotations
quaternion returns a vector4 quaternion given a rotation axis and an angle (in radians).
vector4 quaternion( float, vector );
qrotate returns a vector v rotated by a quaternion.
vector qrotate( vector4, vector );
qmultiply returns the product of two quaternions which, as we've seen, combines their rotations.
vector qmultiply( vector4, vector4 );
VEX for rotating instances
#include <math.h>
float w = 0.01;
int pid = @ptnum;
float tpi = 2.0 * PI;
vector axis;
axis[0] = fit01( rand(0.3*pid+0.3), -1, 1);
axis[1] = fit01( rand(0.3*pid+0.5), -1, 1);
axis[2] = fit01( rand(0.3*pid+0.7), -1, 1);
axis = normalize(axis);
float angle = tpi * rand(0.5*pid);
float dangle = w * @Frame * sqrt(length(@v));
dangle = tpi * frac(dangle/tpi);
angle += dangle;
@orient = quaternion(angle,axis);
To rotate one vector into another, run the following code in 'detail' mode. The cross product creates a vector
orthogonal to both, which is the axis of rotation. The dot product gives the cosine of the angle between the two, so
acos yields the positive rotation angle (from 0 to pi).
vector p0 = point(0,'P',0);
vector p1 = point(1,'P',0);
p0 = normalize(p0);
p1 = normalize(p1);
vector axis = normalize(cross(p1,p0));
float angle = acos(dot(p0,p1));
p@q = quaternion(angle,axis);
This create a quaternion detail attribute 'q', which can be used in a wrangle to rotate p1 to p0.
Some notes on quaternions in VEX
-
Direction of rotation for qrotate(angle,axis) - a rotation of angle around axis corresponds
to the right hand rule
-
If the axis is not normalised - all bets are off. Seems to give some combination
of rotate and scale, but is not intuitive
-
Using @opinput with quaternions - be sure to use p@opinput. Houdini does not implicitly typecast, even if
you're assigning to a vector4. Generally, I find VEX to be a bit sketchy when it comes to implicit typecasting,
so if you find yourself with a weird bug, this is one place to look.