website articles
avoiding trigonometry

Intro


I still think we should use less trigonometry in computer graphics. A good understanding of projections, reflections, and vector operations (as in the true meaning of the dot and the cross/external product) usually comes with a growing feeling of uneasy with the use of trigonometry. More preciselly, while I think that trigonometry is good for inputing data to your algorithm (for the notion of angles is an untuitive to measure of orientation), I feel there's something wrong when I see trigonometry involved in the depths of some 3D rendering algorithm. In fact, I think a kitten is sacrificed somewhere every time there's trigonometry invoved down there. And I'm not that concern with speed or precission really, but with conceptual elegance I guess... Well, let me explain.

I have already discussed in other places before how the dot and cross products encode all the information you need to deal with orientation operations, for these two "orthogonal" operations measure cosines and sines of angles respectivelly. They are equivalent to those in so many ways that it feels as if one could simply stick to the products and get rid of angles and trigonometry altogether. In practice, you can indeed do so much by staying in the simple vector-land of Euclides, without trigonometry, that it makes you wander - are we doing somwthing wrong? Probably yes. However, unfortuntally, even experienced professionals tend to abuse of trigonometrics, and so do things way too complicated, cumbersome, and far from elegant. And perhaps indeed, "wrong".

So, without making the article any more abstract, lets consider one such use case of dot and cross products as a replacement for trigonometry, and see what I'm taking about.


The wrong way to do orientate a space/object


Say you have this function, which computes a matrix that rotates vectors around a normalized axis vector v, by an angle a. Any 3D engine or realtime rendering math library will have one such rutine, which has been most probably blindly copied from another engine, wikipedia or an opengl tutorial... (yes, at this point you should confess it, and depending of your mood, perhaps feel bad about it). The function will look something like this:

mat3x3 rotationAxisAngle( const vec3 & v, float a )
{
    const float si = sinf( a );
    const float co = cosf( a );
    const float ic = 1.0f - co;

    return mat3x3( v.x*v.x*ic + co,       v.y*v.x*ic - si*v.z,    v.z*v.x*ic + si*v.y,
                   v.x*v.y*ic + si*v.z,   v.y*v.y*ic + co,        v.z*v.y*ic - si*v.x,
                   v.x*v.z*ic - si*v.y,   v.y*v.z*ic + si*v.x,    v.z*v.z*ic + co );
}

Imagine you are now somewhere in the internals of your demo or game, perhaps writing some animation module, when you find yourself in need to orientate/rotate an object in a give direction. You want to rotate it such that one of its axis, say the z axis, aligns to a given direction vector d (say, the tangent of an animation path). You decide of course to build the matrix that will encode the transformation by using rotationAxisAngle(). So, you first measure the angle between your object's z axis a?nd the desired orientation vector. Since you are a graphics programmer, you know you can do this by taking the dot product and then extracting the angle with an acos() call... Also, because you know that sometimes acosf() can return weird values if your dot products returns anything outside of the -1..1 range, you decide to clamp its argument accordingly (at this point you might even dare blame the machine's precission for not making the length of your normalized vectors exactly 1). At this stage one kitten has already been murdered. But since you don't know about it, you proceed writing your code. So you next compute the rotation axis, which again you know it's the cross product of your z vector with the desired direction d, for all points in your object will rotate in planes parallel to the one defined by those two vectors. Then you decide to normalize the vector, just in case... (kitten has been resurected, and murdered again). In the end, your code looks like this:

    const vec3   axi = normalize( cross( z, d ) );
    const float  ang = acosf( clamp( dot( z, d ), -1.0f, 1.0f ) );
    const mat3x3 rot = rotationAxisAngle( axi, ang );

To see why this works but it's still wrong in many ways, lets expand the code of rotationAxisAngle() in place and see what's really going on:

    const vec3   axi = normalize( cross( z, d ) );
    const float  ang = acosf( clamp( dot( z, d ), -1.0f, 1.0f ) );
    const float  co = cosf( ang );
    const float  si = sinf( ang );
    const float  ic = 1.0f - co;
    const mat3x3 rot = mat3x3( axi.x*axi.x*ic + co,         axi.y*axi.x*ic - si*axiz,    axi.z*axi.x*ic + si*axi.y,
                               axi.x*axi.y*ic + si*axi.z,   axi.y*axi.y*ic + co,         axi.z*axi.y*ic - si*axi.x,
                               axi.x*axi.z*ic - si*axi.y,   axi.y*axi.z*ic + si*axi.x,   axi.z*axi.z*ic + co );
