VIPER: Volume Invariant Position-based Elastic Rods

06/12/2019 ∙ by Baptiste Angles, et al. ∙ 0

We extend the formulation of position-based rods to include elastic volumetric deformations. We achieve this by introducing an additional degree of freedom per vertex – isotropic scale (and its velocity). Including scale enriches the space of possible deformations, allowing the simulation of volumetric effects, such as a reduction in cross-sectional area when a rod is stretched. We rigorously derive the continuous formulation of its elastic energy potentials, and hence its associated position-based dynamics (PBD) updates to realize this model, enabling the simulation of up to 26000 DOFs at 140 Hz in our GPU implementation. We further show how rods can provide a compact alternative to tetrahedral meshes for the representation of complex muscle deformations, as well as providing a convenient representation for collision detection. This is achieved by modeling a muscle as a bundle of rods, for which we also introduce a technique to automatically convert a muscle surface mesh into a rods-bundle. Finally, we show how rods and/or bundles can be skinned to a surface mesh to drive its deformation, resulting in an alternative to cages for real-time volumetric deformation.



There are no comments yet.


page 1

page 3

page 4

page 5

page 8

page 9

page 10

page 13

This week in AI

Get the week's most popular data science and artificial intelligence research sent straight to your inbox every Saturday.

1. Introduction

In recent years, the computer graphics community has invested exceptional efforts in adapting the (non real-time) physical simulation algorithms at the core of cinematic special effects (e.g.: [Ziva Dynamics, 2018]) to the realm of (real-time) interactive applications (e.g.: games, AR/VR). Many of these advancements have been possible thanks to a new class of physics solver, pioneered by Müller et al. [2007], realized on top of Verlet-class integrators [Bender et al., 2015]. These position-based solvers are capable of elegantly modeling constrained Newtonian dynamics, including rigid-body, cloth, ropes, rods and fluids in a unified framework. A primary example is the NVIDIA FLEX system [Macklin et al., 2014], capable of modeling complex and varied physical phenomena in real-time by leveraging modern GPU hardware.

Rods with volume

Within this technological landscape, of particular relevance to our work is the modeling of elastic “rods” [Pai, 2002; Spillmann and Teschner, 2007; Umetani et al., 2014; Kugelstadt and Schömer, 2016]. These models extend “ropes” by augmenting each segment composing the rod with an orthogonal coordinate frame, hence allowing the modeling of torsion on top of stretching/bending. Our VIPER rods extend these formulations by accounting for volume preservation, a phenomena not modeled by existing position-based rod models. Many interesting phenomena require this constraint (e.g. soft-bodies, fluids). For example, water is the largest constituent of most animal tissues ( in muscles) hence modeling quasi-incompressible phenomena is of critical importance to achieve believable motion. We address this problem by adding a per-vertex scaling degree of freedom – a measure of the local rod cross-section – and optimizing for this quantity within the physics solve. Our rod segments are hybrid surface representations, they are explicitly parameterized by the position and scale/radius of their vertices, but their surface is defined implicitly. This hybrid structure makes them particularly well suited for efficient collision detection/resolution [Green, 2010].

Anatomical modeling

Such a physical model not only satisfies our fixation in efficiently simulating rubber bands, but has immediate applications towards the modeling of muscles. As illustrated in Figure 3, striated skeletal muscles222From now on, we will refer to “striated skeletal muscle” simply as “muscles”. in human bodies can be represented as a collection of fibers surrounded by connective tissue (i.e. fasciae). Simulating these muscle types efficiently is a fundamental problem, as they represent from to of the average human body mass. In this paper, we propose to efficiently model muscles as a structured collection of volume-preserving rods. This new model can also be interpreted as a generalization of the static volumetric primitives in Implicit Skinning [Vaillant et al., 2013], where skin can then be efficiently modeled as a triangular mesh sliding on an implicit iso-surface defined by our fibers. The VIPER primitive is designed to be integrated with other existing components to produce a complete character representation. These include representations of the skeleton, fat, skin, etc.

Volumetric simulation

In the industry, volume-preserving simulation is typically performed by discretizing the interior of the object with tetrahedra or using an approximating cage. However, to the best of our knowledge, even modestly sized models require a lengthy preprocessing (e.g. a large scale eigen-decomposition for computing the deformation modes; see [Barbič and James, 2005]) before real-time simulation becomes possible. For example, while visually striking, computing the simulation in Figure 3 requires a two-pass optimization (\⃝raisebox{-0.6pt}{1} muscle, \⃝raisebox{-0.6pt}{2} skin). \⃝raisebox{-0.6pt}{1} Muscles are discretized with tets sharing vertices, and each of their steps of simulation requires seconds. \⃝raisebox{-0.6pt}{2} Skin ( vertices) is simulated as cloth layered on top of fat (having tets sharing vertices), where each of the necessary substeps of simulation requires seconds of compute. Overall, this cumulates to minutes of compute/frame. Offline simulation can be exploited to learn sub-spaces, which then enables dynamic deformations in real-time; see [Xu and Barbič, 2016]. However, once learnt, the dynamic behavior of the model is “baked”. Clearly, this is an obstacle towards the ultimate goal of truly interactive physics simulation, and, consequently, interactive modeling. The VIPER primitive not only allows the real-time volumetric simulation of complex anatomical structures, but also provides a viable alternative to cages as a compact control structure for soft-body deformations.

Automatic fiber-bundle modeling

Rod-based representations of muscle fascia are not commonly available – typical asset databases contain tetrahedral mesh models instead. While artists sometimes model main characters to the level of interior muscles, this effort is expensive and not justified for background characters. Thus we also introduce a technique to convert existing assets with minimal user intervention. Our solution builds fiber bundles by first creating a set of slices through the muscle then performing an iterative optimization to interpolate these slices with a given number of rods such that they approximate the input surface well.


