Fast corotational simulation for example-driven deformation

Computers & Graphics 40 (2014) 49–57

Contents lists available at ScienceDirect

Computers & Graphics journal homepage: www.elsevier.com/locate/cag

Technical Section

Fast corotational simulation for example-driven deformation$ Chao Song a,1, Hongxin Zhang b,1, Xun Wang a,n, Jianwei Han a, Huiyan Wang a a b

Department of Computer and Information Engineering, Zhejiang Gongshang University, Hangzhou 310018, China State Key Lab of CAD&CG, Zhejiang University, 310058, Hangzhou, China

art ic l e i nf o

a b s t r a c t

Article history: Received 9 June 2013 Received in revised form 24 December 2013 Accepted 12 January 2014 Available online 4 February 2014

We present a fast corotational finite element framework for example-driven deformation of 3-dimensional solids. An example-driven deformation space and an example space energy is constructed by introducing the modified linear Cauchy strain with rotation compensation. During this simulation, our adopted total energy functional is quadratic, and its corresponding optimization can be quickly worked out by solving a linear system. For addressing the possible errors, we propose an effective error-correction algorithm. Some related factors including the parameters and example weights are also discussed. Various experiments are demonstrated to show that the proposed method can achieve high quality results. Moreover, our method can avoid complex non-linear optimization, and it outperforms previous methods in terms of the calculation cost and implementation efficiency. Finally, other acceleration algorithms, such as the embedding technique for handling highly detailed meshes, can be easily integrated into our framework. & 2014 Elsevier Ltd. All rights reserved.

Keywords: Example-driven deformation Corotational finite element method Physical based simulation

1. Introduction Simulation of deforming objects is one of the active topics in computer graphics, and it has been widely applied in the movie industry, digital entertainment and product design. Over the years, researchers have conducted a large amount of excellent work on model deformation simulation, namely, editing, material modeling, user interaction and other related issues. With the repaid development of deformable model acquisition, example-driven deformation techniques [1–3] have received much attention, and their applications have also become the extension of the key frame animation and physical simulation. Most of the approaches in the example-driven deformation domain construct an energy functional of the object and solve the corresponding optimization problems for obtaining the node positions and reconstructing the shapes. In general, the energy functional is composed of two parts, which are the deformation energy and the constraint energy introduced by examples. Moreover, the later part influences the simulation results by constructing an example deformation space. Both parts can be constructed by geometric and physical methods. In comparison, physical

methods can more easily reflect the physical properties of the model materials. Example-driven deformation can produce more diverse effects than conventional key frame techniques in many applications. However, the induced optimization problem always leads to complex non-linear solutions. This limits the simulation response time as well as the solution scale. In this paper, an efficient computing framework for exampledriven deformation is presented. Inspired by existing approaches, we construct an energy functional of the deformable models in a physically based context. For modeling the deformation metric, we apply a revised Cauchy strain for efficient simulation. At the same, the corotational method is integrated to tackle the large deformation problem and also to form the example deformable space. Compared with previous methods, the proposed method is based on linear strain, and the induced optimization can be worked out quickly by simply solving a linear system. To address the possible errors, an effective error-correction algorithm is proposed. Moreover, it is easy for our framework to integrate the existing speed-up algorithms such as the embedding technique for highly detailed meshes.

2. Related work ☆

This article was recommended for publication by M. Botsch. n Corresponding author. Tel.: þ 86 571 88071024x8534; mobile: þ 86 13858009473; fax: þ86 571 28008303. E-mail addresses: [email protected] (C. Song), [email protected] (H. Zhang), [email protected] (X. Wang), [email protected] (J. Han), [email protected] (H. Wang). 1 Chao Song, Hongxin Zhang are the jointed 1st authors. http://dx.doi.org/10.1016/j.cag.2014.01.003 0097-8493 & 2014 Elsevier Ltd. All rights reserved.

We will briefly review related work on physically based animation, shape interpolation, example driven deformation and other relevant topics. In the 1980s, the pioneering work in the field of physically based animation was proposed by Terzopoulos et al. [4]. The main

50

C. Song et al. / Computers & Graphics 40 (2014) 49–57

