Download Geometric Definition of Scalar Multiples & Addition in Computer Graphics and more Papers Computer Graphics in PDF only on Docsity!
Linear Combination of Transformations
Marc Alexa
Interactive Graphics Systems Group, Technische Universit¨at Darmstadt ∗
Abstract
Geometric transformations are most commonly represented as square matrices in computer graphics. Following simple geometric arguments we derive a natural and geometrically meaningful defi- nition of scalar multiples and a commutative addition of transfor- mations based on the matrix representation, given that the matrices have no negative real eigenvalues. Together, these operations allow the linear combination of transformations. This provides the abil- ity to create weighted combination of transformations, interpolate between transformations, and to construct or use arbitrary transfor- mations in a structure similar to a basis of a vector space. These basic techniques are useful for synthesis and analysis of motions or animations. Animations through a set of key transformations are generated using standard techniques such as subdivision curves. For analysis and progressive compression a PCA can be applied to sequences of transformations. We describe an implementation of the techniques that enables an easy-to-use and transparent way of dealing with geometric transformations in graphics software. We compare and relate our approach to other techniques such as matrix decomposition and quaternion interpolation.
CR Categories: G.1.1 [Numerical Analysis]: Interpolation— Spline and piecewise polynomial interpolation; I.3.5 [Computer Graphics]: Computational Geometry and Object Modeling— Geometric Transformations; I.3.7 [Computer Graphics]: Three- Dimensional Graphics and Realism—Animation;
Keywords: transformations, linear space, matrix exponential and logarithm, exponential map
1 Introduction
Geometric transformations are a fundamental concept of computer graphics. Transformations are typically represented as square real matrices and are applied by multiplying the matrix with a coor- dinate vector. Homogeneous coordinates help to represent addi- tive transformations (translations) and multiplicative transforma- tions (rotation, scaling, and shearing) as matrix multiplications. This representation is especially advantageous when several trans- formations have to be composed: Since the matrix product is asso- ciative all transformation matrices are multiplied and the concate- nation of the transformations is represented as a single matrix.
∗email:alexa@gris.informatik.tu-darmstadt.de
Figure 1: A two-dimensional cow space: Two transformations A and B , both of which include a rotation, a uniform scale, and a translation, form a two-dimensional space of transformations. In this space ( 0 , 0 ) is the identical transformation, ( 1 , 0 ) and ( 0 , 1 ) rep- resent the specified transformations A and B.
For the representation of motion it is necessary to interpolate from one given transformation to another. The common way in computer graphics for blending or interpolating transformations is due to the pioneering work of Shoemake [Shoemake 1985; Shoe- make 1991; Shoemake and Duff 1992]. The approach is to decom- pose the matrices into rotation and stretch using the polar decompo- sition and then representing the rotation using quaternions. Quater- nions are interpolated using SLERP and the stretch matrix might be interpolated in matrix space. Note, however, that the quaternion approach has drawbacks. We would expect that “half” of a trans- formation T applied twice would yield T. Yet this is not the case in general because the factorization uses the matrix product, which is not commutative. In addition, this factorization induces an order dependence when handling more than two transformations. Barr et al. [1992], following Gabriel & Kajiya [1985], have for- mulated a definition of splines using variational techniques. This allows one to satisfy additional constraints on the curve. Later, Ra- mamoorthi & Barr [1997] have drastically improved the computa- tional efficiency of the technique by fitting polynomials on the unit quaternion sphere. Kim et al. [1995] provide a general framework for unit quaternion splines. However, compared to the rich tool-box for splines in the euclidean space, quaternions splines are still diffi- cult to compute, both in terms of programming effort as well as in terms of computational effort. We identify as the main problem of matrix or quaternion repre- sentations that the standard operators are not commutative. In this work we will give geometrically meaningful definitions for scalar product and addition of transformations based on the matrix repre- sentation. We motivate the definitions geometrically. The defini-
Figure 2: Defining scalar multiples of transformations. Intuitively, “half” of a given transformation T should be so defined that applying it twice yields T. This behavior is expected for arbitrary parts of transformations. Consequently, scalar multiples are defined as powers of the transformation matrices.
tions lead to the use of an exponential map into the Lie group of geometric transformations. Once this connection is established we compare our definition to other approaches. The implementation of this approach uses a transform object that transparently offers scalar product and addition operators. This gives API users an easy-to- use, intuitive, and flexible tool whenever it is desirable to combine transforms rather than composing them.
2 Related work
Our approach essentially uses interpolation in Lie groups by means of the exponential map [Marthinsen 2000]. Grassia has introduced this idea for 3D graphics to represent the group of rotations [Gras- sia 1998]. The group of rotations SO(3) and the group of rigid body motions SE(3) are commonly used for motion planning in the field of robotics. Park and Ravani compute interpolating splines for a set of rotations in SO(3) [Park and Ravani 1997]. They com- pare the groups SO(3) and SU(2) (the group of unit quaternions) in detail. One main advantage of using SO(3) for interpolation is bi- invariance, e.g. if two sets of rotations are connected with an affine mapping the resulting curves are connected by the same map. In our context, this property is naturally contained as part of linear- ity. Zefran analyzes SE(3) for general problems in motion planning (see [Zefran 1996] and the references therein). The main problem is that the respective spaces have non-Euclidean geometry and one has a choice of several reasonable metrics [do Carmo 1992; Zefran et al. 1996]. Once a metric is defined, variational methods are used to determine an interpolant [Zefran and Kumar 1998]. In our ap- proach we have rather traded the problem of defining the geomet- rically most meaningful metric and solving a variational problem for simplicity, ease-of-use and transparency. In addition, we ex- tend these methods from rotations and rigid body motion to general transformations.
The results of our techniques are on an abstract level identical to those from Geometric Algebra (GA) [Hestenes 1991], a field recently introduced to the graphics community [Naeve and Rock- wood 2001]. Current implementations of GA [Dorst and Mann 2001] use explicit representations of all sub-elements (i.e. points, lines, planes, volumes), which results in R^3 being represented with 8 × 8 matrices. In a sense, our approach could be seen as an alter- native implementation using more complex operations on the ma- trices, however, in smaller dimension.
3 Motivation and definition of scalar mul- tiples of transformations
Suppose that we have some transformation, T , and we want to de- fine a scalar mutliple, α T. What conditions should such a scalar multiple satisfy? Well, in the particular case α = 12 , i.e., ”half” of T , we want the resulting transformation to have the property that when it’s applied twice, the result is the original transformation T , i.e., that (^) ( 1 2
T
T
= T ; (1)
an illustration of our goal is given in Figure 2. We’ll require analogous behavior for one-third of a transforma- tion, one fourth, and so forth. We’ll also want α T to be a contin- uous function of both α and T. Let’s explore what this entails by examining the consequences for some standard transformations: translation, rotation, and scal- ing.
Translation: If T is a translation by some amount v , then clearly translation by α v is a good candidate for α T ; it satisfies the requirements of equation 1 and its analogues, and has the advantage that it’s also a translation.
Rotation: If T is a rotation of angle θ about the axis v , then ro- tation about the axis v by angle αθ is a good candidate for α T , for simialr reasons.
Scaling: Finally, if T is a scaling transformation represented by a scale-matrix with diagonal entries d 1 , d 2 ,... then diag ( d 1 α , d 2 α ,.. .) is a candidate for α T.
In all three cases, we see that for positive integer values of α, our candidate for α T corresponds to T α^ ; the same is true for the matrix representing the transformation. If we had a way to define arbitrary real powers of a matrix, we’d have a general solution to the problem of defining scalar multiples; we’d define α T to be T α^ (where what we mean by this is that α T is the transformation representated by the matrix M α^ , where M is the matrix for T ). Fortunately, for a very wide class of matrices (those with no neg- ative real eigenvalues), there is a consistent definition of M α^ , and computing M α^ is not particularly difficult (see Appendices A and C). Furthermore, it has various familiar properties, the most critical being that M α^ M β^ = M α+β^ = M β^ M α^ (i.e. scalar multiples of the
Figure 3: The addition of transformations. Given two transformations, A and B , applying one after the other (i.e. multiplying the matrices) generally leads to different results depeneding on the order of the operations. By performing n -th parts of the two transformations in turns the difference of the two orders becomes smaller. The limit n → ∞ could be understood as performing both transformations concurrently. This is the intuitive geometric definition of a commutative addition for transformations based on the matrices.
non-linearity and discontinuity between the group of transforma- tions and the space in which we perform our computations (i.e. the corresponding algebra). For example, both operators are in gen- eral not continous in their input, i.e. small changes in one of the transformations might introduce large changes in the result. Fur- ther potential problems and limitations are discussed together with applications in Section 6.
In order to implement this approach, routines for computing ma- trix exponential and logarithm are required. We suggest the meth- ods described in Appendix C because they are stable and the most complex operation they require is matrix inversion, making them easy to integrate in any existing matrix package.
Using an object-oriented programming language with operator overloading it is possible to design a transform object that directly supports the new operations. The important observation is that the logarithm of a matrix has to be computed only once at the instantia- tion of an object. Any subsequent operation is performed in the log- matrix representation of the transformation. Only when the trans- formation has to be sent to the graphics hardware a conversion to original representation (i.e. exponentiation) is necessary.
Our current implementation needs 3 · 10 −^5 sec to construct a transform object, which is essentially the time needed to compute the matrix logarithm. The conversion to standard matrix represen- tation (i.e. exponentiation) requires 3 · 10 −^6 sec. Timings have been acquired on a 1GHz Athlon PC under normal working conditions.
Note that for most applications transform objects are created at the initialization of processes, while the conversion to standard repre- sentation is typically needed in in every frame. However, we have found the 3μs necessary for this conversion to be negligible in prac- tice.
6 Applications & Results
Using the implementation discussed above, several interesting ap- plications are straightforward to implement.
6.1 Smooth animations
A simple animation from a transformation state represented by A to a transformation B is achieved with C ( t ) = ( 1 − t ) A ⊕ t B , t ∈ [ 0 , 1 ]. Using a cubic Bezier curve [Hoschek and Lasser 1993] al- lows one to define tangents in the start and endpoint of the interpo- lation. Using the Bezier representation, tangents are simply defined by supplying two transformations. Tangents could be used to gen- erate e.g. fade-in/fade-out effects for the transformation. Figure 4 shows a linear and cubic interpolation of two given transformations. To generate a smooth interpolant through a number of key frame transformations Ti one can use standard techniques from linear spaces such as splines [Bartels et al. 1985] or subdivision curves
Figure 4: Interpolation sequences between given transformations A and B. The top row shows a simple linear interpolation using the matrix operators defined here, i.e. ( 1 − t ) A ⊕ t B. The bottom row shows a Bezier curve from A to B with additional control transformations. These extra transformations define the tangents in the start and end point of the sequence.
[Zorin and Schr¨oder 1999]. Note that the transparent implemen- tation of the operators allows solving linear systems of equations in transformations using standard linear algebra packages. Using these techniques one can solve for the necessary tangent matrices which define e.g. a cubic spline. However, we find an interpolating subdivision scheme (e.g. the 4pt scheme [Dyn et al. 1987]) partic- ularly appealing because it is simple to implement.
It seems that implementations of quaternion splines or other elaborated techniques are hardly available in common graphics APIs. Note how simple the implementation of interpolating or ap- proximating transformation curves is with the approach presented here. One simply plugs the transform object into existing imple- mentations for splines in Euclidean spaces.
The exponential map, on the other hand, has some drawbacks. Essentially, a straight line in parameter space doesn’t necessarily map to a straight line (i.e. a geodesic) in the space of transforma- tions. This means the linear interpolation between two transfor- mations as defined above could have non-constant speed. Further- more, also spline curves, which could be thought of as approxi- mating straight lines as much as possible, are minimizers of a ge- ometrically doubtful quantity. Nevertheless, we found the results pleasing.
We would also like to point at an interesting difference to quater- nions: The log-matrix representation allows angles of arbitrary de- gree. Computing the logarithm of a rotation by π and then mul- tiplying this log-matrix leads to a representation of rotations more than 2π. While this could be useful in some applications it might be disadvantageous in others. For example, the interpolation between two rotations of ±(π − ε) results in a rotation by almost 2π rather than a rotation by 2ε. However, using the tools presented in the following section this could be easily avoided.
We have compared the computation times of this approach with standard techniques. A SLERP based on quaternions between two rotation matrices is about 10 times faster than our approach. How- ever, this is only true for the linear interpolation between two trans- formations. Quaternion splines are subtantially slower. They typi- cally do not allow interactively adjusting the key transformations.
6.2 Factoring transformations
Transformations form a linear space in the log-matrix representa- tion. This allows us to write any transformation as a kind of “lin- ear combination” of transformations from an arbitrary “basis”. The quotation marks indicate that this “linear combination” takes place in log-space – an associated space in which such combinations make sense. For example, three rotations Rx , Ry , Rz by an angle 0 < φ < π around the canonical axes form a basis for the subspace of rotations. Since they are orthogonal, any transformation T can be factored by computing inner products of the log-representation:
x = 〈log T , log Rx 〉, y = 〈log T , log Ry 〉, z = 〈log T , log Rz 〉, (8) where the inner product is computed entry-wise, i.e. 〈{ ai j }, { bi j }〉 = (^) ∑ ai jbi j. Note that the values x , y , z do not represent Euler angles because the rotations around the axes are performed concurrently and not one after the other. Rather, x , y , z define axis and angle of rotation with ( x , y , z )/||( x , y , z )|| being the axis and ( x + y + z )/φ being the angle. The factors x , y , z could be useful to avoid the interpolation prob- lem mentioned at the end of the last Section. Assuming a represen- tation as above the inner products will lead to ( x , y , z ) ∈ [− r , r ]^3 , where r depends on the angle of rotation in each of Rx , Ry , Rz. Specifically, values − r and r represent the same orientation and one can imagine the angles to form a circle starting in − r and end- ing in r with 0 diametrical to ± r. To interpolate along the shortest path from R 1 to R 2 one chooses for each of the factors x 1 , y 1 , z 1 and x 2 , y 2 , z 2 the shorter path on the circle. Specifically, if the differ- ence between two corresponding factors is larger than | r |, then the shorter interpolation path is via ± r rather than via 0. Clearly, factoring has more applications than analyzing rotations. It could be done with respect to any (application specific) orthogo- nal or non-orthogonal basis. In order to find the representation of a transformation T in an arbitrary transformation basis { Bi } we first compute the inner products of the bases. The matrix
V =
〈log B 1 , log B 1 〉... 〈log B 1 , log Bn 〉 .. .
〈log Bn , log B 1 〉... 〈log Bn , log Bn 〉
describes a mapping from the orthogonal canonical base to the pos- sible skew or deficient one formed by { Bi }. Computing the inverse
Acknowledgements
I would like to thank Roy Mathias for introducing me to ma- trix functions and Reinhard Klein for discussions on transforma- tions and Lie groups. Johannes Behr has helped with coding and provided the implementation of the walking humanoid animation. Wolfgang M¨uller and the anonymous referees have given invaluable advice and helped tremendously in improving this document. This work has been supported by the BMBF grant “OpenSG PLUS”.
References
ALEXA, M., AND M ¨ULLER, W. 2000. Representing animations by principal compo- nents. Computer Graphics Forum 19 , 3 (August), 411–418. ISSN 1067-7055. BARR, A. H., CURRIN, B., GABRIEL, S., AND HUGHES, J. F. 1992. Smooth interpolation of orientations with angular velocity constraints using quaternions. Computer Graphics (Proceedings of SIGGRAPH 92) 26 , 2 (July), 313–320. ISBN 0-201-51585-7. Held in Chicago, Illinois. BARTELS, R. H., BEATTY, J. C., AND BARSKY, B. A. 1985. An introduction to the use of splines in computer graphics. DENMAN, E. D., AND BEAVERS JR., A. N. 1976. The matrix sign function and computations in systems. Appl. Math. Comput. 2 , 63–94. DO CARMO, M. P. 1992. Riemannian Geometry. Birkh¨auser Verlag, Boston. DORST, L., AND MANN, S. 2001. Geometric algebra: a computation framework for geometrical applications. submitted to IEEE Computer Graphics & Applications. available as http://www.cgl.uwaterloo.ca/˜smann/Papers/CGA01.pdf. DYN, N., LEVIN, D., AND GREGORY, J. 1987. A 4-point interpolatory subdivision scheme for curve design. Computer Aided Geometric Design 4 , 4, 257–268. GABRIEL, S. A., AND KAJIYA, J. T. 1985. Spline interpolation in curved space. In SIGGRAPH ’85 State of the Art in Image Synthesis seminar notes. GOLUB, G. H., AND VAN LOAN, C. F. 1989. Matrix Computations , second ed., vol. 3 of Johns Hopkins Series in the Mathematical Sciences. The Johns Hopkins University Press, Baltimore, MD, USA. Second edition. GRASSIA, F. S. 1998. Practical parameterization of rotations using the exponential map. Journal of Graphics Tools 3 , 3, 29–48. ISSN 1086-7651. HESTENES, D. 1991. The design of linear algebra and geometry. Acta Applicandae Mathematicae 23 , 65–93. HIGHAM, N. J. 1997. Stable iterations for the matrix square root. Numerical Algo- rithms 15 , 2, 227–242. HORN, R. A., AND JOHNSON, C. A. 1991. Topics in Matrix Analysis. Cambridge University press, Cambridge. HOSCHEK, J., AND LASSER, D. 1993. Fundamentals of computer aided geometric design. ISBN 1-56881-007-5. JOLLIFFE, I. T. 1986. Principal Component Analysis. Series in Statistics. Springer- Verlag. KENNEY, C., AND LAUB, A. J. 1989. Condition estimates for matrix functions. SIAM Journal on Matrix Analysis and Applications 10 , 2, 191–209. KIM, M.-J., SHIN, S. Y., AND KIM, M.-S. 1995. A general construction scheme for unit quaternion curves with simple high order derivatives. Proceedings of SIG- GRAPH 95 (August), 369–376. ISBN 0-201-84776-0. Held in Los Angeles, Cali- fornia. LENGYEL, J. E. 1999. Compression of time-dependent geometry. 1999 ACM Sympo- sium on Interactive 3D Graphics (April), 89–96. ISBN 1-58113-082-1. MARTHINSEN, A. 2000. Interpolation in Lie groups. SIAM Journal on Numerical Analysis 37 , 1, 269–285. MOLER, C. B., AND LOAN, C. F. V. 1978. Nineteen dubious ways to compute the matrix exponential. SIAM Review 20 , 801–836. MURRAY, R. M., LI, Z., AND SASTRY, S. S. 1994. A Mathematical Introduction to Robotic Manipulation. CRC Press. NAEVE, A., AND ROCKWOOD, A. 2001. Geometric algebra. SIGGRAPH 2001 course #53. PARK, F. C., AND RAVANI, B. 1997. Smooth invariant interpolation of rotations. ACM Transactions on Graphics 16 , 3 (July), 277–295. ISSN 0730-0301. RAMAMOORTHI, R., AND BARR, A. H. 1997. Fast construction of accurate quater- nion splines. Proceedings of SIGGRAPH 97 (August), 287–292. ISBN 0-89791- 896-7. Held in Los Angeles, California. SHOEMAKE, K., AND DUFF, T. 1992. Matrix animation and polar decomposition. Graphics Interface ’92 (May), 258–264. SHOEMAKE, K. 1985. Animating rotation with quaternion curves. Computer Graphics (Proceedings of SIGGRAPH 85) 19 , 3 (July), 245–254. Held in San Francisco, California. SHOEMAKE, K. 1991. Quaternions and 4x4 matrices. Graphics Gems II , 351–354. ISBN 0-12-064481-9. Held in Boston. WEB3D CONSORTIUM. 1999. H-Anim. http://ece.uwaterloo.ca:80/˜h-anim. ZEFRAN, M., AND KUMAR, V. 1998. Rigid body motion interpolation. Computer Aided Design 30 , 3, 179–189. ZEFRAN, M., KUMAR, V., AND CROKE, C. 1996. Choice of riemannian metrics for rigid body kinematics. In ASME 24th Biennial Mechanisms Conference. ZEFRAN, M. 1996. Continuous methods for motion planning. PhD-thesis, U. of Pennsylvania, Philadelphia, PA. ZORIN, D., AND SCHR ODER¨ , P. 1999. Subdivision for modeling and animation. SIGGRAPH 1999 course # 47.
A Existence of matrix roots
In the following we will analyze the conditions for the existence of matrix roots, which are intuitively parts of the transformation such that all parts are identical and their combined application yields the original transformation. We will rather use this intuitive geometric point of view – a formal proof of the claims made here could be found in [Horn and Johnson 1991, Thm. 6.4.14]. First, it is clear that a reflection cannot be split into several equiv- alent parts and, consequently, transformation matrices must have positive determinant. This property is obviously necessary, how- ever, not sufficient. To understand this, we need to analyze the eigenvalues of the transformation matrix as they are representative for the nature of the transformation. Note that the product of all eigenvalues is the determinant and, therefore, has to be real. If all eigenvalues are (real) positive the transformation is a pure scale and taking roots is simple. If the eigenvalues have an imag- inary part the respective transformation has a rotational (or shear) component. Because the product of all eigenvalues is real they form two conjugate groups. These groups stay conjugate when roots of the eigenvalues are taken so that determinant and, thus, the trans- formation matrix is still real. A problem occurs in case of real negative eigenvalues (i.e. the imaginary part is zero), which is why we have excluded these transsformations so far. Taking roots of these values introduces imaginary parts in the determinant. Because the determinant is positive the number of negative eigenvalues has to be even, which allows one to analyze them pairwise. A pair of eigenvalues essen- tially defines a transformation in 2D and since both eigenvalues are real and negative this transformation contains a scale part and ei- ther a rotation by π or two reflections. If both eigenvalues have the same magnitude the transformation is a rotation by π and a uniform scale. Taking roots intuitively means reducing the angle of rota- tion and adjusting the uniform scale. However, if the eigenvalues have different magnitude the corresponding transformation can be seen as two reflections or as a rotation together with a non-uniform scale. It is impossible to split this transformation into equivalent parts, because the non-uniform scale is orientation dependent and the orientation changes due to the rotation. Note that compared to other rotational angles it is not possible to interpret this transfor- mation as a shear. Rephrasing this in terms of eigenvalues: if the imaginary parts have same magnitude their roots could be assigned different signs so that they form a conjugate pair; if they have dif- ferent magnitude this is not possible. Concluding, a transformation is divisible, if the real negative eigenvalues of the matrix representing the transformation have even multiplicity. Assuming a positive determinant, it is not divisible if it has a pair of negative real eigenvalues with different magni- tude. Geometrically, a pair of real negative eigenvalues with differ- ent magnitude indicate a rotation by π together with a non-uniform scale. This is the only type of transformation that cannot be han- dled (without further treatment) with the approach presented here. A rotation by π together with uniform scales as well as other rota- tional angles together with non-uniform scales are permissible. For later use we denote this class of transformation matrices T. Note, however, that for divisible transforms with real negative eigenval- ues there is no preferred choice for the primary roots and, thus, the scalar multiplication operation is not continous for such arguments.
B Matrix exponential and logarithm & Lie products
The connection between the matrix operators defined in Sections 3 and 4 and matrix exponential and logarithm is not quite obvious. Recall the definition of exponential and logarithms from Equa-
tions 3 and 4. One has to be careful when carrying over equivalence transforms for exponentials and logarithms from the scalar case, i.e. eA + B , eAeB , and eBeA^ are generally not the same. However, a suf- ficient condition for the expressions to be the same is that A and B commute (see [Horn and Johnson 1991, Thm. 6.2.38]). This leads to the identities
emA^ = eA +...+ A^ = eA^ · · · eA^ =
eA
) m , m ∈ N.
Assuming that emA^ ∈ T we can take m -the roots on both sides^1 ,
( emA
) 1 / m = eA^ ⇔ e
1 m A^ =
eA
) 1 / m
thus
erA^ =
eA
) r , r ∈ Q. (9)
By definition e log^ A^ = A and log eA^ = A. Setting A = log B in Eq. 9 and assuming the logarithms of both sides exist yields
log
er^ log^ B
= log
e log^ B
) r
r log B = log ( Br ). (10)
This immediately leads to the result for the scalar multiplication given in Equation 5. From this connection of roots and logarithms it is clear that real matrix logarithms exist exactly when real matrix roots exist (see also [Horn and Johnson 1991, Thm. 6.4.15]). As said before, eA + B^ and eAeB^ are generally not the same if A and B do not commute. A way of connecting these expressions in the general case is the Lie product formula (see [Horn and Johnson 1991, Chapter 6.5] for a derivation):
eA + B^ = lim n →∞
e
1 n Ae
1 n B
) n (11)
Applying this to log A , log B instead of A and B leads to
e log^ A +log^ B^ = lim n →∞
e
1 n log^ Ae
1 n log^ B
) n
= lim n →∞
e log^ A
) 1 / n ( e log^ B
) 1 / n ) n
= lim n →∞
A^1 / nB^1 / n
) n , (12)
which leads to the representation of the addition given in Equa- tion 6. The use of the standard matrix addition in the exponent proves that the addition operator is indeed commutative.
C Implementation
The computation of exponential and logarithm of rotations or rigid body motions could be performed using Rodrigues’ formula (see [Murray et al. 1994]). The transformations considered here are more general, however, including (non-uniform) scales. We are unclear whether Rodrigues’ formula generalizes to this group and, therefore, propose an implementation based on matrix series. Note that Rodrigues’ formula is the method of choice if scaling is not needed because it is both faster and more robust. The computation of matrix functions such as the exponential and the logarithm is non-trivial. For example, evaluating Equation 4 for computing the exponential is numerically unstable. The preferred way of computing matrix functions is in many cases to use a Schur
(^1) This depends also on our choice of primary roots, which could be a
problem where the primary matrix root function is discontinous, i.e. for matrices with negative real eigenvalues
decomposition and evaluate the function on the upper triangle ma- trix [Golub and Van Loan 1989]. However, this work is intended for graphics where standard ma- trix packages only offer elementary matrix operations. For this reason, implementations are provided using only matrix inversion, multiplication, and addition. For the sake of completeness the pseudo-code from some of the original publications is repeated here. Moler and van Loan [1978] have investigated several ways of computing the exponential of a matrix A in an iterative way and recommend a Pad´e approximation with scaling. Scaling A leads to smaller eigenvalues, which in turn, speeds up the convergence of iterative solvers. This is the pseudo-code of the procedure (see also [Golub and Van Loan 1989]):
Compute X = eA j = max( 0 , 1 + blog 2 (‖ A ‖)c) A = 2 −^ j^ A D = I ; N = I ; X = I ; c = 1 for k = 1 to q c = c ( q − k + 1 )/( k ( 2 q − k + 1 )) X = AX ; N = N + cX ; D = D + (− 1 ) k^ cX end for X = D −^1 N X = X^2 j
The number of iterations q depends on the desired accuracy, q = 6 has proven to be a good choice for the applications intended here. The logarithm of a matrix A can be computed using a truncated Taylor series. However, convergence is not guaranteed or poor if A is not near the identity matrix. By using the identity log A = 2 k^ log A^1 /^2
k and, thus, repeatedly taking the square root A can be made close enough to identity. Exploiting this equation and using a Taylor approximation has been introduced by Kenney and Laub [1989] and leads to the following algorithm:
Compute X = log A k = 0 while ‖ A − I ‖ > 0. 5 A = A^1 /^2 k = k + 1 end while A = I − A Z = A ; X = A ; i = 1 while ‖ Z ‖ > ε Z = ZA ; i = i + 1 X = X + Z / i end while X = 2 k^ X
However, this algorithms needs to compute square roots of ma- trices. Higham [1997] has compared several iterative methods to compute matrix square roots and generally recommends the fol- lowing simple method due to Denman and Beavers [1976]:
Compute X = A^1 /^2 X = A ; Y = I while ‖ XX − A ‖ > ε iX = X −^1 ; iY = Y −^1 X = 12 ( X + iY ); Y = 12 ( Y + iX ) end while
Note that all while loops in the pseudo codes should terminate after a fixed number of iterations since numerical problems might lead to poor convergence.
