Quaternion Tricks

Trick #1 : Optimizing the Quaternion Rotation Formula

The quaternion rotation formula is given as

(0, v') = R(0, v)R^†

This formula can be expressed in terms of standard 3D vector products such as scaling, dot products, and cross products. This equation involves two quaternion products, each of which requires 16 multiplies and 12 additions. Thus there are 56 operations in all.

We can express the vector cross product in terms of the quaternion product. We know that if you switch the factors in a cross product, the product switches sign. Knowing this, you might accept without proof the following formula.

q_V × p_V = (1/2) (qp - pq)

Indeed if you flip the order of the arguments, the product will change sign.

We will use this relationship, in a rearranged form, on the first two factors in the quaternion rotation formula.

R(0, v) = (0, 2R_V × v) + (0, v)R

Inserting the expression on the right hand side into the quaternion rotation formula yields

(0, v') = (0, 2R_V × v)R^† + (0, v)RR^†

When examining the first term on the right hand side, it can be shown that the scalar portion is always zero. For the second term, we can use the fact that the rotation quaternion R obeys the identity RR^† = 1. Since all 3 terms have zero scalar component, we can now rewrite this expression in terms of vector operations.

v' = v + 2[ R_S( R_V × v ) + R_V × ( R_V × v ) ]

This expression employs 18 multiplications and 12 adds, totalling 30 unique operations in all. This optimisation basically cuts the expense of the quaternion rotation formula in half.

Here is a function that employs this omptimisation

Vector rotateVector(const Vector& v, const Quat& q)
{
     Vector result;
     float x1 = q.y*v.z - q.z*v.y;
     float y1 = q.z*v.x - q.x*v.z;
     float z1 = q.x*v.y - q.y*v.x;

     float x2 = q.w*x1 + q.y*z1 - q.z*y1;
     float y2 = q.w*y1 + q.z*x1 - q.x*z1;
     float z2 = q.w*z1 + q.x*y1 - q.y*x1;

     result.x = v.x + 2.0f*x2;
     result.y = v.y + 2.0f*y2;
     result.z = v.z + 2.0f*z2;

     return result;
}

Let's express the scalar and vector portions of the quaternion in terms of the axis and angle of the rotation

R_S = cos( θ/2 )

R_V = n sin( θ/2 )

Pluging these values into the optimized quaternion rotation formula, and using a little bit of trigonometric magic (double angle identity), the formula reduces to

v' = v + [ sin(θ) (n × v) + ( 1 - cos(θ) ) n × (n × v) ]

Using the vector triple product identity, and the fact the n is normalized, we can reduce it a bit further.

v' = cos(θ)v + sin(θ)(n × v) + (1 - cos(θ)) (n•v) n

This operation has 19 mults and 12 adds, totalling 31 operations. We didn't do anything to optimize the formula, but we were able to transform it into a form involving axis and angle.

Trick #2: Determining a Rotation From an Initial and Final Vector

One of the questions which I am asked most frequently is how to find a rotation which takes some initial vector into some final vector. There are several variants of this problem, including "I am pointing my gun at Bob, and I would like to point my gun at Jim." I have seen this problem solved in many different ways, many times involving inverse trig functions. We are going to solve this problem with quaternions - and it's really easy!

We are going to call the initial vector A, and the final vector B. Both A and B are normalized 3D vectors.

We are going to multiply these two 3D vectors together, and if you remember, this is accomplished with a quaternion multiplication. To reitterate, our input vectors are just standard 3D vectors, but the output of the product is a quaternion, which represents the product of the inputs. Since the ordering matters, we are going to place A on the right, and B on the left.

BA = -A•B - A × B

We can relate the dot and cross products to functions involving the angle between A and B. If A and B are normalized, then their product has the following form

BA = - cos(θ) - n sin(θ)

The angle θ is the angle between the initial and final vector, and the vector n is a normalized vector that is perpendicular to the input vectors. This is the exact form of a rotation quaternion. Don't worry about the minus sign, since we can flip the sign without changing the rotation we are representing. In fact, this almost represents the quaternion that we are looking for. We need a quaternion with a half angle instead of a full angle.

Instead of multiplying the initial vector A by the final vector B, let's instead find a vector that represents a halfway rotation between A and B. This is very easy to do - just find the average of A and B!

H = (1/2) (A + B)

Now, this H is not normalized, and it needs to be. If we are going to normalize it anyway, we don't need to worry about the factor of 1/2. The normalized version of H is

H = (A + B) / A + B

Knowing that the angle between A and H is half of the angle between A and B, we can use our previous result to show that the rotation which takes A into B is given by

R = HA

Here is a function that will give you the rotation quaternion that will rotate some initial vector into some final vector

Quat getRotationQuat(const Vector& from, const Vector& to)
{
     Quat result;

     Vector H = VecAdd(from, to);
     H = VecNormalize(H);

     result.w = VecDot(from, H);
     result.x = from.y*H.z - from.z*H.y;
     result.y = from.z*H.x - from.x*H.z;
     result.z = from.x*H.y - from.y*H.x;
     return result;
}

It doesn't get much easier than that folks.

Trick #3: Tangent Space Compression / Quaternion Maps

Tangent space is a set of 3 vectors - normal, tangent, and bi-normal - which are used in lighting calculations such as bump mapping, or anisotropic shading techniques, like fur or hair shaders.

The tangent space is a representation of the surface of the mesh. The tangent space vectors are stored per-vertex, and interpolated across the polygon for per-pixel operations. If you use floats to represent the components of the tangent space vectors, then you are dedicating 3 x 3 x 4 bytes = 36 bytes per vertex to the storage of tangent space.

The tangent space vectors form an ortho-normal frame. Consider the rotation which will transform the ortho-normal frame of tangent space into the standard frame in world space. If this rotation is expressed as a matrix, then the 3 columns of this matrix represent the 3 tangent space vectors.

If you can express a rotation as a matrix, then you can also express it as a quaternion. This immediately reduces the tangent space representation to 4 components = 16 bytes. When decompressing the vertex, use the quat->matrix conversion, and pull out the 3 columns of the resulting matrix as your tangent space vectors. The quat->matrix conversion costs 13 multiplies and 15 subtractions. It may be worthwhile to determine if it is possible to optimize this conversion for shader operations.

If you want to save more space, you can go even further, if you have the time - computationally speaking. Since rotation quaternions are normalized, there are only 3 degrees of freedom instead of 4. Thus you only really need to store 3 components of the quaternion, and the 4th can be calculated using the formula

w = sqrt(1 - x² - y² - z²)

This brings the tangent space size down to 12 bytes.

Finally, if you are super hard pressed for space, you can take advantage of the fact that the components of a quaternion are within the range [-1,1], which may not require 32 bits to have a sufficient representation. You may possibly use 16 bit floats to represent these components. Doing so will reduce the size further - 6 bytes.

With a 6-byte tangent space representation, it seems conceivable to employ quaternion texture maps, which store the tangent space directly in the texture, rather than using bump map data to perturb the interpolated tangent vectors. You could then apply a texture compression, like 3DC which gets 2:1 compression. Now tangent space only consumes 3-bytes of texture, and 0-bytes on the vertex. If you were previously using bump maps, and you replace them with quaternion maps, then you aren't really increasing the amount of texture memory used, and you have successfully removed all 36 bytes of tangent space data from the vertex.

Of course, doing quaternion map lookups in the pixel shader would also mean doing the quat->matrix conversion in the pixel shader - or, alternatively recasting the pixel shader lighting equation in terms of quaternion algebra, rather than matrix algebra - whichever ends up being cheaper. To gain perspective on the cost of this conversion, you can perform 5 quat->matrix conversions in roughly the same amount of time as a single matrix multiply. A single quat->matrix conversion is about as costly as transforming 3 vectors with a matrix.

The quaternion map would take the place of bump maps altogether, and would provide the added feature of having fine tuned control over tangent space on a per-pixel level. For anisotropic shaders, like hair, fur, and velvet shaders, you could introduce swirls and wiggles in the middle of a polygon.

There are two rather large problems facing the idea of a quaternion map: 1) a lot of people get freaked out when you use the word quaternion. 2) I can't off the top of my head conceive of a way for an artist or content creator to generate such a map.

At any rate, it seems like a very intriguing idea.