goal of this research field is to improve the simulation reality, speed and stability for computer animation. A large number of studies have been performed, including mass–spring, FEM and meshless approaches. Many issues have been thoroughly discussed, including the approximated solution of a control PDE, processing on large deformations and other topics. A more detailed overview can be found in [5–7]. To tackle the nonlinear issue in a large deformation, the corotational method [8], which is also called stiffness warping [9,10], has been widely applied. This category of methods factorizes the large deformation into two parts, which are the linear deformation and the rigid rotation quantities. In each simulating step, the stiffness matrix composes the linear stiffness and rotation compensation. Thus, the treatment for non-linear factors can be simplified and the simulation procedure can be sped up. Recently, the corotational technique has become a mainstream method in many computer graphics applications [11,12]. Shape interpolation is an important issue in geometry processing. Many techniques have been proposed for morphing and shape editing. Alexa et al. [13,14] proposed an array of methods for shape interpolation as rigid as possible. Subsequently, many researchers have conducted a substantial amount of work addressing the rotation issue during interpolation, such as Lipman et al. [15]. Based on the differential coordinates of a surface mesh, Xu et al. [16] proposed a Poisson equation based formulation for shape interpolation. Huang et al. [12] added dynamic effects in the shape interpolation. There is a substantial amount of other work, such as [17,18] which tackled the interpolation issue of spatial rotation. In addition, as a natural extension of shape interpolation, Kilian et al. [19] discussed techniques to construct deformation space from examples. Surface deformation has received continuous attention for many years [20]. The MeshIK [21] is one of the representative methods for creating deformable models based on deformation examples. This method has been extended in later publications [22,1]. Moreover, it is one of the research bases for many techniques of example driven deformation. Frohlich et al. [1] proposed a framework of example-driven deformation based on discrete shells, and the framework can generate static deformation of triangle meshes, combining examples and physical laws. More recently, Martin et al. [2] presented an example-based method for elastic materials that implements dynamic deformations following the physical laws. They adopted nonlinear Green strain to create a deformation metric space and construct the solved energy functional, including the elastic deformation energy and the energy that reflects the examples effects. Then, they solved and obtained the deformable shapes of the elastic models by using the Newton–Raphson optimization, and a performance nickel is the nonlinear optimization. Later, Schumacher et al. [3] implemented a similar method for elastic– plastic deformation. Specifically, as an extension of the method in [3], Coros et al. [23] presented a method for controlling the motions of active deformable characters. Koyama et al. [24] also formulated an analogous concept to generate plausible animation in a shape-matching framework. In this paper, we present an example-driven deformation framework, in which the enhanced Cauchy strain is adopted to quantize the deformation metric and a corotational technique is applied to compensate the artifact induced by linear calculation in large deformations. In the following sections, the formulation of the problem is presented and the applications of static and dynamic deformation are discussed. At the same time, an effective error-correction algorithm is given. Then, we use an embedding technique to speed up the computing procedure. Moreover, several examples are presented to prove the efficiency of the method, and several related factors in deformation calculating is

briefly discussed. Finally, the limitations of the proposed method and our future work are described.

3. Formulation Without loss of generality, we assume that a given solid Ω  R3 is discretized into a linear tetrahedral mesh with n nodes and m elements. Let X, x A R3n denote position vectors that describe the initial and deformed configurations, respectively. x1 ; …; xk A R3n represent the position vectors of k input deformation examples. We will construct an example driven-deformation space, and define the solved energy functional by using the corotational Cauchy strain. 3.1. Corotational Cauchy strain The behavior of deformable solids is mainly modeled as the movements of the inner points in elastic mechanics. Given a solid, the movement of any point p A Ω can be denoted as p : Ω  R-R3 : ðX; 0Þ↦xðX; tÞ: In a linear FEM simulation, the procedure of solid deformation is basically depicted by the deformation gradient, which can be calculated by F ¼ ∂x=∂X; where F is not rotation invariant. In a large deformation, a rigid rotation can lead to the changing of F and a straightforward linear FEM calculation can produce undesirable artifacts due to the nonlinear nature of the rotation transform. Thus, we apply corotational FEM simulation for the rotation compensation. For any one tetrahedron element j, the corotational deformation gradient is calculated as e F^ j ¼ ∂ðR ej > xej Þ=∂X ¼ Rej > Fej

where superscript e indicates that the vector or matrix is calculated from an element. Rej is a block diagonal matrix, each of whose sub-blocks is a mapping element rotation matrix. The element rotation matrix is a 3  3 unitary matrix that can be computed by the polar decomposition of the deformation gradient Fej . Thus, we have a corotational version of the Cauchy strain metric for the element >

εj ¼ ðF^ j þ F^ j Þ=2; which is rotation invariant. Moreover, the corotational Cauchy strain of an element can be represented as

εej ¼ Bej ðR ej > xej  Xej Þ

ð1Þ

where Xej ¼ ½Xej;1> jXej;2> jXej;3> jXej;4>  > and xej ¼ ½xej;1> jxej;2> jxej;3> jxej;4>  > are the initial and deformed coordinates of the nodes in the j-th element, respectively. Bej is a constant matrix related to the geometry of the element, and the details can be found in [25–27]. Furthermore, a 3  3 strain matrix can be denoted as a 6-dimensional vector by using Voigt notation. For any one shape of the given deformable model, a unique element strain vector is defined