Our fundamental contribution is the design of a novel real-time physics engine for soft-body dynamics. Our system presents several sub-contributions:

  • We enrich position-based solvers by introducing a new volume-preserving cosserat rods model and associated constraints.

  • We demonstrate how these primitives, when assembled into fiber-bundles, are effective in efficiently modeling muscles.

  • We propose a technique for conveniently creating fiber-bundles models from existing simulation assets.

  • We introduce the use of VIPER rods as efficient deformation proxies for soft-body deformation.

2. Related Work

We overview the literature from different angles. We recap example-based modeling frameworks that are commonly used in digital production, as well as recent efforts towards the use of simulated anthropomorphic models. We also review methods that attempt to “emulate” them via geometric processes, and finally processes to calibrate a given model to a target.

Example-based deformation

Digital characters are often modeled via their skin (i.e. skinning), with no consideration of underlying volumetric structures, often resulting in non-physically realistic effects such as the candy wrapper problem of (LBS) Linear Blend Skinning [Jacobson et al., 2014, Fig.3]. While these artifacts can be resolved [Kavan et al., 2007; Le and Hodgins, 2016], skinning solutions lack details such as tendons, muscle bulges, wrinkles, and volume preservation. Example-based approaches such as (PSD) Pose-Space Deformation [Lewis et al., 2000; Kurihara and Miyata, 2004] and BlendShapes [Lewis et al., 2014] interpolate artist-sculpted shapes to emulate all these effects. However, producing effects such as skin sliding over underlying structures, or the collision effects visible parts of the body press together, can require considerable skill and weeks to months of sculpting depending on the required quality [Yuen, 2018]. Such data-driven models can be learnt from measurements for both static [Loper et al., 2015], as well as dynamic [Pons-Moll et al., 2015] humans, but they hardly generalize outside of their corresponding training domains.

[width=] /fig.pdf

Figure 4. (top) The parameterization of a VIPER rod, and its discretization. (left) Its rest configuration, and (right) a deformed configuration.

Physically-based anthropomorphic models

Physically-based simulation of characters has a long history [Terzopoulos and Waters, 1990; Scheepers et al., 1997; Sifakis et al., 2005] but, due to its high computational cost, it has only recently began to see practical use. For skeletal muscle deformation, Lee et al. [2010] provides an excellent overview of the field, in regards to which Saito et al. [2015], with its ability to reach near-interactive runtime, can be considered the state-of-the-art. Similarly to blendshape generation, training data can be exploited to generate efficient low-dimensional simulations [Xu and Barbič, 2016; Schumacher et al., 2012; Bouaziz et al., 2014; Ichim et al., 2017], enabling physically-based digital characters in production settings [Clutterbuck and Jacobs, 2010; Ziva Dynamics, 2018]. A shortcoming of these methods is the requirement that the training set encompasses samples of all configurations of the object/character that will be needed. Physically based approaches also permit a decomposition into layers – skin, muscle, fat, and bones – enabling appropriate algorithms to be used for each. Highly relevant to our work is the simulation of skin layered over volumetric primitives pioneered by Li et al. [2013], and its realization in commercial software [Vital Mechanics, 2018], as well as other existing research [Saito and Yuen, 2017], industry [Ziva Dynamics, 2018], and proprietary solutions [Clutterbuck and Jacobs, 2010] to this complementary problem. Recently, Romeo et al. [2018] proposed the use of PBD for muscle simulation, but its s/frame of processing time makes it unsuitable to interactive applications.

Non physically-based anthropomorphic models

Recent efforts have been made to create alternative representations for sub-skin volumetric models (i.e. representations of muscle, fat, and bones). For example, Maya Muscle [Comet, 2011] represents muscles via NURBS that drive the skin via LBS, whose volume is artist-driven in a PSD fashion. However, due to the explicit nature of NURBS, collisions are expensive, hence the performance of the system does not scale well in the complexity of the model. Rather than driving skin via LBS, Implicit Skinning [Vaillant et al., 2013] lets it slide on top of implicit surfaces via optimization. These surfaces are defined by blending components that abstract entire body parts (i.e. union of bone, muscle, and fat). In contrast to our work, note that this model is purely kinematic (i.e. no physics). Another relevant class of methods simplifies anthropomorphic components even further. For example representing an entire arm as two tapered capsules is advantageous for arm vs. cloth collision detection [Muller, 2008, Cloth/Collision]. Sphere-Meshes generalize these representations, and have recently been used to approximate geometry [Thiery et al., 2013], and track its movement in real-time [Tkach et al., 2016].

Calibrating anthropomorphic (volumetric) models

Ali et al. [2013] pioneered the transfer of anatomical structures from a template to a target human via approximate deformation models of soft tissues, and Zhu et al. [2015] calibrated these models from a set of RGBD images capturing a human in motion. Analogously to these methods, physically inspired models [Saito et al., 2015] can be calibrated to a set of 3D surface scans [Kadleček et al., 2016]

. Of particular relevance to our method is the fiber estimation technique pioneered by 

Choi and Blemker [2013] employed in [Saito et al., 2015, Sec. 3.1.1]. While Saito et al. [2015] is interested in deriving the anisotropic deformation frame, we require an explicit decomposition of the muscle in fiber bundles.

Cosserat Rods

Since being introduced to computer graphics by Pai [2002], Cosserat rods have been the subject of many efforts to achieve efficient and accurate simulation of thin materials using a variety of discretizations and energy formulations [Sueda et al., 2011; Spillmann and Teschner, 2007; Lang et al., 2011; Bergou et al., 2008; Bertails et al., 2006; Soler et al., 2018]. Bergou et al. [2010] model volume conservation in rods by updating the radius of segments based on changes in length. In pursuit of real-time performance, rods have been simulated within a PBD framework [Umetani et al., 2014; Kugelstadt and Schömer, 2016], a capability which we extend by further modeling volume preservation. Finally, note rods have been applied to muscle-skeletal simulation [Sueda et al., 2008], although to represent tendons and without consideration of volumetric effects.

3. Generalized Rod Parameterization

