Benjamin Gilles

Guillaume Bousquet

Frame-based Elastic Models

Keywords: I, 3, 5 [Computer Graphics]: Computational Geometry and Object Modeling-Physically-Based Modeling; I, 3, 7 [Computer Graphics]: Three-Dimensional Graphics and Realism-Animation Physically based animation, deformable solids

We present a new type of deformable model which combines the realism of physically based continuum mechanics models and the usability of framebased skinning methods. The degrees of freedom are coordinate frames. In contrast with traditional skinning, frame positions are not scripted but move in reaction to internal body forces. The displacement field is smoothly interpolated using dual quaternion blending. The deformation gradient and its derivatives are computed at each sample point of a deformed object and used in the equations of Lagrangian mechanics to achieve physical realism. This allows easy and very intuitive definition of the degrees of freedom of the deformable object. The meshless discretization allows on-the-fly insertion of frames to create local deformations where needed. We formulate the dynamics of these models in detail and describe some pre-computations that can be used for speed. We show that our method is effective for behaviors ranging from simple unimodal deformations to complex realistic deformations comparable with Finite Element simulations. To encourage its use, the software will be freely available in the simulation platform SOFA.

INTRODUCTION

Deformable models are essential for computer animation, especially for animating characters and soft objects. In current practice, however, the animator has to choose between two very different approaches (see Sec. 2 for a brief review).

One approach is skinning (also known as vertex blending or skeletal subspace deformation). The deformation is kinematically generated by manipulating "bones," i.e., specific coordinate frames. This method is widely used, not only for its simplicity and efficiency, but because it provides natural and intuitive handles for controlling deformation. Skinning generates smooth deformations using a very sparse sampling of the deformation field. Adaptation is simple since frames can be inserted easily to control local features. These interesting features have made it the most widely used method for character animation. However, as a consequence of its purely kinematic nature (i.e., the frame positions need to be scripted), achieving physically realistic dynamic deformation is a major challenge with this approach.

The other approach is physically based deformation, typically using continuum mechanics. This has the significant advantage that physical realism is "baked in" right from the start. Complex animations are generated by numerical integration of discretized differential equations. However, these methods can be expensive and difficult to use. In the popular Finite Element Method (FEM) framework, the degrees of freedom of the discretized model are the vertices of a mesh, which must be constructed for each simulation object. A relatively fine mesh (i.e., a dense sampling of the deformation field) is required to capture common deformations such as torsion, leading to expensive simulations. Mesh adaptation can be difficult due to the topological constraints of the mesh. Particle-based meshless methods have been proposed to address these problems. While they obviate the need to maintain mesh topology, particles can not be placed arbitrarily, since each material point has to be in the range of at least four non-coplanar particles. Therefore, these methods also need a dense cloud of particles not very different from the vertices of an FEM mesh.

In this paper, we propose a new approach that combines the advantages of both skinning and physically based deformation. Instead of the vertices of a mesh, the degrees of freedom are a sparse set of coordinate frames, which parameterize material points using an advanced skinning method called dual quaternion blending B. Gilles et al. [ [START_REF] Kavan | Geometric skinning with approximate dual quaternion blending[END_REF]]. The equations of motion are derived for the moving frames by applying the principles of continuum mechanics across the volume of the deformed object, and solved using classical implicit time integration.

Contributions

Our main contribution is a new approach which unifies skinning and physically based deformation modeling. This allows the creation of efficient physically sound models with very sparse and intuitive sampling, and on-the-fly adaptation to create local deformations where needed.

In contrast with FEM and traditional particle-based meshless methods, there is no constraint on the number of frames influencing each point of the deformable object, because a single frame is sufficient to represent a displacement field including rotations. This also allows very sparse sampling; indeed, some simple deformable models are easily modeled using a couple of frames.

Complex materials can be accurately modeled using maps of distributed material properties such as density and stiffness. The weight functions of the control frames also provide the animator with an additional, very intuitive, method to model the physical behavior required in an animation. For instance, rigid regions are easily obtained by locally associating a large weight to a given frame and small weights to the others. Moreover, we show that it is possible to encode the mass and stiffness matrices in six-dimensional mass-spring systems which are computationally efficient and allow large deformations. Finally, local deformations can be adaptively modeled, using dynamically inserted control frames with appropriate weight functions.

The remainder of the paper is organized as follows. Previous work is summarized in Section 2. The physics of our deformable model is presented in Section 3. Application to object modeling is presented in Section 4. Results are discussed in Section 5 and conclusions are drawn in Section 6.

BACKGROUND

Physically based deformable models have attracted continuous attention in Computer Graphics, since the seminal work of Terzopoulos [START_REF] Terzopoulos | Elastically deformable models[END_REF]. We refer the reader to the excellent survey of [START_REF] Nealen | Physically based deformable models in computer graphics[END_REF]] on this topic. Here, we briefly re-view the main Lagrangian models of deformable objects, followed by a short introduction to skinning.

Mesh-based methods. Early works on deformable models in Computer Graphics have focused on interconnected particles. In mass-spring systems [START_REF] Platt | Animating facial expressions[END_REF], constraints on edge length are enforced to counter stretching. Bending and shear can be controlled using additional springs. More general constraints such as area or volume conservation can be enforced using appropriate energy functions [START_REF] Teschner | A versatile and robust model for geometrically complex deformable solids[END_REF]]. To realistically model volumetric deformable objects, it is necessary to apply continuum mechanics. The spatial derivatives of the displacement field can be computed using finite differences on a regular grid [START_REF] Terzopoulos | Elastically deformable models[END_REF]]. [START_REF] Terzopoulos | Dynamic with constraints for interactive sculpting[END_REF] studied the case of physically deformable NURBS surfaces for shape modeling. Finite elements [START_REF] Bathe | Finite Element Procedures[END_REF][START_REF] Gourret | Simulation of object and human skin formations in a grasping task[END_REF][START_REF] Brien | Graphical models and animation of brittle fracture[END_REF][START_REF] Cotin | Real-time elastic deformations of soft tissues for surgery simulation[END_REF]] allow irregular meshes, which are generally more convenient to sample objects with arbitrary shapes, but may be poorly conditioned. The spatial domain is subdivided into elements such as triangles, hexahedra or more frequently tetrahedra, in which the displacement field is interpolated using shape functions. At each point the strain can be computed using the spatial derivatives of the displacement field. Accurate material models have been implemented from rheological models relating stress and strain in hyperelastic, viscoelastic, inhomogeneous, transversely isotropic and/or quasi-incompressible media [START_REF] Weiss | Finite element implementation of incompressible, transversely isotropic hyperelasticity[END_REF]. For simplicity, linearized strain has been applied assuming small displacements in rotated frames [Müller and Gross 2004]. Precomputed deformations modes have been used to interactively deform large structures [Barbič and James 2005]. Models based on Cosserat points have been proposed for large deformations in thin structures [START_REF] Pai | STRANDS: Interactive simulation of thin solids using Cosserat models[END_REF]] and solids [START_REF] Nadler | A new 3-d finite element for nonlinear elasticity using the theory of a cosserat point[END_REF]. Since robustness problems such as inverted tetrahedra [START_REF] Irving | Tetrahedral and hexahedral invertible finite elements[END_REF]] or hourglass deformation modes in hexahedra [START_REF] Nadler | A new 3-d finite element for nonlinear elasticity using the theory of a cosserat point[END_REF] have been addressed, meshing remain the main issue in finite elements. To reduce computation time, embedding detailed objects in coarse meshes has become popular in computer graphics [Müller and Gross 2004;[START_REF] Sifakis | Arbitrary cutting of deformable tetrahedralized objects[END_REF][START_REF] Nesme | Preserving Topology and Elasticity for Embedded Deformable Models[END_REF]].