ε ¼ ½εe1 εe2 ⋯εem   R6m To quantify the deformation, we introduce U deform ðxÞ to measure the deformation energy between the current shape x and the initial configuration X: Z 1 1 m U deform ¼ ε > Dɛ dV ¼ ∑ εej > Dɛ ej V ej 2 V 2j¼1 where D is the properties matrix of the material, εej denotes the strain of the j-th element, and Vej is the corresponding element volume.

C. Song et al. / Computers & Graphics 40 (2014) 49–57

51

We can further rewrite the above as follows: U deform ¼

1 m ∑ ðR e > xej  Xej Þ > Kej ðR ej > xej  Xej Þ; 2j¼1 j

ð2Þ

where Kej is the j-th element stiffness matrix which is deduced from D and B in the linear FEM frameworks. According to the physical mean, we call the R ej > xej  Xej the unrotated displacement of an element. For a known deformation example i, all of the unrotated displacements fR e1;i> xe1;i  Xe1 ; Re2;i> > e xe2;i  Xe2 ; …; R em;i xm;i  Xem g can be calculated. In addition, only to simplify the notation, in this paper, we express the unrotated displacement of all elements as R Ti xi X for the i-th deformation example. Moreover, we define the following matrix rewriting rule: m

Finally, the interpolated shape xw with respect to w is constructed xw ¼ arg minU space ðx; wÞ x

After performing a partial derivation, the above optimization problem is equivalent to solving the linear system:

∑ Aej Bej ≔AB

j¼1

k

Thus, Eq. (2) can be rewritten as U deform ¼

Fig. 1. Corotational Cauchy strain based interpolation (the far left shape) between a normal and a twist-bend (the far right shape) cube.

RKðR > xw  XÞ ¼ ∑ wi RKðRi> xi  XÞ i¼1

1 > ðR x  XÞ > KðR > x  XÞ 2

Specially, the energy interpolation between two deformation examples is RKR > xw ¼ ð1  wÞRKR 1> x1 þ wRKR 2> x2

3.2. Example-driven deformation space Given k deformation examples xi ði ¼ 1; …; kÞ, we define a deformable space expanded by k unrotated displacement vectors:

where xw is the generated immediate shape and ∑ki ¼ 1 wi ¼ 1 is the default. Fig. 1 shows an interpolation result of a cube model.

fR1> x1  X; R2> x2  X; …; R k> xk  Xg

4. Example-driven deformation

Let w ¼ ½w1 ; …; wk  > be the weight vector of the examples, and lerpðR 1> x1 X; …; R k> xk  X; wÞ express a weighted example, which is the unrotation displacement of a shape for w in the deformable space. Every w maps to an unrotation displacement. To calculate the current shape x, the best estimation must to find to minimizes the elastic potential energy between x and the closest shape in the example space with a given w. Note that the unrotation displacement difference between x and the weight example can be written as

Based on Section 3, we propose a fast example-driven deformation approach that is called corotational Cauchy example-driven deformation (CCED). The example-driven deformation is regarded as a general deformation simulation with an example space energy constraint. In other words, the penalty method is adopted to exert the constraint conditions, and the total energy functional of the example driven deformation can be written as

Lðx; wÞ ¼ R

>

Π ¼ T þ U deform ðxÞ þ λU space ðx; wÞ  W ext k

s:t:

x  X  lerpðR 1> x1 X; …; R k> xk  X; wÞ

∑ wi ¼ 1; wi Z0

ð4Þ

i¼1

where the rotation R in the current step is obtained from the result of the previous step. The above formula is also an abbreviation according to our rewriting rule. In detail, for any one element j in the model, Lej ðxej ; wÞ ¼ R ej > x X  lerpðRej;1> xej;1  Xej ; …; Rej;k> xej;k  Xej ; wÞ. Further, a weighted example is

where T is the kinetic energy, λ is a control parameter, and > W ext ¼ f ext ðx  XÞ is the work of the external force f ext . The CCED approach can be solved as an optimization problem defined on the energy functional Π in Eq. (4). We will introduce the solution in the static and dynamic cases. Moreover, an errorcorrection scheme will be given for the stability and high performance.

k

lerpðR 1> x1 X; …; R k> xk  X; wÞ ¼ ∑ wi ðR i> xi  XÞ

4.1. Static deformation

i¼1

Because the corotational method is based on transforming the examples shapes to the initial coordinate frame, we do not need to nonlinearly interpolate the rotation quantities as in MeshIK [21] for a weighted example. In fact, our experiments show the same results as were obtained by using an exponential matrix, e.g., ∑ki ¼ 1 ½ðexpðwi log R i ÞÞ > wi xi  X. Thus, the corresponding elastic potential energy between the current shape x and a weighted example, which is produced by the deformation examples and is also called the example space energy, can be calculated as 1 1 m U space ðx; wÞ ¼ L > ðx; wÞKLðx; wÞ ¼ ∑ Lej > ðxej ; wÞKej Lej ðxej ; wÞ 2 2j¼1

