OpenGL:Tutorials:Using Quaternions to represent rotation

From GPWiki

Foreword and warning

Quaternions are all the rage these days for (3D) computer games, so this wiki wouldn't be complete without an explanation about them. Unfortunately, I'm not exactly a quaternion-specialist, so there might be errors here. I hope someone with more knowledge on the topic will review this article. Although this article is in the OpenGL-section, the background information is of course true for Direct3D too. As far as I know, D3D also has some convenience functions for Quaternions.

Just what is a quaternion?

A quaternion is an element of a 4 dimensional vector-space. It's defined as w + xi + yj + zk where i, j and k are imaginary numbers. Alternatively, a quaternion is what you get when you add a scalar and a 3d vector. The math behind quaternions is only slightly harder than the math behind vectors, but I'm going to spare you (for the moment).

Sounds scary? Ok, so now you know never to ask that question again...

Fortunately for you, we will only work with a subset of quaternions: unit quaternions. Those are quaternions with length 1. Every unit quaternion represents a 3D rotation, and every 3D rotation has two unit quaternion representations. Unit quaternions are a way to compactly represent 3D rotations while avoiding singularities or discontinuities (e.g. gimbal lock).

Rotation	Quaternions
Identity (no rotation)	1, -1
180 degrees about x-axis	i, -i
180 degrees about y-axis	j, -j
180 degrees about z-axis	k, -k
angle θ, axis (unit vector)

Why use quaternions

Quaternions have some advantages over other representations of rotations.

• Quaternions don't suffer from gimbal lock, unlike Euler angles.
• They can be represented as 4 numbers, in contrast to the 9 numbers of a rotations matrix.
• The conversion to and from axis/angle representation is trivial.
• Smooth interpolation between two quaternions is easy (in contrast to axis/angle or rotation matrices).
• After a lot of calculations on quaternions and matrices, rounding errors accumulate, so you have to normalize quaternions and orthogonalize a rotation matrix, but normalizing a quaternion is a lot less troublesome than orthogonalizing a matrix.
• Similar to rotation matrices, you can just multiply 2 quaternions together to receive a quaternion that represents both rotations.

The only disadvantages of quaternions are:

• They are hard to visualize.
• You have to convert them to get a human-readable representation (Euler angles) or something OpenGL can understand (Matrix).
• Smooth interpolation between quaternions is complicated by the fact that each 3D rotation has two representations.

Why quaternions are neat

Quaternions are neat because the unit quaternions are a double-cover for the set of 3D rotations.

To understand what that means, consider a Mobius strip. To make one, start with a strip of paper and bring the ends together. Give one of the ends a half-twist, and then attach it to the other end. The resulting surface is a Mobius strip.

The Mobius strip has a very important feature: it's non-orientable. When I want to build and render an object, like a sphere or a cylinder, I need to make sure all my vertex normals point the same way (e.g. outward). If I try to render a Mobius strip, I'll get stuck because I won't be able to assign consistent outward-pointing vector normals. That's what it means to be non-orientable.

The solution to this problem is to simply duplicate each vertex. One copy gets an "outward" pointing normal, and the other copy gets an "inward" pointing normal. When I render the Mobius strip, I'll have to draw each polygon twice, once facing "outward", and again facing "inward". Mathematicians would call this a double-cover.

The set of 3D rotations can be thought of as a 3D surface sitting in 4D hyper-space. Like the Mobius strip, this surface is non-orientable. As a result, if I try to represent a 3D rotation using three numbers (e.g. Euler angles), I'll get stuck with either a singularity or a discontinuity.

Just like with the Mobius strip, I can solve the problem with a double-cover. I simply duplicate each possible 3D rotation; one copy gets a "positive" quaternion, and the other copy gets a "negative" quaternion. The resulting representation is equivalent to a 3-sphere in 4D hyper-space, it's orientable, and it completely avoids the problem of singularities and discontinuities.

Some basic quaternion operations

Here are some methods you will regularly need to work with quaternions.

Normalizing a quaternion