Multi-resolution approaches have been proposed [START_REF] Debunne | Dynamic real-time deformations using space and time adaptive sampling[END_REF][START_REF] Grinspun | Charms: A simple framework for adaptive simulation[END_REF]. In recent work, disconnected or arbitrarily-shaped ele-ments [START_REF] Kaufmann | Flexible simulation of deformable models using discontinuous galerkin fem[END_REF][START_REF] Martin | Polyhedral finite elements using harmonic basis functions[END_REF] have been proposed to alleviate the meshing difficulties.

Meshless methods. Meshless methods do not use an underlying embedding structure but unstructured control points. In computer graphics, meshless methods have been first introduced for fluid simulation and then extended to solid mechanics [Müller et al. 2004;[START_REF] Gross | Kinematic control of free rigid bodies using dual quaternions[END_REF]. Each control point has a given influence that generally decreases with the distance to it. Standard approximation or interpolation methods have been investigated for physical simulation such as Shepard functions, radial basis functions and moving least squares (see [START_REF] Fries | Classification and overview of meshfree methods[END_REF] for an extensive review). Despite the added flexibility due to the absence of elements, the approximation function still requires a certain overlap between the influence regions of control points. Particularly the approximation of rotations requires at least four control points, contrary to our method that explicitly use rotations in the degrees of freedom. [START_REF] Martin | Unified simulation of elastic rods, shells, and solids[END_REF] alleviate this problem using an extended moving least squares method based on deformable frames instead of particles. Their main focus is on a novel differential measure of deformation, called elaston, combined with standard RBF weight functions and dense sampling, resulting in large computation times. In contrast, our focus is on fast computation and ease of modeling, using sparse sampling and customized shape functions.

Besides continuum mechanics-based methods, fast algorithms have been developed for video games to simulate quasiisometry [START_REF] Adams | Meshless modeling of deformable shapes and their motion[END_REF]Müller et al. 2005]. They are not able to model real materials, being based on geometry only.

Skinning. Meshless approximation functions estimate a continuous deformation field at locations influenced by at least four control points while finite element functions interpolate deformation inside elements. A third approach, unexploited in physical simulation until now, consists in interpolating rigid transformations directly. In character animation, this is well known as skinning (or vertex blending or skeletal subspace deformation), in which soft-tissues are deformed by blending rigid transforms [START_REF] Magnenat-Thalmann | Joint dependent local deformations for hand animation and object grasping[END_REF]]. Various rigid motion interpolation methods have been proposed and we refer the reader to [START_REF] Kavan | Geometric skinning with approximate dual quaternion blending[END_REF]] for a detailed review. Specifically, dual quaternion blending offers a good approximation of the linear interpolation of screws at a reasonable computational cost [START_REF] Kavan | Skinning with dual quaternions[END_REF]]. It provides a closed-form solution for more than two transforms contrary to screw interpolation that requires an iterative treatment. We believe that this makes it well suited for parameterizing a physically based deformable model; we exploit these features in our model.

Basics and notation.

We now provide background on dual quaternion skinning and introduce the notations used in our derivation of the corresponding 3d elasticity method. Let p = [x ȳ z] T be the initial coordinates of a material point in the world coordinates system, as illustrated in figure 2. To simplify the derivations and without loss of generality, we assume that [x ȳ z] are also the coordinates of the point in the local parameterization of the object. In this paper, we use a horizontal bar on top to denote the value in the reference, undeformed, position. Using skinning, the displaced position p = [x y z] T of the point is computed as

p = Rp + t (1)
where the rotation matrix R and the translation vector t are computed using a weighted sum of the displacements of one or several moving frames influencing the point. In linear blend skinning, R Fig. 2. An object undergoing skinning deformation using two moving frames.

and t are simple weighted sums of each frame's matrix and vector, and well-known artifacts occur with large rotations. These problems can be alleviated using dual quaternion skinning. Dual quaternions are pairs of quaternions with a special algebra which represent rigid transforms when normalized. Most of quaternion properties hold for dual quaternions. We briefly outline the method here, and refer the reader to [START_REF] Kavan | Skinning with dual quaternions[END_REF][START_REF] Kavan | Geometric skinning with approximate dual quaternion blending[END_REF]] for more details. The position of frame i is represented using a dual quaternion written as a 8d vector q i = [q iT 0 q iT ε ] T where q i 0 is a unit quaternion representing rotation and q i ε = (t i∧ • q i 0 )/2 is a quaternion representing translation to the frame origin t i . Operator • denotes standard quaternion product and ∧ denotes conversion from a 3d vector to a quaternion with null scalar part. Blended displacements are computed as normalized weighted sums of dual quaternions:

b ′ = b b = w i Ti q i w i Ti q i (2)
where the 8 × 8 matrix Ti encodes dual quaternion product with q-1 i , the inverse of the frame transform in the reference position. This product represents the transformation applied to the reference configuration, and it is used as a kinematic parameter in [START_REF] Kavan | Geometric skinning with approximate dual quaternion blending[END_REF]]. The matrix and vector of equation ( 1) are then straightforwardly deduced from b ′ (see appendix A). Figure 2 shows some isolines of the parameterization, depicted as initially orthogonal lines within a volume element ∆v centered on point p. The weights w i are functions of p, so the blended displacement {R, t} must be computed at each vertex of the deformed geometry. We represent frame velocity using the 6×1 vector Ω i = [ω T i ṫT i ] T where ṫi and ω i are the linear and angular velocities. Similar to standard quaternions, there is a linear relationship between the velocity vector and the rate of the frame [Han et al. 2008

]: qi = (ω i • q i )/2, where ωi = [ω T i∧ ( ṫi + t i × ω i ) T ∧ ]
T is a dual quaternion representing rigid body velocity. We encode this linear relationship using the 8 × 6 matrix L i such as qi = L i Ω i . We introduce the directional derivative operator ∇ i to denote differentiation with respect to the six independent degrees of freedom of frame i , i.e. for a given N-dimensional vector u(q), we have

u = i ∇ i u Ω i where ∇ i u is a N × 6 matrix.

THE DYNAMICS OF FRAME-BASED CONTINUUM

As shown in the following diagram, we apply a classical methodology in continuum mechanics: from the degrees of freedom (i.e., frame positions and orientations), we interpolate the displacement field in the material (i.e., dual quaternion skinning, equation ( 2)), from which we compute the strain and strain derivatives through spatial differentiation (section 3.3). Material properties are then used to compute the internal elastic energy from the strain (section 3.4). By differentiation with respect to the DOFs, we obtain,

• B. Gilles et al.

for each frame, a generalized force and force Jacobian (section 3.5) that we integrate in time to compute the acceleration, velocity and the new position. The simulation loop is summarized in the following diagram; the symbols are defined in the next section.

Position

Strain Energy

q ∂ → E Material → W ↑ ↓ ∇ Ω ← Ω Mass ← Γ Velocity Acceleration Force
In the remainder of this section, we first set up the dynamics equation, then we derive the local terms of this equation for a continuum animated using dual quaternion skinning.

Differential Equation

Models derived using Lagrangian mechanics are ordinary differential equations which have the following form in generalized coordinates:

M Ω -Γ(q, Ω) = Γ ext (q, Ω) (3)
Here q and Ω are the positions and velocities described in the previous section, Ω denotes accelerations, Γ internal forces, Γ ext external and inertial forces, and M is the mass matrix. Numerous methods have been proposed to solve this ODE, and without loss of generality we will use Implicit Euler integration (e.g., [START_REF] Baraff | Large steps in cloth simulation[END_REF]), which computes velocity updates by solving the following equation:

