Quaternions

Alex Diener — Jul 08, 2010

Quaternions

This tutorial is about quaternions - a way of representing rotations in three-dimensional space. You should be familiar with the Vector tutorial before reading this one.

Rotations

There are a number of different ways of expressing rotations in 3D space: Quaternions, Euler angles, 3x3 matrices, axis/angles, etc. They all have their own advantages and disadvantages. Euler angles suffer from a problem called gimbal lock, which can make certain rotations impossible. Matrices require more storage space in memory, and are less efficient for vector multiplication. Axis/angles work similarly to quaternions, but are less convenient for interpolation and vector multiplication than quaternions are. Quaternions and matrices tend to be the most convenient options in my experience. To learn more about matrices, I suggest you read the Matrix tutorial.

Euler Angles and Gimbal Lock

Euler angles consist of three rotation angles - one for each axis. The reason that they can be problematic is that these rotations are applied one after another, in a specific order (typically X, Y, Z), each one relative to the last. This makes some rotations impossible. Gimbal lock occurs when the second axis (Y, in this case) is rotated by 90° or -90°. When this happens, the third axis (Z) points in a parallel direction to the what the first one (X) was prior to the second rotation (Y), in essence overriding the first rotation and losing one axis of freedom.

Gimbal lock isn't a problem with quaternions, because they perform their rotation on all axes at once, rather than one at a time.

Axis/Angles

Before we continue, it's necessary that I talk about axis/angles briefly. An axis is a direction vector around which to rotate. The angle is the amount of rotation, typically expressed in radians, around that axis. In a right-handed coordinate system, a positive value for the angle rotates counter-clockwise around the axis. If you point the thumb of your right hand in the direction of the axis, a positive rotation angle rotates in the direction that your fingers curl, and a negative angle rotates the opposite direction. The code in this tutorial assumes a right-handed coordinate system.

What's a quaternion?

A quaternion is a 4-dimensional complex number commonly used to represent a rotation in 3-dimensional space. In this tutorial, I'll be using the following struct for quaternions: 
1
2
3
4
5
6
7
 
struct Quaternion {
  float x;
  float y;
  float z;
  float w;
};
typedef struct Quaternion Quaternion;
A quaternion is somewhat similar to an axis/angle. You can think of the x, y, and z components as the axis, and the w component as the angle. Internally, quaternions work very differently from axis/angles; this is simply a convenient analogy. Quaternion internals are complex and very difficult to understand, and will not be discussed in this tutorial. Fortunately, understanding exactly how quaternions do their magic is not a prerequisite to using them. All you need to know is what they can do, how to use them, and how you can implement them in your code.

How do you use quaternions?

A quaternion represents a rotational transformation. What this means is that a quaternion is a "recipe" for a rotation, so to speak. When applied to a vector, the vector is rotated on the axis of the quaternion, by the amount specified in the angle. The "neutral" quaternion that does no rotation at all is known as the identity quaternion. The value of the identity quaternion is {0, 0, 0, 1}.

A quaternion can be converted to and from an axis/angle without any loss of information. A quaternion can be converted to a rotation matrix, but a matrix cannot necessarily be converted back to a quaternion - this is only possible for certain matrices. The following sample code demonstrates these conversions. Note that this code uses the Vector struct and the Vector_normalize function from the Vector tutorial.

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
 
Quaternion Quaternion_fromAxisAngle(Vector axis, float angle) {
  Quaternion quaternion;
  float sinAngle;
  
  angle *= 0.5f;
  Vector_normalize(&axis);
  sinAngle = sin(angle);
  quaternion.x = (axis.x * sinAngle);
  quaternion.y = (axis.y * sinAngle);
  quaternion.z = (axis.z * sinAngle);
  quaternion.w = cos(angle);
  
  return quaternion;
}

void Quaternion_toAxisAngle(Quaternion quaternion, Vector * axis, float * angle) {
  float sinAngle;
  
  Quaternion_normalize(&quat);
  sinAngle = sqrt(1.0f - (quaternion.w * quaternion.w));
  if (fabs(sinAngle) < 0.0005f) sinAngle = 1.0f;
  
  if (axis != NULL) {
    axis->x = (quaternion.x / sinAngle);
    axis->y = (quaternion.y / sinAngle);
    axis->z = (quaternion.z / sinAngle);
  }
  
  if (angle != NULL) {
    *angle = (acos(quaternion.w) * 2.0f);
  }
}