Our physical model of Cosserat rods consists of a smooth parametric curve in 3D space . An orthogonal frame is attached to every point on the curve. The orthogonal frame is a combination of a uniform scale , and an orthonormal matrix ; see Figure 4. Note that our model generalizes classical elastic rods [Spillmann and Teschner, 2007; Kugelstadt and Schömer, 2016], as in those models the scale is kept constant along the curve. As shown in Figure 4, any point in the parametrized volume of the rod can be transformed to the rest configuration by a function :


where is a point on a disc of radius , aligned with the plane, and centered at its center of mass. Similarly, a second function maps the rod from parameterization to its deformed configuration:


4. Variational Implicit Euler Solver

Our solver is based on the variational form of implicit Euler integration [Martin et al., 2011]. The physical model evolves through a discrete set of time samples, with simulation step size . At time the deformed position is defined as and the velocity as . The rest pose is defined as . The mass is assumed to be uniform over the rod. The sum of the external forces is denoted as . We will now drop the indexing whenever possible to improve readability. We consider position dependent internal forces such that the sum of the internal forces is


where is a potential energy function, and is a matrix of stiffness parameters which we assume to be uniform over the rod, and the notation means . We can then write implicit Euler as an optimization describing the compromise between an inertia potential and the elastic potentials:


With we indicate the inertial prediction for , i.e., its next position in absence of internal forces:


where .

[width=] /litem.pdf
Figure 5. Static behavior – We compare the deformation of a standard cosserat rod to our volumetric invariant version. Our additional degrees of freedom allow us to model the buckling (resp. bulging) caused by the stretching (resp. compression) of the rod.

[width=] /ritem.pdf

Figure 6. Dynamic behavior – Our generalized rods do not only capture static volumetric deformations, but the scale’s velocity allows us to model volume dynamics. In this example, note the rod length is unchanged, but our formulation models a volumetric shockwave travelling through the rod.


We discretize the curve in the parametrized domain using a set of points connected using piecewise linear elements of length ; see Figure 4. We can approximate the curve integral by integrating over these piecewise linear elements using the midpoint rule. For the integration we also define a set of midpoints . A point on a midpoint cross section is then parametrized as:


This is similar to the staggered grid discretization of previous work [Spillmann and Teschner, 2007; Grégoire and Schömer, 2006], where the frames are stored at the midpoints. Contrary to previous approaches, our model also has a scale that we store at the endpoints of the linear elements. We can now rewrite Equation 4 as:


5. Elastic Potentials

In this section we detail the elastic potentials used to simulate our volume preserving rods.

Strain potential

We define the strain at a midpoint as


where is a diagonal matrix of stiffness parameters and is the standard basis. and denote the deformation gradients, i.e., the Jacobian matrices of the deformation functions [Sifakis and Barbic, 2012]. As and map to , the Jacobian matrices are . The Jacobians are not rotationally invariant so we rotate them back to the parametrization domain to be able to compare them on a common ground. As explained in Appendix A integrating the strain energy (10) leads to



is the second moment of area of a disc scaled by the stiffness, and the Darboux vector is denoted by

. Note that we retrieved similar energies in previous works [Kugelstadt and Schömer, 2016] augmented by our additional scale degree of freedom. The energies (11) and (12) respectively measure the stretch along the curve and the cross section, while (13) measures the variation of scale across sections, and (14) measures bending/twisting. Interestingly, Equation 13 can also be interpreted as a measure of surface stretch.

We use an additional energy measuring the second order variation of scale complementing (13) with a measure of surface bending


where is the Laplacian operator. Note that this energy can also been seen as an approximation of:


where . This energy compares the Laplacians of the deformation functions giving a second order measure of the deformation.

Volume potential

By denoting the determinant with , and stiffness by , we define the volume preservation at midpoints as:


Note that the determinant is rotationally invariant, hence it is not necessary to rotate the Jacobians. As explained in Appendix A integrating the volume energy (17) leads to:


In Appendix E, we include a comparison between the volume preservation of our method and the inter-step update method of Bergou et al. [2010].

[width=] /fig.jpg

Figure 7. Being built on VIPER primitives, our physics engine can simulate soft-body deformation and dynamic interactions between hundreds of models in real-time. A peculiarity of our engine is that both collision and simulation are executed on the same geometry; see our video in the additional material.

6. Optimization

To approach the non-linear optimization problem in (9), we linearize the inertia and non-linear elastic terms, and then solve the optimization iteratively in a Gauss-Newton fashion. To warm start the optimization, we first compute a prediction step by ignoring the elastic potentials and by solely minimizing the inertia term – this simply provides an initial guess. We then compute a (set of) correction steps that also include the elastic potentials.

Prediction step

By denoting by the angles parametrizing the rotation matrix and is the second moment of area of a disc in world-space, the predictions for the different degrees of freedom of our model are computed as:

(21) (center prediction)
(22) (frame prediction)
(23) (scale prediction)

where is the sum of the external forces which act on the disc, is the sum of the external torques, and is a quantity which can be seen as the counterpart of the total external torque measuring the sum of the external forces projected on the position vectors; see Appendix B for more details. Note that the center and frame predictions are similar to the rigid-body equations of motion for a stretched disc. On top of these equations, we get a scale prediction describing how the scale of the disc is affected by the velocity and the external forces.

Correction steps

We then compute a set of correction steps including both inertial/elastic terms. We define


as the vector containing all the degrees of freedom, and as the vector of Lagrange multipliers. is a block diagonal matrix containing the stiffness parameters multiplied by the length of the piecewise elements, is a block diagonal matrix stacking the inertia weights multiplied by the length of the piecewise elements, and


stacks the potential energy functions. Denoting the iteration number with , the state is then updated as and , where, as derived in Appendix C:


For realtime performance we opt for using an iterative linear system solver such as block Jacobi or Gauss-Seidel. The update for the -th constraint is:


where we dropped the superscripts to improve readability, and is a relaxation parameter. Note that being a block diagonal matrix with block of size at most , can be efficiently computed. Note that Equation 28 is a generalization of the XPBD update [Macklin et al., 2016] derived for our volume preserving rod model.