M -hC -h 2 K δΩ = h (Γ ext + hKΩ) (4)
where h is the time step, K = dΓ dq is the stiffness matrix, and C = dΓ dΩ the damping matrix. In this paper, we use the popular Rayleigh assumption: C = αM + βK.

Frame positions are updated based on the updated velocities:

q i 0 (t + h) = q i 0 (t) + h 2 (ω i∧ • q i 0 (t)) (5)
t i (t + h) = t i (t) + h ṫi ,
followed by a re-normalization of the quaternion part.

Mass

The generalized mass matrix can be computed by assembling 6 × 6 blocks:

M ij = V ρJ T i J j dV (6)
where ρ is the mass density, and the 3×6 Jacobian J i maps a frame velocity vector to a point velocity such as: ṗ = i J i Ω i . The Jacobian is obtained using chain rule differentiation of equation ( 1):

J i = ∇ i p = ∂p ∂b ′ ∂b ′ ∂b ∂b ∂q i ∇ i q i = w i QN Ti L i (7)
Matrix Ti (8×8) and L i (8×6) were introduced in section 2. Matrix N (8×8) is a projection which filters out the changes of norm of the dual quaternion. Matrix Q (3×8) encodes displacement of the given point associated to a change of the normalized dual quaternion b ′ . The expressions of these matrices are given in appendix A.

Strain

In the rest shape, the columns of the deformation gradient dp/dp form an orthonormal reference frame, while in deformed state stretching corresponds to vector length change, and shearing corresponds to angles, as illustrated in figure 2. The deformation gradient can be computed by differentiating equation ( 1):

F = dp(p, b ′ (p)) dp = ∂p ∂p + ∂p ∂b ′ ∂b ′ ∂b ∂b ∂p = R + QNW (8)
Matrix W (8×3) encodes the dependence of the weighted sum b on the material point and depends on the local variations W i = ∂w i /∂p of the weight functions in the material (see appendix A). Three dimensional strain measures are generally based on the deformation gradient. Due to space restrictions, in the remainder of this paper we focus on the popular Green-Lagrange strain tensor, valid for large displacements. The six independent terms of this tensor E = (F T F -I)/2 can be represented in vector form as:

E = [ǫ xx ǫ yy ǫ zz ǫ xy ǫ yz ǫ zx ] T .

Strain energy and stress

The strain energy density W (p), representing the potential elastic energy per unit volume stored at a given point, is:

W = E T S/2
where S is the stress in a vector form. Due to space restrictions, we focus on Hookean elasticity and obtain the popular St. Venant-Kirchhoff material: S = HE, where H is the standard 6 × 6 Hooke's stiffness matrix. The density of strain energy is thus

W = E T HE/2.
Non-linear models could be easily obtained by updating H and its derivatives as a function of E at each time step. For non-uniform materials, H simply depends on p and for anisotropic materials, H takes a generalized form as discussed in [START_REF] Bathe | Finite Element Procedures[END_REF]].

More general models can be derived using the same methodology. In appendix B, we provide formulas for volume preserving materials, and general hyperelastic materials based on the invariants of the deformation tensor.

Force and stiffness matrix

The computation of forces depends on the derivative of the strain with respect to the DOFs:

∇ i E = 1 2 (F T ∇ i F + ∇ i F T F) (9)
This 3 × 3 × 6 third-order tensor is computed by differentiating equation (8) (see appendix A). We write it as a 6 × 6 block B i using the vector form of the strain:

B i = ∇ i E (10)
The generalized force Γ i acting on frame i is the negative of the energy gradient with respect to the degrees of freedom:

Γ i = - V ∇ i W T dV = - V B T i HEdV (11)
The stiffness matrix K used in implicit integration methods (see Eq. 4) is the Jacobian of the force. Each 6 × 6 block K ij encodes the variation of the force on i due to a displacement of j:

K ij = ∇ j Γ i = - V [B T i HB j + (∇ j B T i )HE]dV (12)
The first term encodes the change of intensity of the internal forces, while the second involves a third-order tensor which encodes the change of direction. In our experiments, the latter can be safely ignored, but it may be necessary in case of large angular velocities.

External forces

Since the dynamics equation ( 3) is applied to the moving frames, which are the independent DOFs, each 3d external force f ext applied to a given point p must be mapped to 6d forces applied to the moving frames. Using the Jacobian and applying the power conservation law

f ext T ṗ = i Γ ext T
i Ω i for any Ω i , we find that the external force is dispatched over the frames as:

Γ ext i = J T i (p)f ext (13)

MODELING OBJECTS

The quantities derived in section 3 are analytically defined at any point in space, except for the weights and weight derivatives that are numerically estimated as we will see in section 4.1. In this section, we address weight, sampling and integration issues.

Weight functions

Weight functions encode frame influences and must be defined at each integration point and geometry vertex. In contrast to other meshless methods, our weights do not need to constitute a partition of unity thanks to the intrinsic normalization in equation ( 2), and can be negative (see example with cubic splines in section 5).

When normalized, they are analogous to shape functions in FEM. Moreover, by setting quasi-infinite weights near frames, our deformation field becomes interpolating and only approximating, which is convenient to constrain and interact with the model us-Dirichlet constraints. This is not the case with moving-leastsquares meshless methods [START_REF] Adams | Meshless modeling of deformable shapes and their motion[END_REF][START_REF] Martin | Unified simulation of elastic rods, shells, and solids[END_REF], where penalty terms need to be added. A frame weight function is typically a decreasing function of the distance to the frame origin in the rest configuration. A cutoff value can be set, with the constraint that each point must always be influenced by at least one frame. In figure 3, we show a U -shaped object controlled by two frames 0 and 1, and we display with colors the weight of frame 0 computed with different methods. Euclidean distance fails to capture the geometry of the object because close parts in Euclidean space move similarly (point p is equally influenced by frames 0 and 1 in the example). Though geodesic distances are more accurate, artifacts can occur far from all frames. At p, the influence of frame is not negligible compared with the influence of frame 1. As a result, this point moves along with frame 0 even if frame 1 remains fixed, contrary to what we would expect. To solve this, we use harmonic weight functions [START_REF] Joshi | Harmonic coordinates for character articulation[END_REF]] computed at initialization time. To compute the weight field of a frame, we fix its weight to a maximum value at its origin and to 0 at the origin of the others. To ease the numerical solution of the harmonic equation, our samples are aligned on a voxel grid, and we iteratively convolve the grid with a Gaussian mask until convergence. We assume that boundaries are dissipative, meaning a constant extrapolation of the weights in the Gaussian kernel. This results in the desired weight distribution as shown in the right of Figure 3 where p is only influenced by frame 1. The weight gradients are required in the computation of matrix W involved in the strain formulation. We approximate them using central differences in the weight map. It is straightforward to rigidly attach a part of an object to a given frame, by constraining the weights of all other frames to a null value, as shown in Figure 8. This is more convenient than using high stiffnesses, which are numerically difficult. The transitions to the deformable parts are natu-Fig. 3. Left: Discretization for volume integrals computation. Right: Weight related to frame q 0 using different distance measures (red=max weight and white=min weight) rally obtained thanks to the smoothness of the harmonic functions. Exploring the possibility of tuning the material stiffness using the shape functions is clearly an avenue for future work. This is more difficult using FEM and traditional particle-based approaches because points do not encode rotations. [START_REF] Nesme | Preserving Topology and Elasticity for Embedded Deformable Models[END_REF]] encode the material stiffness within coarse elements using shape functions after a fine-level static analysis, but they use traditional shape functions and material parameters at the fine level.

Volume integrals