Matrix Quaternion_toMatrix(Quaternion quaternion) {
  Matrix matrix;
  
  matrix.m[0]  = (1.0f - (2.0f * ((quaternion.y * quaternion.y) + (quaternion.z * quaternion.z))));
  matrix.m[1]  =         (2.0f * ((quaternion.x * quaternion.y) + (quaternion.z * quaternion.w)));
  matrix.m[2]  =         (2.0f * ((quaternion.x * quaternion.z) - (quaternion.y * quaternion.w)));
  matrix.m[3]  = 0.0f;
  matrix.m[4]  =         (2.0f * ((quaternion.x * quaternion.y) - (quaternion.z * quaternion.w)));
  matrix.m[5]  = (1.0f - (2.0f * ((quaternion.x * quaternion.x) + (quaternion.z * quaternion.z))));
  matrix.m[6]  =         (2.0f * ((quaternion.y * quaternion.z) + (quaternion.x * quaternion.w)));
  matrix.m[7]  = 0.0f;
  matrix.m[8]  =         (2.0f * ((quaternion.x * quaternion.z) + (quaternion.y * quaternion.w)));
  matrix.m[9]  =         (2.0f * ((quaternion.y * quaternion.z) - (quaternion.x * quaternion.w)));
  matrix.m[10] = (1.0f - (2.0f * ((quaternion.x * quaternion.x) + (quaternion.y * quaternion.y))));
  matrix.m[11] = 0.0f;
  matrix.m[12] = 0.0f;
  matrix.m[13] = 0.0f;
  matrix.m[14] = 0.0f;
  matrix.m[15] = 1.0f;
  
  return matrix;
}
If this looks like a bunch of gibberish, don't worry - as I said before, you don't need to understand the internals of quaternions in order to use them. Quaternion_fromAxisAngle takes an axis vector (normalized) and an angle (in radians), and returns a quaternion representing that rotation. Quaternion_toAxisAngle does the inverse. Quaternion_toMatrix fills a 4x4 matrix with the appropriate values to do the same transformation.

You can multiply two quaternions together as shown in the following code. The effect of this is to concatenate the two rotations together into a single quaternion. Quaternion multiplication is not commutative - the transformation of quaternion1 is applied before quaternion2, so to speak. I'm also including a convenience function you can use to rotate a quaternion by an axis/angle, as well as convenience functions for performing these transformations inline. Note that this code uses the Vector_cross and Vector_dot functions from the Vector tutorial.

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
 
void Quaternion_multiply(Quaternion * quaternion1, Quaternion quaternion2) {
  Vector vector1, vector2, cross;
  Quaternion result;
  float angle;
  
  vector1.x = quaternion1->x;
  vector1.y = quaternion1->y;
  vector1.z = quaternion1->z;
  vector2.x = quaternion2.x;
  vector2.y = quaternion2.y;
  vector2.z = quaternion2.z;
  angle = (quaternion1->w * quaternion2.w) - Vector_dot(vector1, vector2);
  
  cross = Vector_cross(vector1, vector2);
  vector1.x *= quaternion2.w;
  vector1.y *= quaternion2.w;
  vector1.z *= quaternion2.w;
  vector2.x *= quaternion1->w;
  vector2.y *= quaternion1->w;
  vector2.z *= quaternion1->w;
  
  quaternion1->x = vector1.x + vector2.x + cross.x;
  quaternion1->y = vector1.y + vector2.y + cross.y;
  quaternion1->z = vector1.z + vector2.z + cross.z;
  quaternion1->w = angle;
}

Quaternion Quaternion_multiplied(Quaternion quaternion1, Quaternion quaternion2) {
  Quaternion_multiply(&quaternion1, quaternion2);
  return quaternion1;
}

void Quaternion_rotate(Quaternion * quaternion, Vector axis, float angle) {
  Quaternion rotationQuaternion;
  
  rotationQuaternion = Quaternion_fromAxisAngle(axis, angle);
  Quaternion_multiply(quaternion, rotationQuaternion);
}

Quaternion Quaternion_rotated(Quaternion quaternion, Vector axis, float angle) {
  Quaternion_rotate(&quaternion, axis, angle);
  return quaternion;
}
Quaternion normalization works the same way as vector normalization, with an extra dimension: 
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
 
void Quaternion_normalize(Quaternion * quaternion) {
  float magnitude;
  
  magnitude = sqrt((quaternion->x * quaternion->x) +
                   (quaternion->y * quaternion->y) +
                   (quaternion->z * quaternion->z) +
                   (quaternion->w * quaternion->w));
  quaternion->x /= magnitude;
  quaternion->y /= magnitude;
  quaternion->z /= magnitude;
  quaternion->w /= magnitude;
}