ð3Þ

Compared with using jLðx; wÞj2F , which is an example energy metric in [2], Eq. (3) is more convenient for setting energy parameter later and has a more explicit physical means.

Static deformation simulation is a mainstream method for shape editing, deformable model generation and modification. The common characteristics of these applications are low or uniform speed. Thus, the kinetic energy can be omitted in the simulation. As discussed in EBEM [2], we also relax the weighting constraints and apply a simplified version of the functional defined in Eq. (4) for performing a static deformation:

Π static ¼ U deform ðxÞ þ λU space ðx; wÞ  W ext

ð5Þ

After expanding the above equation by substituting Eqs. (2) and (3), the static energy functional Π static becomes a quadratic. Minimizing it corresponds to exploring that the gradient, which is subject to both x and w, is zero. In a physical sense, the part of the gradient corresponding to x is simply the forces per node that must sum up to zero, whereas the part corresponding to w seeks that the choice for w is optimal, which is the only choice of w that

52

C. Song et al. / Computers & Graphics 40 (2014) 49–57

will lead to a lower example potential. The optimization can be deduced and eventually converted into a linear system in the form of Ax ¼ b, as follows: " #  " # Axx Axw bx x ¼ ð6Þ Awx Aww bw w where Axx ¼ ð1 þ λÞRKR > > Axw ¼ Awx ¼  λ½RKðR1> x1  XÞ; …; RKðR k> xk  XÞ >

ðAww Þji ¼ λðRj> xj  XÞ > KðR i> xi XÞ bx ¼ ð1 þ λÞRKX þ f ext ðbw Þj ¼  λðRj> xj  XÞ > KX

ðj ¼ 1; …; kÞ

In the above equations, the parameter λ is an important control scalar that indicates the effect of the examples during the deformation. The coefficient matrix A is symmetric and nonnegative definite, and at the same time, the position of non-zero entries retains unchanged and only their values change in each calculation step. The Cholesky decomposition can be performed for fast simulation. Moreover, in Eq. (4), many items can be precomputed at the beginning of the simulation. One static deformation example is illustrated in Fig. 2. 4.2. Dynamic deformation Dynamic deformation can be generated when the kinetic energy is not zero. The technique of dynamic deformations driven by examples can generate vivid deformable model animation, although the generated deformation effects are not better controllable than key control animation.

In Eq. (4), the kinetic energy is T ¼ 12x_ > Mx_ according to the kinetics law, where M is the mass matrix. To solve the optimization problem while considering the stability, we adopt the technique of Newmark time integration, which is a class computing method in kinetics. In brief, vt þ 1 ¼ vt þ ½ð1  δÞat þ δat þ 1 h and 2 at þ 1 ¼ xt þ vt h þ ½ð0:5  αÞat þ αat þ 1 h are adopted, where h is the time step length, superscript t expresses the t-th time step, v ¼ x_ and a ¼ x€ are the velocity and acceleration, respectively. When the parameters δ Z0:5 and α Z0:25ð0:5 þ δÞ2 , the Newmark method is unconditionally stable. After performing time integration, it is necessary to seek that x and w be optimal in the sense that they yield the lowest possible example energy at the end of a time step. Thus, similar to solving Eq. (5), the simulation of dynamic deformation is eventually performed by solving the following linear system at each timestep: t þ1

At þ 1 xt þ 1 ¼ b

:

The above equation can be written as " t þ1 #" # 2 t þ1 3 þ1 Atxw bx Axx xt þ 1 ¼ 4 t þ1 5 t þ1 tþ1 tþ1 w Awx Aww bw

ð7Þ

where Atxxþ 1 ¼

1

M þ ð1 þ λÞRKR >   1 1 t tþ1 t tþ1 t t v bx ¼ bx þð1 þ λÞRKX þ M x þ þ ð2 α  1Þa þ f ext αΔt αΔt 2

αΔt 2

The other coefficients can be calculated by a method similar to Eq. (6), and R is calculated from the result of the last time step in Eq. (7). In our implementation, such as in Fig. 3, M is a lumping mass matrix of FEM.

Fig. 2. Static deformation of a cat. (a) is an initial model, (b) is an example model, (c)–(f) are deformation results. The fixed regions are painted in blue, and the green regions are the transformed parts. (For interpretation of the references to color in this figure caption, the reader is referred to the web version of this paper.)

C. Song et al. / Computers & Graphics 40 (2014) 49–57

2: 3: 4: 5: 6: 7: 8: 9: 10: 11: 12: 13: 14: 15: 16: 17: 18: Fig. 3. Dynamic deformation of a cube under gravity with a single twisted cube example.