To apply the dynamics equation (3) to a given object, we need to integrate the quantities defined in section 3 over the volume of the object. In finite elements, quadrature rules are generally available to compute the integrals in each cell as weighted sums of a small number of function values at precise parametric points. This is not possible in our meshless method. We compute the integral of any function f by regularly discretizing the volume inside the bounding box of the undeformed object:

V f (p)dV ≃ k f (p k )∆v , (14)
where volume ∆v corresponds to the sampling resolution, and the function is null outside the object (see figure 3 left and figure 4). This straightforwardly allows us to model complex objects with non-uniform mass or stiffness, provided that a property map such as a 3d texture or a medical image is available. At each time step, we compute the terms of the system's differential equation (equation ( 3)) by summing the local quantities defined in sections 3.2-3.6 over all the samples, as shown in Algorithm 1.

foreach frame i do Compute L i ; // sec. 2 Γ i = 0, M i = 0; for frame j in 0..i do K ij = K ji = 0; end end foreach volume sample s do compute b ′ s ; // sec. 2 compute R s , t s , Q s , N s , F s , E s ; // sec. 3.3 foreach frame i do compute J si ; // eq. 7 compute B si , Γ i += B T is H s E is ∆v ; // eq. 11 M i += ρ s J T si J si ∆v ; // eq. 6 for frame j in 0..i do K ij -= B T is H s B js ∆v ;
// eq. 12 end end end foreach frame i do for frame j in 0..i-1 do

K ji = K T ij ; end end
Algorithm 1: Computation of forces and system matrices using the samples.

Fast pre-computed models

The computation time of the volume integrals as explained in section 4.2 is proportional to the number of volume samples. A tradeoff between efficiency and speed can be obtained by tuning the resolution of the grid. It may also be possible to pre-compute a reduced set of representative voxels at initialization time in the same spirit as [START_REF] An | Optimizing cubature for efficient integration of subspace deformations[END_REF]]. To speed up the computations even more, we show how to pre-compute the mass and stiffness matrices at initialization time. This frees us from iterating over the points at each time step.

Mass. For simplicity we lump the mass of each frame by neglecting the cross terms in equation ( 6): Mi = Mii . The resulting global mass matrix is block diagonal, which simplifies the time integration step (equation ( 3)). Note that the mass has only an influence on the dynamics of the system, and not on the static solutions (i.e., equilibrium configurations), which makes the simplification hardly noticeable in many applications. We pre-compute the mass matrices in the reference configuration and we update them at each time step. Let R i be the 3 × 3 rotation matrix of frame i with respect to its initial orientation and R i the associated rotation matrix in 6d. The updated mass at time t is then given by Mi as following:

M i = R i Mi R T i , R i = R i 0 0 R i (15)
Stiffness. To avoid stiffness computations at each time step, we transform the initial stiffness matrix into a generalized joint spring network.

Fig. 5. Principle of the pre-computed 6d-spring system. In green: registration for stiffness warping; in black: initial frames; in blue: deformed frames At initialization time, for each frame pair (i, j), we compute the hinge frame m in the middle of the two, as illustrated in Figure 5. In a deformed configuration, the joint is represented by two frames m i = (q i • q-1 i • m) and m j = (q j • q-1 j • m) respectively attached to frames i and j, and the gap in between is a measure of the deformation created by the relative displacement of the moving frames. In the reference state, block Kij encodes the 6d force created at t i , the origin of frame i, as a function of a displacement of frame j expressed at t j ,the origin of frame j. The equivalent joint stiffness Km encodes the 6d force as a function of a 6d displacement, both expressed at t m , the origin of frame m. Computing the joint stiffness based on the off-diagonal block requires simple rigid body mechanics to change the points where the 6d forces and displacements are expressed. Let Θ(j) = ω u(j) denote a small 6d displacement (velocity times unit time step) expressed at point t j . The same displacement expressed at point t m is given by:

Θ(m) = I 0 l j × I Θ(j) = D(l j )Θ(j) , (16)
where l j = (t jt m ), (.) × is the cross product matrix operator. Note that D(a)D(b) = D(a + b) for arbitrary vectors a, b, so multiplying equation 16 with D(-l j ) gives the inverse relation

Θ(j) = D(-l j )Θ(m). Let Γ(m) = τ(m)
f denote a 6d force applied at point t m . The same force expressed at point t i is given by:

Γ(i) = I -l i × 0 I Γ(m) = D(l i ) T Γ(m) (17)
where l i = (t it m ). Based on this, we derive the joint stiffness corresponding to an off-diagonal block as:

Km = D(-li ) T Kij D(-lj ) (18)
The joint stiffness is considered constant in the joint frame m. We project 6d vectors from this frame to the world frame using matrix

R m = R m 0 0 R m
, where R m is the 3 × 3 matrix of frame m. The joint stiffness in world coordinates is thus given by:

K m = R m Km R T m (19)
The corresponding block of the warped stiffness matrix is derived by expressing the force applied to frame i a function of the displacement of frame j:

K ij = D(l i ) T K m D(l j ) (20)
To derive the other blocks, we notice that the displacements of i and j create opposite forces because elastic forces come from relative displacements. Moreover, a joint force on one frame is associated with the opposite force on the other frame, due to Newton's third law. Based on this, and by changing the points where the displacements and forces are expressed, we derive the remaining blocks as:

K ii = -D(l i ) T K T m D(l i ) (21)
K jj = -D(l j ) T K m D(l j ) (22)
K ji = D(l j ) T K T m D(l i ) (23)
Linear and angular momentum are preserved by construction (K ii = -D(t it j ) T K ji ). One can easily verify that, in the reference configuration, K = K. When we have a global rotation, K is just rotated by R m , as expected. We can conclude that the stiffness matrix derived from our joint springs has all the expected properties.

Given frame positions, each joint force is computed based on the gap between the joint frames shown in Figure 5. Let θ i and θ j be the transformations from frames m to m i and m to m j , respectively, expressed as 6d displacements. Since θ i = -θ j by construction, we straightforwardly deduce the forces applied to the frames as:

Γ i = D(l i ) T (K m + K T m )θ j (24) Γ j = -D(l j ) T (K m + K T m )θ j (25)
The net internal forces and stiffness matrices are computed by summing the contribution of all pairs of frames. Only the frames with intersecting weight fields interact. The weight cutoff value can be used to tune the sparsity pattern. The initial 6 × 6 stiffness matrices Km of these generalized joint springs are pre-computed by accumulating the contributions of sample points at arbitrarily fine resolutions. They can thus precisely encode complex material and geometry for small deformations, and they provide acceptable precision for large deformations. The location of the joint can be arbitrary. While the halfway location seems a reasonable choice, the study of its optimal position with respect to the object geometry is kept for future work.

Using this approach, the computation of the matrices using Algorithm 1 is applied only once at initialization. Then at each time step, the much simpler Algorithm 2 is applied.

foreach frame i do Γ i = 0; for frame j in 0..i do K ij = K ji = 0; end end foreach frame i do M i = R i Mi R T i ; // sec. 4.3 for frame j in 0..i-1 do add K ij , K , K ii , jj ; // eq. 20-23 add Γ i , Γ j ;
// eq. 24-25 end end Algorithm 2: Computation of forces and system matrices using the mass-spring system.

Adaptivity

A local deformation can be modeled by inserting a frame at a desired location p * within the deformed object, as shown in Figure 9.