7. Real-time physics engine

We implemented our real-time solver on a GPU by leveraging the Thrust framework provided by the CUDA library, which we execute on a single NVIDIA GTX 1080 graphics card. In our engine, any volumetric object is modeled as a collection of tapered capsule primitives, or “pills” (for brevity), while the floor is represented as a simple halfplane constraint. During simulation, at each time step, we start by first animating kinematic objects (e.g. bones). We then compute the inertial predictions, and perform collision detection (Sec. 7.1) to generate collision resolution constraints (Sec. 7.2). We then solve for all constraints using a Jacobi solver: constraints are computed in parallel, and the resulting positional displacements are averaged out by a reduction. The transformations of VIPERs can then be skinned to any surface mesh model (Sec. 7.4).

Our model provides a viable alternative to cages as a compact control structure for soft-body deformations. We demonstrate this in Figure 7, where we rigged a simple octopus character using rods, and solve for collisions and soft-body deformations in real time. As shown in the accompanying video, our non-optimized prototype achieves real-time performance ( sim, render) for scenes containing up to 100 octopuses, each rigged with 37 pills, whose deformation is skinned to a triangular surface mesh of faces. Note how for robust collision detection we need far fewer pills than the volumetric particles used in Macklin et al. [2014].

[width=] /fig.pdf rest configurationdeformed configuration

Figure 8. (left) An elastic band mesh and its VIPER discretization. (right) The VIPER simulation and their deformation “skinned” to the rest-pose mesh.

7.1. Collision detection

To detect collisions between physical primitives, we adopt the approach presented by Green [2010] popularized in the context of GPU particle fluid simulation. Towards this goal, we approximate each pill by its bounding sphere, and (conservatively) detect collisions to be later handled in the resolution loop detailed in Section 7.2. Collisions are detected with the assistance of a uniform grid with cell width chosen as the diameter of the largest bounding sphere, such that collisions between spheres centered in non-neighboring cells are impossible. Extending [Green, 2010], we also count the spheres in the neighboring cells of each particle, and use this information to construct in parallel a list of all potential collisions. This is in contrast to the original algorithm which for each sphere processes neighbors in series, and therefore degrades to a partially serial algorithm in cases where many particles fall into a single grid cell. Further details regarding this process are provided together with an executable 1D example in the form of a Jupyter notebook in our additional material.

7.2. Collision handling

In a generic physics engine one would implement collision detection/response between any pair of available proxies. For the sake of efficiency, in our framework we only tackle two collision proxies: (kinematic) half-planes and (dynamic) pills. Within the combinatorial set of collision pairs, the main challenge is pill-to-pill collisions. In what follows, we first compute the meta-parameters of the collisions, that are then resolved in the optimization via PBD constraints [Müller et al., 2007].

Collision metadata

Given two pills and , each modeled as a pair of spheres, for example, , the fundamental queries we need to answer are: \⃝raisebox{-0.6pt}{1} is there a collision? \⃝raisebox{-0.6pt}{2} what is the collision point/normal? \⃝raisebox{-0.6pt}{3} what is the inter-penetration amount? As typical in efficient collision resolution, we introduce a single PBD constraint modeling collision forces corresponding to the largest pill-to-pill inter-penetration. In our solution, we leverage the geometric structure of the problem: \⃝raisebox{-0.6pt}{1} a pill can be interpreted as the union of infinitely many spheres whose position and radii are linearly interpolated between its endpoints; \⃝raisebox{-0.6pt}{2} the largest inter-penetration corresponds to the inter-penetration between any pair of spheres, one in pill and one in pill . By first defining the LERP operator as , the interpolated sphere is derived by LERP’ing the endpoint quantities as and . The largest inter-penetration is then given by the solution of the bivariate optimization problem:


Because a closed-form operator providing the barycentric coordinate of the closest-point projection of a point onto the pill is available (see Appendix D) we can further simplify this problem into a scalar optimization problem:


which we solve by Dichotomous Search [Antoniou and Lu, 2007] with a fixed number of iterations (set to 10 in our engine).

Collision constraints

Having detected a collision between two pills and , and having computed (and hence ) by solving (30), we can express a constraint that correlates radii and positions on the two pills in order to resolve the collision in a least squares sense:


7.3. Scale-invariant shape matching – Bundling

To represent more complex geometry than individual rods, such as that in Figure 3, we can gather a collection of rods in a cross-section, and introduce a constraint to explain their joint deformation. We employ the assumption that muscle fibers in a muscle cross-section contract isotropically. Indexing the rods in a cross-section by , our rod deformation model can be expressed as a similarity transform


As illustrated in Figure 10, for each muscle cross-section we define a scale invariant shape-matching energy measuring the deviation of the current rod deformations from a global similarity transform of the rest deformations as


We treat this energy as a hard constraint by finding the optimal and setting as a post-processing step after few iterations. The optimal can be computed following the derivation in [Umeyama, 1991]. The optimal rotation can be found by solving


where . We compute the optimal rotation using the robust approach presented in [Müller et al., 2016]. The optimal scale can be computed as


where adds all entries of the matrix and is the Hadamard product. Finally, the optimal translation can be derived as

[width=] /litem.pdf rest configuration(and constraints)
Figure 9. Scale-invariant shape matching – Shape matching can recover a rigidly transformed configuration, while our model allows for a null-space that includes uniform scale. We employ this model as we work under the assumption that muscle fibers in a muscle cross-section contract isotropically.

[width=] /ritem.pdf

Figure 10. Muscle contraction – (top) Muscle at rest and its excitation (force shortening edge lengths in an area) with traditional shape-matching constraints (middle) vs our novel scaled shape-matching constraints (bottom).

7.4. Skinning VIPERs to surface deformation