Adopting appropriate damping could be a key control factor in the dynamic deformation. A Rayleigh damping can be combined into our simulation framework, and the process for the calculation is straightforward. 4.3. Error correction Our method avoids the calculation for obtaining Hessian matrix. However, it is well-known that the plain corotational method has its disadvantages. There are various studies which have proposed various improvements, such as [28,29]. In fact, the corotational formulation assumes that rotations stay constant in a given time step (or during a calculation step). In reality, this assumption is only an approximation and rotations can change quite significantly, especially in example-driven deformation. The accumulated error could be large because the rotation matrix is based on the result of the previous time(or calculation) step. Two methods can be adopted to address the problem. One method reduces the calculation step (including the time step), and the other performs cross-iteration between x and w in Eqs. (6) or (7). Our cross-iteration method is inspired by block coordinate descent [29]. The error after each step calculation is checked as follows: First, the rotation R ei;c is calculated, which is based on the current step calculation result. Then, adopting the criterion described in [30], we estimate the change between the most recent evaluation steps R ei;l and R ei;c by means of calculating the maximum absolute row sum norm dc ¼ jFei;l  Fei;c j1 , where Fei is the deformation gradient of the element at the previous calculation step or the current calculation step. When dc exceeds a given threshold dmax , which is generally set as dmax ¼ 0:1, we must perform cross-iteration to tackle the errors. Once the crossiteration begins running, the iteration termination criterion based on the variation of w should be combined into the process. In detail, the weights of the current and previous iteration step are denoted as wtc and wtl , respectively. When J wtc  wtl J exceeds a given threshold dwmax , which is set to 0.05 in our experiments, the corresponding object is not in a state of truly stable equilibrium, and the simulation stability will be reduced, especially for the dynamic deformation. The implementation of cross iteration is illustrated in Algorithm 1. Algorithm 1. Error correction for CCED. 1:

Input: x ¼ ½xc ; wc ;

19: 20: 21: 22:

53

bIteration ¼True; while (bIteration) do bIteration¼ False; for each element i do dc ¼ jFei;l  Fei;c j1 ; if (dc 4 dmax ) or J wtc  wtl J 4 dwmax then bIteration ¼True break; end if end for if (bIteration) then wtl ’wc calculate R by xc , update A; solve Ax ¼ b in Eqs. (6) or (7) while fix w, update xc ; calculate R by xc , update A; solve Ax ¼ b in Eqs. (6) or (7), update xc , wc and Fei;c ; Fei;l ’Fei;c ; wtc ’wc end if end while Output: x.

Because the stability change induced by example weight is difficult to estimate, Algorithm 1 is easier to control than step reduction. Fig. 4 demonstrates a set of results for static deformation in which two legs of the dinosaur were fixed and one hand was dragged to the same position. During the calculation, the artifact, which appeared at the beginning of the iteration, is eliminated after 4 iterations according to Algorithm 1. Fig. 5 depicts the static deformation of a horizontal bar under gravity by a one-step calculation. Fig. 5(c) shows that the result driven by one example is an artifact if the error is not corrected. Fig. 5 (d) is the result of the 7-th step, simply using an iteration without cross-processes, and its efficiency is similar to the step reduction method. Fig. 5(e) is the result of the 7-th iterative step by using Algorithm 1. Fig. 5(f) is the variation curve of w with iteration increasing by using the two methods. It is obvious that the equilibrium is achieved faster by using Algorithm 1. Our proposed error correction is able to avoid strange popping artifacts in the case of large times steps or fast rotations. It can be taken as a variant of block coordinate descent on the rotation R and example weight w. In all experiments, we did not find the case of unstableness or non-convergence.

5. Embedding acceleration Mesh embedding techniques can be integrated into our framework for accelerating the simulation. In detail, a highly detailed mesh can be embedded into a coarse deformed model agent. Deformations of the actual model can be calculated according to deformations of the model agent. Therefore the performance of the calculation can be greatly improved. At the same time, this method can be useful for deforming complex non-manifold models. The original solid and its deformed shapes can be surface or volumetric meshes. The embedding model is a tetrahedral mesh in our framework because the mapping between the target and the embedding model can be easily obtained by applying the tetrahedron barycentric coordinate interpolation, or other volumetric coordinate interpolation techniques, such as mean value coordinates [31]. Because all of those methods have an extrapolation ability, our embedding models do not need to strictly wrap the target model. It must be noted that the deformation examples should also be embedded by the homeomorphism embedding model. Given a target model and its deformation examples, automatically generating a homeomorphism embedding model is challenging and we

54

C. Song et al. / Computers & Graphics 40 (2014) 49–57

Fig. 4. Dinosaur: error correction for static deformation. (a) initial model, (b) example model, (c) deformation with error, (d) the results after using the correction algorithm. The blue regions are fixed and the green regions are transformed. (For interpretation of the references to color in this figure caption, the reader is referred to the web version of this paper.)