insert a frame, we need to estimate the associated undeformed position p(p * ). The analytical inversion of the skinning mapping function p(p) is difficult. Fortunately, we can efficiently solve the equation numerically using Newton's method and the deformation gradient F = ∂p/∂p of equation ( 8). The initial p(0) is set to the closest sample point. We then iteratively refine the solution using: p(k+1) = p(k) + F(p (k) ) -1 (p *p(p (k) )). The search typically converges in 3 iterations with a stop criterion of p *p(p (k) ) < 10 -5 . Once the additional frame is inserted, we compute its weight field in the reference state and update the mass-spring network accordingly. This adds a new mass matrix block along with a number of new stiffness blocks. Moreover, all previous mass and stiffness blocks whose spatial domain intersects the inserted weight field must be recomputed. To speed up these computations, the weights of the inserted frame are based on geodesic distances, quickly evaluated in a single propagation pass through the voxel grid up to the desired cutoff value. Each voxel contains the list of frames and weights which influence it, and contributes to the (lumped) mass block using equation ( 6) and to new stiffness blocks using equation (12).

VALIDATION AND RESULTS

Implementation

Our method has been integrated in the simulation platform SOFA (http://www.sofa-framework.org/) and will be freely availin the upcoming release. In principle, any time integration scheme can be used to solve equation ( 3), although we prefer the implicit Euler scheme with a conjugate gradient solver for stability and large time-steps. For implicit schemes, the stiffness matrix K is evaluated to provide force variations given a displacement. The solver returns an updated generalized velocity

Ω i = [ω T i ṫT i ] T for each frame i.
our framework, we decouple physics, collision handling and visualization through three different models. Any type of visual and collision models can be plugged into our method (e.g., discrete/parametric surface/volume, bounding volume hierarchy, etc.) as we provide a volumetric displacement field, as well as mapping functions J i to transfer interaction forces to the frames.

The main simulation loop maps samples to their deformed positions and computes the strain, the mapping function and the strain derivative at their location. For visual model points, only deformed positions are required, while for collision points, the mapping function J i is also needed. To handle collisions and self-collisions, the fast GPU-based detection method presented in [START_REF] Faure | Imagebased collision detection and response between arbitrary volumetric objects[END_REF]] and the method for collision response presented in [START_REF] Allard | Volume contact constraints at arbitrary resolution[END_REF]] are used. Finally, we accumulate, for each frame, sample contributions in terms of mass, force and stiffness (space integration), handle collisions and perform time integration.

For the fast pre-computed version of the simulation (section 4.3), there is no physical sample involved, since the material behavior is encoded in the pre-computed stiffness matrices.

Accuracy

To validate the mechanical accuracy of our method, we model the deformation of a beam using regularly spaced frames along the axis, with piecewise linear weight functions. We apply an extension force to the beam and verify that the force-extension law precisely matches the theoretical St. Venant-Kirchhoff model f = ǫ + 3ǫ 2 /2 + ǫ/2, independently of the number of frames and the volume sample densities. Bending is more complex because it simultaneously involves extension-compression and shear, especially with large displacements as shown in the example in Figure 6. We obtain similar results with FEM, our model being slightly less extended, as illustrated in the left of the figure. This confirms that accurate continuum mechanics can be performed using our model. In the right of Figure 6, we show bars deformed using different numbers of frames. As usual, fewer degrees of freedom result in more stiffness. As with FEM, the weight functions have an important influence on the results. The slight bumps on the interpolated surface are due to the nearly constant curvature between two consecutive frames, resulting from the dual quaternion blending between the two frames. A Catmull-Rom spline interpolation can generate a smoother shape due to the negative weights that couples non-consecutive frames. However, the gradient of the weights is not constant, producing inhomogeneous deformations which is not what we expect with homogeneous materials.

Deformation modeling and adaptivity

The most appealing feature of our method is probably its ability to easily model deformable objects using a reduced number of control primitives. In the video, the T-shaped rubber object shown in Figure 7 (Young's modulus E = 200kP a, Poisson's ratio ν = 0.3) exhibits compression, shear, bending and torsion using only two frames, corresponding to a total of 12 DOF. The same number of only allows to model a single linear tetrahedron in FEM, which can not exhibit torsion and bending! In the example with three frames, the mechanical properties are computed by integration over the volume as described in section 4. The object automatically exhibits an asymmetric stiffness reflecting its asymmetric shape.

The turtle shown in Figure 8 (E = 700kP a, ν = 0.3) illustrates the ability to tune the stiffness using the shape functions. The turtle is animated using seven frames, and the weights associated with the frame located in the shell appear as color maps. In the middle, the weights are computed using standard harmonic coordinates. The shell is not rigid because the weight of the shell frame is lower in the regions near the neck and near the legs, which thus undergo non-rigid motion blending. In the right of the figure, we have imposed a high weight inside the whole shell and let the weight in the remaining volume be computed automatically. This results in a rigid shell, as shown in the video. In the example shown in Figure 1, two frames are sufficient to model a deformable bunny (E = 700kP a, ν = 0.3) . Adding flexibility to an ear is done using an additional frame, which is dramatically simpler than editing the mesh of an FEM model. In the same video, we show that a dynamically inserted frame at the contact point with an object can be used to generate a local deformation, as illustrated in Figure 9. The range of the local deformation can be tuned using the weight function of the inserted frame. Such a high level of adaptivity in a physical model is straightforward with our model, while it is difficult to implement using previous methods.

Performance

The FEM simulation of the bar in Figure 6 runs at 65 fps. While our model with 5 frames runs at 38 fps, the speed increases dramatically to 700 fps using precomputed generalized springs. For a static solution implemented using Newton's method, the convergence to the equilibrium takes 20 times fewer iterations than using FEM, due to a lower number of DOFs.

Table I summarizes the simulation performance for the turtle and bunny models on a 2.33GHz quad-core PC without parallel implementation. The integration of one time step, shown in the column labeled "physical time", takes an average of 7ms. We have tested two resolutions for the computation of volume integrals, yielding minor changes in the deformation. One is suited for interactive simulation while the other is more adapted to off-line rendering. Both are robust to large deformations. Further speedup can be obtained with the fast spring network presented in Section 4.3 whose complexity only depends on the number of interacting frames. However, the overall computation time also depends on the number of vertices of the visual and collision models. As expected, the spring model provides good results in relatively large linear deformations, but is not as robust for large non-linear deformations. This issue will be investigated in future work. The results reported in Table I were obtained using our current implementation, with a lot of room for future optimization. The speed greatly depends on the number of integration points, as shown by the comparison between low resolution and high resolution. At each conjugate gradient iteration used in our ODE solver, the stiffness matrix product requires the mapping of frame displacements to strain changes, then the computation of the corresponding stress changes, and finally the mapping back to frame forces (see eq. 12). The intensive part of this computation is the product with the matrix B, which currently is implemented as a dense matrix even though it is sparse. Another avenue for optimization is to further reduce the number of integration points, which is currently huge compared with the number of independent DOF, even in the low resolution version.

At each iteration, our implicit solver not only computes the force changes due to internal forces, but also due to contact forces. Frame displacements are thus also mapped to displacements of the vertices of the collision model and contact forces are mapped back to the Here also, our implementation does not exploit sparsity; moreover, a further optimization could skip vertices without contact forces. This unnecessary processing of the collision vertices explains why the computation time is higher than expected when using precomputed generalized springs. The model shown in Figure 6 has no collisions and allows an easier comparison with FEM.

As few as 5 frames (30 independent DOFs) allow us to obtain a smooth bent shape similar to the shape obtained using × 10 × 9 FEM nodes (2700 DOFs). Since our method assumes a locally uniform strain around each integration point, we found that integration points were required to obtain the same amount of bending using the same material parameters. We can quantify the difference between two methods by measuring the distance between the endpoints of the centerline, relative to the length of the bar. In this configuration, the difference between FEM and our method is only 2.5%. When the number of integration points are the same in the two methods, the computation time per iteration should be similar, provided that the same level of optimization is applied to each implementation. However, the number of iterations to converge depends on the number of independent DOFs. In this regard our method is significantly better, needing 5 fewer iterations, as can be seen in the accompanying video.