Our VIPER rods can also be employed as a volumetric proxy driving the deformation of a surface mesh; see Figure 1 and Figure 8. In particular, we blend the relative transformations between deformed and rest pose configurations via linear blend skinning. To achieve this, we utilize a modified version of the LBS weight computation by Thiery et al. [2013]. In that paper, weights were computed by fairing an initial assignment where each vertex was fully attached to the nearest element. Instead, we modify this assignment to begin with weights for each vertex that are proportional to the inverse square-distance from the vertex to the surface of each pill. This resolves cases where multiple surfaces are at equal distance, and the nearest neighbor is multiply-defined.

7.5. Energy implementation

The energies derived in Section 5 have been derived using continuous operators. To implement these energies we approximate , , and using finite difference such that


For simplicity of the derivations we used angles to parametrize the rotation . However, our implementation uses quaternions. For small angles the corresponding quaternion is . We also use quaternions to represent rotations and approximate the Darboux vector using


where gives the imaginary part of a quaternion.

8. Anatomical modeling and simulation

We now describe how our rods can be used to model complex anatomical structures such as bones and muscles; see Figure 1. For the former, we use a simplified version of Thiery et al. [2013] where only pill primitives are used – this primarily allows us to reduce the complexity of the collision detection/resolution codebase. We then detail the conversion of digital models of muscles into VIPERs in Section 8.1, and describe a few nuances about the simulation of their motion in Section 8.2.

8.1. Muscles to rods conversion – Viperization

We developed a (weakly-assisted) technique to convert a traditional muscle model into a collection of rods. As outlined in Figure 12, our process involves several phases. We begin by computing a volumetric discretization of the muscle’s surface. We then ask the artist to paint annotations on the surface marking the start (sources) and the end (sinks) of the muscle. A harmonic solve like the one described in Choi and Blemker [2013] is then executed to compute a field that smoothly varies in the range along the muscle’s length. We then sample iso-levels of this field to be surfaces containing the vertices of each of the generated VIPERs. Within each iso-level we require: \⃝raisebox{-0.6pt}{1} VIPERs to have the same radii, and \⃝raisebox{-0.6pt}{2} to be distributed uniformly on the slice (iso-surface). To obtain this, we perform a restricted centroidal Voronoi diagram of points on each slice [Botsch et al., 2010], while simultaneously penalizing the length of each rod. We alternate this variational optimization with a discrete one that re-assigns spheres to different rods in order to minimize the sum of rod lengths. This process, which we refer to as “combing”, starts from one end of the muscle, and executes in order instances of minimum-cost bipartite matching (which we solve in polynomial time via the Hungarian algorithm), where the pairwise costs are the euclidean distances across rod nodes in two adjacent slices. A few examples of the results of this process are illustrated in Figure 1 and Figure 12.

[width=] /litem.pdf
Figure 11.

The mesh-to-VIPER conversion process. Given source/sink constraints, we compute a harmonic function in the volume, and extract a few discrete iso-levels. Within each of these, we execute a restricted CVD to place 5 elements of the same radii on these surfaces. We then execute a combinatorial optimization that connects samples across layers to produce minimal length curves.

[width=] /ritem.pdf pectoralbicepsdeltoid

Figure 12. We illustrate several examples of the VIPERs extracted by our automated process. For the pectoral, our VIPER model employs 8 rods, each discretized by 8 elements. In comparison, the surface meshes by ZIVA contain () on the boundary and its (volumetric) simulation mesh contains () tetrahedral elements. Where , , and .

8.2. Muscle simulation

Once a muscle is viperized as described in Section 8.1, rod centers within the same cross-section are connected via the bundling constraints in Equation 33. Every rod also obeys the constraints described in Section 5. The endpoints of rods are kinematically attached to bones, and intra-muscle collisions are disabled, while inter-muscle collisions are detected and resolved as described respectively in Section 7.1 and Section 7.2. Muscles can also be activated (fiber contractions generating a stronger force, producing a change of shape at constant muscle length and volume) by inserting internal forces, or even simply shortening the length of fibers resulting in the bulging effects illustrated in Figure 10. We also speed-up the solver convergence during fast motion by exploiting the availability of bone transformations. In more detail, each muscle particle has two skinning weights corresponding to the two bones the muscle is attached to. At the beginning of each frame we use the transform of each bone relative to their last frame’s transforms to initialize the displacement of the particle using LBS, and later refine this via simulation.

9. Conclusions & future work

In this paper, we introduced a novel formulation of cosserat rods that considers local volume, and optimizes for its local conservation. The resulting position-based simulation is highly efficient, and is relatively straightforward to implement on graphics hardware. We demonstrated how rod-bundling is a powerful representation for the modeling of volumetric deformation – and in particular for skeletal muscles. Rather than requiring artists to model from scratch, we also introduced an algorithm to procedurally generate VIPERs with minimal user interaction. Finally, by coupling the rod simulation to a surface mesh via skinning, our model can be thought of as a direct alternative to tetrahedral meshes and cages for real-time non-rigid deformation. Most importantly, our generalized rods formulation opens up a number of venues for future work, which we classify in three broad areas, as elaborated below.

Model generalization

While in our rods we discretized the skeletal curve with piecewise linear elements, it would be interesting to investigate whether using continuous curve models such as splines would be tractable – from both mathematical, as well as implementation standpoints. Similarly to [Müller and Chentanez, 2011, 2011], our model could be extended to model anisotropic deformations, i.e. both the rest pose and deformed rod could have a non-circular cross section. Further, while in this paper we treated the modeling of non-circular cross-sections via bundling, the theory of medial axis [Tagliasacchi et al., 2016] tells us how any shape can be approximated via primitives formed by convex-hulls of three-spheres – what Tkach et al. [2016] called “wedges”. Extending our volume-invariant rods models to volume-invariant wedges would provide an elegant generalization of our modeling paradigm.

Anatomical modeling