Fig. 5. A horizontal bar with one end fixed: error correction for static deformation (only a one-step calculation). (a) Deformation under gravity (without example), (b) example model, (c) deformation without error correction (w¼ 0.928732) , (d) the result after 7 simple iterations (w¼ 0.533981), (e) the result after 7 iterations within Algorithm 1 (w¼0.479757), (f) the variation in w with the number of iterations (blue curve: using a simple iteration and red curve: using Algorithm 1). A green semitransparent graph is the initial position of the bar, and the blue regions are the fixed parts. (For interpretation of the references to color in this figure caption, the reader is referred to the web version of this paper.)

are not concerned about it. To merely verify the effectiveness of the method, we first make an embedding model for the input target mesh, and then, we generate deformed meshes for the deformation examples while maintaining the correspondence. An embedded model with a high resolution can be more computationally efficiency when embedding model has a low resolution. However, the deformations of the embedded models follow coarse physics information, so the physical detail could be generally weakened during a direct deformation, as is also mentioned in [2]. Fortunately, deformation examples imply rich information on the physical behavior. Therefore, the exampledriven deformation still can create physically plausible deformation effects, such as the results in Fig. 6.

6. Implementation and discussion Our proposed framework for the example-driven deformation can be easily implemented through expanding the existing FEM

code. In this section, we provide a procedure overview and discuss several important related issues.

6.1. Method overview For a given target model and one or several deformation examples, the calculation procedure can be summarized as follows: Step 1: Pre-computing. We set the target model as the initial configuration X and perform a linear FEM calculation for obtaining linear stiffness matrix Kei of every element. For all of the example models, we calculate the deformation gradient Fei of every element and obtain the element rotation matrix Rei by making a polar decomposition to Fei . Then, the whole matrix R i KR j> and KR i> are computed. At the same time, since the linear stiffness K and all of the rotations Ri are symmetric, we have KRi> ¼ R i K. Step 2: Parameter configuration. In this stage, the parameter λ in the energy functional in Eq. (4) is mainly considered.

C. Song et al. / Computers & Graphics 40 (2014) 49–57

55

Fig. 6. Dynamic deformation of a dinosaur with an embedded cage: where (a) is the initial model, (b) is an embedding cage of the model, (c) is a deformed example of the cage, and (d)–(f) are deformation results. In (d)–(f), the blue regions are fixed, and the green regions are transformed interactively. (For interpretation of the references to color in this figure caption, the reader is referred to the web version of this paper.)

Step 3: Solution. According to applications, one equation is selected and solved with the corotational method from the static (Eq. (6)) and dynamic simulation (Eq. (7)). During computing, we can add position constraints, which is necessary for a static solution. After each calculation step is accomplished, an errorcorrection procedure is performed according to Algorithm 1. In this stage, exerting constraints is as convenient as common deformation simulation without examples. Finally, the calculation results are rendered, similar to any other deformation approaches.

6.2. Discussion Parameter λ: In Eq. (4), if λ is sufficiently large, it is clear that the deformation energy Udeform has minor impacts on the deformation results and can even be omitted. Thus, the problem evolves into simulating deformations that are produced only by example models and external forces. Conversely, if λ is sufficiently small, the deformation examples will hardly affect the calculation result. Therefore, λ should be set to meet specific requirements. During our implementation, the deformation energy and example space energy have the same order of magnitude. The value of λ ¼ λ=ð1 þ λÞ is chosen according to Eq. (6) in the static deformation cases. When λ -1, the deformation energy Udeform has almost no influence on the simulation, and when λ -0, the example space energy Uexample does not affect the simulation. For dynamic deformation simulation, it follows the same pattern. The only difference is that the acceleration change of the object could be taken into consideration for the specific application because example space energy is added. In detail, when the parameter λ is very large, the external force has almost no effects to deformation. As a typical application, in the character or other deformations, including the skeleton (such as the cat in Fig. 1), λ should be sufficiently large unless specialized processing is performed for the skeleton.

