Design Feature: March 2, 1995

Quaternions prove useful for coordinate transformations. The quaternion method is better if one of the coordinate systems keeps moving-which is the usual case in navigation and animated 3-D graphics. You can rotate a quaternion faster than you can rotate a matrix. And, unlike rotation matrices, quaternions never lose their orthogonality, which causes distortion, after many rotations.

Few people use and understand quaternions because most of the literature about quaternions involves obscure mathematical derivations and requires super-imaginary numbers that are the fourth root of -1. I will spare you all that. Instead of talking about how difficult quaternions are to derive, I will show you how easy they are to use.

A quaternion is just an ordered set of four values. You can think of it as a rotational state vector. Simple equations allow you to relate a quaternion to the corresponding rotation matrix, and vice versa.

The quaternion advantage

A quaternion doesn't have one big problem that a rotation matrix has. A rotation matrix is always overspecified. Consider a simple 2-D rotation matrix.

Although the matrix has four elements, it has only one independent variable. The angle [theta] completely determines all four values. You cannot, however, determine [theta] uniquely from just one element of the matrix. That is, if S=0.5, you can't tell if [theta] is 30 or 150°.

But, given just one row or one column of the matrix, you can determine the angle [theta] if both entries are correct. For example, if S=0.5 and C=0.866,[theta] must be 30 rather than 150°.

Suppose, however, an entry is in error, such that S=0.5 and C=0.857. Is the angle 30 or 31°? Representing one angle by two numbers has introduced some uncertainty. Furthermore, two of the four elements in the 2-D rotation matrix are or should be redundant. They provide no additional information when they are accurate, and they provide only confusion when they are inaccurate.

Suppose you want to find the rotation matrix for a 2-D coordinate system after rotating it clockwise 15°, then counterclockwise 5°, then clockwise 20°, and, finally, clockwise 8°. You know that[theta]=-15°+5°-20°-8°=-38°. You can perform that calculation in your head. Therefore, the correct rotation matrix should have Cos(-38°) on the diagonal andąSin(-38°) in the other corners.

But, suppose you actually run these calculations on your computer, computing the rotation matrix for -15°, then multiplying it by the rotation matrix for +5°, then multiplying that product by the rotation matrix for -20°, and, finally, multiplying that result by the rotation matrix for -8°. Again, you should wind up with a rotation matrix with Cos(-38) on the diagonals andąSin(-38) in the corners-but will you? Will A11 be exactly equal to A22? Will A21 exactly equal -A12? Will Square_Root(A11×A11 + A12×A12) exactly equal 1? The answers to these questions all depend on the accuracy of the computer's arithmetic.

Even if the chain of matrix multiplications does yield the right answer, which method would you rather use: mentally adding the angular changes or multiplying the rotation matrices? The mental method involves one addition, three subtractions, one sine, and one cosine. The rotation-matrix method involves four sines, four cosines, and three matrix multiplications. So, in the 2-D case, doing the internal calculations using the first method and then converting the result to a rotation matrix is better.

The 3-D case is similar and entails three independent variables. A rotation matrix maps these three variables to nine matrix elements. The 3-D matrix has added redundancy that gives you nothing but more opportunities for error than does a 2-D matrix. A quaternion, on the other hand, is a rotational state vector that maps the three independent variables into just four elements, greatly reducing the amount of redundancy.

Reducing the redundancy also reduces the number of calculations. To rotate an object, whose orientation a 3-D matrix represents, around a single axis, you must multiply that object's matrix by another rotation matrix.[theta]his procedure requires 27 multiplications, 18 additions, and three assignment operations. To rotate a quaternion about a single axis, you must multiply that quaternion by another quaternion. Quaternion multiplication takes 16 multiplications, 12 additions, and four assignments. Quaternions save 11 multiplications and six additions at the expense of one assignment operation for single-axis rotations.

When working in 3-D, you usually rotate an object about the yaw axis, then the pitch axis, and, finally, the roll axis to establish the object's new orientation. Using quaternions instead of matrices saves 33 multiplications and 18 additions for each three-axis rotation. Because time is usually critical in navigation routines and 3-D animation, saving 33 multiplications per rotation is desirable.

Any four random values—as long as all four aren't zero—create a valid rotation quaternion. This property makes the quaternion useful for generating randomly oriented test cases to test coordinate-conversion algorithms. On the other hand, any nine random values probably don't create a valid rotation matrix because certain elements of the matrix must relate to other elements in a prescribed way. Multiplying an object's orientation matrix by a random matrix not only rotates the object, but also distorts its shape.

This distortion is important because navigation and animation programs multiply the object's orientation matrix every time its position changes. It doesn't take long for the programs to multiply the matrix thousands of times. Truncation and rounding errors add a tiny bit of randomness to each element of the matrix every time the programs multiply it. Soon, the matrix is no longer perfectly orthogonal. The effect is the same as if you multiply the matrix by a random matrix-&$151;distorting the object's shape.