As you have already probably noticed, we are performing an rather imprecise and expensive acos() call, just to undoit immediatelly after by computing a cos() on its return value. So first quesion comes: why not skip the acos/cos chain altogether and save some CPU sycles. Secondly, and more importantly, isn't this observation telling us that we are doing something conceptually wrong and far too complicated, and that there is some sort of simple mathematical principle coming to us, manifesting itself through this expresion simplification?

You might argue that the simplification cannot be done really for one needs the angle anyway in order to feed the sin() function coming right after the cosine. However, this is not true. If you are familiar with the cross product, you might now that the same way dot products encode cosines, cross products encode sines... Most graphics programmers seem to have an intuition for what the dot product does, but a few haven't developed one for the cross product (and only use it to compute normals and rotation axes). Basically, the mathematical principle that is helping us getting rid of the acos/cos pair also tells us that wherever there's a dot product there's probably a cross product around that completes the missing piece of information (the orthogonal piece, the sine).

The proper way to do it


So, yeah, we can extract the sine of the angle between z and d by just looking at the length of their cross product... - remember that z and d are normalized! Which means we can (should!!) rewrite the whole operation this way:
    const vec3   axi = cross( z, d );
    const float  si  = length( axi );
    const float  co  = dot( z, d );
    const mat3x3 rot = rotationAxisCosSin( axi/si, co, si );

and make sure our new rotation matrix building function rotationAxisCosSin() doesn't compute sines and cosines, but takes them as arguments:
mat3x3 rotationAxisCosSin( const vec3 & v, const float co, const float si )
{
    const float ic = 1.0f - co;

    return mat3x3( v.x*v.x*ic + co,       v.y*v.x*ic - si*v.z,    v.z*v.x*ic + si*v.y,
                   v.x*v.y*ic + si*v.z,   v.y*v.y*ic + co,        v.z*v.y*ic - si*v.x,
                   v.x*v.z*ic - si*v.y,   v.y*v.z*ic + si*v.x,    v.z*v.z*ic + co );
}

There's one more thing we can do, and that's to get rid of the normalization/square root, by encapsulating the whole alignment logic in a new function and propagating the 1/si to the matrix:
mat3x3 rotationAlign( const vec3 & d, const vec3 & z )
{
    const vec3  v = cross( z, d );
    const float c = dot( z, d );
    const float k = (1.0f - c)/(1.0-c*c);

    return mat3x3( v.x*v.x*k + c,     v.y*v.x*k - v.z,    v.z*v.x*k + v.y,
                   v.x*v.y*k + v.z,   v.y*v.y*k + c,      v.z*v.y*k - v.x,
                   v.x*v.z*K - v.y,   v.y*v.z*k + v.x,    v.z*v.z*k + c    );
}

Lastly Zoltan Vrana noted that k can be simplified to k = 1/(1+c), which not only is mathematically more elegant, but also moves the singularity in k (and so the whole function) to the case the vectors d and z are parallel (rather than perpendicular), in which case there's no clear rotation that is best anyways. So, the final function looks like this:
mat3x3 rotationAlign( const vec3 & d, const vec3 & z )
{
    const vec3  v = cross( z, d );
    const float c = dot( z, d );
    const float k = 1.0f/(1.0f+c);

    return mat3x3( v.x*v.x*k + c,     v.y*v.x*k - v.z,    v.z*v.x*k + v.y,
                   v.x*v.y*k + v.z,   v.y*v.y*k + c,      v.z*v.y*k - v.x,
                   v.x*v.z*K - v.y,   v.y*v.z*k + v.x,    v.z*v.z*k + c    );
}

Not only we have saved 3 trigonometric functions and avoided an ugly clamp (and a normalization!), but in fact we have conceptually simplified our 3D maths. No trascendental functions, only vectors invoved here. Vectors construct matrices that transforms vectors. And ths is important because, remember, the less trigonometry in your 3D engine, not only the faster and more robust it is, but above all, the more mathematicall elegant (rcoorect!) it is.