-
Notifications
You must be signed in to change notification settings - Fork 0
quaternions
Quaternions happen to be a good way of representating rotations (and orientations).
Quaternions combine real and imaginary numbers, rather like complex numbers, except that there are three distinct imaginaries involved. In math, the three imaginary parts are usually called i, j, and k, but in programming they are more often denoted with x, y, and z. The non-imaginary, real, component, is denoted 'w'.
Despite using the names x, y, and z, these components do not exactly represent Cartesian axes; that is, quaternions are not the same as angle-axis representations. In fact, quaternions are quite weird to understand, not at all intuitive, so why do we use them?
- Compared to Euler angles, they avoid important ambiguities. In particular, they avoid "Gimbal lock". Consider trying to point North while standing at the North pole; it doesn't matter which direction you stand, they are all pointing north. A quaternion can distinctly represent each of these directions, while Euler angles cannot.
- Compared to using rotation matrices, they are smaller (four rather than nine components), reducing memory, bandwidth, and computation costs. Often the mathematics becomes simpler (e.g. inversion).
- They also permit more accurately smoothed rotations, via SLerp.
Conceptually quaternions are more similar to matrices than to vectors: - they can encode a transformation from one coordinate (orientation) frame to another - they can be inverted to define the reverse transformation - they can be applied to other quaternions or to vectors - you apply them by multiplication - the order of operation is important
Of course quaternions can be converted to matrices and vice versa, and quaternions can be constructed from Euler angles and other rotation representations.
Many libraries provide Quaternion types, however not all store the quaternion in the same order. Most use XYZW ordering, but some use WXYZ. Some access components by array notation (q[0]), others by name (q.x); some allow both.
Usually you can simply write conversion functions between these library implementations.
One of the most common tasks: take a given vector and rotate it by a quaternion, which is to say, transform it into the orientation frame of the quaternion. Most libraries provide a function for this operation. E.g. in gl-matrix, the function vec3.transformQuat(v, v, q). However it is also rather simple to implement.
In order to let quaternions and vec3s inter-operate, we turn the vector into a "pure" quaternion, in which the vector gains a .w component of zero: p = vec4(v.x, v.y, v.z, 0.0). Then the transformation is defined as: ((q * p) * 1/q).xyz.
Here 1/q refers to the conjugate of the quaterion. If we can assume that q is a unit quaternion, then the conjugate is obtained by simply negating its x, y, z components.
Note that the inverse transformation (rotating a vector "out" of a particular orientation, or applying the reverse orientation), is thus simply: ((1/q * p) * q).xyz.
// q must be a normalized quaternion
// v is a vec3
// out is the resulting vec3
function quat_rotate(out, q, v) {
// first compute p = q * vec4(v, 0):
const px = q[3]*v[0] + q[1]*v[2] - q[2]*v[1];
const py = q[3]*v[1] + q[2]*v[0] - q[0]*v[2];
const pz = q[3]*v[2] + q[0]*v[1] - q[1]*v[0];
const pw = -q[0]*v[0] - q[1]*v[1] - q[2]*v[2];
// return p * 1/q:
out[0] = px*q[3] - pw*q[0] + pz*q[1] - py*q[2];
out[1] = py*q[3] - pw*q[1] + px*q[2] - pz*q[0];
out[2] = pz*q[3] - pw*q[2] + py*q[0] - px*q[1];
return out;
}
// the inverse rotation:
function quat_unrotate(out, q, v) {
// first compute p = q * vec4(v, 0):
const px = q[3]*v[0] + q[1]*v[2] - q[2]*v[1];
const py = q[3]*v[1] + q[2]*v[0] - q[0]*v[2];
const pz = q[3]*v[2] + q[0]*v[1] - q[1]*v[0];
const pw = -q[0]*v[0] - q[1]*v[1] - q[2]*v[2];
// return p * 1/q:
out[0] = px*q[3] - pw*q[0] + pz*q[1] - py*q[2];
out[1] = py*q[3] - pw*q[1] + px*q[2] - pz*q[0];
out[2] = pz*q[3] - pw*q[2] + py*q[0] - px*q[1];
return out;
}// q must be a normalized quaternion
template<typename T, glm::precision P>
glm::tvec3<T, P> quat_rotate(glm::quat const & q, glm::tvec3<T, P> & v) {
glm::tvec4<T, P> p(
q.w*v.x + q.y*v.z - q.z*v.y, // x
q.w*v.y + q.z*v.x - q.x*v.z, // y
q.w*v.z + q.x*v.y - q.y*v.x, // z
-q.x*v.x - q.y*v.y - q.z*v.z // w
);
return glm::tvec3<T, P>(
p.x*q.w - p.w*q.x + p.z*q.y - p.y*q.z, // x
p.y*q.w - p.w*q.y + p.x*q.z - p.z*q.x, // y
p.z*q.w - p.w*q.z + p.y*q.x - p.x*q.y // z
);
}As with matrices, when using quaternions to represent orientations only, it is essential to maintain them as unit quaternions; that is, their magnitude should equal 1. If this is not maintained they may have unexpected and confusing effects beyond rotation.
Fortunately this normalization is easy to compute: simply compute the magnitude (the square root of the sum of the squares of the components), then divide each component by the magnitude. You may want to perform normalization periodically, perhaps once per frame.
function quat_normalize(out, q) {
const mag = Math.sqrt(q[0]*q[0] + q[1]*q[1] + q[2]*q[2] + q[3]*q[3]);
const div = 1/mag;
out[0] = q[0] * div;
out[1] = q[1] * div;
out[2] = q[2] * div;
out[3] = q[3] * div;
return out;
}Quite often we want a unit vector in the "up", "right" etc. directions of an object's orientation. To do so, simply let the quaternion 'q' transform a unit axis vector:
// get a unit X axis vector from quaternion 'q':
vec3.transformQuat(axis_x, vec3.fromValues(1, 0, 0), q);
// or, using the code above,
quat_rotate(axis_x, vec3.fromValues(1, 0, 0), q);However since the vec3 values are all zeroes or ones, the math can be significantly reduced. And since these operations are commonly needed, it is worth implementing them as functions:
function quat_ux(out, q) {
out[0] = 1 - 2 * ((q[1] * q[1]) + (q[2] * q[2]));
out[0] = 2 * ((q[0] * q[1]) + (q[3] * q[2]));
out[0] = 2 * ((q[0] * q[2]) - (q[3] * q[1]));
return out;
}
function quat_uy(out, q) {
out[0] = 2 * ((q[0] * q.y) - (q.w * q.z));
out[0] = 1 - 2 * ((q[0] * q[0]) + (q.z * q.z));
out[0] = 2 * ((q.y * q.z) + (q.w * q[0]));
return out;
}
function quat_uz(out, q) {
out[0] = 2 * ((q[0] * q.z) + (q.w * q.y));
out[0] = 2 * ((q.y * q.z) - (q.w * q[0]));
out[0] = 1 - 2 * ((q[0] * q[0]) + (q.y * q.y));
return out;
}In C++ (using the GLM library):
template<typename T, glm::precision P>
inline glm::tvec3<T, P> quat_ux(glm::tquat<T, P> const & q) {
return glm::tvec3<T, P>(
T(1) - T(2) * ((q.y * q.y) + (q.z * q.z)),
T(2) * ((q.x * q.y) + (q.w * q.z)),
T(2) * ((q.x * q.z) - (q.w * q.y))
);
}
template<typename T, glm::precision P>
inline glm::tvec3<T, P> quat_uy(glm::tquat<T, P> const & q) {
return glm::tvec3<T, P>(
T(2) * ((q.x * q.y) - (q.w * q.z)),
T(1) - T(2) * ((q.x * q.x) + (q.z * q.z)),
T(2) * ((q.y * q.z) + (q.w * q.x))
);
}
template<typename T, glm::precision P>
inline glm::tvec3<T, P> quat_uz(glm::tquat<T, P> const & q) {
return glm::tvec3<T, P>(
T(2) * ((q.x * q.z) + (q.w * q.y)),
T(2) * ((q.y * q.z) - (q.w * q.x)),
T(1) - T(2) * ((q.x * q.x) + (q.y * q.y))
);
}One way to describe a quaternion is as the rotation from one direction vector to another. For example, what is the rotation I need to apply to change my current heading to point toward my home?
// create unit quat as the rotation from vec3 u to vec3 v
// and store in 'out'
// u and v need not be unit vectors
function rotationBetweenVectors(out, u, v) {
let uv = vec3.dot(u, v);
let uu = vec3.dot(u, u);
let vv = vec3.dot(v, v);
let norm = Math.sqrt(uu * vv);
// real & imaginary components of result:
let real = norm + uv;
let imag = vec3.create();
// handle special case:
if (real < 1.e-6 * norm) {
// bases are opposing, so no single perpendicular basis is preferable
// (cross product of u and v is ill-defined)
// just make up any orthogonal axis to u:
real = 0.;
if (Math.abs(u[0] > Math.abs(u[2]) {
vec3.set(imag, -u[1], u[0], 0);
} else {
vec3.set(imag, 0, -u[2], u[1]);
}
} else {
// axis of rotation is simply cross product of u,v:
vec3.cross(imag, u, v);
}
return quat.normalize(out, quat.set(out, imag.x, imag.y, imag.z, real));
}We can perform linear interpolation on quaternions (by linearly-interpolating each component), however it does not usually give what we want. One way to imagine what a quaternion represents is to think of it as a person standing on a particular point the surface of the Earth, facing a particular direction. Given two such quaternions, the linear interpolation between them draws a direct line through the Earth's core. Not only does this change the distance from the Earth's centre during the interpolation, it also changes the speed with respect to the surface.
In contrast SLerp, which stands for "Spherical Linear Interpolation", traces the shortest path between the two quaternions along the surface of the Earth, which not only maintains constant distance from the centre, but also a constant "speed" (a constant ratio between the interpolation factor and the distance traveled along the surface). When used for animating objects, cameras, etc. it leads to much smoother and more graceful behaviour.
The average of two quaternions can be done simply using SLerp. For more than two quaternions, here's a solution:
// note: this assumes all the quaterions in list are unit quaternions,
// i.e. they are already normalized.
let qsum = quat.create();
for (let q in list) {
// there is an ambiguity in quaternions that q and -q represent the same orientation.
// consider similarity to running total and negate if too different
if (quat.dot(q, qsum) < 0) {
quat.scale(q, q, -1);
}
// now add each x, y, z, w component
quat.add(qsum, qsum, q);
}
// since we want a unit vector output anyway,
// just normalize to get average:
return quat.normalize(qsum, qsum);