Quaternion Quaternion_normalized(Quaternion quaternion) {
  Quaternion_normalize(&quaternion);
  return quaternion;
}
This function returns the inverse of a quaternion: The quaternion that performs the opposite transformation of the quaternion passed into the function. 
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
 
void Quaternion_invert(Quaternion * quaternion) {
  float length;
  
  length = 1.0f / ((quaternion->x * quaternion->x) +
                   (quaternion->y * quaternion->y) +
                   (quaternion->z * quaternion->z) +
                   (quaternion->w * quaternion->w));
  quaternion->x *= -length;
  quaternion->y *= -length;
  quaternion->z *= -length;
  quaternion->w *= length;
}

Quaternion Quaternion_inverted(Quaternion quaternion) {
  Quaternion_invert(&quaternion);
  return quaternion;
}
You can multiply a vector by a quaternion (without having to go to an intermediate matrix representation) as shown in the following code. For those of you who aren't yet fluent in this sort of terminology, "multiplying" a vector by a quaternion essentially means performing the quaternion's rotation on the vector. 
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
 
Vector Quaternion_multiplyVector(Quaternion quaternion, Vector vector) {
  Quaternion vectorQuaternion, inverseQuaternion, resultQuaternion;
  Vector resultVector;
  
  vectorQuaternion.x = vector.x;
  vectorQuaternion.y = vector.y;
  vectorQuaternion.z = vector.z;
  vectorQuaternion.w = 0.0f;
  
  inverseQuaternion = Quaternion_inverted(quaternion);
  resultQuaternion = Quaternion_multiplied(vectorQuaternion, inverseQuaternion);
  resultQuaternion = Quaternion_multiplied(quaternion, resultQuaternion);
  
  resultVector.x = resultQuaternion.x;
  resultVector.y = resultQuaternion.y;
  resultVector.z = resultQuaternion.z;
  
  return resultVector;
}
And finally, there's SLERP. SLERP stands for spherical linear interpolation. As the name suggests, it allows you to interpolate between two quaternions. The alpha parameter specifies the point of interpolation: 0 for start, 1 for end, 0.5 for the midpoint between start and end, etc. So, if start is a quaternion that rotates 90° on the X axis, end is a quaternion that rotates 180° on the X axis, and alpha is 0.5, then the function will return a quaternion that rotates 135° on the X axis. 
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
 
#define SLERP_TO_LERP_SWITCH_THRESHOLD 0.01f

Quaternion Quaternion_slerp(Quaternion start, Quaternion end, float alpha) {
  float startWeight, endWeight, difference;
  Quaternion result;
  
  difference = ((start.x * end.x) + (start.y * end.y) + (start.z * end.z) + (start.w * end.w));
  if ((1.0f - fabs(difference)) > SLERP_TO_LERP_SWITCH_THRESHOLD) {
    float theta, oneOverSinTheta;
    
    theta = acos(fabs(difference));
    oneOverSinTheta = (1.0f / sin(theta));
    startWeight = (sin(theta * (1.0f - alpha)) * oneOverSinTheta);
    endWeight = (sin(theta * alpha) * oneOverSinTheta);
    if (difference < 0.0f) {
      startWeight = -startWeight;
    }
  } else {
    startWeight = (1.0f - alpha);
    endWeight = alpha;
  }
  result.x = (start.x * startWeight) + (end.x * endWeight);
  result.y = (start.y * startWeight) + (end.y * endWeight);
  result.z = (start.z * startWeight) + (end.z * endWeight);
  result.w = (start.w * startWeight) + (end.w * endWeight);
  Quaternion_normalize(&result);
  
  return result;
}
Those are the basics in a nutshell. There are some advanced topics I didn't cover, such as interpolation between more than two quaternions at once, but the basics go a long way. Contact me if you have any questions, or ask on the message board. Also, be sure to check out the Quaternion playground companion application at the bottom of this page.

Code Implementation

Quaternion.h : Download (2.3 KB)

Quaternion.c : Download (7.2 KB)

Dependencies:

Matrix.h : Download (2.8 KB)

Matrix.c : Download (9.8 KB)

Vector.h : Download (1.5 KB)

Vector.c : Download (2.4 KB)

Companion application:

uDevGames 2011

Convergence — Best Gameplay
Kung Fu Killforce — Best Overall Game, Best Audio, Best Presentation
Flying Sweeden — Best Graphics, Most Original
Time Goat — Best Story

Site


Submit News, Article

Contact Us

Forum Threads

Visit the Forums