Game Development Reference
while the scaling will displace the vector even more from the shell. So the com-
bined operation acts as if we had applied S twice:
A pq A pq =( RS ) T ( RS )= S T R T RS = S T S = S 2 .
So, unfortunately, we have to take the square root of this matrix equation to
obtain S . As this is a common problem in mathematics and physics, this problem
has been addressed a lot and there are good numerical methods to calculate this
quantity. The usual approach is to diagonalize the matrix S 2 :
S 2 = V diag( λ ) V T ,
where λ are the eigenvalues of the matrix S 2 .
A very good overview, as well as some state-of-the-art algorithms for the di-
agonalization of 3
3 matrices, has been given by [Kopp 06]. Once the matrix is
diagonalized, we can take the square root of the diagonal entries.
S = V diag( √ λ ) V T
to obtain the matrix S .
There is also what is called the Denman-Beavers square root iteration —this
works without diagonalization. It is easy to implement and very robust, although
not as efficient (see [Denman and Beavers 76]).
We will use the Jacobi algorithm here, which is the oldest but is also a very
robust algorithm. It starts off with the identity matrix for V and applies so-called
Jacobi sweeps on it (see [Kopp 06]).
Since the rotation matrices we are dealing with are “almost diagonal” already,
it will take only one to two Jacobi sweeps on the average for each vertex. Since
this operation is done very often, we should think about caching the matrix V
from the previous time step instead of starting off with the identity matrix at every
time step. This induces further memory usage but reduces computation time.
Using S , we can now calculate the rotational part:
R = A pq S − 1 .
I want to thank Ury Zhilinsky for his input and support during my work on this
chapter. The MD5 model used was built upon polygonal data from the Stanford
3D Scanning Repository.