For several reasons, the problem is worse in a navigation routine than in a 3-D animation program. First, because navigation routines run for hours or days longer than animation routines, errors have more time to build. Second, the distortion of the "shape" of the matrix distorts the guidance commands,causing coupling. Such coupling could cause starboard-yaw commands to have the side effect of producing small counterclockwise roll commands. Third, people's lives depend on navigation routines more often than they do on animation routines. Therefore, keeping the orientation matrix orthogonal is important.

You can keep the matrix orthogonal by using double-precision arithmetic or by periodically normalizing the matrix. Unfortunately, both of these solutions take time. But, unlike matrices, quaternions never become distorted. By the way, you can easily fix a matrix that has become distorted and needs to be normalized by converting it to a quaternion and then converting the quaternion back to a matrix. I guarantee that the resulting matrix will be orthogonal.

Quaternions will never completely replace rotation matrices. The rotation matrix is unsurpassed for rotating vectors. Just multiply the matrix by the vector and get the vector result. But this facility with vectors doesn't mean you have to use the rotation matrix for all the intermediate transformation calculations. Computing a series of 2-D rotations as the sum of angular rotations and converting the final result to a matrix is better; computing a series of 3-D rotations using a quaternion is better, too. Do all the rotation calculations on a quaternion; then, convert the quaternion to a rotation matrix.

Quaternion/matrix conversions

The conversions between quaternions and rotation matrices are simple. Consider quaternion Q and an equivalent rotation matrix, M.

Most mathematicians would tell you to use the following equations to convert Q to M:

M11 = Q12 + Q22 - Q32 - Q42C

M12 = 2*(Q2*Q Q1*Q4)
	 
M13 = 2*(Q2*Q4 + Q1*Q3)	

M21 = 2*(Q2*Q3 + Q1*Q4)
 
M22 = Q12 - Q22 + Q32 - Q42

M23 = 2*(Q3*Q4-Q1*Q2)

M31 = 2*(Q2*Q4-Q1*Q3)
	
M32 = 2*(Q3*Q4 + Q1*Q2)	

M33 = Q12 - Q22 - Q32 + Q42

These equations assume that quaternion Q has been normalized. If Q is not normalized, you must divide each of the elements by Q1² + Q2² + Q3² + Q4².

These equations aren't difficult to compute as they stand, but clever collection of terms in five temporary variables (TX, TY, TZ, TQ and TK) can further reduce computation time. Algorithm 1 is a better way to convert a quaternion to a matrix.

Algorithm 1
TX	=	Q2*Q2

TY	=	Q3*Q3
	  
TZ	=	Q4*Q4
	  
TQ	=	TY+TZ
	  
if (TQ + TX + Q1*Q1) is not 0 then 

TK = 2 / (TQ + TX + Q1*Q1)

else

TK	=	0

M11	=	 1 - TK*TQ 

M22	=	 1 - TK*(TX + TZ)

M33	=	 1 - TK*(TX + TY)

TX	=	TK*Q2
 	 
TY	=	TK*Q3
	  
TQ	=	(TK*Q4)*Q1 

TK	=	TX*Q3

M12	=	TK - TQ
	  
M21	=	TK + TQ

TQ	=	TY*Q1

TK	=	TX*Q4

M13	=	TK+TQ
 	
M31	=	TK-TQ

TQ	=	TX*Q1
	  
TK	=	TY*Q4

M23	=	TK - TQ
	  
M32	=	TK + TQ

You can convert M to Q in any of four easy ways. Some ways work better than others for certain situations. The conversion is much like computing an arctangent. You can write a series that computes the arctangent of Y/X. The series converges rapidly when X is much larger than Y, but it takes a long time to converge when X is small, especially when X is 0. Therefore, good arctangent routines take advantage of the fact that arctangent (Y/X)=90°- arctangent(X/Y). Good arctangent routines use the alternate form whenever Y>X. Good arctangent routines also use tricks to fold angles from the second, third, and fourth quadrants into the first quadrant before computing the result, and the routines then adjust the result based on the quadrant.

Algorithm 2 selects one of the four matrix-to-quaternion conversion routines, depending on the values on the diagonal of the rotation matrix.

Algorithm 2
if M11 => 0 and (M22 + M33)=> 0 then 

Q1	=	1 + M11 + M22 + M33

Q2	=	M32 - M23 

Q3	=	M13 - M31
	 	
Q4	=	M21 - M12


if M11 =>0 and (M22 + M33) 0 then

Q1	=	M32 - M23 
	
Q2	=	1 + M11 - M22 - M33

Q3	=	M21 + M12  
	
Q4	=	M13 + M31


if M11 0 and (M22 - M33) =>0 then

Q1	=	M13 - M31
	