We also compared the methods at coarser resolutions, as shown in Figure 10. We used the same number of FEM nodes as integration points in the frame-based method, as with the previous comparison. Three frames were used with 2 × 2 × 2 points, and five frames with 4 × 4 × 4 points. The frame-based model becomes softer as the resolution decreases, while the FEM model becomes stiffer. We notice that the FEM converges faster than our method on the coarse example. However, our method is closer to the reference result. In the 4 × 4 × 4 example, the distance between the endpoint of the centerline in our model from the endpoint in the reference FEM is 13%, while it is 25% in the FEM. In the 2 × 2 × 2 example, the distance is 15% in our model, while it is 42% in the FEM. Our method clearly outperforms the FEM when fine meshes are required to preserve smoothness while the spatial deformation frequencies are sufficiently low so that it can be captured by a reduced number of DOFs. It can be seen as a new type of reduced model, more local and more intuitive than the reduced models produced by traditional modal analysis.

CONCLUSION

We have presented a new type of deformable model using continuum mechanics applied to objects undergoing skinning deformation fields. Our approach allows the creation of sparse meshless models with arbitrary constitutive laws, and we have demonstrated it using St. Venant-Kirchhoff materials. The models are robust to large displacements and deformations. We have shown that the behavior of objects with complex materials and geometries can be simulated using a small number of moving frames. Moreover, the deformation can be encoded in a computationally efficient way using six-dimensional mass-spring systems precomputed from the mass and stiffness matrices. Compared with FEM, adaptivity is easier because no volumetric mesh is used, and mixing rigid and deformable frames is straightforward using appropriate shape functions. Sampling is easier than with traditional particle-based meshless methods because there is no constraint on the number and on the placement of the frames. Moreover, the ability to use constrained weight fields, rather than radial basis functions, simplifies the design and simulation of objects composed of rigid and deformable parts.

Currently, the range of the local deformations obtained by inserting frames at contact points depends on the distance to the other frames. This limitation will be addressed in future work. More generally, the optimal placement of frames and the design of weight functions will have to be investigated, as well as automatic insertion and deletion of frames to handle large deformations or even topological changes. Hardware implementations would certainly accelerate the mapping of the visual and collision models, the computation of the volume integrals and allow better on-the-fly resampling. The relation between stiffness and weight functions could be exploited such as in [START_REF] Nesme | Preserving Topology and Elasticity for Embedded Deformable Models[END_REF]] to simulate fine distributions of heterogeneous materials using coarse frame sampling. Another interesting study will be the comparison, in the context of physically-based simulation, of various blending methods developed in geometric modeling as the linear combination of matrices (i.e., linear blend skinning) or the linear combination of the log-quaternion representations et al. 2007].

q i ε = [a i 0 b i 0 c i 0 w i 0 a i ε b i ε c i ε w i ε ] T
where q i 0 is the real part representing rotations and q iT ε = (t i∧ • q i 0 )/2 is the dual part representing translations to the frame origin t i = [t i x t i y t i z ] T and t i∧ = [t i x t i y t i z 0] T . Relative transformations from the rest configuration:

q i • q-1 i = Ti q i , Ti = H -[(q i 0 ) -1 ] 0 H -[(q i ε ) -1 ] H -[(q i 0 ) -1 ] H -(q) =    w c -b a -c w a b b -a w c -a -b -c w    , q =    a b c w   
where H -is the Hamilton operator that encodes right quaternion products.

The linearly blended dual quaternion is:

b = w i Ti q i = Q 0 Q ε Normalized dual quaternion: b ′ = b b = 1 Q 0 Q 0 Q ε - Q ε .Q 0 Q 0 3 0 Q 0 = [q T 0 q T ε ] T = [a 0 b 0 c 0 w 0 a ε b ε c ε w ε ] T
The rigid transformation equivalent to b ′ , that maps reference material point positions to current positions through frame displace-ments, is [START_REF] Kavan | Geometric skinning with approximate dual quaternion blending[END_REF]

]: R(b ′ ) =   1 -2b 2 0 -2c 2 0 2a 0 b 0 -2w 0 c 0 2a 0 c 0 + 2w 0 b 0 2a 0 b 0 + 2w 0 c 0 1 -2a 2 0 -2c 2 0 2b 0 c 0 -2w 0 a 0 2a 0 c 0 -2w 0 b 0 2b 0 c 0 + 2w 0 a 0 1 -2a 2 0 -2b 2 0   and t(b ′ ) = 2   -w ε a 0 + a ε w 0 -b ε c 0 + c ε b 0 -w ε b 0 + a ε c 0 + b ε w 0 -c ε a 0 -w ε c 0 -a ε b 0 + b ε a 0 + c ε w 0  
There is a linear relationship between velocity and rate of the frame: qi = (ω i • q i )/2 [Han et al. 2008]. It can be represented by a 8 × 6 matrix L i such as qi = L i Ω i and we have:

L i (q i ) = ∇ i q i = 1 2 L i0 0 L iε L i0 , L i0 =    w i 0 c i 0 -b i 0 -c i 0 w i 0 a i 0 b i 0 -a i 0 w i 0 -a i 0 -b i 0 -c i 0    L iε =   w i ε + c i 0 t i z + b i 0 t i y c i ε -w i 0 t i z -b i 0 t i x -b i ε + w i 0 t i y -c i 0 t i x -c i ε + w i 0 t i z -a i 0 t i y w i ε + c i 0 t i z + a i 0 t i x a i ε -w i 0 t i x -c i 0 t i y b i ε -w i 0 t i y -a i 0 t i z -a i ε + w i 0 t i x -b i 0 t i z w i ε + b i 0 t i y + a i 0 t i x -a i ε -b i 0 t i z + c i 0 t i y -b i ε + a i 0 t i z -c i 0 t i x -c i ε -a i 0 t i y + b i 0 t i x  
To compute the Jacobian (equation ( 7)) and the deformation gradient (equation ( 8)), we need the first-order derivatives Q, N and W that are given by (see also the long equation ( 26)) :

N(b, b ′ ) = ∂b ′ ∂b = N 0 0 N ε N 0
with N 0 = -1 Q 0 (q 0 q T 0 -I) and N ε = -1 Q 0 (q 0 q T ε + q ε q T 0 ) -Qε.Q 0 Q 0 3 (Iq 0 q T 0 )

W(p) = ∂b ∂p = Ti q i W i with W i = ∂w i ∂p W i represents the weight gradients and are numerically computed in the voxel grid by central differences.

The derivative of the deformation gradient is needed to compute the derivative of the strain (equation ( 9)). Applying the derivation chain rule and using Einstein summation convention, we obtain a third-order tensor:

∇ i F abc = w i [ ∂R ∂b ′ abd N de Ti ef L i f c ] + Q ad N de [ ∂W ∂q i ebf L i f c ] +w i [ ∂Q ∂b ′ ade N ef Ti f g L i gc ]N dh W hb + w i Q ad [ ∂N ∂b def Ti f g L i gc ]W eb
Notations can be simplified by expressing each × 3 submatrix ∇ i F k = ∂F/∂Ω ik using the operators ∆R, ∆W, ∆Q and ∆N:

∇ i F k = w i ∆R([N Ti L i ] k ) + QN Ti ∆W([L i ] k ) +w i ∆Q([N Ti L i ] k )NW + w i Q∆N([ Ti L i ] k )W
Each operator produces a matrix from a × 1 column vector