// normalising a quaternion works similar to a vector. This method will not do anything
// if the quaternion is close enough to being unit-length. define TOLERANCE as something
// small like 0.00001f to get accurate results
void Quaternion::normalise()
{
	// Don't normalize if we don't have to
	float mag2 = w * w + x * x + y * y + z * z;
	if (fabs(mag2 - 1.0f) > TOLERANCE) {
		float mag = sqrt(mag2);
		w /= mag;
		x /= mag;
		y /= mag;
		z /= mag;
	}
}

The complex conjugate of a quaternion

// We need to get the inverse of a quaternion to properly apply a quaternion-rotation to a vector
// The conjugate of a quaternion is the same as the inverse, as long as the quaternion is unit-length
Quaternion Quaternion::getConjugate()
{
	return Quaternion(-x, -y, -z, w);
}

Multiplying quaternions

To multiply two quaternions, write each one as the sum of a scalar and a vector. The product of and is where

// Multiplying q1 with q2 applies the rotation q2 to q1
Quaternion Quaternion::operator* (const Quaternion &rq) const
{
	// the constructor takes its arguments as (x, y, z, w)
	return Quaternion(w * rq.x + x * rq.w + y * rq.z - z * rq.y,
	                  w * rq.y + y * rq.w + z * rq.x - x * rq.z,
	                  w * rq.z + z * rq.w + x * rq.y - y * rq.x,
	                  w * rq.w - x * rq.x - y * rq.y - z * rq.z);
}

Please note: Quaternion-multiplication is NOT commutative. Thus q1 * q2 is not the same as q2 * q1. This is pretty obvious actually: As I explained, quaternions represent rotations and multiplying them "concatenates" the rotations. Now take you hand and hold it parallel to the floor so your hand points away from you. Rotate it 90° around the x-axis so it is pointing upward. Now rotate it 90° clockwise around its local y-axis (the one coming out of the back of your hand). Your hand should now be pointing to your right, with you looking at the back of your hand. Now invert the rotations: Rotate your hand around the y-axis so its facing right with the back of the hand facing upwards. Now rotate around the x axis and your hand is pointing up, back of hand facing your left. See, the order in which you apply rotations matters. Ok, ok, you probably knew that...

Rotating vectors

To apply a quaternion-rotation to a vector, you need to multiply the vector by the quaternion and its conjugate.

// Multiplying a quaternion q with a vector v applies the q-rotation to v
Vector3 Quaternion::operator* (const Vector3 &vec) const
{
	Vector3 vn(vec);
	vn.normalise();
 
	Quaternion vecQuat, resQuat;
	vecQuat.x = vn.x;
	vecQuat.y = vn.y;
	vecQuat.z = vn.z;
	vecQuat.w = 0.0f;
 
	resQuat = vecQuat * getConjugate();
	resQuat = *this * resQuat;
 
	return (Vector3(resQuat.x, resQuat.y, resQuat.z));
}

How to convert to/from quaternions

In the following, I will present the methods necessary to convert all kind of rotation-representations to and from quaternions. I'll not show how to derive them because, well, who cares? (oh, and because I don't know how)

Quaternion from axis-angle

To rotate through an angle θ, about the axis (unit vector) , use:

// Convert from Axis Angle
void Quaternion::FromAxis(const Vector3 &v, float angle)
{
	float sinAngle;
	angle *= 0.5f;
	Vector3 vn(v);
	vn.normalise();
 
	sinAngle = sin(angle);
 
	x = (vn.x * sinAngle);
	y = (vn.y * sinAngle);
	z = (vn.z * sinAngle);
	w = cos(angle);
}

Quaternion from Euler angles

// Convert from Euler Angles
void Quaternion::FromEuler(float pitch, float yaw, float roll)
{
	// Basically we create 3 Quaternions, one for pitch, one for yaw, one for roll
	// and multiply those together.
	// the calculation below does the same, just shorter
 
	float p = pitch * PIOVER180 / 2.0;
	float y = yaw * PIOVER180 / 2.0;
	float r = roll * PIOVER180 / 2.0;
 
	float sinp = sin(p);
	float siny = sin(y);
	float sinr = sin(r);
	float cosp = cos(p);
	float cosy = cos(y);
	float cosr = cos(r);
 
	x = sinr * cosp * cosy - cosr * sinp * siny;
	y = cosr * sinp * cosy + sinr * cosp * siny;
	z = cosr * cosp * siny - sinr * sinp * cosy;
	w = cosr * cosp * cosy + sinr * sinp * siny;
 
	normalise();
}