Q2	=	M21 + M12

Q3	=	1 - M11 + M22 - M33

Q4	=	M32 + M23


if M11 0 and (M22 - M33) 0 then

Q1	=	M21 - M12
	 	
Q2	=	M13 + M31

Q3	=	M32 + M23
	 	
Q4	=	1 - M11 - M22 + M33

Quaternion multiplication

Just as you can multiply two rotation matrices to get another rotation matrix, you can multiply quaternions to compute the effect of a series of rotations. Like matrix multiplication, quaternion multiplication is not commutative because the order of rotation matters. To see why the order matters, imagine an airplane pointing north. If you rotate it 90° to its right, it will be pointing east. If you then pitch it up 90°, it will be standing on its tail. But if you take an airplane pointing north, pitch it up 90°, and then rotate it to its right 90°, it will wind up standing on its starboard wing tip.

Let's define three quaternions to have elements Q, L, and R:

        Q1	       L1	        R1 

    Q = Q2         L = L2           R = R2 	

        Q3	       L3	        R3 

        Q4	       L4	        R4

To compute the quaternion product Q=L*R, use Algorithm 3:

Algorithm 3

Q1	=	L1 * R1 -  L2 * R2 -  L3 * R3 -  L4 * R4 

Q2	=	L2 * R1 + L1 * R2 -  L4 * R3 + L3 * R4 

Q3	=	L3 * R1 + L4 * R2 + L1 * R3 -  L2 * R4 

Q4	=	L4 * R1 -  L3 * R2 + L2 * R3 + L1 * R4 

Rotating quaternions

A quaternion rotates an object an angle [theta] about an axis. Let Ai + Bj + Ck be a unit vector outward from the origin representing an axis of revolution. Let [theta] represent the (right-handed) angle of rotation about that axis. Algorithm 4 shows how to find the quaternion that produces this rotation.

Algorithm 4

Q1	=	 Cos([theta] / 2)

Q2	=	A * Sin([theta] / 2)

Q3	=	B  * Sin([theta] / 2)
	 
Q4	=	C  * Sin([theta] / 2)

Algorithm 4a

To rotate in the yaw direction:

Q1 = Cos([theta] / 2)

Q2 = 0

Q3 = 0
	  
Q4 = -Sin([theta] / 2)


Algorithm 4b

To rotate in the pitch direction:

Q1 = Cos([theta] / 2)

Q2 = Sin([theta] / 2)
	  
Q3 = 0 Q4 = 0


Algorithm 4c

To rotate in the roll direction:

Q1 = Cos([theta] / 2)

Q2 = 0

Q3 = Sin([theta] / 2)

Q4 = 0

Rotation rates

Calculations often yield the rotation rates about the three axes. Although the order of rotation is important in theory, the rotation order does not matter much in practice if the magnitude of the rotation is small. Therefore, if rates are in radians per second, and DELTA_T is a small value (seconds), Algorithm 5can compute the new orientation of the quaternion, Q.

Algorithm 5

Y	=	YAW_RATE   *  DELTA_T	

P	=	PITCH_RATE *  DELTA_T	

R	=	ROLL_RATE   *  DELTA_T


T	=	1 - (R*R + P*P + Y*Y) / 12

T1	=	Q1*T - (R*Q2 + P*Q3 + Y*Q4) / 2 

T2	=	Q2*T + (R*Q1 - P*Q4 + Y*Q3) / 2 

T3	=	Q3*T + (R*Q4 + P*Q1 - Y*Q2) / 2


Q4 = Q4*T + (-R*Q3 + P*Q2 + Y*Q1) / 2

Q3 = T3	

Q2 = T2	

Q1 = T1

This second-order approximation is accurate to 1 sec of error when the values of Y, P, and R are less than 0.17 radians (10°). When you write the code, multiply by the reciprocals of the constants because multiplication is faster than division.

If a quaternion, Q, represents the current orientation of an object and if you want small angles Y, P, and R to rotate the object (in radians), computing the quaternion's new orientation takes only 24 multiplications, 15 additions, and no sines or cosines. That calculation is much faster than the 81 multiplications and 54 additions multiplying an orientation matrix by three rotation matrices would take-and that isn't counting the three sines and three cosines you would have to compute to create the three rotation matrices.

The derivations of these simple algorithms are complex. Those complex, inscrutable derivations tend to scare people away from quaternions. Fortunately, you don't need to derive the equations every time you use these algorithms. Just pretend they are matrix operations optimized for speed, with the fortunate side effect of keeping the matrix orthogonal after every multiplication.

Do-While Jones is a software engineer who has worked for the defense industry of a free-world nation since 1971. He has obtained a patent for a radar signal-processing algorithm and has an EE degree. The author of the book Ada in Action, he has published more than 40 articles in popular computer magazines. You may contact him at do_while@ridgecrest.ca.us on the Internet.