As highlighted by our supplementary video, rod primitives can be exploited for efficient approximate modeling and simulation of complex structures. Nonetheless, the dynamics of the human body are the result of the complex interplay between muscle, fat, and deformation of skin. Enriching our model to also account for these factors would be an interesting extension. For example, rather than driving the muscle surface via skinning as in Section 7.4, one could represent a muscle as a controllable implicit blend [Angles et al., 2017], approximate fat as an elastic offset between muscles and skin, and simulate skin as an elastic surface whose vertices lie on a (potentially) user-controlled iso-level of the implicit function. Further, while artistic editing of physically driven anatomical systems can be difficult due to the complexity of simulation, our framework could enable interactive modeling, similar to what ZBrush/ZSphere currently provides for authoring static geometry. By extending the works in [Tkach et al., 2016; Tkach et al., 2017], an efficient anatomical model for a particular user could also be constructed by fitting to RGBD data.


Our solver has not yet directly leveraged the multi-resolution structure of rod geometry. More specifically, the curve parameterization of rods offers a domain over which designing prolongation/restriction operators needed for a geometric multi-grid implementation becomes straightforward. An orthogonal dimension for optimization would be to consider the existence of multi-resolution structures within the cross-sectional domain; see Figure 3. This could be exploited in offering multi-scale interaction for artists in editing our rod models, as well as producing LOD models for efficient simulation. Finally, while we employed out-of-the-box geometry processing tools to convert a triangular mesh into a rod model, we believe fitting a fiber-bundle model to a given solid could be achieved without the (often finicky) conversion to tet mesh, but rather as a direct optimization over fiber placements.


  • [1]
  • Ali et al. [2013] Dicko Hamadi Ali, Tiantian Liu, Benjamin Gilles, Ladislav Kavan, François Faure, Olivier Palombi, and Marie-Paule Cani. 2013. Anatomy transfer. ACM TOG (2013).
  • Angles et al. [2017] Baptiste Angles, Marco Tarini, Loic Barthe, Brian Wyvill, and Andrea Tagliasacchi. 2017. Sketch-Based Implicit Blending. ACM TOG (Proc. SIGGRAPH Asia) (2017).
  • Antoniou and Lu [2007] Andreas Antoniou and Wu-Sheng Lu. 2007. Practical Optimization: Algorithms and Engineering Applications.
  • Barbič and James [2005] Jernej Barbič and Doug L. James. 2005. Real-Time Subspace Integration for St. Venant-Kirchhoff Deformable Models. ACM TOG (2005).
  • Bender et al. [2015] Jan Bender, Matthias Müller, and Miles Macklin. 2015. Position-Based Simulation Methods in Computer Graphics.. In Proc. Eurographics (Technical Course Notes).
  • Bergou et al. [2010] Miklós Bergou, Basile Audoly, Etienne Vouga, Max Wardetzky, and Eitan Grinspun. 2010. Discrete viscous threads. In ACM TOG.
  • Bergou et al. [2008] Miklós Bergou, Max Wardetzky, Stephen Robinson, Basile Audoly, and Eitan Grinspun. 2008. Discrete elastic rods. In ACM TOG.
  • Bertails et al. [2006] Florence Bertails, Basile Audoly, Marie-Paule Cani, Bernard Querleux, Frédéric Leroy, and Jean-Luc Lévêque. 2006. Super-helices for predicting the dynamics of natural hair. In ACM TOG.
  • Botsch et al. [2010] Mario Botsch, Leif Kobbelt, Mark Pauly, Pierre Alliez, and Bruno Lévy. 2010. Polygon mesh processing. AK Peters/CRC Press.
  • Bouaziz et al. [2014] Sofien Bouaziz, Sebastian Martin, Tiantian Liu, Ladislav Kavan, and Mark Pauly. 2014. Projective Dynamics: Fusing Constraint Projections for Fast Simulation. ACM TOG (2014).
  • Choi and Blemker [2013] Hon Fai Choi and Silvia S Blemker. 2013. Skeletal muscle fascicle arrangements can be reconstructed using a laplacian vector field simulation. PloS one (2013).
  • Clutterbuck and Jacobs [2010] Simon Clutterbuck and James Jacobs. 2010. A Physically Based Approach to Virtual Character Deformations. In ACM SIGGRAPH Talk sessions.
  • Comet [2011] Michael Comet. 2011. Maya Muscle. (Accessed on Aug. 8th, 2018).
  • Gould [1986] Nicholas Ian Mark Gould. 1986. On the accurate determination of search directions for simple differentiable penalty functions. IMA J. Numer. Anal. (1986).
  • Green [2010] Simon Green. 2010. Particle simulation using CUDA. NVIDIA whitepaper (2010).
  • Grégoire and Schömer [2006] Mireille Grégoire and Elmar Schömer. 2006. Interactive Simulation of One-dimensional Flexible Parts. In Proc. of ACM Symposium on Solid and Physical Modeling.
  • Ichim et al. [2017] Alexandru-Eugen Ichim, Petr Kadleček, Ladislav Kavan, and Mark Pauly. 2017. Phace: Physics-based face modeling and animation. ACM TOG (2017).
  • Jacobson et al. [2014] Alec Jacobson, Zhigang Deng, Ladislav Kavan, and J.P. Lewis. 2014. Skinning: Real-time Shape Deformation. SIGGRAPH Course,
  • Kadleček et al. [2016] Petr Kadleček, Alexandru-Eugen Ichim, Tiantian Liu, Jaroslav Křivánek, and Ladislav Kavan. 2016. Reconstructing personalized anatomical models for physics-based body animation. ACM TOG (Proc. SIGGRAPH Asia) (2016).
  • Kavan et al. [2007] Ladislav Kavan, Steven Collins, Jiří Žára, and Carol O’Sullivan. 2007. Skinning with Dual Quaternions. In Proceedings of the 2007 Symposium on Interactive 3D Graphics and Games.
  • Kugelstadt and Schömer [2016] Tassilo Kugelstadt and Elmar Schömer. 2016. Position and orientation based Cosserat rods.. In Proc. SCA.
  • Kugelstadt and Schömer [2016] T. Kugelstadt and E. Schömer. 2016. Position and Orientation Based Cosserat Rods. In Proc. SCA.
  • Kurihara and Miyata [2004] Tsuneya Kurihara and Natsuki Miyata. 2004. Modeling Deformable Human Hands from Medical Images. In Proceedings of the 2004 ACM SIGGRAPH Symposium on Computer Animation (SCA-04).
  • Lang et al. [2011] Holger Lang, Joachim Linn, and Martin Arnold. 2011. Multi-body dynamics simulation of geometrically exact Cosserat rods. Multibody System Dynamics (2011).
  • Le and Hodgins [2016] Binh Huy Le and Jessica K. Hodgins. 2016. Real-time Skeletal Skinning with Optimized Centers of Rotation. ACM Trans. Graph. (2016).
  • Lee et al. [2010] Dongwoon Lee, Michael Glueck, Azam Khan, Eugene Fiume, and Ken Jackson. 2010. A survey of modeling and simulation of skeletal muscle. ACM TOG (2010).
  • Lewis et al. [2014] J.P. Lewis, Ken Anjyo, Taehyun Rhee, Mengjie Zhang, Fred Pighin, and Zhigang Deng. 2014. STAR: Practice and Theory of Blendshape Facial Models. In Eurographics.
  • Lewis et al. [2000] J. P. Lewis, Matt Cordner, and Nickson Fong. 2000. Pose Space Deformation: A Unified Approach to Shape Interpolation and Skeleton-Driven Deformation. In Proc. ACM SIGGRAPH.
  • Li et al. [2013] Duo Li, Shinjiro Sueda, Debanga R Neog, and Dinesh K Pai. 2013. Thin Skin Elastodynamics. ACM TOG (Proc. SIGGRAPH) (2013).
  • Loper et al. [2015] Matthew Loper, Naureen Mahmood, Javier Romero, Gerard Pons-Moll, and Michael J Black. 2015. SMPL: A skinned multi-person linear model. ACM TOG (2015).
  • Macklin et al. [2016] Miles Macklin, Matthias Müller, and Nuttapong Chentanez. 2016. XPBD: Position-based Simulation of Compliant Constrained Dynamics. In Proc. of the International Conference on Motion in Games.
  • Macklin et al. [2014] Miles Macklin, Matthias Müller, Nuttapong Chentanez, and Tae-Yong Kim. 2014. Unified particle physics for real-time applications. ACM TOG (Proc. SIGGRAPH) (2014).
  • Martin et al. [2011] Sebastian Martin, Bernhard Thomaszewski, Eitan Grinspun, and Markus Gross. 2011. Example-based Elastic Materials. ACM TOG (2011).
  • Muller [2008] Matthias Muller. 2008. NVIDIA PhysX SDK 3.4.0 Documentation. (Accessed on Aug. 9th, 2018).
  • Müller et al. [2016] Matthias Müller, Jan Bender, Nuttapong Chentanez, and Miles Macklin. 2016. A Robust Method to Extract the Rotational Part of Deformations. In Proceedings of the 9th International Conference on Motion in Games (MIG ’16).
  • Müller and Chentanez [2011] Matthias Müller and Nuttapong Chentanez. 2011. Adding Physics to Animated Characters with Oriented Particles.
  • Müller and Chentanez [2011] Matthias Müller and Nuttapong Chentanez. 2011. Solid simulation with oriented particles. ACM TOG (2011).
  • Müller et al. [2007] Matthias Müller, Bruno Heidelberger, Marcus Hennix, and John Ratcliff. 2007. Position based dynamics. Journal of Visual Communication and Image Representation (2007).
  • Pai [2002] Dinesh K Pai. 2002. Strands: Interactive simulation of thin solids using cosserat models. In Computer Graphics Forum.
  • Pons-Moll et al. [2015] Gerard Pons-Moll, Javier Romero, Naureen Mahmood, and Michael J. Black. 2015. Dyna: A Model of Dynamic Human Shape in Motion. ACM TOG (Proc. SIGGRAPH) (2015).
  • Romeo et al. [2018] Marco Romeo, Carlos Monteagudo, and Daniel Sánchez-Quirós. 2018. Muscle Simulation with Extended Position Based Dynamics. In Spanish Computer Graphics Conference (CEIG).
  • Saito and Yuen [2017] Jun Saito and Simon Yuen. 2017. Efficient and Robust Skin Slide Simulation. In Proceedings of the ACM SIGGRAPH Digital Production Symposium (DigiPro ’17).
  • Saito et al. [2015] Shunsuke Saito, Zi-Ye Zhou, and Ladislav Kavan. 2015. Computational bodybuilding: Anatomically-based modeling of human bodies. ACM TOG (Proc. SIGGRAPH) (2015).
  • Scheepers et al. [1997] Ferdi Scheepers, Richard E. Parent, Wayne E. Carlson, and Stephen F. May. 1997. Anatomy-Based Modeling of the Human Musculature (SIGGRAPH ’97).
  • Schumacher et al. [2012] Christian Schumacher, Bernhard Thomaszewski, Stelian Coros, Sebastian Martin, Robert Sumner, and Markus Gross. 2012. Efficient simulation of example-based materials. In Proceedings of the 11th ACM SIGGRAPH/Eurographics conference on Computer Animation.
  • Sifakis and Barbic [2012] Eftychios Sifakis and Jernej Barbic. 2012. FEM Simulation of 3D Deformable Solids: A Practitioner’s Guide to Theory, Discretization and Model Reduction. In ACM SIGGRAPH 2012 Courses (SIGGRAPH ’12).
  • Sifakis et al. [2005] Eftychios Sifakis, Igor Neverov, and Ronald Fedkiw. 2005. Automatic Determination of Facial Muscle Activations from Sparse Motion Capture Marker Data. In ACM TOG (Proc. SIGGRAPH).
  • Smith et al. [2018] Breannan Smith, Fernando De Goes, and Theodore Kim. 2018. Stable Neo-Hookean Flesh Simulation. ACM TOG 37, 2 (2018), 12.
  • Soler et al. [2018] Carlota Soler, Tobias Martin, and Olga Sorkine-Hornung. 2018. Cosserat Rods with Projective Dynamics. In Computer Graphics Forum.
  • Spillmann and Teschner [2007] J. Spillmann and M. Teschner. 2007. CORDE: Cosserat Rod Elements for the Dynamic Simulation of One-Dimensional Elastic Objects. In Proc. SCA.
  • Sueda et al. [2011] Shinjiro Sueda, Garrett L Jones, David IW Levin, and Dinesh K Pai. 2011. Large-scale dynamic simulation of highly constrained strands. ACM TOG (2011).
  • Sueda et al. [2008] Shinjiro Sueda, Andrew Kaufman, and Dinesh K Pai. 2008. Musculotendon simulation for hand animation. ACM TOG (2008).
  • Tagliasacchi et al. [2016] Andrea Tagliasacchi, Thomas Delame, Michela Spagnuolo, Nina Amenta, and Alexandru Telea. 2016. 3D Skeletons: A State-of-the-Art Report. Proc. Eurographics (State of the Art Reports) (2016).
  • Terzopoulos and Waters [1990] Demetri Terzopoulos and Keith Waters. 1990. Physically-based Facial Modeling, Analysis, and Animation. Journal of Visualization and Computer Animation (1990).
  • Thiery et al. [2013] Jean-Marc Thiery, Émilie Guy, and Tamy Boubekeur. 2013. Sphere-Meshes: shape approximation using spherical quadric error metrics. ACM TOG (2013).
  • Tkach et al. [2016] Anastasia Tkach, Mark Pauly, and Andrea Tagliasacchi. 2016. Sphere-Meshes for Real-Time Hand Modeling and Tracking. ACM TOG (Proc. SIGGRAPH Asia) (2016).
  • Tkach et al. [2017] Anastasia Tkach, Andrea Tagliasacchi, Edoardo Remelli, Mark Pauly, and Andrew Fitzgibbon. 2017. Online Generative Model Personalization for Hand Tracking. ACM TOG (Proc. SIGGRAPH Asia) (2017).
  • Umetani et al. [2014] Nobuyuki Umetani, Ryan Schmidt, and Jos Stam. 2014. Position-based elastic rods. In Proc. SCA.
  • Umeyama [1991] Shinji Umeyama. 1991. Least-Squares Estimation of Transformation Parameters Between Two Point Patterns. IEEE Trans. Pattern Anal. Mach. Intell. (1991).
  • Vaillant et al. [2013] Rodolphe Vaillant, Loïc Barthe, Gaël Guennebaud, Marie-Paule Cani, Damien Rohmer, Brian Wyvill, Olivier Gourmel, and Mathias Paulin. 2013. Implicit Skinning: Real-time Skin Deformation with Contact Modeling. ACM TOG (Proc. SIGGRAPH) (2013).
  • Vital Mechanics [2018] Vital Mechanics 2018.
  • Xu and Barbič [2016] Hongyi Xu and Jernej Barbič. 2016. Pose-space Subspace Dynamics. ACM TOG (2016).
  • Yuen [2018] Simon Yuen. 2018. Personal communication. Head of Creatures, Method Studios.
  • Zhu et al. [2015] Lifeng Zhu, Xiaoyan Hu, and Ladislav Kavan. 2015. Adaptable anatomical models for realistic bone motion reconstruction. CGF (Proc. EuroGraphics) (2015).
  • Ziva Dynamics [2018] Ziva Dynamics 2018. Ziva Dynamics.