Quaternion to Matrix

// Convert to Matrix
Matrix4 Quaternion::getMatrix() const
{
	float x2 = x * x;
	float y2 = y * y;
	float z2 = z * z;
	float xy = x * y;
	float xz = x * z;
	float yz = y * z;
	float wx = w * x;
	float wy = w * y;
	float wz = w * z;
 
	// This calculation would be a lot more complicated for non-unit length quaternions
	// Note: The constructor of Matrix4 expects the Matrix in column-major format like expected by
	//   OpenGL
	return Matrix4( 1.0f - 2.0f * (y2 + z2), 2.0f * (xy - wz), 2.0f * (xz + wy), 0.0f,
			2.0f * (xy + wz), 1.0f - 2.0f * (x2 + z2), 2.0f * (yz - wx), 0.0f,
			2.0f * (xz - wy), 2.0f * (yz + wx), 1.0f - 2.0f * (x2 + y2), 0.0f,
			0.0f, 0.0f, 0.0f, 1.0f)
}

Quaternion to axis-angle

Given a quaternion , the (non-normalized) rotation axis is simply , provided that an axis exists. For very small rotations, gets close to the zero vector, so when we compute the normalized rotation axis, the calculation may blow up. In particular, the identity rotation has , so the rotation axis is undefined.

To find the angle of rotation, note that w = cos(θ / 2) and .

// Convert to Axis/Angles
void Quaternion::getAxisAngle(Vector3 *axis, float *angle)
{
	float scale = sqrt(x * x + y * y + z * z);
	axis->x = x / scale;
	axis->y = y / scale;
	axis->z = z / scale;
	*angle = acos(w) * 2.0f;
}

Example

Ok, with the above Quaternion class, It's very simple to create a camera class that has one such Quaternion to represent its orientation:

void Camera::movex(float xmmod)
{
	pos += rotation * Vector3(xmmod, 0.0f, 0.0f);
}
 
void Camera::movey(float ymmod)
{
	pos.y -= ymmod;
}
 
void Camera::movez(float zmmod)
{
	pos += rotation * Vector3(0.0f, 0.0f, -zmmod);
}
 
void Camera::rotatex(float xrmod)
{
	Quaternion nrot(Vector3(1.0f, 0.0f, 0.0f), xrmod * PIOVER180);
	rotation = rotation * nrot;
}
 
void Camera::rotatey(float yrmod)
{
	Quaternion nrot(Vector3(0.0f, 1.0f, 0.0f), yrmod * PIOVER180);
	rotation = nrot * rotation;
}
 
void Camera::tick(float seconds)
{
	if (xrot != 0.0f) rotatex(xrot * seconds * rotspeed);
	if (yrot != 0.0f) rotatey(yrot * seconds * rotspeed);
 
	if (xmov != 0.0f) movex(xmov * seconds * movespeed);
	if (ymov != 0.0f) movey(ymov * seconds * movespeed);
	if (zmov != 0.0f) movez(zmov * seconds * movespeed);
}

In this code, xrot, yrot, xmov, ymov and zmod are floats representing how fast the player wants to rotate/move around/on this axis. "seconds" is the time passed since the last call to tick. rotspeed and movespeed represent how fast the camera can rotate or move. PIOVER180 is defined as pi/180, so multiplying with it converts from degrees to radians.

You might be wondering why in rotatex we multiply "rotation * nrot" and in rotatey "nrot * rotation". As I said, multiplication is not commutative. The first rotates the existing quaternion around x (looking up and down), the second rotates an upward-quaternion around the existing rotation. This way, we look left/right around the global y-axis, while rotation up/down is around the local x-axis. This is the behaviour you have in a 3D shooter. Try to change the order of rotations to see what happens.