The weighting constraints in Equation (4): Strictly speaking, the constraints which are ∑ki ¼ 1 wi ¼ 1 and wi 4 0 should be added in Eq. (4) to form interpolation presentations. However, these constraints restrict the representative ability of the example space. Relaxing the constraints can be allowed and is valuable for some real applications in our framework. The positive constraint wi Z0 is not always necessary for deformations because sometimes it is reasonable to obtain extrapolations of the deformed examples. A similar situation is also mentioned in [2]. For example, a bridge model is demonstrated in Fig. 7 with two deformed examples. The weight value of the 1st example (Fig. 7(b)) is shown in Fig. 8, in which the abscissa axis denotes the frame number and the vertical axis denotes the weight value. In Fig. 7(f), the weight w is negative (w ¼ 0:32), and the result is acceptable. The normalization constraint ∑ki ¼ 1 wi ¼ 1 in Eq. (4) can be exerted by adding a penalty ηð∑i wi  1Þ2 to the functional Π, where η is the penalty parameter. The constraint can make the strong conjunction among examples in calculation. For example, in Fig. 9, we constrain the value of weights with w1 þw2 ¼ 1 so as to create a smooth transition between the two deformation example states. Conversely, unless the constraint ∑i wi ¼ 1 is exerted, in theory, the probability of which some deformation examples do not affect the calculation results would increase when more than 2 examples are included in the calculation. 6.3. Performance In general, the performance of our method is close to the stiffness warping approach [9]. Compared with their simulation, our method only increases the scale of linear system with the number of given examples as well as an additional processing of error correction. Other additional calculation is due to the precomputer step. Because the number of examples is limited, the computational cost of each calculation step mainly depends on the number of nodes and tetrahedra in the simulating object.

56

C. Song et al. / Computers & Graphics 40 (2014) 49–57

Fig. 7. Dynamics deformation of a bridge model with two deformed examples. Images (d)–(f) are the results generated by hitting the top of the nearer tower.

Fig. 8. The curve of weight w1 of the bridge example. The blue curve shows the wight of the first example during deformation in Fig. 7. The framerate of deformation is 25 fps. (For interpretation of the references to color in this figure caption, the reader is referred to the web version of this paper.)

We implemented all aforementioned methods in Microsoft windows 7 and Visual Studio 2008. The Intel MKL library for matrix calculation and the PARDISO API for solving linear systems were applied. The performance of our results is listed in Table 1. All data are collected in a common laptop with an Intel i5 2540M CPU and 3GB memory. For a comparison, we implemented the elastic deformation based on the method in [3] (also called by ESEBM) by using linear tetrahedral element and Newton–Raphson optimization. We respectively simulated the dynamic deformation of cube, bridge and bar models under the same constraint condition and physical parameters as Figs. 3, 7 and 9. The results that represent our method have similar effects, and the average simulation speed of our method is about 6.3, 6.2 and 6.5 times as fast as that of ESEBM. In addition, by using different resolution (555 and 10 002 DOFs), same physical parameters and λ (λ ¼100), the deformation of the dinosaur (the same as Fig. 6) is compared. On average, our single step speed is about 5.2 faster than that of ESEBM.

7. Conclusions

Fig. 9. A bar pendulum with two examples. Examples are S-shaped and U-shaped. From right to left, images show several deformed shapes at different simulation time during a half pendulum periodicity.

A corotational simulation framework for example-driven deformation is presented in this paper. By constructing an exampledriven deformation space with corotational technique, we propose an efficient simulation solver of the deformation energy functional based on the linear Cauchy strain metric. To ensure realistic results and simulation stability, we propose an effective algorithm for error correction. The proposed scheme enriches the corotational methods system. However, our method applies to the large deformation with large rigid rotation and small strain, similar to

C. Song et al. / Computers & Graphics 40 (2014) 49–57

57

Table 1 Performance of results (single calculation step). Model

Cat Cube Dinosaur Dinosaur Bridge Bar

Figure

Fig. Fig. Fig. Fig. Fig. Fig.

2 3 4 6 7 9

Tetrahedra

10 212 15 360 5249 379 (cage) 14 946 450

DOF

8580 9147 1493 10 002 9244 208

Examples

2 1 1 1 2 2

λ

100 20 100 100 10 10

classical rotational method. For the deformation with a large strain, the existing methods are still recommended to ensure a realistic simulation. In the future, we will extend our work and explore a combination with the modal analysis technique for higher performance. Moreover, we will also research example-driven deformation for special materials such as plastic and anisotropic materials in our framework. Acknowledgments This work was supported by in part the Natural Science Foundation of China (No. 61170098, 61100137, 61002003), Natural Science Foundation of Zhejiang Province (No. Y1091084, Z12F020017, Z1111051, Y12F020131), and the Scientific Research Fund of Zhejiang Provincial Education Department (No. Y200804713). Appendix A. Supplementary material Supplementary data associated with this article can be found in the online version of http://dx.doi.org/10.1016/j.cag.2014.01.003. References [1] Fröhlich S, Botsch M. Example-driven deformations based on discrete shells. Comput Graph Forum 2011;30(8):2246–57. [2] Martin S, Thomaszewski B, Grinspun E, Gross MH. Example-based elastic materials. ACM Trans Graph 2011;30(4):72. [3] Schumacher C, Thomaszewski B, Coros S, Martin S, Sumner RW, Gross MH. Efficient simulation of example-based materials. In: Symposium on computer animation; 2012. p. 1–8. [4] Terzopoulos D, Platt J, Barr A, Fleischer K. Elastically deformable models. In: Proceedings of ACM SIGGRAPH; 1987. p. 205–14. [5] Gibson SFF, Mirtich B. A survey of deformable modeling in computer graphics. Technical Report TR-97-19, MERL, Cambridge, MA; 1997. [6] Witkin A. Physically based modeling: principles and practice – constrained dynamics. In: Computer graphics; 1997. p. 11–21. [7] Nealen A, Müller M, Keiser R, Boxerman E, Carlson M. Physically based deformable models in computer graphics. Comput Graph Forum 2006;25 (4):809–36. [8] Hauth M, Straßer W. Corotational simulation of deformable solids. In: WSCG; 2004. p. 137–44.

