From the proceedings of ACM SIGGRAPH 2002
Interactive Skeleton-Driven Dynamic Deformations
Steve Capell Seth Green Brian Curless Tom Duchamp Zoran Popovic ́ University of Washington
Abstract
This paper presents a framework for the skeleton-driven animation of elastically deformable characters. A character is embedded in a coarse volumetric control lattice, which provides the structure needed to apply the finite element method. To incorporate skele- tal controls, we introduce line constraints along the bones of sim- ple skeletons. The bones are made to coincide with edges of the control lattice, which enables us to apply the constraints efficiently using algebraic methods. To accelerate computation, we associate regions of the volumetric mesh with particular bones and perform locally linearized simulations, which are blended at each time step. We define a hierarchical basis on the control lattice, so for detailed interactions the simulation can adapt the level of detail. We demon- strate the ability to animate complex models using simple skeletons and coarse volumetric meshes in a manner that simulates secondary motions at interactive rates.
CR Categories: I.3.5 [Computer Graphics]: Computational Ge- ometry and Object Modeling—Physically Based Modeling I.3.7 [Computer Graphics]: Three-Dimensional Graphics and Realism— Animation;
Keywords: animation, deformation, physically-based animation, physically-based modeling
1 Introduction
Physical simulation is central to the process of creating realis- tic character animations. In the film industry, animators require detailed control of the motion of their characters, but creating physically-based secondary motions is difficult and time consum- ing to do by hand. Recently, techniques have been developed for automatically simulating these secondary motions. These methods are built atop skin, muscle, and bone models and can generate de- tailed, dynamic motions. However, constructing these models is time consuming, and the simulations are computationally expen- sive.
By contrast, in video game or virtual reality applications where interactivity is critical, character animation is built atop much sim- pler models. The shapes are composed of convenient primitives and are controlled by line segment based skeletons. Deformations of body parts are purely kinematically driven, using, e.g., blended coordinate frames. Incorporating realistic physically-based dynam- ics using the film industry’s approach is currently impractical.
Copyright ⃝c 2002 by the Association for Computing Machinery, Inc. Permission to make digital or hard copies of part or all of this work for personal or classroom use is granted without fee provided that copies are not made or distributed for commercial advantage and that copies bear this notice and the full citation on the first page. Copyrights for components of this work owned by others than ACM must be honored. Abstracting with credit is permitted. To copy otherwise, to republish, to post on servers, or to redistribute to lists, re- quires prior specific permission and/or a fee. Request permission from Permis- sions Dept, ACM Inc., fax +1-212-869-0481 or e-mail permissions@acm.org. ⃝c 2002 ACM 1-58113-521-1/02/0007 $5.00
In this paper, we attempt to bring dynamic simulation into the realm of real-time, skeleton-driven animation. The challenge is to find the right combination of physical principles, geometric mod- eling, computational tools, and simplifying assumptions that yield compelling animations at interactive rates.
Our approach is based on the equations of motion of elastic solids, simulated in a finite element setting. The volumetric finite element mesh need only be specified coarsely, subject to the re- quirement that it encompass the geometric model on which simula- tion will be performed. This last requirement is necessary in order to ensure complete integration over the interior of the object. In fact, as long as the interior of the object is well-defined, simulation of its elastic deformation is possible regardless of the the surface representation or complexity.
The volumetric mesh we choose is not restricted to a regular grid; rather, it is comprised of elements such as tetrahedra and hexahe- dra. This flexibility permits construction of meshes that conform better to the surface of the object, improving simulation quality. In addition, to support adaptive level of detail during simulation, we construct a hierarchical basis, which allows detail to be introduced or removed as needed.
Since our ultimate goal is simulation of skeletally controlled characters, our framework supports line constraints, where lines correspond to bones. In order to incorporate these constraints eas- ily, we require the volumetric mesh to contain edges coincident with the bones. Finally, to achieve interactive rates, we linearize the equations of motion, solve them over volumetric regions associated with each bone, and blend the deformations where regions overlap.
We believe that this work makes a number of contributions. Our crafting of the function space in order to make constraint han- dling easier is, to our knowledge, novel. We introduce blended local linearization of nonlinear equations, in the context of de- formable animated characters. We generalize a method of solving constraints using linear subspace projection. We also introduce a constraint that allows one-dimensional bones to behave as three- dimensional bones. Finally, we believe our most important contri- bution is putting together a collection of techniques that allows us to interactively animate arbitrary shapes with skeletal controls while generating realistic dynamic deformations.
2 Related Work
Probably the most common technique for deforming articulated characters is to define the position of the surface geometry as a function of an underlying skeletal structure or set of control pa- rameters. Recent advances in this area can be found in the work of Lewis et al. [2000], Singh and Kokkevis [2000], and Sloan et al. [2001]. Our work builds on the notion of skeletal control, but within a physically-based framework.
In the late 1980’s, Terzopoulos et al. pioneered the field of physically-based deformable models for computer graphics. Using Lagrangian equations of motion and finite differences they simu- lated elastic [1987] and inelastic [1988] behaviors, combined with a rigid body motion term to compensate for instabilities with stiff bodies [1988].
1
From the proceedings of ACM SIGGRAPH 2002
Much of the research that followed sought to add more sophisti- cated constraint solvers, accelerate the solutions under a variety of approximations, and add stability to permit larger timesteps. Platt and Barr [1988] improved on existing constraint methods (e.g., the penalty method) by introducing reaction constraint and augmented Lagrangian constraint approaches. Applying these constraints to already complex simulations, however, was not a step toward inter- activity.
Pentland and Williams [1989] simplified the problem by solv- ing for the vibrational modes of a body and keeping only the lower frequency modes and obtained realtime simulations of physically plausible deforming bodies. Witkin and Welch [1990] used La- grangian dynamics to solve for low-order polynomial, global de- formations, coupled with constraints enforced through Lagrange multipliers. Baraff and Witkin [1992] later extended this method for better handling of non-penetration constraints. Finally, Metaxas and Terzopoulos [1992] combined Witkin and Welch’s global de- formation framework with local finite element surface deformations and Lagrange multiplier constraints to animate superquadric sur- faces. In each of these examples, while achieving interactive rates, the deformations were substantial approximations to detailed volu- metric deformations and were not demonstrated on complex shapes.
To accelerate computations, some hierarchical methods have also been employed. Terzopoulos et al. [1988] use a multigrid solver for a surface-based inelastic simulation. The approach of Metaxas and Terzopoulos [1992] is analogous to a two level simu- lation that uses global deformations at the coarse level and finite ele- ments for finer surface deformations. More recently, Debunne et al. [1999] built an octtree of particles that interact according to Lame ́’s equation, resulting in interactive simulations. The particles are sim- ulated using an explicit Euler solver that steps each particle adap- tively in time and at differing spatial resolutions. To animate a sur- face, the particles are linked to each surface point by a weighting scheme. Their approach was recently extended to use unstructured tetrahedral hierarchies [Debunne et al. 2001].
To add stability to computations, implicit solvers have proven to be quite effective. Terzopoulos et al. [1987; 1988] used semi- implicit solvers in their initial work. Baraff and Witkin [1998] used an implicit scheme to permit large timesteps in notoriously unstable cloth simulations. Desbrun et al. [1999], also working with cloth models, showed that the implicit solution method acts as a filter that stabilizes stiff systems. They also add a rotation term to preserve angular momentum and a correction term after each time step to simulate nonlinear elasticity.
Free-form deformation (FFD), introduced by Sederberg and Parry [1986], is also closely related to our work. FFD involves embedding an object in a domain that is more easily parameterized than the object itself. The main advantages of FFD are that arbitrary objects can be easily deformed and the space of deformations can be crafted independently of the representation and resolution of the object. Since its introduction, the flexibility of FFD has been im- proved by introducing lattices of arbitrary topology [MacCracken and Joy 1996], and dynamic free-form deformation has been in- troduced to apply FFD to animation [Faloutsos et al. 1997]. Our framework builds on FFD by also embedding the object in a coarse control lattice. But unlike the work of Faloutsos et al., where a diag- onal stiffness matrix was used, we use the principles of continuum elasticity to compute the dynamics of the object being deformed.
Another approach to fast, physically-based deformations is to solve quasi-static solutions, i.e., compute the equilibrium state of the system given forces and constraints, and then animate by ad- justing the forces and constraints over time. Gourret et al. [1989] explored such a technique for volumetric finite elements, and James and Pai [1999] developed an interactive boundary element solution under the assumption of constant material properties inside the vol- ume. Other quasi-static approaches have also been favored for sur-
S
K=K0
Ω
(a)
(b)
2
Figure 1: (a) An object Ω ⊂ R3 instrumented with a skeletal complex S and a control latticeK=K0. FunctionsonΩaredefinedintermsofK,whichformsaneighborhood of the object. A trilinear basis function is associated with each vertex of K (includ- ing the red vertices on the skeleton). (b) The control lattice subdivided once to form K1 (which includes the black, red, and green vertices and edges). Additional basis functions are associated with newly introduced vertices (shown in green).
gical planning and simulation [Bro-Nielsen and Cotin 1996; Koch et al. 1996; Roth et al. 1998]. Quasi-static solutions, however, are approximations that do not capture the true dynamics of motion.
In some cases, interactivity has been achieved without resort- ing to significant simplifications in the dynamic model. Recently, Picinbono et al. [2000] described an interactive surgery simulation using a nonlinear finite element method. They were able to achieve interactive rates for a virtual liver composed of about 2000 tetrahe- dra. However, more complex objects such as the ones used in our work would require many more tetrahedra.
There are also a number of papers that approach the problem of deforming humans and animals using anatomical modeling [Aubel and Thalmann 2000; Wilhelms and Gelder 1997]. They differ from our approach in that we are more concerned with interactivity and the appearance of realism than actual anatomical modeling.
Another possibility for deformation is to perform thin-shell com- putations. Recent work by Cirak and Ortiz [2001] demonstrates the use of subdivision elements for computing the dynamics of thin shells. However, thin shells are insufficient for modeling the inte- rior of solid objects.
In the next section we introduce the basic mathematical and physical formulation that underlies our framework. We then de- scribe our simulation and control methodology (Section 4), discuss results (Section 5), and conclude (Section 6).
3 Formulation
Each object that we wish to animate is represented as a domain Ω ⊂ R3. We make no assumptions about Ω other than we know its interior. We instrument each object with a skeleton, i.e. an anno- tated transformation hierarchy suitable for animation. We refer to each transformation in the hierarchy, as well as its associated origin point, as a joint. The skeleton defines a graph S, whose vertices correspond to joints and whose edges correspond to line segments between two joints that share a parent-child relationship. Joints and bones are located interior to the object, so the graph S is a piecewise linear subset S ⊂ Ω, that we call a skeletal complex (see Figure 1).
Motion of the object is represented by a time dependent function p:Ω×R→R3 : (x,t)→p(x,t). (1)
Let
pS :S×R→R3 (2)
denote the restriction of the map to S. Rigidity of the bones implies that pS is an isometry on each edge of S. In particular, pS is a piecewise linear function on S.
Our goal is to solve for the dynamic motion of the object given the motion of the skeleton. Since we model the object as an elastic
body, the function p(x, t) is then the solution of a system of partial differential equations, subject to the constraint p(x, t) = pS (x, t) for x ∈ S.
To solve the system numerically, we apply the finite element method (see, e.g., [Prenter 1975]). We separate the map p(x, t) into a constant rest state r(x) and a dynamic displacement d(x, t), each of which is represented as a finite sum. The rest state of the object is given by the identity map r : Ω → R3, which has the expansion:
r(x)=raφa(x)=raφa(x)=x (3) a
where the functions φa(x) are elements of a finite basis B, and ra ∈ R3. As demonstrated in the above equation, we use the Einstein summation convention throughout this paper: whenever a term contains the same index as both a subscript and a superscript, the term implies a summation over the range of that index. The displacement is expanded similarly:
d(x, t) = qa(t) φa(x) . (4)
where qa(t) ∈ R3 are the dynamically evolving coefficients that determine the deformation of the object over time. The state of the system is simply the sum of the rest state and the displacement:
p(x, t) = (ra + qa(t))φa(x) (5)
We represent the state of the body at time t as a column vector of generalized coordinates q = q(t) whose a-th component is the coefficient qa(t) in Equation (4), and we model the dynamics of the body as a system of second order ordinary differential equations. The system is obtained by applying the finite element method to the Lagrangian formulation of the equations of elasticity (see, e.g., [Shabana 1998]). In the remainder of the section we describe the basis B and formulate the finite element problem.
3.1 The Hierarchical Basis
In order to allow our simulations to adapt to local conditions, we employ a hierarchical basis. Such bases are well established as use- ful tools for numerical computation (see, e.g., [Bank 1996]). Our construction mirrors that of what are referred to as lazy wavelets in [Stollnitz et al. 1996]. The basis B is defined in terms of re- peated subdivision of a control lattice surrounding the object (see Figure 1). It is desirable that the control lattice conform to the shape of the object while being as coarse as possible. An advantage of us- ing an unstructured lattice instead of a regular grid to define the de- formation function space is that the lattice can be tailored to fit the object. More precisely a control lattice K is a finite union K = ∪iCi of convex cells Ci satisfying the following conditions:
(i) For all i, j, i ̸= j the intersection Ci,j = Ci ∪ Cj is either empty or a face, edge, or vertex of both Ci and Cj.
(ii) TheedgesofSareedgesofcellsofK.
(iii) ThedomainΩiscontainedintheinteriorofK.
(iv) For all i, each vertex of Ci has valence 3 (within Ci).
Condition (iv) still allows a variety of cell shapes including hexa- hedra, tetrahedra and triangular prisms.
We now show how to construct a collection of functions B = {φa} on K whose restriction to Ω is a linearly independent set of continuous functions on Ω. Let V0 ⊂ V1 ⊂ V2 ⊂ . . . be the nested sequence of function spaces described in appendix A.1. The set VJ consists of the piecewise trilinear functions on the complex KJ obtained from K by J hexahedral subdivisions. For each vertex a of K, positioned at xa, let φa denote the unique function in V0 such that φa(xa) = 1 and φa(xb) = 0 for b ̸= a a vertex of K. Include φa
in B if the restriction of φa to Ω is non-zero. Proceed inductively as follows. Let a be a vertex of KJ+1 that is not a vertex of KJ and let φa be the unique function in VJ+1 such that φa(a) = 1 and that vanishes at all other vertices of KJ+1. Include φa in B if its restriction to Ω is non-zero and if its restriction to S is zero.
Although the elements of B are defined on all of K, we are only interested in their values on Ω; we will, therefore, interpret B as collection of functions on Ω. One can show that the set B is linearly independent set of functions, which we call the hierarchical basis.
By construction φa ∈ B for each joint vertex a ∈ S, and the re- striction of φa to S is linear on each bone of S. Moreover if φa ∈ B, for a not a joint vertex, then φa vanishes identically on S. Conse- quently, the function pS (x, t) can be written in the form
pS(x, t) = (ra + qa(t))φa(x) . (6) a∈S
and because φa(a) = 1, the vector (ra + qa(t)) is the location of the joint vertex a at time t.
3.2 Equations of Motion
By virtue of Equation (5), we can express the kinetic energy T and elastic potential energy V as functions of q ̇ and q, respectively, where q ̇ denotes the time derivative of q. The equations of motion are then the Euler-Lagrange equations
d ∂T(q ̇) + ∂V(q) + Qext − μq ̇ = 0 (7) dt ∂ q ̇ ∂ q
The kinetic energy of a moving body is a generalization of the
12 familiar mv :
2
T=1 ρ(x)p ̇·p ̇dΩ=1Mabq ̇a·q ̇b (8) 2Ω2
where ρ(x) is the mass density of the body, and Mab =
From the proceedings of ACM SIGGRAPH 2002
∂ T /∂ q ̇ and ∂ V /∂ q denote gradients with respect to q ̇ and
where
q, respectively. The term Qext is a generalized force arising from external body forces, such as gravity. The last term is a generalized dissipative force, added to simulate the effect of friction. We will now derive each of the first three terms of Equation 7, ultimately yielding a system of ODEs to be solved in generalized coordinates.
Ω
ρ φa φb dΩ. Equation (8) yields the formula d ∂T=Mq ̈.
ab
(9)
is called the mass
The elastic potential energy of a body captures the amount of work required to deform the body from the rest state into the current configuration. It is expressed in terms of the strain tensor and stress tensor. Strain is the degree of metric distortion of the body. A standard measure of strain is Green’s strain tensor:
dt ∂q ̇
The matrix M composed of the elements M matrix. We discuss its computation in Section 3.3.
∂di ∂dj eij = +
+δkl
∂dk ∂dl ∂xi ∂xj
(10)
∂xj ∂xi
The diagonal terms of the strain tensor represent normal deforma- tions while those off the diagonal capture shearing. Forces acting on the interior of a continuum appear in the form of the stress ten- sor, which is defined in terms of strain:
τij = 2G
ν tr(e)δij + eij (11) 1 − 2ν
3
where tr(e) = δijeij. The constant G, called the shear modulus or modulus of rigidity, determines how hard the body resists deforma- tion. The coefficient ν, called Poisson’s ratio, determines the extent to which strains in one direction are related to those perpendicular to it. This gives a measure of the degree to which the body preserves volume. The elastic potential energy V(q), which is analogous to the familiar definition of work as force times distance, is given by the formula
V=G ν tr2(e)+δijδkleikejldΩ (12)
By combining Equations 4, 10, and 12 we can express the elastic potential V and its derivatives (with respect to q) as polynomial functions of q. The coefficients of these polynomials are integrals that can be precomputed. Details are described in appendix A.2.
The matrix S = ∂2V is referred to as the stiffness matrix. ∂q∂q
To add realism, we include the force of gravity in our formula- tion. Gravity is an example of a body force that affects all points inside the body. We treat gravity as a constant acceleration field specified by the vector g. The gravitational potential energy is then the integral
Vg = ρg·p= ρφag·qa. (13) ΩΩ
The generalized gravitational force is the gradient
From the proceedings of ACM SIGGRAPH 2002
Ω
1 − 2ν
∂Vg Qga==ρφag (14)
The above force can be interpreted as the familiar mg except that the mass term represents all of the mass associated with a particular basis function.
3.3 Numerical Integration
In order to compute the gravity terms and the mass and stiffness matrices we precompute the integrals in equations (8), (14), and (29). The integration is done numerically using the following steps:
1. Subdivide K to the desired level for numerical integration.
2. Compute the values of the basis functions at each vertex.
3. Tetrahedralize the domain. After subdividing once, the do- main is composed of only hexahedral cells. We then divide each of these cells into tetrahedra in order to approximate functions on the domain as piecewise linear.
4. Compute the integrals over each domain tetrahedron using piecewise linear approximations to the basis functions. If all four vertices of a tetrahedron fall outside the surface of the object, its contribution to the integrals is neglected.
With the integrals computed, equation (7) can now be solved using a nonlinear Newton-Raphson solver.
4 Skeletal Simulation
The fully nonlinear elastic formulation described in the previous section is computationally expensive, and does not take into consid- eration the skeleton. In this section we introduce a set of techniques, tailored for fast skeleton-driven animation, that approximate the nonlinear dynamics.
Figure 2: The upper left image shows an input model instrumented with a skeleton and local coordinate systems. The upper right image shows the model embedded in (half of) a control lattice. The lower left image shows how the skeleton coincides with edges and vertices of the control lattice. The lower right image shows the entire control lattice, as well as the division of the object into regions for local linearization. Each region is associated with one of the local coordinate systems in the upper left image. Note the color blending where regions overlap.
4.1 Instrumentation
Prior to simulation, a model must be instrumented with a skeleton and control lattice. Although recent work by Teichmann and Teller addresses automated skeleton construction [1998], we currently let the animator specify the skeleton in order to achieve the desired level of control. We have implemented a simple system that allows a skeleton to be constructed manually in just a few minutes. The user creates a joint by clicking on the object with the mouse. If the ray through the mouse point (from the camera projection center) intersects the object at least twice, a joint is placed midway be- tween the first two intersections. This positioning scheme produces joints that are centrally located inside the object. Two joints can be selected to define a bone, and with the selection of a root joint, a transformation hierarchy can be created automatically.
We currently use a constructive procedure that allows the user to build the control lattice interactively by adding cells incremen- tally and repositioning the control vertices as needed. Several hours are required for an experienced user to create a moderately complex control lattice. The abundance of volumetric meshing schemes sug- gests that automatic creation of the control lattice is possible, and we hope to address this problem in the future. Figure 2 shows the skeleton and control mesh for a kangaroo model.
4.2 Solving the System
Due to the computational expense of solving the full nonlinear equations of elasticity, we seek simplifications that make the equa- tions easier to solve. One possibility is to linearize the equations of motion at the beginning of each timestep as was done by Baraff
∂qa Ω
4
and Witkin in their work on cloth simulation [1998]. In our experi- ence, simulations using this method are essentially indistinguish- able from results obtained using a nonlinear implicit method to solve the system, as long as the timestep is not so large as to al- low radical shape change during a single step. After applying their implicit solver to our formulation, the resulting equations are:
∆q = h(q ̇+∆v) (15)
α
x, y, and z components of the 3-vectors ac into the n-vectors a , where α ∈ {x, y, z}, define the matrix Cac = φa (xc ), and separate ∆v into its x, y, and z components ∆vα, equation (20) becomes:
CT∆vα =aα, α∈{x,y,z} (21)
So each constraint requires that ∆v be constant along three par- ticular directions in R3n. Maintaining the constraints involves the following steps:
1. At the beginning of each timestep, ∆v is initialized so that
equation (21) holds. This is accomplished by computing the
QR-decomposition of C and transforming equation (21) into
RT bα = aα, ∆vα = Qbα, from which ∆v can be easily com-
(M−hμI+h2S)∆v
= hμq ̇ − ∂V −Qext −hSq ̇(16) ∂q
where h is the timestep, μ is the damping coefficient, I is the identity matrix, ∆v is the change in the velocity q ̇ during the timestep, ∆q is the change in q during the timestep, M is the mass matrix, and S is the stiffness matrix. All quantities are evaluated at the beginning of the timestep. Equation (16) is a sparse linear system that can be solved for ∆v using a Conjugate Gradients (CG) solver. Then ∆v is substituted into equation (15) to obtain ∆q.
4.3 Bone Constraints
In our framework the skeleton is controlled directly by keyframe data or some other source external to the dynamic simulation. From the viewpoint of the simulation, the skeleton is simply a compli- cated constraint. Because we have restricted the bones to lie along edges in the control lattice, and the basis is interpolating, it is es- pecially easy to handle the bone constraints algebraically. Each control point that lies on a bone corresponds to a component of ∆v that is known a priori, rather than having to be computed. Simpli- fying equation (16) to the form A∆v = b, we can sort the variables into known (∆vk and bk) and unknown quantities (∆vu and bu) and form the following system:
A11 A12 ∆vk =bu (17) A21A22 ∆vu bk
The reason that some components of the vector b are now unknown arises from the fact that the external forces required to enforce the bone constraints are unknown, and they appear on the r.h.s. of equa- tion (16). In order to solve for ∆vu we simply solve the system:
A22∆vu = bk − A21∆vk (18)
The advantage of this approach is that adding skeletal constraints actually reduces the computational cost by shrinking the system that must be solved.
4.4 Linear Subspace Constraints
Because we would like our objects to interact with other objects, position constraints are also important. The framework of Baraff and Witkin [1998] provides an elegant solution for particle systems. During each internal step of a CG solver, they project out certain components of ∆v corresponding to constrained particles. Here we show that this technique can be extended to include position constraints at any point in a continuous body. Position constraints in our framework are of the form:
dc(t) = qaφa(xc) (19) which simply says that the displacement at xc conforms to some
known function dc. Evaluating equation (15) at xc results in: ∆vaφa(xc) = dc(t + h) − dc(t) − q ̇aφa(xc) (20)
h
The r.h.s. of the above equation is simply a constant ac that can be computed at the beginning of each timestep. If we accumulate the
puted. Although QR-decomposition of an n × m matrix re- 2
From the proceedings of ACM SIGGRAPH 2002
quires O(nm ) time, in our case the number of constraints m
is typically small, so the computational cost is low.
2. Each column c of C has an associated projection matrix P = I − ccT /cT c, which, when applied to a vector, eliminates the component in the direction of c. These projectors are applied during CG such that incremental updates to ∆v are orthogonal to the vectors c, ensuring that equation (21) remains true (for
details see [Baraff and Witkin 1998]).
In our current framework, conflicting constraints can be detected during QR-decomposition and removed. In the future we hope to augment this method to solve over-constrained systems more ele- gantly, as was done for FFD by Hsu et al. [1992].
4.5 Blended Local Linearization
A major bottleneck in our system is the computation of the stiff- ness matrix at the beginning of each timestep (the elastic potential is a quartic function of q). A well-known simplification is to lin- earize the strain tensor by dropping the last term in equation (10), which results in a quadratic elastic potential and thus a constant stiffness matrix (which is composed of the first three addends in equation (28)). As compared to other simplifications such as us- ing a mass-spring-based elastic potential, linearization of strain has the advantage that it is a very good approximation, but only when the deformation is small; for large deformations, severe distortions occur.
A notable case for the linear strain model is when the object undergoes a large rigid rotation, coupled with a small deformation. While the elastic potential based on nonlinear strain does not penal- ize rotations, the linear strain model does, while failing to penalize certain shearing deformations. Terzopoulos et al. [1988] addressed this case by modeling the deformation relative to a frame of refer- ence that follows the gross motion of the object. Since the relative deformation is assumed to be small, the linearized strain is a rea- sonable approximation. This approach is common practice in the engineering literature, such as in the textbook of Shabana [1998], in which multibody systems composed of interconnected parts are considered. In such systems, the deformation of each part can be measured from a local reference configuration that factors in the ro- tation of the part. As long as the deformation of each part is small relative to its rotated reference configuration, the linear strain model is a good approximation.
To apply these ideas to articulated characters, we first recognize that the soft tissues of vertebrates do not typically undergo large deformations relative to nearby bones. Based on this assertion, our approach is to divide the object into regions, each of which can be simulated using the linear strain model.
The user divides the object into regions by assigning weights to the control vertices, forming a partition of unity over the object. A piece of the object can belong to a single region or can be di- vided fractionally among several regions. We encode the weights
5
From the proceedings of ACM SIGGRAPH 2002 ii
for region i in a diagonal square matrix W , where Waa is the weight associated with vertex a in region i. The lower right image in Fig- ure 2 shows a partitioned object, colored according to the region assignments. Our current system requires that the user select indi- vidual weights for each control vertex, but a more intuitive painting interface would be straightforward to implement. It would also be helpful to automate the task of region assignment (recent work by Li et al. [2001] may be adaptable to our problem domain).
i
From the region assignments we form a cell complex K corre-
sponding to region i:
Ki={C∈K:∀va∈v(C),Wi >0} (22)
aa
where v(C) is the set of control vertices on cell C. Each region has an associated function space:
Bi ={φa|Ki :φa ∈B, φa|Ki ̸=0} (23) where φa|Ki denotes the restriction of φa to Ki. We define a rect-
i
angular matrix Q to select the basis functions that have nonzero restrictions to region i. The element Qiab = 1 if and only if φa ∈ B corresponds to φb ∈ Bi. The pseudocode for taking a single simu- lation step is:
Figure 3: The left image shows the kangaroo at rest. Brown spheres represent active basis functions. The cyan sphere represents a position constraint. On the right, the position constraint has been moved causing adaptation of the basis. The red spheres in the right image represent newly introduced detail coefficients.
denote this region Ωβ ⊂ Ω. The following potential describes the constraint:
1
2 Ωβ
The above potential is quadratic, so its Hessian is simply a constant
that can be added to the stiffness matrix: 1 [ri, qi, q ̇i] ⇐ [Qir, Qiq, Qiq ̇] ∂2U
U =
d · ddΩ (25)
foreach region i do foreach a do
=I φaφb (26) aaaa Ωβ
2 qi :=qi −Ti(ri)+ri
∂qa∂qb
end
3 Construct Ai and bi from equation (18)
4 Solve Ai∆vi = bi
end
5 ∆v⇐iWiQiT∆vi
6 q ̇⇐q ̇+∆v
7 q ⇐ q + h q ̇
Line 1 extracts the regional variables from the global system.
Line 2 converts qi so that it corresponds not to displacement from the rest state, but to displacement from the rest state transformed according to the transformation of the bone coordinate system. The homogeneous transformation Ti, extracted from the current config- uration of the skeleton, represents the transformation of the bone from its rest position to its current position. But it is not enough to simply transform qi, because the transformation itself must be sub- tracted from qi. A displacement field dT that transforms the object x according to the transformation T(x), has the following form:
dT +x=T(x) (24)
It is from the above expression that line 2 is derived. Line 3 builds
the linear system required to solve for the local equations of mo-
tion, including the extraction of bone constraints, and line 4 solves
the linear system using CG. Line 5 merges the solutions from each
i
region, each weighted according to the user-assigned weights in W . Finally, the state of the global system is updated in lines 6 and 7.
4.6 Twist Constraint
In natural creatures with three-dimensional bones, the flesh cannot twist (i.e. rotate) around the axis of the bone without causing the flesh to deform. Such deformations are resisted by emergent elas- tic forces, so the twisting is limited. But flesh can rotate about a line constraint without deforming. To avoid such unnaturally free movement, we introduce a soft constraint to penalize all displace- ment (not just deformation) within a fixed radius of the bones. We
where I is a 3 × 3 identity matrix. The above constraint must be computed relative to the rigidly transformed bone, which fits well into our local computation framework.
4.7 Adaptation
Because we use a hierarchical basis, our simulator can add detail where needed. We apply the simple heuristic that detail is more helpful where there are large deformations (similar to, e.g., [De- bunne et al. 2001]). If the object is sufficiently deformed over the support of a particular basis function, then all of the basis functions in the next finer level with support overlapping the area of high dis- tortion are introduced into the simulation. Likewise, basis functions are removed when there is little deformation in their support. Each level of the basis has an associated threshold for determining when to refine and another for determining when to coarsen. As noted in [Debunne et al. 2001], a lower threshold is required for coarsen- ing than refining in order to prevent the simulation from oscillating between levels of resolution.
Regardless of the criteria employed, adapting the basis is straightforward in our framework. Most of this simplicity comes from refining the basis, not the geometry, as was done by Gortler and Cohen [1995] and recently generalized by Grinspun et al. [2002]. For some fixed number of basis levels we precompute the mass and stiffness matrices and store them in a sparse data struc- ture. Adapting the basis simply corresponds to extracting and re- linquishing certain components from these matrices, which can be done very quickly. The resultant subsets of the basis are linearly independent regardless of which basis functions we choose. Fig- ure 3 shows adaptation of the kangaroo model. For more details regarding our adaptation methodology see [Capell et al. 2002].
5 Results
The accompanying video shows the results of applying our frame- work to two triangle meshes that we acquired from the Internet.
6
Figure 4: Frames from an interactive animation. There is no noticeable warping due to strain linearization, and the different materials (e.g., ears, horns) behave distinctly.
Figure 5: On the left is the global linear solution, which shows significant warping when the cow turns its head to one side. In the center is the fully nonlinear solution. On the right is the blended local linear solution, which shows no noticeable warping of the head. A slight protrusion can be seen in the neck of the right image due to region blending.
The control mesh for the kangaroo model has 448 cells and 177 vertices; the cow control mesh has 572 cells and 214 vertices. On a 1 Ghz PC, both the cow and kangaroo animated at about 100 Hz using only the coarse basis functions, which is clearly within range for interactive applications (with adaptation, simulation time varies depending on the degree of adaptation required). Figure 4 shows frames of an animation of the cow model (using the coarse basis), which demonstrates the ability of our system to handle variable ma- terial properties; the ears flop around realistically while the horns stay rigid. This feature is possible to do interactively because the control mesh can be carefully crafted to respect material bound- aries, and because our computation of the stiffness matrix takes variable material properties into account.
For our datasets, the blended local linear and global linear so- lutions required about the same amount of computation time. Yet the blended local linear solution produced much more pleasing re- sults, as demonstrated in Figure 5. The blended local linear solution looks similar to the fully nonlinear solution, while the global linear solution is badly warped.
6 Conclusion
We have introduced a method for interactive simulation of de- formable bodies controlled by an underlying skeleton. By choosing a volumetric mesh that aligns with the bones, we are able to meet the bone constraints rapidly. We extend a fast constraint solver that works directly within an iterative solver. We also introduce a twist constraint that mimics the effects of three-dimensional bones when only one-dimensional bones are being modeled. Our method per- forms with the speed of simple linear-strain models of elasticity, but does not suffer from distortions arising from global linearity.
There are many avenues for future work. We would like to au- tomatically generate skeletons and especially control lattices, the latter being the most labor intensive aspect of our framework. Our assumptions about small deformations break down near the joints.
It may be possible to address this problem by using nonlinear elas- ticity near the joints. The deformations near joints might also be improved by specifically tailoring adaptation to the problem. Fi- nally, it would be convenient to include dynamic, not just fully con- strained, bones.
Acknowledgments The authors would like thank Chris Twigg, Mira Dontcheva and Samantha Michel for their instrumental work in creating and converting Maya skeletal animations, and Shawn Bonham and Sean Smith for additional help. This work was sup- ported by the Animation Research Labs, Microsoft Research, NSF grants DMS-9803226 and CCR-0092970, and an Intel equipment donation.
A Appendix
A.1 Review of Trilinear Functions
A trilinear function on the standard unit cube C3 = {x = (x, y, z) : 0 ≤ x, y, z ≤ 1} is a function of the form
From the proceedings of ACM SIGGRAPH 2002
7
(i) (ii) (iii)
The value of f at the centroid of C 3
3
3
is the average of its values
f(x,y,z) = a0 +a1x+a2y+a3z+a4xy+a5xz+a6yz+a7xyz. 3
The function f is determined by its values at the vertices of C : let φˆ (s) denote the hat function
1−|s| for|s|≤1
φˆ ( s ) = .
0 for |s| > 1. and let φ0(x, y, z) = φˆ(x)φˆ(y)φˆ(z). Then
f(x) =
0≤i,j,k≤1
fi,j,kφ0(x−i,y−j,z−k)
fi,j,k
The value of f at the midpoint of an edge of C3 is the average of its values at the endpoints of the edge.
The value of f at the centroid of a face of C is the average of its values at the corners of the face.
= f (i, j, k). It is easy to
where
satisfy the following interpolation or hexahedral subdivision rules:
at the eight vertices of C .
If we subdivide the unit cube into 8 sub-cubes in the standard
way, we can use these subdivision rules to determine the value of f
at the vertices of each sub-cube. Repeatedly subdividing and apply-
ing the subdivision rules yields the value of f at each diadic point
(i/2J,j/2J,k/2J) of C3. Because the diadic points are dense in C3,
the subdivision rules completely determine f from its values at the 3
vertices of C . More generally, starting with values of a function at the vertices of the standard cubic tiling of R3 and applying the subdivision rules to each cubic cell determines a piecewise trilinear function on R3.
We can generalize this construction to define piecewise trilinear functions on any control lattice in which the vertices of each 3-cell of K have valence 3. Starting with the values of f at the vertices of K, we infer its values at the centroid of every edge, face and 3-cell of K. This gives values of f at every vertex of the refined complex K1 obtained by subdivision (see [MacCracken and Joy 1996] for details). Because the vertices of each 3-cell of K have valence 3, the subdivided complex K1 has only hexahedral cells, so after one subdivision, the subdivision process behaves just as for cubes in R3.
There is a corresponding nested sequence of function spaces V0 ⊂V1 ⊂V2 ⊂…
check
that trilinear
functions
defined on K. To define VJ , subdivide J-times to obtain the complex KJ and specify values at each vertex of KJ. The subdivision rules then determine a function on all of K. Thus, each function in VJ, for J = 0,1,2…, is determined by its values at the vertices of K.
A.2 Derivatives of Elastic Potential
The gradient and Hessian of V from equation (12) are: ca ac ac
∂V ∂qc
∂2V ∂qc∂qf
FALOUTSOS, P., VAN DE PANNE, M., AND TERZOPOULOS, D. 1997. Dynamic free- form deformations for animation synthesis. IEEE Transactions on Visualization and Computer Graphics 3, 3 (July–Sept.), 201–214.
GORTLER, S. J., AND COHEN, M. F. 1995. Hierarchical and variational geometric modeling with wavelets. Symposium on Interactive 3D Graphics, 35–42.
GOURRET, J.-P., THALMANN, N. M., AND THALMANN, D. 1989. Simulation of object and human skin deformations in a grasping task. Computer Graphics (Pro- ceedings of SIGGRAPH 89) 23, 3 (July), 21–30.
GRINSPUN, E., KRYSL, P., AND SCHRO ̈ DER, P. 2002. Charms: A simple framework for adaptive simulation. To appear in the Proceedings of SIGGRAPH 2002.
HSU, W. M., HUGHES, J. F., AND KAUFMAN, H. 1992. Direct manipulation of free-form deformations. Computer Graphics (Proceedings of SIGGRAPH 92) 26, 2 (July), 177–184.
JAMES, D. L., AND PAI, D. K. 1999. Artdefo – accurate real time deformable objects. Proceedings of SIGGRAPH 99 (August), 65–72.
KOCH, R. M., GROSS, M. H., CARLS, F. R., VON BU ̈REN, D. F., FANKHAUSER, G., AND PARISH, Y. 1996. Simulating facial surgery using finite element methods. Proceedings of SIGGRAPH 96 (August), 421–428.
LEWIS, J. P., CORDNER, M., AND FONG, N. 2000. Pose space deformation: A unified approach to shape interpolation and skeleton-driven deformation. In Pro- ceedings of SIGGRAPH 2000, 165–172.
LI, X., WOON, T. W., TAN, T. S., AND HUANG, Z. 2001. Decomposing polygon meshes for interactive applications. In ACM Symposium on Interactive 3D Graph- ics, 35–42.
MACCRACKEN, R., AND JOY, K. I. 1996. Free-form deformations with lattices of arbitrary topology. Computer Graphics (Proceedings of SIGGRAPH 96) 30, 181– 188.
METAXAS, D., AND TERZOPOULOS, D. 1992. Dynamic deformation of solid prim- itives with constraints. Computer Graphics (Proceedings of SIGGRAPH 92) 26, 2 (July), 309–312.
PENTLAND, A., AND WILLIAMS, J. 1989. Good vibrations: Modal dynamics for graphics and animation. Computer Graphics (Proceedings of SIGGRAPH 89) 23, 3 (July), 215–222.
PICINBONO, G., DELINGETTE, H., AND AYACHE, N. 2000. Real-time large dis- placement elasticity for surgery simulation: Non-linear tensor-mass model. In Pro- ceedings of the Third International Conference on Medical Robotics, Imaging and Computer Assisted Surgery: MICCAI 2000, 643–652.
PLATT, J. C., AND BARR, A. H. 1988. Constraint methods for flexible models. Computer Graphics (Proceedings of SIGGRAPH 88) 22, 4 (August), 279–288.
PRENTER, P. M. 1975. Splines and Variational Methods. John Wiley and Sons. ROTH, S. H. M., GROSS, M. H., TURELLO, S., AND CARLS, F. R. 1998. A
bernstein-be ́zier based approach to soft tissue simulation. Computer Graphics Fo- rum 17, 3, 285–294.
SEDERBERG, T. W., AND PARRY, S. R. 1986. Free-form deformation of solid geo- metric models. Computer Graphics (Proceedings of SIGGRAPH 86) 20, 4 (Aug.), 151–160.
SHABANA, A. 1998. Dynamics of Multibody Systems. Cambridge University Press. SINGH, K., AND KOKKEVIS, E. 2000. Skinning characters using Surface-Oriented
Free-Form deformations. In Proceedings of the Graphics Interface 2000, 35–42. SLOAN, P.-P. J., ROSE, C. F., AND COHEN, M. F. 2001. Shape by example. In
Symposium on Interactive 3D Graphics, 135–144.
STOLLNITZ, E. J., DEROSE, T. D., AND SALESIN, D. H. 1996. Wavelets for Com-
puter Graphics: Theory and Applications. Morgan Kaufmann, San Francisco, CA. TEICHMANN, M., AND TELLER, S. 1998. Assisted articulation of closed polygonal
models. In Computer Animation and Simulation ’98, 87–101.
TERZOPOULOS, D., AND FLEISCHER, K. 1988. Modeling inelastic deformation: Vis- coelasticity, plasticity, fracture. Computer Graphics (Proceedings of SIGGRAPH 88) 22, 4 (August), 269–278.
TERZOPOULOS, D., AND WITKIN, A. 1988. Physically based models with rigid and deformable components. IEEE Computer Graphics and Applications 8, 6 (Nov.), 41–51.
TERZOPOULOS, D., PLATT, J., BARR, A., AND FLEISCHER, K. 1987. Elastically deformable models. Computer Graphics (Proceedings of SIGGRAPH 87) 21, 4 (July), 205–214.
WILHELMS, J., AND GELDER, A. V. 1997. Anatomically based modeling. In Pro- ceedings of SIGGRAPH 97, 173–180.
WITKIN, A., AND WELCH, W. 1990. Fast animation and control of nonrigid struc- tures. Computer Graphics (Proceedings of SIGGRAPH 90) 24, 4 (August), 243– 252.
2A1 qa + A2 qa + B qa
+2qd qa · Cadc + (qa · qb) Ccab 11
2Acf +Afc +IBfc +2Iqa ·Cafc 121
+2qd ⊗Cfdc +2Ccfb ⊗qb +I qd ·Cfcd 112
acf cfd caf +qa⊗C2 +I qd·C2 +qa⊗C2
fac bfc abcf +C2 ⊗qa +C2 ⊗qb +I(qa ·qb)D1
+2(qd ⊗qa)Dafcd +I(qd ·qe)Dfdce 12
+(qa ⊗qd)Dadcf +(qa ⊗qe)Dafce 22
=acdcadbac(27)
From the proceedings of ACM SIGGRAPH 2002
+qa qd ·C2 +qa qd ·C2
+qd (qa ·qb)Dabcd +qa (qd ·qe)Dadce
=
(28)
+(qa ·qb)C2 12
where I is a 3 × 3 identity matrix and Aab= 4Gν∂φa⊗∂φbdΩ
Ω 1−2ν ∂x ∂x 2 Ω ∂x ∂x
1
Aab= 4G∂φa ⊗∂φbdΩ
ab ∂φa ∂φb
B = 4G · dΩ
Ω ∂x ∂x
abc 4Gν∂φa∂φb ∂φc
C1 =Ω1−2ν∂x ∂x·∂x dΩ abc ∂φa ∂φb ∂φc
(29)
Note that Aab is a 3 × 3 matrix, Cabc is a 3-vector, and Bab and Dabcd iii
are scalar quantities.
References
AUBEL, A., AND THALMANN, D. 2000. Realistic deformation of human body shapes. In Proceedings of Computer Animation and Simulation 2000, 125–135.
BANK, R. E. 1996. Hierarchical bases and the finite element method, vol. 5 of Acta Numerica. Cambridge University Press, Cambridge, 1–43.
BARAFF, D., AND WITKIN, A. 1992. Dynamic simulation of non-penetrating flexible bodies. Computer Graphics (Proceedings of SIGGRAPH 92) 26, 2, 303–308.
BARAFF, D., AND WITKIN, A. 1998. Large steps in cloth simulation. In Proceedings of SIGGRAPH 98, 43–54.
BRO-NIELSEN, M., AND COTIN, S. 1996. Real-time volumetric deformable models for surgery simulation using finite elements and condensation. Computer Graphics Forum (Proceedings of Eurographics ’96) 15, 3, 57–66.
CAPELL, S., GREEN, S., CURLESS, B., DUCHAMP, T., AND POPOVIC ́, Z. 2002. A multiresolution framework for dynamic deformations. University of Washington, Department of Computer Science and Engineering, Technical Report 02-04-02.
1
CIRAK, F., AND ORTIZ, M. 2001. Fully c -conforming subdivision elements for finite
deformation thin-shell analysis. International Journal for Numerical Methods in Engineering 51, 7 (July), 813–833.
DEBUNNE, G., DESBRUN, M., BARR, A., AND CANI, M.-P. 1999. Interactive multiresolution animation of deformable models. Eurographics Workshop on Ani- mation and Simulation.
DEBUNNE, G., DESBRUN, M., CANI, M.-P., AND BARR, A. H. 2001. Dynamic real-time deformations using space & time adaptive sampling. In Proceedings of SIGGRAPH 2001, 31–36.
DESBRUN, M., SCHRO ̈ DER, P., AND BARR, A. 1999. Interactive animation of struc- tured deformable objects. Graphics Interface ’99 (June), 1–8.
C2 =Ω4G∂x ∂x·∂x dΩ
abcd 4Gν ∂φa ∂φb ∂φc ∂φd
D1 = Ω 1−2ν ∂x · ∂x ∂x · ∂x 2 Ω ∂x ∂x ∂x ∂x
dΩ Dabcd = 4G∂φa · ∂φb ∂φc · ∂φd dΩ
8