Appendix A Elastic Potentials


As defined in Section 5 we can write the strain energy as


where we drop the subscript from for brevity. This equation can be extended as




We can derive the gradient operator for each part of this summation leading to


where is the Darboux vector; and analogous expressions can be derived for . We can now observe that (41) can be broken up into the sum of distinct energies


as the cross terms evaluate to . After integrating over the disc we can reformulate (49) as


where is the second moment of the area of a disc scaled by the stiffness.


As defined in Section 5 we can write the volume energy as


Based on the strain derivation, the determinant can be computed as


We can now integrate over the disc leading to


Appendix B Prediction Step

As described in (9), the inertial potential over the disc is of the form


We aim at finding the unknowns that minimize (64) giving us the prediction update. We denote by the angles parametrizing the rotation matrix. Because the deformation function is non linear, i.e., due to the rotational degrees of freedom, we rely on a Gauss-Newton iterative scheme for the minimization. We linearize (64) at leading to


where k is the iteration number, and . At each iteration we minimize (65), and then apply the update , where we initialize . The matrix can be written as , where

is a cross product skew-symmetric matrix. The vector

is defined as . To compute the prediction updates we will use a single iteration of Gauss-Newton so as . We will now drop the superscripts and the subscripts to improve readability. As Equation 65 is quadratic we can be find its minimum by solving the normal equation


