by James McLellan
Love them? Hate them? Quaternions are very useful in graphics and physical simulation applications for their ability to accurately express rotations about 3-axes. OpenGL, OpenSceneGraph, Bullet and many other libraries natively support some kind of quaternion math, so that the average programmer does not have to spend one moment longer than necessary thinking about them.
But if you're caught in an application where a library with a handy quaternion isn't available, you might be stuck. There are many web sites that tell you how to do quaternion math. Some get it right. Some don't.
So, in case you are ever that person that needs a 3D rotation right now without having a handy library to fall back on, here's how to code a quaternion.
This Excel file contains a simple construction of a quaternion matrix that you can quickly use to test.
Download:
quaternions.xls
To make your own quaternion transformation:
w = cos(phi) * cos(theta) * cos(psi) + sin(phi) * sin(theta) * sin(psi)
If you are visually identifying good quaternion formations, "w" should equal 1.0 for no-rotation conditions, and should equal the cosine of the angle for single-axis rotations. For example:
<x, y, z, w> = <0.0, 0.0, 0.0, 1.0> This is a no-rotation condition
<x, y, z, w> = <0.5, 0.0, 0.0, 0.86> This is a 30 degree rotation along the x-axis. The "w" value would be the same if it was along the y-axis or z-axis instead.
w = cos(phi) * cos(theta) * cos(psi) + sin(phi) * sin(theta) * sin(psi);
x = sin(phi) * cos(theta) * cos(psi) - cos(phi) * sin(theta) * sin(psi);
y = cos(phi) * sin(theta) * cos(psi) + sin(phi) * cos(theta) * sin(psi);
z = cos(phi) * cos(theta) * sin(psi) - sin(phi) * sin(theta) * cos(psi);
To use your quaternion, you'll construct a 3x3 matrix that you can multiply with a point to rotate it.
q[0][0] = 1-2 * (y * y + z * z);
q[0][1] = 2 * (x * y - w * z);
q[0][2] = 2 * (w * y + x * z);
q[1][0] = 2 * (x * y + w * z);
q[1][1] = 1 - 2 * (x * x + z * z);
q[1][2] = 2 * (y * z - w * x);
q[2][0] = 2 * (x * z - w * y);
q[2][1] = 2 * (w * x + y * z);
q[2][2] = 1 - 2 * (x * x + y * y);
So, if you wanted to rotate a wind vector that was coming in at x = -20.0, y = 10.0, z = 5.0 it would look like:
a[0] = -20.0;
a[1] = 10.0;
a[2] = 5.0;
r[0] = q[0][0] * a[0] + q[1][0] * a[0] + q[2][0] * a[0]; // the rotated x-value
r[1] = q[0][1] * a[1] + q[1][1] * a[1] + q[2][1] * a[1]; // the rotated y-value
r[2] = q[0][2] * a[2] + q[1][2] * a[2] + q[2][2] * a[2]; // the rotated z-value
If theta=30, phi=0, psi=0; the result would be r[0] = result.x = -14.8, r[1] = result.y = 10, r[2] = result.z =+14.33.
There are also instances where you want to stack rotations. For example, say you are rotating a child object with a parent, grandparent and you want to include those rotations.
Although addition might seem like the first solution, it's not. What you want, instead, is to multiply the two quaternions. A good walkthrough of the derivation is available on Wikipedia. I'll skip straight to the answer:
a1 = quaternion1.w;
a2 = quaternion2.w;
b1 = quaternion1.x;
b2 = quaternion2.x;
c1 = quaternion1.y;
c2 = quaternion2.y;
d1 = quaternion1.z;
d2 = quaternion2.z;
result.w = a1*a2 - b1*b2 - c1*c2 - d1*d2;
result.x = a1*b2 + a2*b1 - c1*d2 - d2*c1;
result.y = a1*c2 + a2*c1 + b1*d2 - d1*b2;
result.z = a1*d2 + d2*d1 - b1*c2 + b2*c1;
You can also test it out using the Excel spreadsheet attached. Try stacking a quaternion rotation theta=45 (other angles zero) to a quaternion of psi=90 (other angles zero). It should equal a single quaternion rotation of theta=45 and psi=90.