Inigo Quilez   ::     ::  

Intro


I think we should use less trigonometry in computer graphics. As my understanding of projections, reflections, and vector operations improved over time, a growing unease grew as well every time I saw trigonometry at the core of 3D algorithms. These days I'm at a point where if I see some asin, acos, atan, sin, cos or tan in the middle of a 3D algorithm, I'll assume the code was written by an inexperienced programmer and needs review.

Now, don't get me wrong. Trigonometry is convenient and necessary for data input, for feeding the larger algorithm. What's wrong is when angles and trigonometry suddenly emerge deep in the internals of a 3D engine or algorithm out of nowhere, where only vectors and geometry were involved as inputs. And each time such a wrong happens, a kitten is being killed. Not a kitten that cared about performance or precision necessarily, but a kitten that appreciated elegance and truth. Let me explain.

I've already discussed in other places before how the dot and cross products encode all the information you need to deal with orientation related operations. These two products capture all the information that you've learnt to associate with the cosine and sine of angles, and more. And since they operate on vectors, which are the native way we represent objects, there's no reason to not stay in vector-land all the way and stay away from the more abstract (less dimensional) world of trigonometry and its multiple identities. So, let's make this claim of mine more concrete with an example (we could the same exercise with many others, this is just one).


The wrong way to do orientate a space/object


Say you have this function, which computes a matrix that rotates vectors around a given axis v by an amount of a radians. Many 3D engines, renderers or math libraries will have one such routine. 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 ); }

So now imagine you are somewhere in the internals of your project, writing some code that needs to set the orientation of an object in a given direction. For example, you are aligning a spaceship to an animation path, by making sure the spaceship's z axis aligns with the path's tangent or direction vector d. So, you decide you want to use rotationAxisAngle() function for that. Then, you first measure the angle between the spaceship's z axis and the desired orientation vector d. Since you are a graphics programmer, you know you can do this by taking the dot product of the two vectors and then extracting the angle the form from it with an acos() call. Great. Also, because you know that sometimes acos() can return nonsense if its argument lays outside of the -1..1 range, you decide to clamp the dot product to -1..1, "because floating point". Fine. But at this stage one kitten has already been murdered. But since you don't know that yet, you proceed to writing the rest of the code. So, next you compute the rotation axis, which you know it's the cross product of your z and d vectors, for all points in the spaceship need to rotate in planes parallel to that spanned by these two vectors. Then you decide to normalize the rotation axis vector, just because you like normalizing things. And so, 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 );

This does work. But I claim it's wrong. To see why, let's expand the code of rotationAxisAngle() in place and see the whole picture of what's mathematically 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 immediately after by computing a cos() on its return value. This begs the question: why not skip the acos() and cos() chain altogether. Isn't this simplification more than just an optimization? Could it be a deeper mathematical structure or relationship manifesting itself? Some simplifications and symmetrys in mathematical expressions are certainly accidental and meaningless, but more often than not they are the conequence of deep truths, so let's keep exploring see if we can learn anything here.

Even though we could potentially simplify cos(acos(x)), it seems we need to compute the angle anyways for our sin() function coming right after the cosine. However... if you are familiar with the cross product, you might now that cross products encode sines just like dot products encode cosines. It might not be surprising though you never realized this, and that's fine, that's why we are here. So this is the point - wherever there's a dot product, chances are there's a cross product lying around that completes the picture of what's geometrically being described. For any parallel thing you can measure (dot products, cosines) there's a perpendicular component that you can measure too (cross product, sines). So, check this out:

The proper way to do it


We can extract the sine of the angle between z and d by just looking at the length of their cross product. Indeed, z and d are normalized so the length of their cross product is the sine of the angle they form times one and times one again! Which means we can and should rewrite our operation this way:

mat3x3 rotationAlign( const vec3 & d, const vec3 & z ) { const vec3 v = cross( z, d ); const float si = length( v ); const float co = dot( z, d ); 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 ); }

Good, we are getting somewhere now. Let's keep going and also get rid of the normalization/square root by 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.0f-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, we can simplify k to k = 1/(1+c) which is more elegant and also moves the two singularities in k and so the whole function (d and z are parallel) to a single one (when d and z are the same vector, in which case there's no clear rotation that is best anyways (Lastly Zoltan Vrana for the tip). 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 ); }

Conclusion


So look at that. There aren't any trigonometric or inverse trigonometric functions, nor clamp or normalizations. This is great from a stability and performance point of view. And also, we have solved a vector problem using only vectors, which feels like it should have been the right thing right from the beginning.

Now, I'm not going to claim that in this particular example it's clearly intuitive what that matrix means. We can see that it is

v * v^T ( c, -v.z, v.y) M = ------- + ( v.z, c, -v.x) 1+c (-v.y, v.x, c)

which has something very quaterniony to it if you ask me, with that "real" component c and that vector component v, but that's for another day. The point is, almost always you can perform the same exercise of simplifying trigonometric calls and slowly untangle a simpler vector form expression for the same problem. Also, oftentimes, instead of discovering a vector form that way, you can actually rethink the whole problem from scratch as a vector problem from the beginning. You need to have some intuition of dot products, cross products, projections, etc for that, but once you got it you'll find that most trigonometric identities out there are just saying something about simple vector configurations, and you'll never have to memorize one again. But that also is for another day.

I'll just conclude saying that part of the problem I find in real life is due to poorly designed APIs, like rotationAxisAngle(v,a) or WebXR stereo projection information, which work on angles instead of actual vectors. Sometimes the cost is a performance and stability cost like we saw, other times the cost takes the form of unnecessary confusion (are the angles in radians or degrees, and they negative or positive in the inner side of the frustum, etc etc). Either way, let's strive to do more vectors operations and less trigonometry!