Max iterations in Algorithm 1

5 6 3 3 4 2

Average time cost (ms) Update A=At

Linear system

Whole step

49 68 27 2 67 2

103 320 63 3 494 3

168 434 182 19 581 6

[9] Müller M, Dorsey J, McMillan L, Jagnow R, Cutler B. Stable real-time deformations. In: Proceedings of ACM SIGGRAPH symposium on computer animation. ACM New York, NY, USA; 2002. p. 49–54. [10] Müller M, Gross M. Interactive virtual materials. In: Proceedings of graphics interface 2004; 2004. p. 239–46. [11] Irving G, Teran J, Fedkiw R. Tetrahedral and hexahedral invertible finite elements. Graph Models 2006;68(2):66–89. [12] Huang J, Tong Y, Zhou K, Bao H, Desbrun M. Interactive shape interpolation through controllable dynamic deformation. IEEE Trans Vis Comput Graph 2011;17(7):983–92. [13] Alexa M, Cohen-Or D, Levin D. As-rigid-as-possible shape interpolation. In: SIGGRAPH; 2000. p. 157–64. [14] Alexa M. Linear combination of transformations. In: ACM SIGGRAPH; 2002. p. 380–7. [15] Lipman Y, Sorkine O, Levin D, Cohen-Or D. Linear rotation-invariant coordinates for meshes. ACM Trans Graph 2005;24(3):479–87. [16] Xu D, Zhang H, Wang Q, Bao H. Poisson shape interpolation. Graph Models 2006;68(3):268–81. [17] Park FC, Ravani B. Smooth invariant interpolation of rotations. ACM Trans Graph 1997;16(3):277–95. [18] Winkler T, Drieseberg J, Alexa M, Hormann K. Multi-scale geometry interpolation. Comput Graph Forum 2010;29(2):309–18. [19] Kilian M, Mitra NJ, Pottmann H. Geometric modeling in shape space. ACM Trans Graph 2007;26(3):64. [20] Botsch M, Sorkine O. On linear variational surface deformation methods. IEEE Trans Vis Comput Graph 2008;14(1):213–30. [21] Sumner RW, Zwicker M, Gotsman C, Popovic J. Mesh-based inverse kinematics. ACM Trans Graph 2005;24(3):488–95. [22] Der KG, Sumner RW, Popovic J. Inverse kinematics for reduced deformable models. ACM Trans Graph 2006;25(3):1174–9. [23] Coros S, Martin S, Thomaszewski B, Schumacher C, Sumner RW, Gross MH. Deformable objects alive!. ACM Trans Graph 2012;31(4):69. [24] Koyama Y, Takayama K, Umetani N, Igarashi T. Real-time example-based elastic deformation. In: Symposium on computer animation; 2012. p. 19–24. [25] Zienkiewicz OC. The finite element method. third ed. UK: McGraw-Hill (U.K.) Ltd.; 1977. [26] Taylor G, Bailey C, Cross M. A vertex-based finite volume method applied to non-linear material problem in computational solid mechanics. Int J Numer Methods Eng 2003;56(1):507–29. [27] Bathe KJ. Finite element procedures in engineering analysis. Prentice-Hall, Inc., New Jersey, US; 1982. [28] Liu T, Bargteil AW, O'Brien JF, Kavan L. Fast simulation of mass–spring systems. ACM Trans Graph 2013;32(6) 209:1–7. Proceedings of ACM SIGGRAPH Asia 2013, Hong Kong. [29] Sorkine O, Alexa M. As-rigid-as-possible surface modeling. In: Proceedings of the fifth Eurographics symposium on geometry processing, SGP '07, Aire-laVille, Switzerland. Switzerland: Eurographics Association; 2007. p. 109–16, ISBN 978-3-905673-46-3. [30] Mezger J, Thomaszewski B, Pabst S, Straßer W. Interactive physically-based shape editing. In: Symposium on solid and physical modeling; 2008. p. 79–89. [31] Ju T, Schaefer S, Warren JD. Mean value coordinates for closed triangular meshes. ACM Trans Graph 2005;24(3):561–6.