V = [V T 0 V T ε ] T = [v 0 v 1 v 2 v 3 v 4 v 5 v 6 v 7 ]
T (for instance the k th column of N Ti L i ). They are obtained by differentiating R, W, Q and N (see also long equations ( 27) and ( 28

∆N 0 = - 1 Q 0 2 [q 0 V T 0 + V 0 q T 0 + V 0 .q 0 (I -3q 0 q T 0 )] ∆N ε = - 1 Q 0 2 [q 0 V T ε + V ε q T 0 + V ε .q 0 (I -3q 0 q T 0 )] ∆N c = - 1 Q 0 2 [q ε V T 0 + V 0 q T ε + V 0 .q ε (I -3q 0 q T 0 )
-3V 0 .q 0 (q ε q T 0 + q 0 q T ε )] -2

Q ε .Q 0 Q 0 2 ∆N 0

B. EXTENSION TO OTHER MATERIALS

Volume preservation: For volume preserving materials, we need to compute the unit volume change with respect to the degrees of freedom. The unit volume is the determinant of the deformation gradient: J = det(F). Its derivative is obtained from already defined matrices:

∇ i J k = ∂J ∂Ω ik = det(F) lj F -T lj ∇ i F ljk
elements have a null strain and strain energy, resulting in a stable but undesired configuration. To circumvent this, an energy term with bulk modulus k is generally added to penalize volume change from the initial configuration. Possible expressions for the energy density, force and stiffness are:

W k 2 (J-1) 2 , Γ i = -k(J-1)∇ i J T , K ij ≈ -k∇ i J T ∇ j J
Isotropic hyperelastic material: a widely spread approach in mechanical engineering is to express the energy density W in terms of the three invariants of the right Cauchy-Green deformation tensor C = F T F. They are function of the strain and volume:

I1(C) = T r(C) = 2ε xx + 2ε yy + 2ε zz + 3 I2(C) = (T r(C) 2 -T r(C 2 ))/2 = 4(ε xx + ε yy + ε zz + ε xx ε yy + ε xx ε zz +ε yy ε zz -ε 2 xy + ε 2 xz + ε 2 yz ) + 3 I3(C) = det(C) = J 2
Their derivatives are obtained using the rows of B:

∇ i I1 = 2(B i1 + B i2 + B i3 ) ∇ i I2 = 2∇ i I1 -8(ε xy B i4 + ε xz B i5 + ε yz B i6 ) + 4(ε xx (B i2 + B i3 ) + ε yy (B i1 + B i3 ) + ε zz (B i1 + B i2 )) ∇ i I3 = 2J∇ i J
The relationship ∂W/∂I i is given by the material model (e.g., Mooney-Rivlin). Plugging this relationship into the following equations, we can compute forces and stiffness matrices as:

Γ i = -[ ∂W ∂I1 ∇ i I1 + ∂W ∂I2 ∇ i I2 + ∂W ∂I3 ∇ i I3] T K ij ≈ -[∇ i I1 T ∇ j ∂W ∂I1 + ∇ i I2 T ∇ j ∂W ∂I2 + ∇ i I3 T ∇ j ∂W ∂I3 ]
For incompressible materials, one can use reduced invariants instead. They are obtained using previously introduced variables:

Ī1 = J -2/3 I1 (29) Ī2 = J -4/3 I2 Ī3 = J

Fig. 1 .

 1 Fig. 1. Frame-based model features. Physical realism with very sparse sampling: two frames are sufficient to model a dynamically deformable bunny (Left). Intuitive animator control: The head and the ears are almost rigid (Middle), but by inserting an additional frame (Right), the ear becomes flexible.

Fig. 4 .

 4 Fig. 4. Sample points in a deformed configuration with hue mapping of the strain (hi res. bunny example: 15000 samples, 2 frames)

•B.

  Gilles et al.

Fig. 6 .

 6 Fig. 6. Left: comparison of our model (green) with 5 frames along the main axis and 10 × 10 × 9 samples in the volume, with a St-Venant-Kirchhoff model (red) with 10 × 10 × 9 points, with same loads and material parameters. Right: comparison of beams with 5, 3 and 2 frames in our model, from front to back.

Fig. 7 .

 7 Fig. 7. An object with asymmetric stiffness automatically computed based on its asymmetric shape.

Fig. 8 .

 8 Fig. 8. Modeling stiffness using weight functions. Top: (Left) The deformation frames. (Middle) Shell weights automatically computed. (Right) Weights set to model a rigid shell. Bottom: the turtle with a rigid shell.

Fig. 9 .

 9 Fig. 9. More or less global deformation produced by dynamically inserting a frame with two different weight functions

Fig. 10 .

 10 Fig. 10. Comparison of FEM simulation (in blue) and frame-based simulation (in green) for 4 × 4 × 4 (left) and 2 × 2 × 2 (right) FEM nodes and integration points. The reference high-resolution FEM is shown in red. The integration points are shown as disks. The green curve is the outline of a visual model bound to the frames.

•B.

  Gilles et al.

  a 0 x + b 0 y + c 0 z cεb 0 x + a 0 y + w 0 z -bεc 0 xw 0 y + a 0 z aε+ w 0 xc 0 y + b 0 z w 0 -c 0 b 0 -a 0 -cε + b 0 xa 0 yw 0 z -wε + a 0 x + b 0 y + c 0 z aε + w 0 xc 0 y + b 0 z bε + c 0 x + w 0 ya 0 z c w -a 0 -b 0 bε + c 0 x + w 0 ya 0 z -aεw 0 x + c 0 yb 0 z -wε + a 0 x + b 0 y + c 0 z cεb 0 x + a 0 y + w 0 z -b 0 a 0 w 0 b 0 v 1 + c 0 v 2 ) -c 0 v 3 + b 0 v 0 + a 0 v 1w 0 v 2 b 0 v 3 + c 0 v 0 + w 0 v 1 + a 0 v 2 c 0 v 3 + b 0 v 0 + a 0 v 1 + w 0 v 2 -2(a 0 v 0 + c 0 v 2 ) -a 0 v 3w 0 v 0 + c 0 v 1 + b 0 v 2 -b 0 v 3 + c 0 v 0w 0 v 1 + a 0 v 2 a 0 v 3 + w 0 v 0 + c 0 v 1 + b 0 v 2 -2(a 0 v 0 + b 0 v 1 ) yv 1 + zv 2v 7 zv 3 + yv 0xv 1 + v 6 -yv 3 + zv 0xv 2v 5 xv 3 + zv 1yv 2 + v 4 v 3 -v 2 v 1 -v 0 -zv 3yv 0 + xv 1v 6 xv 0 + yv 1 + zv 2v 7 xv 3 + zv 1yv 2 + v 4 yv 3zv 0 + xv 2 + v 5 v 2 v 3 -v 0 -v 1 yv 3zv 0 + xv 2 + v 5 -xv 3zv 1 + yv 2v 4 xv 0 + yv 1 + zv 2v 7 zv 3 + yv 0xv 1 + v 6 -v 1 v 0 v 3

Table I .

 I Summary of simulation performance

	Model	Visual	Collision	Physical	Precomp.	fps
		time (ms)	time (ms)	time (ms)	time (s)
		(# points)	(# points)	(# points/frame)
	Turtle	4	40	106	5	7
	(low res.)	(6642)	(6642)	(2313 / 7)
	Turtle	4	40	648	110	1.4
	(hi res.)	(6642)	(6642)	(15452 / 7)
	Turtle	4	40	10	110	19
	(springs)	(6642)	(6642)	(0 / 7)
	Bunny	10	7	51	1.5	15
	(low res.)	(34000)	(2500)	(2323 / 2)
	Bunny	10	7	250	34	3.7
	(hi res.)	(34000)	(2500)	(15526 / 2)
	Bunny	10	7	10	34	37
	(springs)	(34000)	(2500)	(0 / 2)

ACM Transactions on Graphics, Vol. , No. , Article , Publication date: 2010.

ACKNOWLEDGMENTS

This work is funded in part by the French National Research Agency (ANR), the Canadian Institutes of Health Research, Canada Research Chairs Program, NSERC, Peter Wall Institute for Advanced Studies, MITACS, and European project "Passport for Liver Surgery" (FP7, ICT-2007.5.3). We would like to thank Florent Falipou and Michaël Adam for preparing the examples and the video.

APPENDIX A. DETAILED MATRIX EXPRESSIONS

Dual quaternion corresponding to the current frame position:


<ca>raphaël Laurent
email: raphael.laurent@gipsa-lab.grenoble-inp.fr

Clément Moulin-Frier
email: clement.moulin-frier@inria.fr

Pierre Bessière
email: pierre.bessiere@college-de-france.fr

Jean-Luc Schwartz
email: jean-luc.schwartz@gipsa-lab.grenoble-inp.fr

Julien Diard
email: julien.diard@upmf-grenoble.fr

<CT>Integrate, yes, but what and how? A computational approach of perceptuo-motor fusion in speech perception

We consider a computational model comparing the possible roles of "association" and "simulation" in phonetic decoding, demonstrating that these two routes can contain similar information in some "perfect" communication situations and highlighting situations where their decoding performance differs. We conclude that optimal decoding should involve some sort of fusion of association and simulation in the human brain.

<C-Text begins>

In their target article, Pickering & Garrod (P&G) propose an ambitious model of language perception and production. It is centered on three main ingredients. First, it considers the complete hierarchy of layers of language processing, from message to semantics to syntax to phonology and finally, to speech. Second, it features predictive forward models, so that temporally extended sequences, such as whole sentences and dialogues, can be processed. Third, it features dual processing routes, the "association" route and "simulation" route, so that auditory and motor knowledge can be involved simultaneously, rejecting the classic dichotomy between perception and action processes.

In this commentary, we set aside the temporal and hierarchical aspects, and focus on the domain of speech perception and production, where sequences are typically short (e.g., syllable perception and production), and processing limited to phonological decoding.

Even in this more restricted field, the age-old debate between purely motor-based accounts and purely sensory-based accounts of perception and production now appears to be a false dilemma [START_REF] Schwartz | The perception-for-actioncontrol theory (PACT): A perceptuo-motor theory of speech perception[END_REF]. Indeed, neurophysiological and behavioral evidence strongly suggests a dual route account of information processing in the central nervous system, with both a direct, associative route and an indirect, simulation route.

The target article amply documents the evidence, we do not repeat examples here.

In our view, the debate is now shifting toward the issue of the functional role of each route and their integration. That is to say, a central question of the debate is to ask what is integrated and how integration proceeds in the human brain.

We would argue that conceptual models such as proposed in the target article would unfortunately have a difficult time bringing light to these questions. To support this argument, we consider the question of perceptual decoding of phonetic units, for which we have developed a computational framework [START_REF] Moulin-Frier | Adverse conditions improve distinguishability of auditory, motor and perceptuo-motor theories of speech perception: An exploratory Bayesian modeling study[END_REF]) based on Bayesian programming (Bessière et al. 2008;[START_REF] Colas | Common Bayesian models for common cognitive issues[END_REF][START_REF] Rl Lebeltel | Bayesian robot programming[END_REF].

With this framework, various models of speech perception can be simulated and quantitatively compared. One model is purely auditory, exploiting what P&G call "association." A second model is purely motor, exploiting what they call "simulation." A third one is sensory-motor, integrating the association and simulation processes.

All of these models can then be implemented and compared in various experimental configurations. Three major results emerge from such comparisons.

<<NL>(1) Under some hypotheses, with perfectly identified communication noise and no difference between motor repertoires of the speaker and the listener (i.e., when conditions for speech communication are "perfect"), motor and auditory theories are indistinguishable. Therefore, the "association" and "simulation" routes provide exactly the same information in these perfect communication conditions. The reason is that, in our learning scenario, the auditory classifier is learned by association from data obtained through a motor production process, and possesses enough mathematical power of expression.

This casts an interesting light on the question of what information is encoded in the association and simulation routes: Labeling a box as an "association" route, in a conceptual model, is not enough to be certain that it is different, from an information processing point of view, from another box of the model. Computational descriptions however, by virtue of rigorous mathematical notation, have to be precisely defined, and their content can be systematically assessed. This also explains why behavioral evidence has historically not been able to discriminate between motor and auditory theories of perception and production: They are sometimes simply indistinguishable. Unfortunately, we believe this difficulty was not avoided in the target article, in particular when P&G detail experimental evidence for their model (e.g., target article, sect. 3.2.1, para. 7, "these four studies support forward modeling, but they do not discriminate between prediction-by-simulation and prediction-by-association"; and sect. 3.2.3, para. 6, "all of these findings provide support for the model of prediction-by-simulation [...]. Of course, comprehenders may also perform prediction-by-association […].").

(2) In the general case where "perfect conditions" for communication are not met, mathematical comparison of the models emphasizes the respective roles of motor and auditory knowledge in various conditions of speech perception in adverse conditions. Therefore, the information provided by the "association" and "simulation" routes is more or less distinct and prominent depending on the communication conditions. In other words, this demonstrates that adverse conditions provide leverage for discriminating hypotheses about the perceptual and motor processes involved. This is convergent with recent findings from neuroimaging and transcranial magnetic stimulation (TMS) studies [START_REF] D'ausilio | The role of the motor system in discriminating degraded speech sounds[END_REF][START_REF] Meister | The essential role of premotor cortex in speech perception[END_REF][START_REF] Zekveld | Top-down and bottom-up processes in speech comprehension[END_REF]), as well as computational studies [START_REF] Castellini | The use of phonetic motor invariants can improve automatic phoneme discrimination[END_REF]).

(3) In any case, sensory-motor fusion provides better perceptual performance than pure auditory or motor processes. Therefore, complementarities of information provided by the "association" and "simulation" routes could be efficiently exploited in the framework of integrative theories such as those hinted at in the discussion of the target article. It is now obvious in the field of audiovisual perception that auditory and visual cues are complementary, with a great deal of work already done on sensor fusion. In our opinion, comparable work can now be done on how to integrate auditory and motor processes in speech perception. In this view, the proposal by P&G that "comprehenders emphasize whichever route is likely to be more accurate" (sect. 4, para. 6) can be regarded as a first candidate model, which would have to be made mathematically precise and compared with alternative explanations, possibly driven by neuroanatomical findings (e.g., both auditory and motor processes are performed automatically in parallel and compete, or they both bring information in an ongoing fusion process, etc.).

<NL ends>

An obvious challenge, of course, is to bridge the gap between computational approaches such as ours, which are usually restricted to isolated syllable production and perception, and conceptual models as proposed in the target article, that tackle continuous flows of speech and consider semantic, syntactic and phonology layers of processing.

However, in our view, the main challenge for future studies is first to assess what kind of information is present in "association" and "simulation" routes, and second, to better understand how computational fusion models, describing the integration of these two routes, can account for experimental neurocognitive data.

  <C-Text ends> <RFT>References[Raphaël Laurent, Clément Moulin-Frier, Pierre Bessière, Jean-Luc Schwartz and Julien Diard][RL] <Refs> Bessière, P., Laugier, C. & Siegwart, R. ed. (2008) Probabilistic reasoning and decision making in sensory-motor systems, volume 46 of Springer tracts in advanced robotics. Springer-Verlag[RL]