Interestingly, the left hand side can be simplified to a block diagonal matrix of the form


by noticing that . Moreover, as the center of mass of the disc is placed at the origin and . We can now integrate the diagonal elements leading to


where is the second moment of area of a disc in world-space. The right hand side can be simplified as


where is the sum of the external forces which act on the disc, is the sum of the external torques and is a quantity which can be seen as the counterpart of the external torque by measuring the external forces applied along the center direction.

Center update

By solving the linear system (66) for we find the prediction update for the center


As we integrate on a disc located at a midpoint this update is valid for the center of the disc located at this point. We approximate the update over the end points by using the same update rule.

Scale update

Similarly, the scale update can be computed solving the linear system for leading to


We also approximate the update over the end points using the same update rule

Frame update

The frame update can be computed by solving the linear system for leading to


Appendix C Correction Step

Inertia approximation

From the derivation in Appendix B we can obtain an approximation of the inertia term as




are the inertial predictions for the different degrees of freedom. The variational form of implicit Euler (9) can then be written in the form


where and are vectors containing all the degree of freedoms and their predictions, is a block diagonal matrix stacking the stiffness parameters scaled by the length of the piecewise elements, is a block diagonal matrix stacking the inertia weights scaled by the length of the piecewise elements, and stacks the potential energy functions.

Variational Solver

To solve this optimization we can linearize the elastic potentials and write an iterative Gauss-Newton optimization