Geometric Skinning with Approximate Dual Quaternion Blending
L a d i s l a v K a v a n ∗ 1 S t e v e n C o l l i n s 1 J i ˇr ́ı Zˇ a ́ r a 2 C a r o l O ’ S u l l i v a n 1 1Trinity College Dublin, 2Czech Technical University in Prague
84.9 FPS 197.4 FPS 55.1 FPS
Log-matrix Blending Dual Quaternions Spherical Blend Skinning Dual Quaternions
Figure 1: A comparison of dual quaternion skinning with previous methods: log-matrix blending [Cordier and Magnenat-Thalmann 2005] and s p h e r i c a l b l e n d s k i n n i n g [ K a v a n a n d Zˇ a ́ r a 2 0 0 5 ] . T h e p r o p o s e d a p p r o a c h n o t o n l y e l i m i n a t e s a r t i f a c t s , b u t i s a l s o m u c h e a s i e r t o i m p l e m e n t and more than twice as fast.
Abstract
Skinning of skeletally deformable models is extensively used for real-time animation of characters, creatures and similar objects. The standard solution, linear blend skinning, has some serious drawbacks that require artist intervention. Therefore, a number of alternatives have been proposed in recent years. All of them suc- cessfully combat some of the artifacts, but none challenge the sim- plicity and efficiency of linear blend skinning. As a result, linear blend skinning is still the number one choice for the majority of de- velopers. In this paper, we present a novel skinning algorithm based on linear combination of dual quaternions. Even though our pro- posed method is approximate, it does not exhibit any of the artifacts inherent in previous methods and still permits an efficient GPU im- plementation. Upgrading an existing animation system from linear to dual quaternion skinning is very easy and has a relatively minor impact on run-time performance.
CR Categories: I.3.5 [Computer Graphics]: Computational Geometry and Object Modeling – Geometric Transformations— [I.3.7]: Computer Graphics—Three-Dimensional Graphics and Re- alism – Animation
Keywords: skinning,rigidtransformations,blending,dualquater- nions, linear combinations
1 Introduction
Skinning and skeletal animation is the technology behind charac- ter animation in many applications. In some situations, physically accurate skin deformation, which supports muscle bulging and dy- namic effects, is desirable. In other situations, however, a fast algo-
∗e-mail: kavanl@cs.tcd.ie
rithm capable of skinning multiple models interactively is needed, for example in videogames and crowd simulations.
The standard algorithm for low-cost skinning is known by many names: linear blend skinning, vertex blending, skeletal subspace deformation or enveloping. It is sometimes used not only for skin deformation (as the name suggests) but also to animate other deforming elements, for example cloth, because it is consider- ably faster than physically based cloth simulation [Cordier and Magnenat-Thalmann 2005]. The basic principle is that skinning transformations are represented by matrices and blended linearly. It is very well known that the direct linear combination of matri- ces is a troublesome way of blending transformations. This pro- duces artifacts in the deformed skin, even if we restrict the skinning transformations to rigid ones (i.e., composition of a rotation and translation). In spite of these shortcomings, linear blending is still a very popular skinning method, but perhaps only because there is no simple alternative.
Recent previous work suggests converting rigid transformation ma- trices to (quaternion,translation) pairs and blending them instead of their matrix equivalents [Hejl 2004; Kavan and Zˇa ́ra 2005]. This works, but at a cost: Hejl’s algorithm [2004] imposes con- straints on the model’s rigging (specifically, a vertex can only be influenced by neighbouring bones, otherwise artifacts can occur). This could be inconvenient, because linear blending has no such constraints (and it is exploited for many 3D models). Kavan and Zˇ a ́ra’s method [2005] does not have this restriction, but uses a com- plex and computationally expensive Singular Value Decomposition scheme. Obviously, the simplicity of linear blend skinning is lost in both cases.
The representation of rigid transformations by matrices or (quaternion,translation) pairs illustrates just two possible pa- rameterizations of SE(3), i.e., the group of rigid transforma- tions. Nothing prevents us from blending, for example, 3-tuples (axis,angle,translation) or pairs (axis,translation · sin(angle)). Even if we restrict ourselves to blending via linear combinations (motivated by efficiency and simplicity of implementation), we can construct infinitely many different blending methods just by con- sidering different parameterizations of rigid transformations.
122 FPS
1
A theoretically optimal rigid transformation blending method has been proposed previously [Govindu 2004]. The algorithm pos- sesses all desired mathematical properties that guarantee correct skinning, but unfortunately it is iterative and thus prohibitively slow for most real-time applications. Therefore, we propose instead a closed-form approximation, based on dual quaternions – a general- ization of regular quaternions first proposed in the nineteenth cen- tury [Clifford 1882].
These concepts can be elegantly illustrated in 2D Euclidean space. Assume we have 3 points p1,p2 and p3 lying on a spherical arc (representing rotation about the origin), such as in Figure 2. Direct averaging of point coordinates (left) produces pavg that no longer lies on the arc. This is the reason for artifacts in linear blend skin- ning – in extreme cases the average can even coincide with the arc’s center. In 2D Euclidean space, this can be easily amended by aver- aging angles corresponding to the points p1,p2 and p3 (see Figure 2 right).
pp2 pp2 p’
1 1avg
pavg
p3 p3
Figure 2: Averaging points (left) versus averaging angles (right).
The principle of our approach is the same as when blending nor- mals, e.g., as in Phong shading (see Figure 3). When blending normals n1 and n2 with weights 0.3 and 0.7, we first blend them linearly, producing vector nb, and subsequently normalize it to the resulting normal nfinal. Even though this is not equivalent to the theoretically perfect intrinsic blending, the approximation is often sufficient (particularly when n1 and n2 are close) and the algorithm is very fast. Dual quaternions permit us to apply the same trick on SE (3).
n1 nfinal nb
n2
Figure 3: Standard normal blending trick: vectors n1 and n2 are first blended linearly giving nb and then re-normalized, giving n f inal .
Moreover, the mathematical properties of dual quaternions ensure that none of the skin collapsing effects (see Section 2.2) exhibited by previous techniques will manifest themselves. Blending of dual quaternions can be elegantly computed in a vertex shader with com- plexity comparable to standard linear blending. Dual quaternions are more memory efficient, requiring only 8 floats per transfor- mation (essentially, two regular quaternions), instead of the 12 re- quired by matrices. In an existing application, it is extremely easy to replace a linear blend skinning implementation by a dual quater- nion one. All that is necessary is a slight modification of the vertex shader and conversion of the matrices to dual quaternions before passing them to the shader. The model files as well as the internal
data structures do not need any change at all – the only difference is in the transformation blending.
This paper extends the work presented in [Kavan et al. 2007]. We have added a more complete dual quaternion tutorial (Appendix A), providing proofs of all important statements and highlighting the connection to spatial kinematics. Addressing practical issues, we propose a more efficient vertex shader (Section 4), discuss quater- nion antipodality issues (crucial for robust implementation, Sec- tion 4.1) and propose a method to integrate scale/shear joint trans- formations (Section 4.2).
Conventions. We denote scalars by lower-case letters, vectors, complex numbers and quaternions in bold and matrices by capital letters. Dual quantities are distinguished from non-dual by a caret; for example aˆ denotes a dual number and qˆ a dual quaternion. The i-th component of vector v is written as vi, thus also v = (v1,…,vn). The dot product of vectors v and w is denoted as ⟨v,w⟩ and ∥v∥ is the usual vector norm. The cross product is denoted as v × w and takes precedence over vector addition, as is usual.
2 Related Work
The first reference to dual quaternions (historically called bi- quaternions) appears in [Clifford 1882]. More recent publications study dual quaternions from the viewpoint of theoretical kinematics [Bottema and Roth 1979; McCarthy 1990]. To date, dual quater- nions have been applied mainly in computer vision and robotics [Daniilidis 1999; Perez and McCarthy 2004], but only rarely in computer graphics [Luciano and Banerjee 2000].
Geometric algebras. Dual quaternions are a special case of the more general concept of geometric algebras. These algebras naturally contain not only vectors and quaternions, but also k- dimensional subspaces [Wareham et al. 2005]. This leads to very elegant and dimension-independent expressions of geometric prop- erties, but sometimes unfortunately also to an increase in time and memory complexity of the resulting implementation [Fontijne and Dorst 2003]. Dual quaternions, in turn, are not so general, but are more compact and faster to manipulate. For computer graph- ics practitioners, the big advantage of dual quaternions is that they are based on regular quaternions – a well-known tool in computer graphics [Shoemake 1985].
Blending vs interpolation. A vast amount of literature has been devoted to the problem of transformation interpolation [Barr et al. 1992; Juttler 1994; Marthinsen 1999; Belta and Kumar 2002; Hofer and Pottmann 2004; Li and Hao 2006; Wang et al. 2008]. This is not surprising, because the construction of interpolation curves for given key transformations (e.g., camera orientations) is a funda- mental problem in computer animation. Unfortunately, in skinning, we face a different problem: the blending of rigid transformations, i.e. their weighted average (confusingly, in some literature this is also called interpolation). The weighted averages can be used to construct interpolation curves (see [Buss and Fillmore 2001]) but not vice versa.
2.1 Skinning
Historically, the idea of skin deformation by an underlying skeleton is credited to [Magnenat-Thalmann et al. 1988]. Since then, sev- eral different approaches to skeletal animation have emerged. Our approach – dual quaternion skinning – falls into the category of ge- ometric methods. For completeness, however, we also survey other related techniques.
Physically based methods. A logical approach to character anima- tion is to simulate the internal structure of the body: bones, muscles
2
and fat tissues. This can work either with explicit anatomy knowl- edge [Scheepers et al. 1997; Aubel and Thalmann 2000; Teran et al. 2005], or without [Capell et al. 2002; Guo and Wong 2005; Pratscher et al. 2005]. Physically based methods generally obtain a high level of realism (delivering also dynamic effects and muscle bulges), but at high computational costs.
Capturing real subjects. Several methods successfully exploit modern motion capture and/or 3D scanning devices to capture skin deformation of real people [Allen et al. 2002; Anguelov et al. 2005; Park and Hodgins 2006; Allen et al. 2006]. These approaches are highly accurate, but require expensive hardware and are of course limited to existing subjects only.
Example based techniques. Multiple input meshes can be used both to resolve the artifacts of linear blending and to add addi- tional effects like muscle bulging. Example based methods use either direct interpolation between example meshes [Lewis et al. 2000; Sloan et al. 2001], approximation by principal components of example deformations [Kry et al. 2002] or fit the linear blending parameters to match the provided examples [Mohr and Gleicher 2003]. A more accurate (yet more complex) example interpola- tion method has been proposed [Kurihara and Miyata 2004] and augmented with an innovative GPU approach [Rhee et al. 2006]. Recently, example based skinning methods have been improved by using rotational instead of linear regression [Wang et al. 2007; We- ber et al. 2007]. Generally, this class of algorithms offers a level of realism limited only by the number of input examples. However, the production of examples can be costly, requiring a lot of artist labour.
Another possible way to overcome the limitations of linear blending with the aid of examples is to use more than one weight per matrix, resulting in a method called Multi-Weight Enveloping [Wang and Phillips 2002]. This idea has recently been refined by Merry et al., whose system is called Animation Space [2006]. The great advan- tage of Animation Space is that it is a linear framework, i.e., that blending is done in a linear space (albeit multi-dimensional) and yet it still significantly outperforms both linear blend skinning as well as multi-weight enveloping in terms of deformation quality [Jacka et al. 2007].
Geometric methods. In this case, only one input mesh is provided (designed in a reference pose). The skeleton-to-skin binding is de- fined in a direct, geometrical way. The most popular method, estab- lished with linear blend skinning, is to bind each vertex to one or more joints. In the latter case, the weight (amount of influence) of all influencing matrices must be specified. This weighting, as well as skeleton fitting, is typically done manually. However, an auto- matic procedure has been described recently [Baran and Popovic ́ 2007], thus simplifying the rigging process considerably.
Advanced blending methods, e.g., direct quaternion blending [Hejl 2004], log-matrix blending [Cordier and Magnenat-Thalmann 2005] and spherical blending [Kavan and Zˇ a ́ra 2005] use the same rigging structure as linear blend skinning. Even though these tech- niques remove some of the artifacts, they still fall short of deliver- ing natural skin deformation in all postures (see Figure 1 and Sec- tion 2.2).
Alternative rigging. Some researchers propose combatting skin- ning artifacts by implementing a different rigging method, for ex- ample based on swept surfaces [Hyun et al. 2005] or auxiliary curved skeletons [Yang et al. 2006; Forstmann and Ohya 2006; Forstmann et al. 2007]. In some cases, this also allows advanced effects to be animated, such as muscle bulging and skin creasing. The disadvantages include complexity of the GPU implementation (even though the latest method from Forstmann et al. [2007] is al- most as fast as linear blend skinning) and inconsistency with the
established rigging pipeline: new rigging tools and data formats are needed. In this paper, we argue that the problems of linear blending do not stem from incorrect or insufficient rigging, but from incor- rect blending. With dual quaternion skinning, it is therefore not necessary to either change the rigging structures or to update exist- ing 3D models.
Another important class of methods is based on generalization of barycentric coordinates [Ju et al. 2005; Joshi et al. 2007]. This enables the user to deform a character using a simpler mesh (defor- mation cage). This approach is very appealing mainly in off-line production, where high quality and direct control of the deforma- tions are more important than run-time efficiency.
2.2 Geometric Skinning
In this section, we elaborate on geometric skinning methods with the rigging structure adopted from linear blend skinning (which is the de facto standard in the videogames industry). A 3D object conforming to this standard consists of skin, a skeleton and ver- tex weights. The skin is a 3D triangular mesh with no assumed topology or connectivity and the skeleton is a rooted tree (both are designed in a reference pose). The nodes of the skeleton repre- sent joints and the edges can be interpreted as bones. However, each bone can be easily identified by its origin, so the difference between joints and bones is rather moot in our case (and in the lit- erature, these terms are often used interchangeably; we will use the term joint). The transformations relating joints in the hierarchy are assumed to be rigid (until Section 4.2, where we propose a method to integrate scale/shear joint transformations). The vertex weights describe the skin-to-skeleton binding, i.e., the amount of influence of individual joints on each vertex.
Let us assume that there are p joints in our model. In the rest- pose, each joint has an associated local coordinate system. The transformation from the rest-pose of joint j ∈ {1, . . . , p} to its ac- tual position in the animated posture can be expressed by a rigid transformation matrix – let us denote this matrix as Cj ∈ SE(3).
We assume that vertex v is attached to joints j1,…, jn with weights w = (w1,…,wn). The indices j1,…, jn are integers referring to the joints that influence a given vertex – in other words, they are in- dices into the array of joints. There is usually a fixed upper bound on n (the number of influencing joints), typically 4, due to graphics hardware considerations. The weights are normally assumed to be convex, i.e., ∑ni=1 wi = 1 and wi ≥ 0. However, this non-negativity is not exploited in our algorithms (analogously to linear blend skin- ning), so artists can feel free to experiment with negative vertex weights. The weight wi represents the influence of joint ji on ver- tex v.
The vertex position in the mesh deformed by linear blend skinning is then computed as
n
v′ = ∑ wiC ji v (1)
i=1
that is, transforming vertex v by all influencing joint transforma- tions Cji and taking a weighted average. This is reminiscent of Figure 2 left, which suggests why artifacts such as the “candy- wrapper” occur with linear blend skinning. To explain this arti- fact, consider a very simple arm rig with only two joints: j1 cor- responding to the shoulder and j2 to the elbow joint (see Figure 4 left). Vertex v in the figure is equally influenced by both joints, i.e, w1 = w2 = 0.5, in order to achieve smooth skinning. Let us fur- ther assume that the arm is animated by twisting joint j2 by 180 degrees around the x-axis. The joint transformations therefore can
3
be written as
Cj1=01,Cj2= 0 1
We assume that the rotational blending is a 2D interpolation of the angle (see Figure 2 right) and that interpolation of the translation vectors is linear. For illustrative purposes only, Figure 5 left shows the resulting transformation Cblend applied to the whole mesh in- stead of individual vertices. As the translational part of both Cj1
I 0 Rx(180◦) 0
where I denotes the 3 × 3 identity matrix and Rx denotes rotation
and C j2
is zero, the blended matrix is simply R(αt) 0
01
about the x-axis. We see that averaging C j1 v and C j2 v produces vertex v′ exactly at the position of joint j2, i.e., the skin collapses to a single point. Examples of this effect with a realistic 3D model are shown in Figure 14 left.
Cblend(t)= z
v Cj1v
Let us examine what happens if we choose a different coordinate system in which to express our transformations, for example that associated with the shoulder joint (see Figure 5 right). If u is the translation between the shoulder and the elbow joint and Tu is the corresponding translation matrix, we can express the transforma- tions with respect to the shoulder joint as
j1
j2
v’ Cj2v
Figure 4: A typical “candy-wrapper” artifact of linear blend skin- ning. On the left is a reference pose and on the right an animated one.
To gain a better insight into linear blend skinning, we can re-write Equation (1) using distributivity of matrix-vector multiplication
I 0 j1 uj1u 01
nn ∑wiCjiv= ∑wiCji v
i=1 i=1
(2)
j2 uj2u 0 1
Note that while the matrix C′j2 differs from Cj2 , they both represent the same transformation (as is obvious from the bottom row of Fig- ure 5). If we now interpolate between C′j1 and C′j2 using the same method as before, we obtain
01
This is unfortunate because Cb′ lend (t) no longer represents the same transformation as Cblend(t) for 0 < t < 1. To see this, compare the transformations of the elbow joint position by Cblend(t) and Cb′ lend (t). While the former leaves the elbow joint fixed, the lat- ter produces the following trajectory
u′(t)=Cb′lend(t) u1 =Rz(αt)u−tRz(α)u+tu
For example, for α = 120◦ and u = (2,0,0) (as in Figure 5) we obtain u′(0.5) = (2.5, √3/2, 0). This represents an unwanted drift
away from the desired elbow position u, as illustrated in Figure 5 right. This has, of course, catastrophical consequences for skinning. In practice, the situation is even worse, as the origin (i.e., the default center of rotation) is usually even further away (typically near the character’s center of mass). See Figure 6 for an example with a 3D character model.
2.3 Center of Rotation Selection
The solution of the problem described in the previous section is straightforward, i.e., it is sufficient to set the rotation center to be in the appropriate joint. This is the basic idea of Hejl’s method [2004], which assumes that the rotation center of vertex v is fixed and coin- cides with the joint nearest to vertex v. For blending the rotational component, Hejl proposes to apply a linear combination of regular quaternions. This trick is of a similar nature to that of normal aver- aging (see Figure 3), but is instead performed on the unit quaternion hypersphere. Even though this is just an approximation of rigorous spherical averages [Buss and Fillmore 2001], it is sufficiently ac- curate for skinning and can be efficiently implemented on graphics hardware.
In cases similar to that of the elbow (two bones connected by a joint), this works perfectly. However, for more complex joint in- fluences, e.g., when more than two influencing joints are involved,
The right-hand side of Equation (2) shows that skinning can also be computed by blending the transformations C ji first and then apply- ing the result to v, thereby obtaining the result in one step.
The transformation blending used in the right-hand side of Equa- tion (2) is a direct linear combination of matrices. It is well known that this method is troublesome, because the blended ma- trix ∑ni=1 wiC ji is not necessarily a rigid transformation, even if all Cji arerigid(i.e.,thesetoforthonormalmatricesisnotclosedunder addition). This is another way of explaining the linear blend skin- ning problems – undesired scaling creeps into our transformation matrix. In the extreme case shown in Figure 4, the blended matrix can even become rank deficient.
An obvious idea to improve skinning quality is to replace the linear blending of matrices in Equation (2) with a more sophisticated ma- trix blending technique. In particular, if the new blending technique will always produce a rigid transformation, we can expect that no candy-wrapper-like artifacts will occur. An ideal way would be to apply one of the SE(3) intrinsic methods [Govindu 2004; Kavan et al. 2006]. This removes the artifacts of linear blend skinning but at a significant cost. The intrinsic methods are iterative and each it- eration requires non-trivial computations (for details please refer to Appendix B.1). As a result, this approach has not become popular in real-time skinning.
Ideally, we would like to achieve high-quality skin deformation but with computational complexity comparable to that of linear blend skinning. The applied matrix blending method does not have to be exact, but should correctly handle both the rotational and transla- tional parts of the transformations Cj1 ,...,Cjn . Handling the trans- lational component is non-trivial, because it depends on the chosen coordinate system (i.e., center of rotation). Let us illustrate what happens if we simply interpolate the translation vectors linearly. For our demonstration, we use the example of a human arm bending at the elbow. We set the shoulder transformation C j1 to the identity, while C j2 corresponds to the elbow bend. With respect to the coor- dinate system of the elbow joint, the transformation matrices have a simple form
I 0 , Cj2 = Rz(α) 0 0101
R (αt) t(u−R (α)u)
C′ =TCT−1=
C′ = T C T−1 = Rz(α) u−Rz(α)u
Cb′lend(t)= z
z
Cj1 =
4
t=0
t = 0.25
t = 0.5
t = 0.75
t=1
Figure 5: Interpolation of rigid transformations with rotation center in the elbow (left) and in the shoulder joint (right). The shoulder is not a suitable rotation center because it leads to an undesirable drift of the interpolated transformation.
selection of the center of rotation is not so simple. With character models, such difficult situations usually occur around the arm-pit or dorsum. To illustrate the problem more clearly, we use the example of a skeletally animated piece of cloth.
Let us assume that we have a fully outstretched piece of cloth influ- enced by two joints, as in Figure 7 left. Vertices close to the bone emanating from j1 will follow transformation C j1 , vertices close to the bone emanating from j2 will follow transformation C j2 and ver- tices in between will blend between these two transformations (to simulate the stretching effect). A problem occurs when switching the center of rotation from j1 to j2. For example, vertices v1 and v2 in the figure are close to each other and therefore they will have similar vertex weights, i.e., approximately 0.5 for both j1 and j2. Therefore, both vertices will be rotated by approximately the same
Figure 6: Artifacts produced by blending rotations with respect to the origin (left) are even worse than those of linear blend skinning (right).
amount (i.e., half of the rotation present in C j1 ). However, vertex v1 is slightly closer to j1 and vertex v2 is slightly closer to j2, which means that v1 will rotate about j1, while v2 will rotate about j2. This causes a problem depicted in Figure 7 right. Since the rota- tion centers of v1 and v2 differ considerably, the deformed vertices v′1,v′2 will be quite far apart. This is unfortunate, as nearby vertices v1,v2 in the reference skin should be mapped to nearby vertices v′1,v′2 in the deformed skin. Violation of this condition manifests itself as cracks, see Figure 15 left.
j1
j2
v1'
v1 v2' v2
Figure 7: Problems with rotating vertices around their nearest joints. On the left is a reference pose and on the right an animated one.
A method to cope with these problems has been proposed by Kavan and Zˇa ́ra [2005]. This approach, called spherical blend skinning, also works by blending translations and rotations (represented by quaternions) independently. However, the center of rotation is se- lected in a more sophisticated way. In particular, the center of rota- tion is not fixed, but is computed at run-time using the actual joint transformations. To summarize, spherical blend skinning defines the center of rotation for vertex v influenced by joints j1 , . . . , jn as the point r which minimizes
∑ ∥Cjar−Cjbr∥ 1≤a
3 Rigid Transformation Blending
From the discussion of previous geometric skinning methods, we see what an ideal rigid transformation blending method for skinning should look like. In particular, it should properly blend the centers of rotation (as with log-matrix blending), but also take advantage of quaternions to blend the rotational parts. A straightforward idea would be to simply blend the rotation centers linearly and couple them with rotations computed using quaternions. Unfortunately, as demonstrated below, this does not work. In this section, we will therefore discuss how the rotation centers should be blended prop- erly.
Before we start, however, we need to prepare a formula expressing the center of rotation from a given rotation and translation (which will be useful not only in the following section, but also in Ap- pendix A.3).
Lemma 1. If a 2D rigid transformation is given by translation t ∈ R2 and angle of rotation α, then its center of rotation r is given as
1 α
r = 2 t+z×t cotan 2 (4)
where z = (0, 0, 1) is the z-axis.
Proof. The formula can be derived elegantly using complex num- bers – let us therefore assume that t is a complex number. The center of rotation being sought, r ∈ C, is the stationary point of our rigid transformation, i.e.,
r = t + eiα r
From this equation we can express the center of rotation as follows:
t 1−e−iα 1−e−iα 1−cosα +isinα r = 1−eiα · 1−e−iα = 2(1−cosα)t = 2(1−cosα) t
Using one of the well-known trigonometric identities sinα =cotanα
and the fact that multiplication by a complex unit is equivalent to a cross product with the z-axis, we obtain Equation (4).
r2
j1
v1′ v’ 2
v1 r v2
1
j3 j2
Figure 8: Problem with rotation centers in spherical blend skinning.
On the left is a reference pose and on the right an animated one.
In this example, we modify the rigging from Figure 7 by adding one more joint, j3. Let us assume v2 is influenced by all three joints (with j3 having only a very small weight, e.g., w3 = 0.01), while v1 is too far from j3 and thus is influenced only by j1 and j2. Even though vertex v2 is only negligibly influenced by joint j3, the algorithm computing rotation centers considers all influencing joints as equal (because the rotation center will be reused at other vertices). Therefore, the center of rotation r1 used to rotate vertex v1 will be different from the center of rotation r2 used to rotate v2, and we run into similar problems as before, though typically not so pronounced (see also Figure 17 left).
Intuitively, it is necessary to compute a proper center of rotation for every vertex individually, taking the joint weights into account. With spherical blend skinning, this is unfortunately not possible in real-time, as it would require execution of the SVD algorithm once per vertex. However, an interesting alternative would be to blend the rotation centers themselves. This idea has been investigated by Alexa [2002], using matrix exponentials and logarithms. The main principle is to linearly combine matrix logarithms instead of the matrices themselves (therefore the method is also known as log- matrix blending). The relationship between matrix logarithms and rotation centers might not be immediately obvious. However, note that the logarithm of a matrix M ∈ SE(3) can be written as (see [Murray et al. 1994])
⎛ 0 −θa θa m ⎞ 321
l o g M = ⎜⎝ θ a 3 0 − θ a 1 m 2 −θa2 θa1 0 m3
0000
⎟⎠
( 3 )
where θ is the angle of rotation, a = (a1,a2,a3), ∥a∥ = 1, is the di- rection of the rotation axis, m = (m1,m2,m3) = θr×a+da, while r is the center of rotation and d is the pitch (ratio of translational and rotational velocity). Note that the center of rotation occurs in the cross product r × a, which simply expresses the fact that any point r + t a is also a valid rotation center (where t ∈ R is arbitrary). All rotation centers therefore form a line in space, known also as the screw axis, see Appendix A.3.
Linear combination of matrix logarithms therefore involves linear combination of rotation centers. Log-matrix blending thus avoids the center of rotation problem inherent in quaternion-based meth- ods, making it an appealing technique for skinning [Cordier and Magnenat-Thalmann 2005]. Unfortunately, log-matrix blending has another shortcoming, first pointed out in [Bloom et al. 2004]. The problem is that the matrix logarithm represents rotations in the
1−cosα 2
6
Note that the same formula applies in 3D, just using an arbitrary axis of rotation instead of the z-axis. This is because the translation component parallel to the rotation axis does not affect the center of rotation (see Appendix A.3).
3.1 2D Case
We will present a simple SE(2) example which will illustrate the issues with blending centers of rotation. Let us assume we have two 2D rigid transformations M1,M2 ∈ SE(2), where M1 is given by angle α1 and translation vector (2cosα1,2sinα1), and M2 by angle α2 and translation vector (2cosα2,2sinα2), see Figure 9.
r1
1
M2
2 M1
r1
r2 2
(-1,0)
1
resulting trajectory
(0,0) (2,0)
r2
(-1,0)
2 1
Figure 10: Trajectory resulting from linear interpolation of rotation centers.
Let us therefore derive the desired non-linear interpolant. First, however, a remark regarding interpolation of the angle: in 2D, it would be natural to simply interpolate the rotation angle lin- early. Unfortunately, the corresponding 3D counterpart would in- volve proper spherical averages [Buss and Fillmore 2001], which is undesirable, because for fast skinning we would prefer linear quaternion blending. Therefore, to make our discussion relevant, we will consider the angle of rotation to be interpolated using the 2D version of linear quaternion blending. This amounts to interpo- lation between (cosα1,sinα1) and (cosα2,sinα2) along a straight line followed by projection back to the spherical arc, as shown in Figure 3. Moreover, to be fully compatible with quaternions, we will work with half of the rotation angle in trigonometric functions (noting that quaternions employ half of the angle of rotation).
Specifically, if we denote pi = (cos αi , sin αi ), the interpolated an-
gle α (t ) will be given by
(0,0)
(2,0)
Figure 9: Interpolation between transformations M1 and M2 ∈ SE(2) with rotation centers r1 and r2. Note that the angle of ro- tation is the same with respect to any center.
First of all, we need to express the rotation centers r1 and r2. In fact, we can easily find the rotation centers for all transformations from the family
⎛ cosα −sinα 2cosα ⎞ π ⎝ sinα cosα 2sinα ⎠, α ∈ 0,
22
α(t) (1−t)p1 +tp2 0012 12
α(t)
cos 2 ,sin 2 = ∥(1−t)p +tp ∥ (7)
If we plug t = (2cosα,2sinα) into Equation (4), we obtain, after some simplifications,
sinα α r=−1,1−cosα=−1,cotan2 (5)
According to Equation (5), all we need, to compute the center of rotation corresponding to α(t), is sinα(t) and cosα(t). This can be retrieved from Equation (7) using the well known identities
2 2 cos α(t) − sin α(t)
The rotation center corresponding to Mi therefore is αi
α(t) α(t)
cosα(t) =
sinα(t) = 2sin 2 cos 2
22
ri =(ri,x, ri,y)= −1, cotan 2 , i=1,2 (6) Now, we can turn our attention to interpolation between M1 and
The actual derivations are somewhat lengthy, therefore we em- ployed Maple [Char et al. 1983] to perform all of the substitutions and simplifications. According to the listing (see Figure 11), the resulting formula for the y-coordinate of the rotation center is given by
M . Obviously, M can be obtained from M by composing M 2211
with origin-centered rotation with angle α2 − α1 . Therefore, the natural way to interpolate between M1 and M2 is along the spherical arc, depicted in Figure 9. An important question is what happens with the center of rotation of the interpolated transformation M(t). Obviously, the rotation center of M(t) will lie on the line segment determined by r1 and r2. However, from Equation (6), we see that the interpolation of rotation centers will not be linear. Let us test how far from the correct solution the linear combination of rota- tion centers is. The resulting trajectory for α1 = 30 and α2 = 60 degrees is shown in Figure 10. We see that it is quite far from the desired spherical arc, and therefore we can conclude that non-linear interpolation of rotation centers is indeed necessary.
(1−t)cos α1 +t cos α2 22
ry(t)= α1 α2 (8) (1−t)sin 2 +tsin 2
(for the x-coordinate we of course have rx(t) = −1). Equation (8) can be re-written to
ry(t) =
α1 α1 α2 α2 (1−t)sin 2 cotan 2 +tsin 2 cotan 2
α1 α2 (1−t)sin 2 +tsin 2
7
> p1x := cos(alpha1/2): p1y := sin(alpha1/2): > p2x := cos(alpha2/2): p2y := sin(alpha2/2):
> ptx := (1-t)*p1x + t*p2x: pty := (1-t)*p1y + t*p2y:
> norm := sqrt(ptx^2 + pty^2):
> cosAlphaThalf := ((1-t)*p1x + t*p2x) / norm:
> sinAlphaThalf := ((1-t)*p1y + t*p2y) / norm:
> cosAlphaT := simplify(cosAlphaThalf^2 – sinAlphaThalf^2);
cosAlphaT := −(2cos(1 α1)2 −4cos(1 α1)2 t+ 22
2cos(1 α1)t cos(1 α2)+2cos(1 α1)2 t2− 222
2cos(1 α1)t2 cos(1 α2)+2t2 cos(1 α2)2 −1+2t− 222
2sin(1 α1)t sin(1 α2)−2t2 +2sin(1 α1)t2 sin(1 α2)) 2222
(−2cos(1 α1)t cos(1 α2)+2cos(1 α1)t2 cos(1 α2)−1+ 2222
2t −2sin(1 α1)t sin(1 α2)−2t2 +2sin(1 α1)t2 sin(1 α2)) 2222
> sinAlphaT := simplify(2*sinAlphaThalf* cosAlphaThalf);
sinAlphaT := −2((−sin(1 α1)+sin(1 α1)t −t sin(1 α2)) 222
(−cos(1 α1)+cos(1 α1)t −t cos(1 α2))) 222
(−2cos(1 α1)t cos(1 α2)+2cos(1 α1)t2 cos(1 α2)−1+ 2222
2t −2sin(1 α1)t sin(1 α2)−2t2 +2sin(1 α1)t2 sin(1 α2)) 2222
>
ryT := simplify(sinAlphaT/(1-cosAlphaT));
ryT := −((−sin(1 α1)+sin(1 α1)t −t sin(1 α2)) 222
(−cos(1 α1)+cos(1 α1)t −t cos(1 α2))) 222
(−1+2t −2sin(1 α1)t sin(1 α2)−2t2+ 22
2sin(1 α1)t2 sin(1 α2)+cos(1 α1)2 −2cos(1 α1)2 t+ 2222
cos(1 α1)2 t2 +t2 cos(1 α2)2) 22
> simplify(denom(ryT)+(-sin(1/2*alpha1)+ sin(1/2*alpha1)*t-t*sin(1/2*alpha2))^2);
0
> ryTfinal := (cos(1/2*alpha1)-cos(1/2*alpha1)*t+ t*cos(1/2*alpha2)) / (sin(1/2*alpha1)- sin(1/2*alpha1)*t+t*sin(1/2*alpha2));
>
cos(1 α1)−cos(1 α1)t +t cos(1 α2) ryTfinal := 2 2 2
sin(1 α1)−sin(1 α1)t +t sin(1 α2) 222
simplify(ryTfinal – ryT);
0
This allows us to conclude that when the angle of rotation is inter- polated using linear quaternion blending, the corresponding center of rotation r(t) is interpolated non-linearly according to the formula
(1−t)sin α1 r1 +t sin α2 r2
r(t)= 2 2 (9)
(1−t)sinα1 +tsinα2 22
This is, in fact, linear interpolation weighted by the factor sin αi . 2
This suggests that we could store the rotation center in the form of
ri sin αi instead of just ri. Actually, if we add a unit quaternion to 2
represent the rotation, we obtain a representation of SE(2) that can be interpolated linearly without introducing discrepancy between the rotation and its center.
Figure 11: Center of rotation derivation in Maple α1 α2
Have we discovered anything new? The answer is no. In fact, this representation of SE(2) was first discovered in the 19th cen- tury [Clifford 1882] and is known today as planar dual quaternions [McCarthy 1990]. In particular, any planar dual quaternion qˆ has the form
qˆ = c o s θ + s i n θ ( k + ε i r y − ε j r x ) ( 1 0 ) 22
where θ is the angle of rotation and (rx,ry) is the center of rotation. Note that i, j,k are the usual (Hamilton’s) quaternion units, while ε is the dual unit, i.e., a number with property ε2 = 0 (the intuition is that ε is so small that its square vanishes completely). In the following, we will assume that the reader is familiar with basic dual quaternion operations and their properties – otherwise please see the enclosed tutorial (Appendix A).
3.2 3D Case
Using dual quaternions, the interpolation derived in Section 3.1 can be written very concisely. We call this method Dual quaternion Linear Blending (DLB) and define it as follows
DLB(t;qˆ1,qˆ2) = (1−t)qˆ1 +tqˆ2 ∥(1−t)qˆ1 +tqˆ2∥
where qˆ1,qˆ2 are unit dual quaternions representing the input trans- formations. Note that if qˆ1,qˆ2 are planar dual quaternions, the cor- responding rotation centers exactly obey Equation (9). However, DLB is defined for general dual quaternions and therefore can be applied to arbitrary 3D rigid transformations.
What is the geometrical interpretation of the 3D DLB? Recall that a unit dual quaternion is a special representation of the screw pa- rameters. In particular, as shown in Appendix A.3, any unit dual quaternion qˆ can be written as
qˆ =cosθ0 +s sinθ0 +ε s sinθ0 −θε sinθ0 +s θε cosθ0 202ε222022
where θ0 is the angle of rotation, s0 is the direction of the axis of rotation, θε is the amount of translation along the rotation axis and sε = r × s0 is the moment of the rotation axis (where r is the center of rotation). Again, as in Equation (3), we see that the center of rotation occurs in the cross product with the direction of the axis of rotation. Note that if s0 = (0, 0, 1), the expression above reduces to Equation (10) (i.e., the 2D case is nothing more than the 3D case restricted to rotations about the z-axis).
The interpretation of the terms occuring in qˆ is as follows. Obvi- ously, the non-dual part, i.e., cos θ0 + s0 sin θ0 is simply the reg-
(1−t)sin 2 r1,y +t sin 2 r2,y
2 2
ular quaternion representing the rotational component. Regarding
= α1 α2 (1−t)sin +t sin
the dual part, we see that the moment sε is again multiplied by sin θ0 , as was the case in 2D (Section 3.1). The term θε sin θ0 is
22222 8
the “padding” which normalises qˆ (it is analogous to the cos θ0 ele- 2
ment of a regular quaternion). This is not surprising because a rigid
transformation has only 6 degrees of freedom, whereas there are 8
elements in a dual quaternion. The last term, i.e., s0 θε represents 2
translation along the axis of rotation. We see that in this term, the
weighting factor happens to be cos θ0 (instead of sin θ0 ). 22
If required, DLB can be computed even without using dual quater-
nions. The principle is to convert the input matrices to the screw pa-
rameters (i.e., the screw axis, the angle of rotation and the amount
of translation), and then linearly blend these parameters multiplied
by the appropriate weighting functions (i.e., either sin θ0 or cos θ0 ). 22
However, the practical value of this approach is questionable, be- cause it is obviously less efficient than blending dual quaternions. In fact, an immediate optimization would be to precompute the trigonometric functions and their products, which would lead to nothing but disguised dual quaternion elements. For readers more interested in the constructive approach to dual quaternions and the related kinematic issues, we refer them to [McCarthy 1990; Bot- tema and Roth 1979].
A big advantage of DLB is that it works for more than two rigid transformations. If these transformations are expressed as unit dual quaternions qˆ1,…,qˆn with convex weights w = (w1,…,wn), the generalized DLB is simply
while left-invariance requires
∀T ∈ SE(3) : TΦ(t;N1,N2) = Φ(t;TN1,TN2)
If Φ satisfies both right- and left-invariance, we say it is bi- invariant. It is obvious that bi-invariance implies coordinate in- variance, requiring
∀T ∈ SE(3) : TΦ(t;N ,N )T−1 = Φ(t;TN T−1,TN T−1) 1212
If the interpolation Φ is defined via linear combination, there is an immediate connection between distributivity and coordinate invari- ance. This is the case of DLB, defined in Section 3.2. It is easy to see that bi-invariance of DLB follows directly from distributivity of dual quaternions, as shown in the following lemma.
Lemma 2. For any unit dual quaternions pˆ,qˆ1,qˆ2 and any inter- polation parameter t ∈ [0, 1], both of the following equations are true
DLB(t;qˆ1pˆ,qˆ2pˆ) = DLB(t;qˆ1,qˆ2)pˆ DLB(t;pˆqˆ1,pˆqˆ2) = pˆDLB(t;qˆ1,qˆ2)
Proof. For proof of the right-invariance it is sufficient to use the right-distributivity of dual quaternions which implies ∥(1−t)qˆ1pˆ + tqˆ2pˆ∥ = ∥((1 − t)qˆ1 + tqˆ2)pˆ∥ = ∥(1 − t)qˆ1 + tqˆ2∥∥pˆ∥ = ∥(1 − t)qˆ1 +tqˆ2∥ (as pˆ is a unit dual quaternion). We can thus write
DLB(w;qˆ1,…,qˆn) =
w qˆ +…+w qˆ 1 1 n n
∥w1qˆ1 +…+wnqˆn∥
(11)
This is very useful for skinning, where we often need to blend more than two joint transformations. However, note that DLB is only an approximation of rigorously defined weighted averages. The per- fectly correct blending algorithm can be obtained by generalizing the methods presented in [Buss and Fillmore 2001], see Section 3.4 and Appendix B.1.
3.3 Distributivity and Coordinate Invariance
In Section 3.1, we have shown that the requirement of correct han- dling of rotation centers in SE(2) leads to planar dual quaternions. This opens a question as to whether this is just a coincidence, or rather a consequence of some more fundamental property. In the following, we argue that the latter is the case and that the crucial property of dual quaternions is called distributivity. In the language of geometry, it corresponds to coordinate invariance.
Let us start by recalling the definitions. Distributivity (of multipli- cation with respect to addition) requires that for any dual quater- nions qˆ1,qˆ2,qˆ3
DLB(t;qˆ1pˆ,qˆ2pˆ)
= (1−t)qˆ1pˆ +tqˆ2pˆ = ∥(1−t)qˆ1pˆ +tqˆ2pˆ∥
= (1−t)qˆ1 +tqˆ2 pˆ = DLB(t;qˆ1,qˆ2)pˆ ∥(1−t)qˆ1 +tqˆ2∥
(qˆ1 +qˆ2)qˆ3 qˆ3(qˆ1 +qˆ2)
= qˆ1qˆ3 +qˆ2qˆ3 = qˆ3qˆ1 +qˆ3qˆ2
D L B ( t ; 1 , qˆ 2 qˆ 1 ) qˆ 1
because 1 corresponds to the identity and qˆ 2 qˆ ∗1 to the relative trans-
Note that because of non-commutativity, we do indeed need both of these equations, sometimes referred to as right and left distribu- tivity.
Coordinate invariance appears in many contexts in computer graph- ics. In this paper, we are concerned with coordinate invariance of transformation blending. Let us denote a general interpolation be- tween two arbitrary transformations N1 , N2 ∈ SE (3) as Φ(t ; N1 , N2 ), where t ∈ [0, 1] is the interpolation parameter. Right-invariance re- quires Φ to satisfy
∀T ∈ SE(3) : Φ(t;N1,N2)T = Φ(t;N1T,N2T)
9
Proof of the left-invariance is a direct analogy of the proof above (using left-distributivity of dual quaternions).
It may seem that bi-invariance and distributivity are just theoretical properties without any practical implications. However, this is not the case, as we now illustrate on the example from Figure 9. Let us assume that unit planar dual quaternions qˆ1,qˆ2 correspond to our transformations M1,M2 ∈ SE(2). As discussed in Section 3.1, the natural interpolation between M1 and M2 is to consider the relative motion between M1 and M2, interpolate it, then compose the re- sult with M1. In our example, that relative transformation was the origin-centered rotation with angle α2 −α1. In the language of dual quaternions, this can be written as
∗
formation (i.e., in our case, pure rotation). However, thanks to bi- invariance, we have
DLB(t;1,qˆ2qˆ∗1)qˆ1 = DLB(t;qˆ1,qˆ2) (12)
We have re-derived the result from Section 3.1, i.e., that linear blending of dual quaternions represents a plausible way to interpo- late between M1 and M2. Another interpretation of Equation (12) is that we do not have to evaluate the relative motion (i.e., qˆ2qˆ∗1) ex- plicitly, because linear blending of dual quaternions automatically blends only that relative component.
Bi-invariance can also explain the problems with linear blending of (quaternion,translation) pairs encountered in [Hejl 2004; Ka- van and Zˇa ́ra 2005]. To see this, let us formally define an al- gebra over (quaternion,translation) pairs with multiplication ⊗
corresponding to composition of transformations, i.e., for two (quaternion,translation) pairs (q0,t0), (q1,t1), this is defined as
∗ (q0,t0)⊗(q1,t1) = (q0q1,t0 +q0t1q0)
interpreting the 3D vector t1 as a quaternion with zero scalar part (as is usual). The right distributivity would require that
(q0 +q1,t0 +t1)⊗(q2,t2) = (q0,t0)⊗(q2,t2)+(q1,t1)⊗(q2,t2) Expanding the right hand side yields
(q0,t0)⊗(q2,t2)+(q1,t1)⊗(q2,t2) = ∗ ∗ = (q0q2 +q1q2,t0 +t1 +q0t2q0 +q1t2q1)
while the left hand side expands to (q0 +q1,t0 +t1)⊗(q2,t2) =
= (q0q2 +q1q2,t0 +t1 +(q0 +q1)t2(q0 +q1)∗)
It can be shown that the term (q0 +q1)t2(q0 +q1)∗ is not equivalent to q0t2q∗0 +q1t2q∗1, which precludes right distributivity. This sim- ply reflects the fact that blending (quaternion,translation) pairs ro- tates about the origin. Note that this might be advantageous in some settings (e.g., in rigid body physics, where we can define the origin to coincide with the center of mass). However, as demonstrated in Section 2.2, this is a complication in skinning.
As a final remark, we note that both linear blending of matrices and linear blending of matrix logarithms are coordinate invariant. The former follows from the distributivity of matrix multiplication, while the latter follows from the properties of the matrix exponen- tial and logarithm (see [Moakher 2002])
∀N,T ∈SE(3): exp(TNT−1)=Texp(N)T−1 ∀N,T ∈ SE(3) : log(TNT−1) = T log(N)T−1
This explains why we did not encounter any rotation-center issues with linear blend skinning or with log-matrix blending.
3.4 Accuracy of Dual Quaternion Linear Blending
As argued in Section 3.3, DLB is a plausible method to interpolate rigid transformations. However, it is not perfect, as it is not a group- intrinsic method (i.e., it involves the “normal-interpolation” trick shown in Figure 3). In this section we discuss whether this will introduce artifacts when employing DLB in skinning. For clarity, we will start with the case of two transformations. The first step is to establish the perfectly correct blending method, i.e., one that respects the geometry of the underlying group (in our case SE(3)).
In the case of SO(3), the theoretically perfect solution is Spheri- cal Linear Interpolation (SLERP) [Shoemake 1985]. Recall that for two unit quaternions q1, q2 with parameter t ∈ [0, 1], the formula is SLERP(t;q1,q2)=(q2q∗1)tq1 (assumingthat⟨q1,q2⟩≥0;thiscan always be enforced by negating one of the quaternions). With the aid of dual quaternions, SLERP can be easily generalized to SE(3). We call the resulting method Screw Linear Interpolation (ScLERP), for reasons that will soon become clear. For any two unit dual quaternions qˆ1,qˆ2, it is given as ScLERP(t;qˆ1,qˆ2) = (qˆ2qˆ∗1)tqˆ1. What is its geometric interpretation?
can immediately observe two important properties: the axis nˆ of the screw motion is constant (independent of t), and the angle of rota- tion t α0 , as well as the amount of translation t αε , vary linearly with respect to the interpolation parameter t. This means that ScLERP is a constant speed and shortest path interpolation, as could have been expected because it is a generalization of SLERP. As discussed in Section 3.3, another crucial property we expect from ScLERP is coordinate invariance, which we prove in the following lemma.
Lemma3. ScLERPisbi-invariant,thatisforanyunitdualquater- nions pˆ,qˆ1,qˆ2 and any interpolation parameter t ∈ [0,1], both of the following equations are true
ScLERP(t;qˆ pˆ,qˆ pˆ)=ScLERP(t;qˆ ,qˆ )pˆ 12 12
ScLERP(t;pˆqˆ1,pˆqˆ2) = pˆScLERP(t;qˆ1,qˆ2) Proof. Theright-invarianceiseasytoshow,because
ScLERP(t;qˆ1pˆ,qˆ2pˆ)
= (qˆ2pˆpˆ∗qˆ∗1)tqˆ1pˆ = (qˆ2qˆ∗1)tqˆ1pˆ = = S c L E R P ( t ; qˆ 1 , qˆ 2 ) pˆ
Proving the left-invariance is a little more tricky: ScLERP(t;pˆqˆ1,pˆqˆ2) = (pˆqˆ2qˆ∗1pˆ∗)tpˆqˆ1. It is now sufficient to show that (pˆqˆ2qˆ∗1pˆ∗)t = pˆ(qˆ2qˆ∗1)tpˆ∗, because this gives u s ( pˆ qˆ 2 qˆ ∗1 pˆ ∗ ) t pˆ qˆ 1 = pˆ ( qˆ 2 qˆ ∗1 ) t pˆ ∗ pˆ qˆ 1 = pˆ S c L E R P ( t ; qˆ 1 , qˆ 2 ) . H o w e v e r , t h e p o w e r c a n b e w r i t t e n a s ( pˆ qˆ 2 qˆ ∗1 pˆ ∗ ) t = exp(t log(pˆqˆ2qˆ∗1pˆ∗)). Thanks to Lemma 14, we can derive
exp(t log(pˆqˆ2qˆ∗1pˆ∗)) which concludes the proof.
21 21
= exp(tpˆ log(qˆ2qˆ∗1)pˆ∗) =
= pˆ exp(t log(qˆ qˆ∗))pˆ∗ = pˆ(qˆ qˆ∗)tpˆ∗
We therefore see that ScLERP is an interpolation with the same behavior as SLERP, but generalized to SE(3). Now, we can use ScLERP as the gold standard to compare DLB against. As shown in Lemma 2 and Lemma 3, both methods are bi-invariant. Note that from this follows the most important feature, i.e., that the error will not be dependent on the choice of the coordinate systems (otherwise the error could be unbounded). In fact, we can exploit this common property of ScLERP and DLB to simplify their comparison.
Specifically, instead of comparing DLB(t;qˆ1,qˆ2) directly with
ScLERP(t;qˆ1,qˆ2), we rewrite them as DLB(t;1,qˆ2qˆ∗)qˆ1 and ∗1
ScLERP(t;1,qˆ2qˆ1)qˆ1, which is correct because of right-invariance. Since qˆ 1 is the same in both expressions, it is sufficient to compare just DLB(t;1,qˆ2qˆ∗1) with ScLERP(t;1,qˆ2qˆ∗1), which is an easier problem. As qˆ2qˆ∗1 is a unit dual quaternion, it can be written as
qˆ2qˆ∗ = cos αˆ +nˆ sin αˆ . This enables us to derive 122
1−t +tqˆ qˆ∗ 1−t +t cos(αˆ )+nˆt sin(αˆ ) DLB(t;1,qˆ2qˆ∗1)= 21= 2 2
∥1−t +tqˆ2qˆ∗1∥ ScLERP(t;1,qˆ2qˆ∗1) = (qˆ2qˆ∗1)t = cos
αˆ t 2
∥1−t +tqˆ2qˆ∗1∥ αˆ
+nˆ sin t 2
from which we see that both DLB and ScLERP use the same, con-
We can see that qˆ 2 qˆ ∗1 is a unit dual quaternion, which represents the relative motion between qˆ and qˆ . The power can be written as
stant screw axis nˆ .
Therefore, the one difference between DLB and ScLERP can only be in the motion along the screw axis, i.e., in the angle of rota- tion and amount of translation. Since DLB(t;1,qˆ2qˆ∗1) is a unit dual quaternion, we can also write it in the form
∗t αˆ 1αˆ 2
(qˆ2qˆ1) = cos(t 2 )+nˆ sin(t 2 ) for some dual angle αˆ and dual vec-
tor nˆ (see Appendix A.4). The dual vector nˆ represents the axis of the screw motion and the dual angle t αˆ = t α0 + ε t αε contains both
ˆˆ
222
the angle of rotation (t α0 ) and the amount of translation (t αε ). We
∗ βt βt DLB(t;1,qˆ2qˆ1)=cos 2 +nˆsin 2
10
By considering only the scalar part of this equation, we see that ˆ αˆ
according to Equation (18). Therefore,
∥ bˆ ∥ ∥ bˆ ∥ ∥ b 0 ∥
bˆ ′0 So, the rotational part of bˆ ′ is b′ = b0
cosβt =1−t+tcos(2) 2 ∥1−t+tqˆ2qˆ∗∥
1
(13)
bˆ 1 b b b⟨b,b⟩ bˆ′= =(b0+εbε) =0+ε ε−00ε
∥ b 0 ∥ 3
∥ b 0 ∥
and the translation is given
It is possible to compute an upper bound on the difference between DLB and ScLERP by expressing the dual angle βˆt from Equa-
tion (13) and comparing it with ScLERP’s dual angle αˆ t. This is not difficult but requires a lengthy mathematical analysis: we therefore employed Maple to carry out the computations (see Appendix B). The result is that the angles of rotation in DLB and ScLERP always differ by less than 8.15 degrees (which is in accordance with the results reported in [Kavan and Zˇa ́ra 2005] for the case of regular quaternions). The amount of translation always differs by less than 15.1% of the translation present in qˆ2qˆ∗1. Note that these results are upper bounds – in practice, the difference is usually much smaller. Therefore, in applications such as skinning, the difference between DLB and ScLERP will most likely not be visible at all, although this would need to be verified with a perceptual study.
The same is true for n > 2, even though the error analysis in this more general case is not so simple. The problem is that the gener- alization of SLERP for n > 2 leads to spherical averages proposed in [Buss and Fillmore 2001]. Therefore, the gold standard of SE(3) blending for more than two transformations involves generalizing Buss and Filmore’s algorithms to unit dual quaternions (which is non-trivial because the set of unit dual quaternions is not a hyper- sphere). We discuss this in more detail in Appendix B.1.
4 Implementation Notes
Equation (11) is the key to fast and plausible skinning. In this sec- tion, we discuss the implementation issues on a typical heteroge- neous architecture consisting of a CPU and GPU. The first step is to convert the skinning matrices C1,…,Cp (where p is the total num- ber of joints) to dual quaternions qˆ 1 , . . . , qˆ p , unless our application works with them already. This will typically be done on a CPU and does not take long, because the conversion to a dual quaternion in- volves just one quaternion multiplication and the number of joints p is usually quite small. Let us recall a concise multiplication for- mula for regular quaternions. If two regular quaternions q1 , q2 are written using their scalar and vector parts, i.e., q1 = a1 + r1 and q2 = a2 + r2 , then their multiplication can be expressed as
q1q2 = a1a2 −⟨r1,r2⟩+a1r2 +a2r1 +r1 ×r2 (14) The resulting dual quaternions qˆ1,…,qˆp are then sent to the GPU
as uniform parameters, each represented by a 2 × 4 matrix.
The skin deformation, i.e., the DLB itself and the vertex and nor- mal transformations take place in the vertex shader. In contrast to our previous approach [Kavan et al. 2007], we do not convert the resulting unit dual quaternion to an homogeneous matrix, but ap- ply it directly to transform each vertex and normal (which results in shorter vertex shader code). The first step of DLB(w; qˆ j1 , . . . , qˆ jn ), i.e., the computation of the linear combination bˆ = ∑ni=1 wiqˆ ji , is straightforward, as it consists only of a per-component linear com- bination. However, the subsequent normalization, i.e., the compu- tation of bˆ ′ = bˆ /∥bˆ ∥ can be optimized as follows.
If the non-dual and dual parts of bˆ are b0 and bε , then the norm is ∥bˆ ∥ = ∥b0 ∥ + ε ⟨b0 ,bε ⟩ , according to Equation (22). The inverse is
bˆ ′ε
0 ∥b0∥
by the vector part of 2b′ε (b′0 )∗ . The latter expands to
b b ⟨b ,b ⟩ b∗ bεb∗ ⟨b ,b ⟩ 2b′ε(b′0)∗=2 ε−00ε 0=2 0−0ε ∥b0∥ ∥b0∥3 ∥b0∥ ∥b0∥2 ∥b0∥2
Since the scalar part of 2b′ε (b′0 )∗ = 0 (because bˆ ′ is unit), it means
no need to evaluate ⟨b0 , bε ⟩/∥b0 ∥2 , because its pur- pose is only to cancel out the scalar part of bε b∗0 /∥b0 ∥2 . Therefore, we can compute the translational part of the matrix M just by com- puting the vector part of 2bεb∗0/∥b0∥2 (and its scalar part can be safely ignored). This means that instead of computing the full dual quaternion normalization, all we need to compute is c0 = b0/∥b0∥ and cε = bε /∥b0 ∥ and retrieve the translation as the vector part of 2cε c∗0, which can be done efficiently using Equation (14).
The regular quaternion c0 is used to rotate the input vertex v and the normal vn. The key to a fast GPU implementation is the following classical formula [Shreiner et al. 2007], which describes how to efficiently express quaternion rotation in terms of cross products.
Lemma4. Letq=a+rbeaunitregularquaternionwithscalar part a and vector part r. Rotation of a vector (v0 , v1 , v2 ) repre- sentedbytheregularquaternionv=v0i+v1j+v2kcanbecom- puted as
v′ = v+2r×(r×v+av) (15) where v′ is the vector v rotated by q.
Proof. The proof consists of re-arranging the well-known expres- sion v′ = qvq∗. Using Equation (14), we can expand
that there is
v′
= qvq∗ = (a+r)v(a−r) = (−⟨r,v⟩+av+r×v)(a−r) = = −a⟨r,v⟩+a⟨v,r⟩+⟨r,v⟩r+a2v+a(r×v)−a(v×r)−
(r×v)×r =
= ⟨r,v⟩r+a2v+2a(r×v)+r×(r×v)
given by
∥b0 ∥
1 = 1 −ε⟨b0,bε⟩ ∥ bˆ ∥ ∥ b 0 ∥ ∥ b 0 ∥ 3
Recall Lagrange’s formula
r×(r×v) = r⟨r,v⟩−v⟨r,r⟩ which, added to the previous equation, results in
v′ = ⟨r,v⟩r+a2v+2a(r×v)+2r×(r×v)−r⟨r,v⟩+v⟨r,r⟩ = = v(a2 +∥r∥2)+2a(r×v)+2r×(r×v)
from which Equation (15) readily follows.
Note that the shader compiler could, in theory, perform the opti- mizations of qvq∗ itself. However, according to our experiments, shader compilers produce suboptimal code when compared to im- plementation of Equation (15) (which translates to efficient code easily as modern GPUs provide fast cross product operations). An optimized dual quaternion skinning implementation can thus be summarized as Algorithm 1.
The resulting algorithm is encouragingly simple. However, ac- tual vertex shader code would add some lighting computations and transformation by the model-view-projection matrix.
11
Algorithm 1
Input: dual quaternions qˆ 1 , . . . , qˆ p (uniform parameters) vertex position v and normal vn
joints indices j1,…, jn and weights w1,…,wn Output: transformed vertex position v′ and normal v′n
bˆ = w 1 qˆ j 1 + . . . + w n qˆ j n
// denote the non-dual part of bˆ as b0 and the dual one as bε
c0 = b0/∥b0∥
cε =bε/∥b0∥
// denote the scalar part of c0 as a0 and its vector part as d0
//′ denote the scalar part of cε as aε and its vector part as dε
v′ = v + 2d0 × (d0 × v + a0 v) + 2(a0 dε − aε d0 + d0 × dε )
vn = vn +2d0 ×(d0 ×vn +a0vn)
// note that v′n must be transformed by the inverse transpose ma- trix
4.1 Coping with Antipodality
Conversion of homogeneous matrices to dual quaternions is straightforward, up to the point where an appropriate sign for the resulting dual quaternion must be chosen. Due to dual quaternion antipodality, both qˆ i and −qˆ i represent the same rigid transforma- tion, but their interpolation can be different (see Figure 22). The problem is exactly the same as for regular quaternions and occurs with all previous quaternion-based methods [Hejl 2004; Kavan and Zˇa ́ra 2005]. In the following, we restrict ourselves to the case of regular quaternions (as generalization to dual quaternions is triv- ial).
In the case of only two regular quaternions q j1 , q j2 (i.e., in the case of classical SLERP [Shoemake 1985]), the situation is simple: we take either q j2 or −q j2 , whichever gives a non-negative dot product with q j1 . In the general case (n > 2), we can by analogy require a non-negative dot product for all pairs of the resulting quaternions. More formally, we are looking for n numbers s1,…,sn ∈ {0,1} such that the following holds
∀h,i ∈ {1,…,n} : ⟨(−1)sh qjh ,(−1)si qji ⟩ ≥ 0 (16)
A sequence of s1 , . . . , sn will be called valid if it satisfies Equa- tion (16). With dual quaternions, the only difference is that the dot product in Equation (16) is taken only between the non-dual parts of qˆ jh , qˆ ji (to keep the notation simpler, we thus discuss only the
However, as we show below, the algorithm actually always finds a valid sequence s1,…,sn, if one exists.
Lemma 5. If there exists a valid sequence s1,…,sn ∈ {0,1} for regular quaternions q j1 , . . . , q jn , then Algorithm 2 finds one of the valid sequences.
Proof. First of all, note that if s′1,…,s′n is valid, then so is 1− s′1 , . . . , 1 − s′n , as negating both quaternions does not change the sign of the dot product in Equation (16). Therefore, if a valid se- quence exists then we can assume that the sequence may be given as r1,…,rn ∈ {0,1} such that r1 = 0. We show that our algorithm will return si = ri for i = 1,…,n. This is obviously true for i = 1, ass1=r1=0.Fori>1,weknowthat
⟨q j1 , (−1)ri q ji ⟩ ≥ 0
because the sequence r1,…,rn is valid. The case of ri = 0 means that ⟨qj1,qji⟩≥0 and therefore our algorithm sets also si =0. In the case of ri = 1, we have ⟨q j1 , q ji ⟩ < 0, so the ’else’ branch in the if-statement is taken and the algorithm sets si = 1.
In extreme situations, it is possible that no valid sequence exists. This means that for any choice of signs si, there will always exist a pair of rotations that will be interpolated along the longer path, see Figure 22. Luckily, such situations almost never occur in skinning.
Even though our antipodality resolution algorithm is quite simple and fast, it needs to be executed in the vertex shader, because it depends on the joint set j1,..., jn influencing a particular vertex. If a shorter vertex shader code is required, it is possible to resort to the following approximate method, which does not guarantee a valid sequence even if one exists. However, it allows the signs to be pre-computed before sending dual quaternions to the GPU, so the vertex shader will not have to care about antipodality. The idea is simple: we assume that each joint in the actual skeleton posture has rotated by less than 180 degrees with respect to its parent (i.e., has taken the shorter path). This assumption is practically always valid with character models, as joint rotations larger than 180 degrees are not normally possible. Therefore, we can set s1 = 0 and select the remaining signs s2,...,sn so that
∀i ∈ {2,...,n} : ⟨(−1)sπ(i) qπ(i),(−1)si qi⟩ ≥ 0
where π(i) denotes the parent joint of joint i. An algorithm to find these signs is a simple modification of Algorithm 2. Even though, theoretically, this algorithm does not necessarily produce valid signs for every vertex, we have not encountered any problems with our test data.
4.2 Non-rigid Joint Transformations
Dual quaternions cannot represent non-rigid transformations, such as scale and shear. This means that dual quaternion skinning, unlike linear blend skinning, is restricted only to rotating and/or translat- ing joints. This does not seem prohibitive at first, given that any physical joints can also only rotate and/or translate. However, vir- tual characters can benefit from allowing non-rigid transformations. For example, non-uniform scaling can be applied to vary a charac- ter’s proportions without having to modify the 3D model itself – a useful feature, for example in crowd animation. Some artists also use non-uniformly scaling joints to imitate muscle bulging and sim- ilar effects (which can even be done automatically as discussed by Mohr and Gleicher [2003]).
A natural way to support non-rigid transformations would be via higher-dimensional geometric algebras, which generalize dual
regular quaternion case). A very simple algorithm follows.
Algorithm 2
Input: Regular quaternions qj1,...,qjn
Output: Sequence s1,...,sn ∈ {0,1}
s1 = 0
for i = 2 to n do
if ⟨q j1 , q ji ⟩ ≥ 0 then si = 0
else
si = 1 end if
end for
to find
s1 , . . . , sn
Note that Algorithm 2 is already commonly known, for example in the gaming community. In the light of the more sophisticated techniques of [Park et al. 2002; Johnson 2003], we initially con- jectured that this algorithm provides just an approximate solution.
12
first phase second phase
Figure 12: Adding scale transformations: in the first phase, we re-scale the mesh in the rest-pose. In the second phase, rigid joint transforma- tions are applied using dual quaternion skinning.
quaternions. For example, 5D conformal geometric algebra [Ware- ham et al. 2005] supports rigid transformations along with a dila- tion (uniform scale). Unfortunately, to support the whole range of non-rigid transformations, including non-uniform scale and shear, one would have to move to even higher dimensional geometric alge- bras, suggesting a great cost (as an n-dimensional geometric algebra requires 2n coordinates). Even though this approach would be the- oretically interesting, for practical purposes we propose a simpler way of incorporating non-rigid transformations.
The idea is to separate the joint transformations into non-rigid and rigid parts and perform skinning in two phases. In the first phase, we deal with the non-rigid portion of the transformations, inflating the rest-pose mesh as required. In the second phase, we apply the rigid part to bend the mesh to the final shape. As no rotations are involved in the first pass, we can safely apply linear blending of transformation matrices [Shoemake and Duff 1992]. As only rigid transformations are involved in the second pass, we can apply dual quaternion skinning as described in Section 4. Composition of both steps yields the desired result, see Figure 12. Note that Kurihara and Nishita [2007] proposed a similar technique independently.
To decompose each skinning matrix Ci , i = 1, . . . , p into a non-rigid and rigid part, we could apply polar decomposition [Shoemake and Duff 1992]. However, a more efficient solution is possible when considering how the matrices Ci are formed from individual joint transformations. Note that this would not work if no skeleton is present [James and Twigg 2005]. However, most character anima- tion systems use skeletal animation. In the following we therefore assume we have a skeleton. Let us denote the transformation from the model coordinate system (frame) to joint i’s frame in the rest- pose as Ai (absolute matrix) and the transformation from the model frame to joint i in the animated skeleton as Fi (final matrix). Then the skinning transformation Ci (composed matrix) can be written as Ci = Fi(Ai)−1.
The absolute matrix Ai is simply a concatenation of individual joint transformations, i.e.,
method works in this case as well, we did not find it practical, be- cause the scale/shear transformation Si affects all descendants of joint i. For example, when scaling a spine joint (e.g., to obtain a larger belly), vertices of the arms get scaled in the same direc- tion, which is typically not desirable. Therefore, we propose the scale/shear transformation Si to be considered local, i.e., to affect only joint i’s transformation. In our opinion, this results in more intuitive editing, as one usually wants to enlarge/thicken individual joints rather than the whole sub-tree.
According to the two-phase skinning paradigm introduced above,
we define different matrices for each phase: A′i,Fi′,Ci′ for phase
one, and A′′,F′′,C′′ for phase two. Concerning the first phase, ′ iii
Ai=R1...Rπ(i)Riasbefore,but
F i ′ = R ′1 . . . R ′π ( i ) R ′i S i
where R′k , k = 1, . . . , i, is matrix Rk with its translational part scaled by Sπ(k) (rotational part being kept intact). This is to account for bone elongations caused by the non-rigid transformation present in the parent joint (note that R′1 = R1). The composed matrix for the first phase, Ci′ = Fi′(A′i)−1, thus performs local scaling of joint i in the rest-pose (see Figure 12).
Concerning the second phase, the absolute matrix of joint i is A′′ = R′ ...R′ R′ = F′S−1
Ai =R1...Rπ(i)Ri
where π(i) denotes the parent joint of joint i and R is a “relative
It is simple to modify the vertex shader code from Section 4 in order to implement this method. Matrices C′′,i = 1,..., p are converted
i
matrix” describing the transformation from joint π(i) to i in the rest-pose skeleton. Matrix R1 describes the root joint transforma- tion with respect to the model coordinate system (often R1 = I, as the root is usually placed at the origin). Note that matrices Ri are
i
to dual quaternions qˆi, which are passed to the GPU along with
assumed to be rigid.
The final matrix Fi is formed in a similar way, but accounting for
the joint transformations Ti ∈ SE(3)
Fi = R1T1 ...Rπ(i)Tπ(i)RiTi
In order to enable non-rigid joint transformations, we could simply multiply Ti with a scale/shear matrix Si. Even though our proposed
matrices Ci′ , i = 1, . . . , p. In the vertex shader, for a given vertex v associated with joints j1,..., jn with weights w1,...,wn, we first transform the vertex by ∑ni=1 w ji C′ji and then by the blended dual quaternion as in Section 4. We deal with vertex normals in a similar way, just applying the inverse transpose of ∑ni=1 w ji C′ji . The results of this method are shown in Figure 13, while the performance issues are discussed in the next section.
5 Results and Comparison
In our experiments, we use a human model with 5002 vertices, 9253 triangles and 54 joints. First, we consider only rigid joint
13
i 1 π(i)i ii
i.e., the final matrix of the first phase without the scale/shear factor Si. The final matrix of the second phase deals with the joint rigid transformations Ti
F′′=R′T ...R′ T R′T i 1 1 π(i) π(i) i i
and thus differs from the original final matrix Fi just by employing
the elongated bone transformations R′i. Therefore, the composed
matrix of the second phase, C′′ = F′′(A′′)−1, is rigid. iii
Figure 13: Proportions of the original model (left) are changed by applying non-uniform scaling to the character’s arms (right). Our technique allows scaling transformations to be combined with dual quaternion skinning (top), thus achieving more realistic deforma- tions than with linear blend skinning (bottom).
Figure 14: Comparison of linear (left) and dual quaternion (right) blending. Dual quaternions preserve rigidity of input transforma- tions and therefore avoid skin collapsing artifacts.
transformations and compare the proposed dual quaternion skin- ning with linear blending, direct quaternion blending [Hejl 2004], log-matrix blending [Cordier and Magnenat-Thalmann 2005] and spherical blend skinning [Kavan and Zˇa ́ra 2005]. As discussed in Section 2.2, some artifacts are better visualized on a simple model of cloth (6000 vertices, 12000 triangles and 49 joints). Note that the only variable in our experiments is the transformation blend- ing – the input data (model files and postures) are always the same. The visual results confirm that our DLB method is indeed free of all the artifacts exhibited by previous methods (see Figures 14, 15, 16 and 17).
In order to compare computational performance, we have imple- mented both CPU and GPU versions of dual quaternion skinning.
Figure 15: Comparison of direct quaternion blending (left) and dual quaternion (right) blending. Only the latter delivers smooth defor- mation.
The average performance of the CPU implementations is reported in Figure 18 and the number of instructions of our vertex shaders can be found in Figure 19. Note that our implementation of log-matrix blending uses an optimization for rigid transformations based on the Rodrigues formula, as suggested in [Alexa 2002]. Our vertex shader assumes n = 4 (which is common due to graphics hardware considerations) and does not perform any optimizations if there are fewer than 4 influencing joints.
From the measurements, we see that dual quaternion, linear and direct quaternion blending [Hejl 2004] have quite similar perfor- mance both on the GPU and CPU. Although our algorithm is slightly slower than both linear and direct quaternion blending, we believe that this is not a high price to pay for the elimination of arti- facts. When compared to log-matrix and spherical blending, we see that dual quaternion skinning is more than twice as fast (and also much easier to implement).
Figure 16: Comparison of log-matrix (left) and dual quaternion (right) blending. The shortest-path property of dual quaternion blending guarantees natural skin deformations.
14
90 degrees 179 degrees 181 degrees 270 degrees
Figure 20: Example of a skin flipping artifact caused by the shortest path property. When rotating from 179 to 181 degrees, skin discontinu- ously changes its shape, because the other rotation direction becomes shorter.
Figure 17: Comparison of spherical (left) and dual quaternion (right) blending. Dual quaternions do not need to cluster vertices and therefore naturally avoid artifacts.
12 [ms] 8
4 0
Figure 18: Average CPU performance for skin deformation of the human model in milliseconds (Pentium 4, 3.4 GHz): LBS – linear blend skinning, QB – direct quaternion blending, Log – log-matrix blending, SBS – spherical blend skinning, DLB – dual quaternion linear blending (pre-computed antipodality).
Linear blend skinning supports non-rigid transformations at no ex- tra cost, which is a big advantage. Previous techniques, i.e., direct quaternion blending, log-matrix blending and spherical blend skin- ning, considered only rigid transformations. Our method to add scale/shear transformations presented in Section 4.2 requires an ex- tra set of 3 × 4 matrices to be passed to the GPU (along with dual quaternions), which means that we need in total 20 scalars per joint
60 [# of instructions] 50
40
30
20 10 0
Figure 19: Number of vertex shader instructions: LBS – linear blend skinning, QB – direct quaternion blending (pre-computed an- tipodality), DLB – dual quaternion linear blending (pre-computed antipodality), QB′ and DLB′ – same as before but with antipodality resolved in the vertex shader.
(compare with the 12 required by linear blend skinning and the 8 required by standard dual quaternion skinning).
The two-phase skinning also has an impact on the run-time perfor- mance, requiring an extra 29 shader instructions to blend matrices C1′ ,...,Cn′ (for n = 4). This means that we end up with about twice as many instructions as required for linear blend skinning. This may or may not be an issue, depending on where the bottlenecks of our application are and what the target hardware is.
The comparison of dual quaternion skinning to the Animation Space method [Merry et al. 2006] is somewhat troublesome. Ac- cording to our discussion with the author of Animation Space [Merry 2007], it is likely that no objective experimental comparison is possible, because the modelling methods are different (in partic- ular, Animation Space makes use of additional examples). There- fore, it would be hard to judge if a particular result was due to the properties of the algorithm itself or simply how much effort went into the modelling.
6 Conclusions and Future Work
In this paper, we propose a novel skinning method based on the blending of dual quaternions. Our method is efficient, simple to implement and does not require the modification of existing mod- els or authoring tools (advanced methods exploiting the linearity of linear blend skinning, such as [Mohr and Gleicher 2003], could be adapted). We therefore believe that it provides a highly practi- cal alternative to the popular but inaccurate linear blend skinning method.
However, the proposed algorithm has some limitations. When non- rigid joint transformations are required, our method requires longer shader code and more memory than linear blend skinning. Another
15
52
DLB'
42
QB'
43
DLB
36
LBS
33
QB
11.78
Log
10.83
SBS
4.66 QB
5.07 DLB
3.67 LBS
potential shortcoming of dual quaternion skinning could be a “flip- ping artifact”, which occurs with joint rotations of more than 180 degrees, see Figure 20. This is a corollary of the shortest path prop- erty: when the other path becomes shorter, the skin changes its shape discontinuously. An alternative would be to support multi- ple revolutions and count the number of twists. Even though this would produce continuous deformations, it is doubtful whether the final animation would look more natural. In most situations, how- ever, this is not an issue as extreme rotations leading to flipping artifacts are usually prevented via joint constraints.
Taking a broader perspective, blending of dual quaternions has po- tentially many more applications, not limited to skinning. More generally, it is a method for rigid transformation blending with in- teresting properties. Therefore, it could potentially be useful in contexts such as motion blending, analysis or compression [Alexa 2002]. The investigation of other applications of dual quaternions in computer graphics promises to be an interesting area for future work.
6.1 Implementation
In order to facilitate further research and experiments, we have provided some on-line resources: http://isg.cs.tcd.ie/projects/DualQuaternions/ Specifically, C source code converting between dual quaternions and (quaternion,translation) pairs is available, along with Cg vertex shaders that implement our skinning methods.
7 Acknowledgements
We wish to thank the anonymous reviewers for their valuable com- ments, Bruce Merry for discussion on Animation Space and Carlo H. Se ́quin for early insights into the topic. We acknowledge the support of the Higher Education Authority of Ireland and Science Foundation Ireland. This work has been partly supported by the Ministry of Education of the Czech Republic under the research programs LC-06008 (Center for Computer Graphics) and MSM 6840770014.
A Dual Quaternion Tutorial
Dual quaternions are not a standard tool in computer graphics, in contrast to regular quaternions. Basic information about them can be found in the literature [Bottema and Roth 1979; McCarthy 1990], which provides a good theoretical background, but without information about applications in computer graphics. In order to bridge this gap, in this appendix we provide a brief tutorial summa- rizing the properties which are most important to computer graph- ics practitioners. We assume that the reader is already familiar with regular quaternions, otherwise see for example [Dam et al. 1998; Hanson 2006]. Dual quaternions can be considered as quaternions whose elements are dual numbers. Let us therefore start our discus- sion with this simpler algebra.
A.1 Dual Numbers
The algebra of dual numbers, denoted as Rˆ, is similar to complex numbers: any dual number aˆ can be written as aˆ = a0 + ε aε , where a0 is the non-dual part, aε the dual part and ε is a dual unit sat- isfying ε2 = 0. The dual conjugate is analogous to the complex conjugate: aˆ = a0 − ε aε . Multiplication of two dual numbers is given as
(a0+εaε)(b0+εbε)=a0b0+ε(a0bε +aεb0) (17)
16
The following lemma establishes that the dual conjugate behaves like the complex conjugate with respect to multiplication.
Lemma 6.
For any aˆ,bˆ ∈ Rˆ, it is true that aˆbˆ = aˆbˆ.
Proof.
aˆbˆ = a0b0+ε(a0bε+aεb0)=
= ab−ε(ab+ab)=(a−εa)(b−εb)=aˆbˆ 000εε00ε0ε
The following lemmas present the formulas for dual quaternion in- version and square root.
Lemma 7. The inverse of a dual number a +εa ∈ Rˆ, where a ̸= 0ε0
0, is given as
1
a 0 + ε a ε
= 1 −εaε (18) a 0 a 20
Proof. In order to find the inverse of a dual number a0 + ε aε , we have to solve for b0 , bε in the following equation
(a0 + ε aε )(b0 + ε bε ) = 1
Using Equation (17), this reduces to the following two equations:
a0b0 = 1 a0bε+aεb0 = 0
Therefore, b0 = 1/a0 and bε = −aε /a20 , provided that a0 ̸= 0.
Note that purely dual numbers, that is dual numbers with a0 = 0, do not have an inverse. This is a fundamental difference from complex numbers, because every non-zero complex number has an inverse.
Lemma 8. The square root of a dual number a0 +εaε ∈ Rˆ, such thata0>0isgivenas
√ √ aε a0+εaε = a0+ε2√a
0
(19)
Proof. Tofindthesquarerootofadualnumbera0+εaε,wehave to solve the following equation
(b +εb )2=a +εa 0ε0ε
Writing out the non-dual and dual components, we obtain b2 = a , 2b b = a
000εε
from which the formula for dual square root readily follows (pro-
vided that a0 > 0).
Finally, note that the Taylor series of a function of dual argument
reduces to a finite sum
f(a0+εaε)= f(a0)+εaε f′(a0)
because the higher powers of ε are zero. For example, dual sine and cosine functions are therefore given as
sin(a0+εaε) cos(a0+εaε)
= sin(a0)+εaε cos(a0) (20) = cos(a0)−εaε sin(a0) (21)
A.2 Dual Quaternions
A dual quaternion qˆ can be written as qˆ = wˆ +ixˆ+ jyˆ+kzˆ, where wˆ is the scalar part (dual number), (xˆ, yˆ, zˆ) is the vector part (dual vec- tor), and i, j,k are the usual quaternion units. The dual unit ε com- mutes with quaternion units, for example iε = εi. Just like ordinary quaternions, dual quaternions are also associative, distributive, but not commutative. A dual quaternion can also be considered as an 8-tuple of real numbers, or as the sum of two ordinary quaternions, qˆ = q0 + ε qε . Conjugation of a dual quaternion is defined using classical quaternion conjugation: qˆ ∗ = q∗0 + ε q∗ε . However, dual
number conjugation also applies to dual quaternions: qˆ = q0 − ε qε , so we end up with two different conjugations (and we will actually need both of them). It is not difficult to verify that the order in which the conjugations are performed does not matter
∗∗∗∗∗ ∗∗ qˆ =q0+εqε =q0−εqε =(q0−εqε) =qˆ
√√ The norm of a dual quaternion is defined as ∥qˆ∥ = qˆ∗qˆ = qˆqˆ∗.
It is possible to simplify this expression, as is shown in the follow- ing lemma.
Lemma 9. For any dual quaternion qˆ = q0 +εqε , such that q0 ̸= 0, the norm can be written as
∥qˆ∥=∥q0∥+ε⟨q0,qε⟩ (22) ∥q0∥
Proof. Letusfirstexpand
qˆ ∗ qˆ = ( q ∗0 + ε q ∗ε ) ( q 0 + ε q ε ) = q ∗0 q 0 + ε ( q ∗ε q 0 + q ∗0 q ε ) =
= ∥q0∥2 +2ε⟨q0,qε⟩
By taking the dual square root (Equation (19)), we obtain the final
formula.
Unit dual quaternions are those satisfying ∥qˆ∥ = 1. According to the previous lemma, a dual quaternion qˆ is unit if and only if ∥q0∥ = 1 and ⟨q0,qε⟩ = 0. We denote the set of unit dual quater- nions as Qˆ 1 . Geometrically, Qˆ 1 is a 6-dimensional manifold em- bedded in 8-dimensional Euclidean space (called an image-space of dual quaternions [McCarthy 1990]).
Dual quaternions inherit many properties of regular quaternions. The following lemma states that conjugation of dual quaternions also swaps the order of multiplication (as is the case with regular quaternions).
Proof. Using Lemma 10 and the fact that a dual number commutes with a dual quaternion, we can compute
∥pˆqˆ∥2 = (pˆqˆ)∗(pˆqˆ) = qˆ∗pˆ∗pˆqˆ = ∥pˆ∥2qˆ∗qˆ = ∥pˆ∥2∥qˆ∥2 Taking the dual square root on both sides concludes the proof.
The inverse of a dual quaternion is defined only when q0 ̸= 0. In
this case, we have
− 1 qˆ ∗ qˆ =∥qˆ∥2
Note that the inverse of a unit dual quaternion is just conjugation.
Let us now turn our attention to the representation of rigid trans- formations using dual quaternions. As expected, unit dual quater- nions naturally represent 3D rotation, when the dual part qε = 0. If we have a 3D vector (v0 , v1 , v2 ), we define the associated unit dual quaternion as vˆ = 1 + ε(v0i + v1 j + v2k). The rotation of vec- tor (v0,v1,v2) by a dual quaternion qˆ can then be written as qˆvˆqˆ∗. This can be easily verified, because if qε = 0 then qˆ = q0 and qˆ vˆ qˆ ∗ simplifies to
q (1+ε(v i+v j+v k))q∗ =1+εq (v i+v j+v k)q∗ 0012000120
where q0(v0i + v1 j + v2k)q∗0 is the familiar formula for rotation by a regular quaternion.
What is interesting is that dual quaternion multiplication can also represent 3D translation. A unit dual quaternion ˆt, defined as
ˆt = 1 + ε ( t 0 i + t 1 j + t 2 k ) 2
corresponds to translation by vector (t0,t1,t2) (note that dual quater- nions work with half of the translation vector, analogous to classical quaternions, which work with half of the angle of rotation). To see this, let us expand
ˆt vˆ ˆt ∗
= ˆt(1+ε(v0i+v1 j+v2k)) t
ε 1+ 2(t0i+t1 j+t2k) = t t
= ˆt1+ε v0+0 i+v1+1 j+v2+2 k = 222
Lemma10. Foranydualquaternionspˆ,qˆ,itistruethat ( pˆ qˆ ) ∗ = qˆ ∗ pˆ ∗
Proof. Thisiseasytoverifybydirectcomputation
(pˆqˆ)∗ = (p0q0 +ε(pεq0 +p0qε))∗ =
= (p0q0)∗+ε((pεq0)∗+(p0qε)∗)= = q∗0p∗0+ε(q∗0p∗ε +q∗εp∗0)=qˆ∗pˆ∗
( 2 3 )
= 1+ε((v0 +t0)i+(v1 +t1)j+(v2 +t2)k),
which shows that the unit dual quaternion ˆt performs translation by
(t0,t1,t2).
General rigid transformation is a composition of rotation and trans- lation. Therefore, let us first apply the regular quaternion q0 (rep- resenting the rotational component) and then the dual quaternion ˆt (representing the translational component)
ˆt ( q 0 vˆ q 0 ∗ ) ˆt ∗ = ( ˆt q 0 ) vˆ ( q 0 ∗ ˆt ∗ ) = ( ˆt q 0 ) vˆ ( q 0 ∗ ˆt ∗ ) = ( ˆt q 0 ) vˆ ( ˆt q 0 ) ∗ Therefore, we can see that ˆtq0 is a dual quaternion that performs
rigid transformation. Expanding, we obtain
εε
ˆtq0= 1+2(t0i+t1j+t2k) q0=q0+2(t0i+t1j+t2k)q0 (25)
Now we can state the following:
Lemma 12. Every rigid transformation can be represented by a unit dual quaternion, and conversely, every unit dual quaternion represents a rigid transformation.
Proof. The first part of the statement follows from Equation (25), if we can show that ˆtq0 is unit. However, both ˆt and q0 are unit and therefore, using Equation (24), ∥ˆtq0 ∥ = ∥ˆt∥∥q0 ∥ = 1.
Another important property (again analogous to regular quater- nions) is so called multiplicative property of the dual quaternion norm.
Lemma11. Foranydualquaternionspˆ,qˆ,itistruethat ∥pˆqˆ∥ = ∥pˆ∥∥qˆ∥
(24)
17
To prove the second part, we consider a unit dual quaternion pˆ = p0 +εpε. We need to find (t0,t1,t2) and q0 so that
q +ε(t i+t j+t k)q =p +εp 0201200ε
1 Obviously, q0 = p0 and (t0,t1,t2) is given by the equation 2 (t0i +
t1 j + t2 k)p0 = pε . This equation can be solved if we show that the scalar part of pεp∗0 is zero. However, the scalar part of pεp∗0 is equivalent to ⟨pε,p0⟩, which is zero because pˆ is unit.
Let us assume that we already have a routine for conversion be- tween a 3 × 3 rotation matrix and a unit quaternion, as well as a routine for quaternion multiplication. Equation (25) then shows how to convert a 4 × 4 rigid transformation matrix to a unit dual quaternion. The opposite conversion, from a unit dual quaternion q0 + ε qε to a matrix is also straightforward. The rotation is just a matrix representation of q0 and the translation is given by the vector part of 2qε q∗0. Note that the conversion routines in the C language are available on our website (see Section 6.1).
Dual quaternion conjugations can be interpreted as follows. Unit dual quaternion conjugate qˆ∗ corresponds to the inverse transfor- mation of qˆ, as can be seen from Equation (25). The same equa- tion reveals that dual-number-like conjugation, qˆ, inverts transla- tion only, leaving rotation intact. Both conjugations applied to- gether, qˆ ∗ , correspond to the inverse rigid transformation but with the original translational part.
A.3 Connection to Spatial Kinematics
The geometric interpretation of regular quaternions follows from
the formula q = cos θ + s sin θ , which contains the axis of rotation 22
s and the angle of rotation θ . The following lemma generalizes that to dual quaternions. In spite of the lengthier proof, it is important because it reveals the geometrical interpretation of dual quaternions and their connection to spatial kinematics (Chasles’ theorem). Lemma 13. Let θˆ ∈ Rˆ and sˆ ∈ Qˆ1 with zero scalar part. Then
To verify the second condition, we compute
⟨q ,q ⟩= s sinθ0,s θε cosθ0 +s sinθ0 −θε sinθ0 cosθ0 0ε02022ε2222
= θε sinθ0 cosθ0 −θε sinθ0 cosθ0 =0 222222
because⟨s0,s0⟩=1and⟨s0,sε⟩=0.
To prove the second part of the statement, let us assume that qˆ ∈ Qˆ 1
isfixed. Thetaskistofindθˆ=θ0+εθε andsˆ=s0+εsε sothat
Equation (26) holds. According to Lemma 12, we know that qˆ
canbewrittenasqˆ=q0+εtq0. Asq0=cosθ0 +s0sinθ0,the 222
unknowns θ0 and s0 are given by converting q0 to an axis-angle representation. Now we must find θε and sε , which are given by
1tq =s θε cosθ0 −θε sinθ0 +s sinθ0 (27) 2002222ε2
Let us expand the left-hand side, i.e., write out the quaternion prod- uct using Equation (14)
11θθ tq0=tcos0+s0sin0 =
2222
θˆ θˆ qˆ = c o s + sˆ s i n
= 1tcos θ0 − 1 sin θ0 ⟨t,s0⟩+ 1(t×s0)sin θ0 (28) 222222
We will consider three cases: (i) q0 = 1, (ii) q0 = −1, (iii) q0 ̸= ±1. In the first case, we can set θ0 = 0. The Equations (27) and (28) thus reduce to
θε s0 = 1t 22
Therefore, it is sufficient to set θ = ∥t∥ and s = t/∥t∥, while s ε0ε
can be, for example, the zero vector. The second case is quite sim- ilar, except that we assume θ0 = 2π . The choice of θε , s0 and sε is the same as in the previous case.
In the last case, when q0 ̸= ±1, we know that θ0 ̸= 2kπ for any integer k. By comparing the scalar parts of Equations (27) and (28), we obtain
−θε sinθ0 =−1sinθ0⟨t,s0⟩ 2222
( 2 6 ) is a unit dual quaternion. Conversely, for every qˆ ∈ Qˆ1, there exists
22
θˆ ∈ Rˆ and sˆ ∈ Qˆ with zero scalar part such that Equation (26)
Since sin θ0 ̸= 0, this simplifies to θε = ⟨t, s0 ⟩, yielding the formula 12
holds.
Proof. First, let us show that qˆ is always unit. Therefore, we ex-
pand Equation (26) to give
qˆ = cos θ0 +εθε +(s0 +εsε)sin θ0 +εθε 22
and using Equations (20) and (21)
for θε . Comparison of the vector parts of Equations (27) and (28) gives
θε θ0 θ0 1 θ0 1 θ0 s0 2 cos 2 +sε sin 2 = 2tcos 2 + 2(t×s0)sin 2
This allows us to express sε explicitly
sε sin θ0 = 1 (t−s0θε)cos θ0 + 1(t×s0)sin θ0
22222 qˆ =cosθ0 −εθε sinθ0 +(s0+εsε) sinθ0 +εθε cosθ0 sε = 1(t−s0θε)cotanθ0 +1(t×s0)
222222222
isolating the non-dual and dual parts
θ0 θ0 θεθ0θεθ0 θ0 qˆ=cos2+s0sin2+ε s02cos2−2sin2+sεsin2
The expression of sε can be simplified, if we recall that θε = ⟨t,s0⟩ and, using Lagrange’s equation, t − s0 ⟨t, s0 ⟩ = (s0 × t) × s0 . This enables us to write
1θ
sε = (s0×t)cotan 0 +t ×s0 (29)
q0 qε 22
Toshowthatqˆ=q0+εqε isunit,wehavetoshowthat∥q0∥=1 and ⟨q0 , qε ⟩ = 0. Quaternion q0 is obviously unit (because s0 is).
From this expression we see that obviously ⟨s0 , sε ⟩ = 0 and there- fore, sˆ0 ∈ Qˆ1, as required.
18
We can conclude that dual quaternions are a special representation of the screw parameters – one with advantageous algebraical prop- erties.
A.4 Exponential and Logarithm
Equation (26) can be used to define the exponential of any dual quaternion pˆ with zero scalar part and non-zero non-dual part. Let
usintroducesˆ=pˆ/∥pˆ∥andθˆ=2∥pˆ∥,sothatpˆ=sˆθˆ andsˆisaunit 1θ2
Intuitively speaking, Lemma 13 states that any unit dual quaternion is composed of parameters θ0,θε ,s0 and sε . As we can see from the proof, θ0 is the angle of rotation, and unit vector s0 represents the direction of the axis of rotation (in the degenerate case with θ0 = 0 or θ0 = 2π, s0 represents the direction of the translation vector). Furthermore, θε = ⟨t,s0⟩ is the amount of translation along vector s0. The only slightly less intuitive variable is sε. However, if we recall Equation (4), we observe that the term
dual quaternion with zero scalar part. Then, the exponential can be
(s0×t)cotan 0 +t 22
r=
that occurs in Equation (29) in the form sε = r × s0 is the center of rotation. Therefore, we can conclude that the variables θ0 , θε , s0 and sε are parameters of the associated screw motion.
There is a close connection with a classical result of spatial kine- matics known as Chasles’ theorem [Murray et al. 1994]. Chasle’s theorem states that any rigid transformation can be described by a screw, i.e., rotation about an axis followed by translation in the direction of that axis. For example, in Figure 21 top, we have a teapot which is first rotated about axis s0 , ∥s0 ∥ = 1, and then trans- lated by vector t. The translation vector t can be decomposed to t0 = t∥ + t⊥ , where t∥ is the component parallel to s0 and t⊥ is the component orthogonal to s0 (formally, t∥ = s0 ⟨s0 , t⟩, t⊥ = t − t∥ ). The component t⊥ lies in the plane with normal s0, and therefore can be recovered by selecting the proper center of rotation r in that plane. If we shift the axis of rotation to point r, as shown in Fig- ure 21 bottom, we obtain the corresponding screw motion (i.e., with translational component reduced to t∥).
defined as follows
exppˆ=exp sˆθˆ =cosθˆ+sˆsinθˆ 222
As shown in Lemma 13, exp pˆ is a unit dual quaternion. The inverse mapping, from unit dual quaternions to dual quaternions with zero scalar part, is denoted as log. The following lemma establishes that dual quaternion exponential and logarithm behave in a similar way to their matrix counterparts.
Lemma14. Letqˆ=cosθˆ+sˆsinθˆ,whereqˆ,sˆ∈Qˆ1andsˆhaszero 22
scalar part. Then for any mˆ ∈ Qˆ1, both of the following equations
are true
ˆ ˆ mˆ sˆ θ mˆ ∗ = mˆ e x p sˆ θ
e x p
22 log(mˆ qˆmˆ ∗) = mˆ log(qˆ)mˆ ∗
mˆ ∗
( 3 0 )
(31) Proof. Since the scalar part of sˆ is zero, the same is true for the
scalar part of mˆ sˆ θˆ mˆ ∗ , as can be shown by direct computation. This 2
means that the exp on the left hand side of the first equation is well defined and, according to its definition, we can write
e x p mˆ sˆ θˆ mˆ ∗ = c o s θˆ + mˆ sˆ s i n θˆ mˆ ∗ = mˆ c o s θˆ + sˆ s i n θˆ mˆ ∗ 22222
because a dual number always commutes with a dual quaternion and mˆ mˆ ∗ = 1. This proves Equation (30). The proof of Equa- tion (31) is similar
t|| st
0 ∗ θˆ θˆ ∗ θˆ ∗ θˆ mˆ qˆ mˆ = mˆ c o s 2 + sˆ s i n 2 mˆ = c o s 2 + mˆ sˆ mˆ s i n 2
t T h e r e f o r e l o g ( mˆ qˆ mˆ ∗ ) = mˆ sˆ θˆ mˆ ∗ = mˆ l o g ( qˆ ) mˆ ∗ . ^2
screw axis
q0
r
Figure 21: Conversion of a rotation about axis s0 , ∥s0 ∥ = 1, fol- lowed by translation t (top) to the corresponding screw (bottom).
A power of a unit dual quaternion qˆ is then defined naturally:
qˆuˆ = exp(uˆlogqˆ) = cos uˆθˆ +sˆsin uˆθˆ 22
Dual quaternions exhibit so-called antipodality, i.e., the fact that both qˆ and −qˆ represent the same rigid transformation
(−qˆ)vˆ(−qˆ)∗ = (−qˆ)vˆ(−qˆ∗) = qˆvˆqˆ∗
The mapping between SE(3) and Qˆ1 is thus one to two. Note that this is equivalent to the antipodality of regular quaternions. How- ever, even though both q0 and −q0 represent the same rotation, the powers qt0 and (−q0)t are different: one corresponds to clockwise and the other to counterclockwise rotation, see Figure 22. There- fore, when converting matrices to dual quaternions, we must choose an appropriate sign. This depends on each particular application; our method to resolve antipodality in skinning is discussed in Sec- tion 4.1.
s0
19
qt,0
As mentioned in Section 3.4, the situation with more than two rigid transformations is considerably more complex. In particular, it is necessary to generalize Buss and Fillmore’s spherical averages [2001] to the space of unit dual quaternions (which is not a hy- persphere). Such an algorithm has been proposed in [Kavan et al. 2006] and called Dual quaternion Iterative Blending (DIB) – see Algorithm 3.
Algorithm 3 (DIB)
Input: Unit dual quaternions qˆ1,…,qˆn,
convex weights w = (w1,…,wn), desired precision p Output: Blended unit dual quaternion bˆ
bˆ = DLB(w;qˆ1,…,qˆn) repeat
xˆ = ∑ ni = 1 w i l o g ( bˆ ∗ qˆ i )
bˆ = bˆ e x p ( xˆ ) u n t i l ∥ xˆ ∥ < p return bˆ
An intuitive explanation of this algorithm is shown in Figure 24. In the first step, the input dual quaternions are left-multiplied by bˆ∗, which maps the initial estimate bˆ onto the identity. The logarithm mapping then transforms bˆ ∗ qˆ 1 , bˆ ∗ qˆ 2 into the tangent space of Qˆ 1 at the identity, giving xˆ1 = log(bˆ∗qˆ1), xˆ2 = log(bˆ∗qˆ2). The blended value xˆ = w1 xˆ 1 + w2 xˆ 2 is computed and projected back by the ex- ponential mapping. Finally, multiplication bˆ exp(xˆ) yields the unit dual quaternion closer to the exact solution.
It is easy to see that DIB is bi-invariant, which follows from the bi- invariance of DLB (Lemma 2) and the properties of dual quaternion exponential and logarithm (Lemma 14). Note that DIB is indeed a generalization of ScLERP. In fact, it can be shown that for n = 2, DIB terminates in just a single iteration [Kavan et al. 2006]. For n > 2, more iterations are generally required, but in practice, the al- gorithm converges quite quickly. However, a mathematical analysis of DIB (along the lines of [Buss and Fillmore 2001]) is a non-trivial task that has not been addressed in the literature so far.
It is unfortunately not possible to repeat the error analysis from the previous section for DIB, because we do not have any closed-form expression of DIB. It might still be possible to establish a (poten- tially non-tight) upper bound. However, that would require a deeper analysis of the geometry of Qˆ 1 (the case of n = 2 corresponds only to a simple subspace of Qˆ1, i.e., a screw motion). Note that even an empirical comparison of DLB and DIB (i.e., by sampling all pos- sible rigid transformations) is very challenging for n > 2, because of the phenomenon known as the curse of dimensionality (i.e., the
q1 ~ (-q )1 0~0
Figure 22: Dual quaternions inherit the antipodality of classical quaternions. In this example, transformation of the teapot by qt0 for t ∈ [0, 1] produces a counterclockwise rotation (longer trajectory), while transformation by (−q0)t leads to a clockwise one (shorter trajectory).
B Difference between DLB and ScLERP
As stated in Section 3.4, the problem of comparing DLB and ScLERP requires derivation of the dual angle βˆt from Equation (13)
and then comparing it with αˆ t for t ∈ [0, 1]. For brevity, we will use shorthand C for cos( α0 ) and S for sin( α0 ). Using this convention,
22 Equations (20) and (21) take the following form
αˆ α αˆ α =C−ε εS, sin =S+ε εC
cos
Now, we can expand the term in the denominator of Equation (13)
2222
αˆ αˆ 1−t+tqˆ2qˆ∗1 = 1−t+tcos 2 +nˆtsin 2 =
αε αε = 1−t+tC−εt2S+(n0+εnε)t S+ε2C =
αε αε = 1−t+tC+n0tS+εt − S+nεS+ n0C
2 2 r0 r
The newly introduced quaternions r0 and rε satisfy
which gives βˆ t
1 − t + t c o s ( αˆ )
2 =
1 − t + t C − ε t α ε S
2 =
1
∥1−t +tqˆ2qˆ∗1∥
= 1 −ε⟨r0,rε⟩ ∥r0∥ ∥r0∥3
ε
∥r0∥ = (1−t+tC)2+t2S2
⟨r0,rε⟩ = 1−t+tC+n0tS,−tαε S+tnεS+tαε n0C 22
= (t−1−tC)tαεS+αεt2CS=(t−1)tαεS 222
because nˆ = n0 + ε nε is a unit dual quaternion (with zero scalar part). Therefore, the denominator of our equation can be written as
∗ ⟨r0,rε⟩ ∥1−t+tqˆ2qˆ1∥=∥r0+εrε∥=∥r0∥+ε ∥r ∥
and its inverse
0
cos = 2
∥1−t +tqˆ2qˆ∗1∥
∥1−t +tqˆ2qˆ∗1∥
=
1−t+tC tαεS (1−t+tC)⟨r0,rε⟩ ∥r0∥ −ε 2∥r0∥ + ∥r0∥3
20
^^^^ ^^ x11=b*bx2
> r0 := sqrt((1-t+t*cos(alpha_0/2))^2 + t*t*sin(alpha_0/2)^2 ):
> f0 := (t*cos(alpha_0/2) + 1-t)/r0:
> fe := – t*alpha_e/2*sin(alpha_0/2)/r0 – (1-t + t*cos(alpha_0/2))*(t-1)*t*alpha_e/2* sin(alpha_0/2)/r0^3:
> ang := 2*arccos(f0):
> plot( subs(alpha_0 = Pi, t ->
ang(t)), t =
0..1, y = -0.1..3.14);
3 2.5 2 y 1.5
1 0.5
0
0.2 0.4
0.6
0.8 1
Pi, t ->
Pi, t ->
> anglediff := ang – alpha_0*t:
> evalf(minimize(subs(alpha_0 =
anglediff(t)), t = 0..1));
−.1422292715
> evalf(maximize(subs(alpha_0 =
anglediff(t)), t = 0..1));
.1422292755
> pitch := simplify(-2*fe/sin(ang/2));
pitch:=t2alpha esin(1alpha 0)(−%1−t+t%1) 2
((−1 +2t −2t %1−2t2 +2t2 %1)√1−2t +2t %1+2t2 −2t2 %1
t
%1 := cos( 1 alpha 0) 2
> plot( subs(alpha_0 = Pi, alpha_e=1, t -> pitch(t)), t = 0..1, y = -1..1);
t2(−1+%12) ) −1 + 2 t − 2 t %1 − 2 t 2 + 2 t 2 %1
1 0.8 0.6
y
0.4
0.2 0
0.2 0.4
0.6 0.8 1
> pitchdiff := pitch – alpha_e*t:
> evalf(minimize(subs(alpha_0 = Pi, alpha_e=1,
t -> pitchdiff(t)), t = 0..1));
−.1501415529
> evalf(maximize(subs(alpha_0 = Pi, alpha_e=1,
t -> pitchdiff(t)), t = 0..1));
.1501415529
t
^b Q1 1
q1 ^^
initial estimate
b*q2
b*q1 q2
2
log
^^
^
^^ ^^ x11^xx2
^
b exp(x)
^
3
q1
^ exp(x)
exp
improved estimate
q2
Figure 23: Upper bound of the difference between ScLERP and DLB for n = 2, computed using Maple.
Figure 24: Illustration of one iteration of the DIB algorithm.
fact that the number of required samples grows exponentially with dimension).
In order to determine the error practically, we therefore compare only the results of applying DLB and DIB in skinning. In particular, we compute the maximum difference of vertex positions obtained by DLB and DIB. As expected, we found that the deviation between DLB and DIB is maximal in situations with large joint rotations, such as in Figure 25. However, even in that case, the differences between DLB (Figure 25a) and DIB (Figure 25b) are not noticeable. Therefore, we have visualized the error by color in Figure 25(c, d). Note that in none of our test animations did the error exceed 0.29 units (considering the proportions of our character, the units correspond approximately to inches).
a) b)
c) d)
Figure 25: Comparison of DLB and DIB applied in skinning. The difference between DLB (a) and DIB (b) is imperceptible even for large joint rotations and therefore we visualize it using color in (c) and (d).
21
On the other hand, the difference between computation times of DLB and DIB is quite pronounced. Our CPU implementation takes 5.07ms with DLB and 15.33ms with DIB (with accuracy p = 10−5 ). Therefore, we conclude that the DIB algorithm, albeit theoretically perfect, is not very practical in skinning – the significantly increased cost does not yield significant improvement in deformation quality.
References
ALEXA, M. 2002. Linear combination of transformations. In
SIGGRAPH ’02: Proceedings of the 29th annual conference on Computer graphics and interactive techniques, ACM Press, 380–387.
ALLEN, B., CURLESS, B., AND POPOVIC ́ , Z. 2002. Articulated body deformation from range scan data. In SIGGRAPH ’02: Proceedings of the 29th annual conference on Computer graph- ics and interactive techniques, ACM Press, New York, NY, USA, 612–619.
ALLEN, B., CURLESS, B., POPOVIC ́, Z., AND HERTZMANN, A. 2006. Learning a correlated model of identity and pose- dependent body shape variation for real-time synthesis. In Pro- ceedings of the 2006 ACM SIGGRAPH/Eurographics sympo- sium on Computer animation, Eurographics Association, Aire- la-Ville, Switzerland, 147–156.
ANGUELOV, D., SRINIVASAN, P., KOLLER, D., THRUN, S., RODGERS, J., AND DAVIS, J. 2005. SCAPE: shape completion and animation of people. ACM Trans. Graph. 24, 3, 408–416.
AUBEL, A., AND THALMANN, D. 2000. Realistic deformation of human body shapes. Proc. Computer Animation and Simulation 2000, 125–135.
BARAN, I., AND POPOVIC ́, J. 2007. Automatic rigging and ani- mation of 3D characters. ACM Trans. Graph. 26, 3, 72.
BARR, A. H., CURRIN, B., GABRIEL, S., AND HUGHES, J. F. 1992. Smooth interpolation of orientations with angular velocity constraints using quaternions. ACM Trans. Graph., 313–320.
BELTA, C., AND KUMAR, V. 2002. An SVD-based projec- tion method for interpolation on SE(3). IEEE Transactions on Robotics and Automation 18, 3, 334–345.
BLOOM, C., BLOW, J., AND MURATORI, C., 2004. Errors and omissions in Marc Alexa’s Linear combination of trans- formations. http://www.cbloom.com/3d/techdocs/lcot_ errors.pdf.
BOTTEMA, O., AND ROTH, B. 1979. Theoretical kinematics. North-Holland Publishing Company, Amsterdam, New York, Oxford.
BUSS, S. R., AND FILLMORE, J. P. 2001. Spherical averages and applications to spherical splines and interpolation. ACM Trans. Graph. 20, 2, 95–126.
CAPELL, S., GREEN, S., CURLESS, B., DUCHAMP, T., AND POPOVIC, Z. 2002. Interactive skeleton-driven dynamic defor- mations. In SIGGRAPH ’02: Proceedings of the 29th annual conference on Computer graphics and interactive techniques, ACM Press, New York, NY, USA, 586–593.
CHAR, B., GEDDES, K., AND GONNET, G. 1983. The Maple symbolic computation system. j-SIGSAM 17, 3–4 (Aug./Nov.), 31–42.
CLIFFORD, W. 1882. Mathematical Papers. London, Macmillan.
CORDIER, F., AND MAGNENAT-THALMANN, N. 2005. A data- driven approach for real-time clothes simulation. Computer Graphics Forum 24, 2, 173–183.
DAM, E., KOCH, M., AND LILLHOLM, M., 1998. Quaternions, interpolation and animation. Technical Report DIKU-TR-98/5, University of Copenhagen.
DANIILIDIS, K. 1999. Hand-eye calibration using dual quater- nions. International Journal of Robotics Research 18, 286–298.
FONTIJNE, D., AND DORST, L. 2003. Modeling 3D euclidean geometry. IEEE Comput. Graph. Appl. 23, 2, 68–78.
FORSTMANN, S., AND OHYA, J. 2006. Fast skeletal animation by skinned arc-spline based deformation. In EG 2006 Short Papers, 1–4.
FORSTMANN, S., OHYA, J., KROHN-GRIMBERGHE, A., AND M C D O U G A L L , R . 2007. Deformation styles for spline-based skeletal animation. In SCA ’07: Proceedings of the 2007 ACM SIGGRAPH/Eurographics symposium on Computer animation, Eurographics Association, Aire-la-Ville, Switzerland, 141–150.
GOVINDU, V. M. 2004. Lie-algebraic averaging for globally con- sistent motion estimation. In CVPR (1), 684–691.
GUO, Z., AND WONG, K. C. 2005. Skinning with deformable chunks. Computer Graphics Forum 24, 3, 373–381.
HANSON, A. J. 2006. Visualizing Quaternions. Morgan Kaufmann Publishers Inc.
HEJL, J., 2004. Hardware skinning with quaternions. Game Pro- gramming Gems 4, Charles River Media, 487–495.
HOFER, M., AND POTTMANN, H. 2004. Energy-minimizing splines in manifolds. ACM Trans. Graph. 23, 3, 284–293.
HYUN, D.-E., YOON, S.-H., CHANG, J.-W., SEONG, J.-K., KIM, M.-S., AND JU ̈ TTLER, B. 2005. Sweep-based human deformation. The Visual Computer 21, 8-10, 542–550.
JACKA, D., REID, A., MERRY, B., AND GAIN, J. 2007. A com- parison of linear skinning techniques for character animation. In AFRIGRAPH ’07: Proceedings of the 5th international confer- ence on Computer graphics, virtual reality, visualisation and in- teraction in Africa, ACM, New York, NY, USA, 177–186.
JAMES, D. L., AND TWIGG, C. D. 2005. Skinning mesh anima- tions. ACM Trans. Graph. 24, 3, 399–407.
JOHNSON, M. P. 2003. Exploiting Quaternions to Support Expres- sive Interactive Character Motion. PhD thesis, MIT.
JOSHI, P., MEYER, M., DEROSE, T., GREEN, B., AND SANOCKI, T. 2007. Harmonic coordinates for character ar- ticulation. In SIGGRAPH ’07: ACM SIGGRAPH 2007 papers, ACM, New York, NY, USA, 71.
JU, T., SCHAEFER, S., AND WARREN, J. 2005. Mean value coor- dinates for closed triangular meshes. In SIGGRAPH ’05: ACM SIGGRAPH 2005 Papers, ACM, New York, NY, USA, 561–566.
J U T T L E R , B . 1994. Visualization of moving objects using dual quaternion curves. Computers & Graphics 18, 3, 315–326.
KAVAN, L., AND Zˇ A ́ RA, J. 2005. Spherical blend skinning: A real-time deformation of articulated models. In 2005 ACM SIG- GRAPH Symposium on Interactive 3D Graphics and Games, ACM Press, 9–16.
22
KAVAN, L., COLLINS, S., O’SULLIVAN, C., AND Zˇ A ́ RA, J., 2006. Dual quaternions for rigid transformation blending. Tech- nical report TCD-CS-2006-46, Trinity College Dublin.
KAVAN, L., COLLINS, S., Zˇ A ́ RA, J., AND O’SULLIVAN, C. 2007. Skinning with dual quaternions. In 2007 ACM SIGGRAPH Sym- posium on Interactive 3D Graphics and Games, ACM Press, 39– 46.
KRY, P. G., JAMES, D. L., AND PAI, D. K. 2002. Eigenskin: real time large deformation character skinning in hardware. In Proceedings of the 2002 ACM SIGGRAPH/Eurographics sym- posium on Computer animation, ACM Press, 153–159.
KURIHARA, T., AND MIYATA, N. 2004. Modeling deformable hu- man hands from medical images. In SCA ’04: Proceedings of the 2004 ACM SIGGRAPH/Eurographics symposium on Computer animation, ACM Press, New York, NY, USA, 355–363.
KURIHARA, T., AND NISHITA, T. 2007. Dual-quaternion skinning with non-rigid transformatio. In SCA ’07: Posters, Eurographics Association, Aire-la-Ville, Switzerland, 18–19.
LEWIS, J. P., CORDNER, M., AND FONG, N. 2000. Pose space deformation: a unified approach to shape interpolation and skeleton-driven deformation. In Proceedings of the 27th annual conference on Computer graphics and interactive tech- niques, ACM Press/Addison-Wesley Publishing Co., 165–172.
LI, J., AND HAO, P. 2006. Smooth interpolation on homogeneous matrix groups for computer animation. Journal of Zhejiang Uni- versity 7, 7, 1168–1177.
LUCIANO, C., AND BANERJEE, P. 2000. Avatar kinematics mod- eling for telecollaborative virtual environments. In WSC ’00: Proceedings of the 32nd conference on Winter simulation, So- ciety for Computer Simulation International, San Diego, CA, USA, 1533–1538.
MAGNENAT-THALMANN, N., LAPERRIE` RE, R., AND THAL- MANN, D. 1988. Joint-dependent local deformations for hand animation and object grasping. In Proceedings on Graphics in- terface ’88, Canadian Information Processing Society, 26–33.
MARTHINSEN, A. 1999. Interpolation in Lie groups. SIAM J. Numer. Anal. 37, 1, 269–285.
MCCARTHY, J. M. 1990. Introduction to theoretical kinematics. MIT Press, Cambridge, MA, USA.
MERRY, B., MARAIS, P., AND GAIN, J. 2006. Animation space: A truly linear framework for character animation. ACM Trans. Graph. 25, 4, 1400–1423.
MERRY, B., 2007. Personal communication.
MOAKHER, M. 2002. Means and averaging in the group of rota- tions. SIAM Journal on Matrix Analysis and Applications 24, 1, 1–16.
MOHR, A., AND GLEICHER, M. 2003. Building efficient, accurate character skins from examples. ACM Trans. Graph. 22, 3, 562– 568.
MURRAY, R. M., SASTRY, S. S., AND ZEXIANG, L. 1994. A Mathematical Introduction to Robotic Manipulation. CRC Press, Inc., Boca Raton, FL, USA, 413–414.
PARK, S. I., AND HODGINS, J. K. 2006. Capturing and animating skin deformation in human motion. ACM Trans. Graph. 25, 3, 881–889.
PARK, S. I., SHIN, H. J., AND SHIN, S. Y. 2002. On-line locomo- tion generation based on motion blending. In Proceedings of the 2002 ACM SIGGRAPH/Eurographics symposium on Computer animation, ACM Press, 105–111.
PEREZ, A., AND MCCARTHY, J. M. 2004. Dual quaternion syn- thesis of constrained robotic systems. Journal of Mechanical Design 126, 425–435.
PRATSCHER, M., COLEMAN, P., LASZLO, J., AND SINGH, K. 2005. Outside-in anatomy based character rigging. In SCA ’05: Proceedings of the 2005 ACM SIGGRAPH/Eurographics sympo- sium on Computer animation, ACM Press, New York, NY, USA, 329–338.
RHEE, T., LEWIS, J., AND NEUMANN, U. 2006. Real- time weighted pose-space deformation on the GPU. Computer Graphics Forum 25, 3, 439–448.
SCHEEPERS, F., PARENT, R. E., CARLSON, W. E., AND MAY, S. F. 1997. Anatomy-based modeling of the human muscula- ture. In SIGGRAPH ’97: Proceedings of the 24th annual con- ference on Computer graphics and interactive techniques, ACM Press/Addison-Wesley Publishing Co., New York, NY, USA, 163–172.
SHOEMAKE, K., AND DUFF, T. 1992. Matrix animation and polar decomposition. In GI ’92, Morgan Kaufmann Publishers Inc., San Francisco, CA, USA, 258–264.
SHOEMAKE, K. 1985. Animating rotation with quaternion curves. In Proceedings of the 12th annual conference on Computer graphics and interactive techniques, ACM Press, 245–254.
SHREINER, D., WOO, M., NEIDER, J., AND DAVIS, T. 2007.
OpenGL Programming Guide: The Official Guide to Learning OpenGL, Version 2.1. Addison-Wesley Longman Publishing Co., Inc., Boston, MA, USA.
SLOAN, P.-P. J., ROSE, III, C. F., AND COHEN, M. F. 2001. Shape by example. In Proceedings of the 2001 symposium on Interactive 3D graphics, ACM Press, 135–143.
TERAN, J., SIFAKIS, E., BLEMKER, S. S., NG-THOW-HING, V., LAU, C., AND FEDKIW, R. 2005. Creating and simulating skeletal muscle from the visible human data set. IEEE Transac- tions on Visualization and Computer Graphics 11, 3, 317–328.
WANG, X. C., AND PHILLIPS, C. 2002. Multi-weight enveloping: least-squares approximation techniques for skin animation. In Proceedings of the 2002 ACM SIGGRAPH/Eurographics sym- posium on Computer animation, ACM Press, 129–138.
WANG, R. Y., PULLI, K., AND POPOVIC ́, J. 2007. Real-time enveloping with rotational regression. ACM Trans. Graph. 26, 3, 73.
WANG, W., JU ̈ TTLER, B., ZHENG, D., AND LIU, Y. 2008. Com- putation of rotation minimizing frames. ACM Trans. Graph. 27, 1, 1–18.
WAREHAM, R., CAMERON, J., AND LASENBY, J. 2005. Appli- cations of conformal geometric algebra in computer vision and graphics. Lecture Notes in Computer Science 3519, 329–349.
WEBER, O., SORKINE, O., LIPMAN, Y., AND GOTSMAN, C. 2007. Context-aware skeletal shape deformation. Computer Graphics Forum (Proceedings of Eurographics) 26, 3.
YANG, X., SOMASEKHARAN, A., AND ZHANG, J. J. 2006. Curve skeleton skinning for human and creature characters: Research articles. Comput. Animat. Virtual Worlds 17, 3-4, 281–292.
23