Abstract:
Animating strands (such as long hair), for movies, videos, etc. is accomplished using computer graphics by use of differential algebraic equations. Each strand is subject to simulation by defining its motion path, then evaluating dynamic forces acting on the strand. Collision detection with any objects is performed, and collision response forces are evaluated. Then for each frame a differential algebraic equations solver is invoked to simulate the strands.

Description:
FIELD OF THE INVENTION 
       [0001]    This disclosure is directed to simulating the motion of thin and long objects (“strands” such as hair) in the field of computer graphics and more specifically for animation. 
       BACKGROUND OF THE INVENTION 
       [0002]    The technical problem of simulating by computer graphics the dynamics of strands in terms of the forward dynamics problem for rigid multi-body serial chains is well known; see S. Hadap,  Hair Simulation,  (Ph.D. Thesis, Dept. of Informatics, Univ. of Geneva, Switzerland, 2003). One method uses recursive algorithms for solving the forward dynamics problem of rigid multi-body serial chains. Such algorithms are well known; see R. Featherstone,  Robot Dynamics Algorithms,  (Kluwer Academic Publishers, 1987). Recursive and residual algorithms for the simulation of multi-body systems are also known; see J. I. Rodriguez, J. M. Jimenez, F. J. Funes, J. G. de Jalon, “Recursive and Residual Algorithms for the Efficient Numerical Integration of Multi-Body Systems,”  Multibody System Dynamics,  11(4), (Springer Netherlands, May 2004): 295-320. 
       SUMMARY 
       [0003]    Examples of objects whose motion can be animated well using the present method and apparatus include long hair, wavy or curly hair, braids, animal ears, tails, foliage, and jewelry. The method is implemented in computer software, of the type that animation artists use in production of computer animated motion pictures or videos. 
         [0004]    The present method advantageously simulates accurately the physical effects of inertia, centrifugal and Coriolis forces, gravity, wind or other force fields on strands. The method handles correctly non-zero bend and twist in the rest shape of the strands, as well as changes in bend and twist during the simulation. In many systems for simulating strands the joint constraints between the strand elements (“links”) are satisfied only approximately, a procedure which is prone to numerical errors. The present method enforces these constraints exactly. The method uses a formulation based on the minimal set of coordinates needed to specify the strands&#39; configuration. This simplifies the task of specifying the rest shape and ensures optimal use of computing resources. The method employs algorithms for kinematics and dynamics computations that are optimal; the number of operations grows linearly in the number of links. The method correctly computes collisions between strands and scene geometry. The method simulates correctly the effects of constraints and collisions, including static and dynamic friction associated with collisions. The method uses analytical methods for simulating the effects of constraints. Using the present method, it is possible to art-direct intuitively and effectively the physically based simulation of the animated objects. 
         [0005]    In the present method, the equations of motion for strands are formulated as a system of differential algebraic equations (DAE). Analytical constraints are computed as external forces, and thus the user is able to impose analytical constraints in the framework of DAE. An approximate representation of static and dynamic friction as unilateral constraints is also provided. For solving the system of DAEs corresponding to strand dynamics, the DAE solver software DASPK well known in the computer field is used in one embodiment. 
     
    
     
       BRIEF DESCRIPTION OF THE DRAWINGS 
         [0006]      FIG. 1  is a flowchart of the present method. 
           [0007]      FIG. 2  shows a strand as a multi-body system. 
           [0008]      FIG. 3  shows constraints acting on the  FIG. 2  system. 
           [0009]      FIG. 4  shows collisions with an obstacle. 
           [0010]      FIG. 5  shows an apparatus for performing the  FIG. 1  method. 
       
    
    
     DETAILED DESCRIPTION 
       [0011]    The present method uses a formulation of strand dynamics in terms of Differential Algebraic Equations (“DAE”). DAE are well known; see K. E. Brenan, S. L. Campbell, L. R. Petzold,  Numerical Solution of Initial - Value Problems in Differential - Algebraic Equations,  (SIAM, 1996). The motion equations corresponding to strands are highly non-linear and have large numerical stiffness. Non-linear behavior is due to Coriolis forces and the non-linear elastic model used to limit high deformations. Numerical stiffness in the simulation of strand dynamics is caused by small rotational inertia due to thin geometry, large bending and torsional stiffness-to-mass ratio, and non-straight rest shape. 
         [0012]    The DAE approach is used here to improve the numerical stability of the strand system in the presence of non-linearity and numerical stiffness. Improved numerical stability allows larger time steps in the simulation, with a corresponding improvement of the system&#39;s interactive response time.  FIG. 1  thus shows the present method  12  in flowchart form, as would be carried out by a computer program (software) in the context of animation of strands. In the first block  16 , a conventional data structure is provided representing the strands. For instance, the strands are represented as instances of a class (‘class’ is used here in the sense of an object-oriented programming language construct, see for example B. Stroustroup “The C++ Programing Language”, Addison-Wesley 2001). The data that is being grouped by the class that represents a strand includes:
       position, velocity, acceleration, orientation, angular velocity and angular acceleration of the strand as a whole.   position, velocity, acceleration, orientation, angular velocity and angular acceleration of an animated strand that is used as reference for art-direction.   parameters used to implement the elasticity model of the strand: torsional stiffness and anisotropic bend stiffness coefficients.   a vector of instances of a class representing the segments that make up the strand.   a vector of data structures representing joints between segments.   a vector of data structures representing joint velocities.
 
Next at  20 , motion paths in space are defined for the base link (see  FIG. 2 ) for each strand of interest. Next at  22 , possible dynamic forces are computed (see  FIG. 3 ). Then at  26 , collisions with any obstacles are detected (see  FIG. 4 ) Then at  30 , the response forces from the detected collisions are computed. Then at  34 , the DAE solver is invoked to solve the resulting strand simulation for the current animation frame. At  38 , it is determined if any additional frames are to be simulated (animated). Detail of each block of  FIG. 1  is provided hereinafter.
       
 
         [0019]    DAEs are a well known general form of differential equations, where the equations are expressed in implicit rather than explicit form: 
         [0000]        F ( y,{dot over (y)},t )=0   (1) 
         [0020]    Compare Eq. (1) with a differential equation given in the more familiar explicit form: 
         [0000]        {umlaut over (q)}=M   −1   G ( q, {dot over (q)}, t )   (2) 
         [0000]    where the highest order derivatives are separated on the left hand side of the equations. Because of this freedom, DAEs form a much larger class of differential equations than the class of explicit differential equations. 
         [0021]    The dynamic state of a strand is completely described for purposes of animation by the state vector: 
         [0000]    
       
         
           
             
               
                 
                   y 
                   = 
                   
                     [ 
                     
                       
                         
                           q 
                         
                       
                       
                         
                           
                             q 
                             . 
                           
                         
                       
                     
                     ] 
                   
                 
               
               
                 
                   ( 
                   3 
                   ) 
                 
               
             
           
         
       
     
         [0000]    where q are is a vector of joint variables and {dot over (q)} is a vector of the joint velocities. Both vectors have size n, where n is the total number of degrees of freedom allowed by the joints which couple adjacent strand links. In general, the total number of degrees of freedom n is different from the number of links m, because joints may have more than one degree of freedom. 
         [0022]    The DAE solver advances the solution from state [y, {dot over (y)}], to [y, {dot over (y)}] t+1 . The well known DASPK software is used as the DAE solver in one embodiment. The DAE solver requires a definition of the residue function: 
         [0000]      rs(t)≡F(y(t), {dot over (y)}(t), t)   (4) 
         [0023]    There are 2n equations in the DAE system describing the dynamics of a strand: the first n are the force balance equations for the n links, see Equations (8) and (10) below. The last n are consistency conditions between 
         [0000]    
       
         
           
             
               y 
               = 
               
                 
                   
                     [ 
                     
                       
                         
                           q 
                         
                       
                       
                         
                           
                             q 
                             . 
                           
                         
                       
                     
                     ] 
                   
                    
                   
                       
                   
                    
                   and 
                    
                   
                       
                   
                    
                   
                     y 
                     . 
                   
                 
                 = 
                 
                   [ 
                   
                     
                       
                         
                           q 
                           . 
                         
                       
                     
                     
                       
                         
                           q 
                           ¨ 
                         
                       
                     
                   
                   ] 
                 
               
             
             , 
           
         
       
     
         [0000]    due to {dot over (q)} appearing in both. 
         [0024]    A recursive formulation of the kinematic and dynamic equations is used here. This is well known; see R. Cottle, J. S. Pang, R. Stone,  The Linear Complementarity Problem,  (Academic Press, 1992). The recursive formulation has the advantages that it achieves linear time in n, while using reduced coordinates and thus satisfying exactly all the constraints between the bodies making up the strands. 
         [0025]    In the recursive formulation, the equations of motion of the system are expressed only implicitly as recurrence relations between physical quantities. The computation of the residual function rs(t) has two recursive steps: forward kinematics: given the velocity and acceleration of the root link v 0 , a 0  and the state vector [q, {dot over (q)}] determine the velocities and accelerations of all links; and forward dynamics: given all the above plus the inertia tensors of the links and external forces determine the residual (unbalanced) forces for all links. 
         [0026]    To describe the dynamics of the multi-body system, the present method uses spatial algebra and the associated well known notation, see R. Featherstone,  Robot Dynamics Algorithms,  (Kluwer Academic Publishers, 1987). Consider a serial branched multi-body system (a strand) having m links connected by m joints as shown in  FIG. 2 . There are no loops in this multi-body system. The base link is denoted by link 0  and is typically constrained to a motion path (see  FIG. 1 ). This means that the velocity and acceleration v 0 , a 0  of the base link are specified by the animation system instead of being a result of the dynamics simulation. The other links are numbered in a first manner, i.e. ensuring that the parent of a link has always a smaller index than the link itself. 
         [0027]    The kinematics equations in spatial notation between link j  and its parent link i  are: 
         [0000]        {circumflex over (v)}   j   ={circumflex over (X)}   j   {circumflex over (v)}   i   +ŝ   j   {dot over (q)}   j    (5) 
         [0000]        â   j   ={circumflex over (X)}   j   â   i   +ŝ   j   {umlaut over (q)}   j   +{circumflex over (v)}   j   {circumflex over (×)}ŝ   j   {dot over (q)}   j    (6) 
         [0028]    The spatial velocity {circumflex over (v)} j , spatial acceleration â j  and spatial inertia matrix Î j  are all defined in the local frame {circumflex over (F)} j  that is rigidly attached to the principal inertia axis of link j . The joint variables q j  define the spatial transformation {circumflex over (X)} j  that may be used to transform spatial quantities defined in the parent&#39;s frame {circumflex over (F)} i  to the frame {circumflex over (F)} j . The spatial matrix of the joint axis connecting link j  to its parent link is ŝ j . The computation of the forward kinematics of the system amounts to applying Equations (5-6) for j=1 . . . m while starting from the given velocity {circumflex over (v)} 0  and acceleration â 0  of the root link. 
         [0029]    In the DAE framework, computing forward dynamics means computing a scalar residue value for each degree of freedom of a joint. These residue values are associated with the links&#39; accelerations. First, the force balance condition for each link are used to obtain residual forces. This computation is again recursive and, due to data dependencies, starts from the outermost links and proceeds towards the base link. 
         [0030]    Consider one of the outermost links and denote it by link m . The force balance equation for link m  is: 
         [0000]        {circumflex over (r)}ŝ   m   =Î   m   â   m   +{circumflex over (v)}   m   {circumflex over (×)}Î   m   {circumflex over (v)}   m   −{circumflex over (f)}   m    (7) 
         [0000]    where {circumflex over (f)} m  is the external spatial force acting on link m  and {circumflex over (r)}ŝ m  is the spatial force residue. 
         [0031]    The residue for link m  is obtained by projecting the force residue {circumflex over (r)}ŝ m  onto the associated joint&#39;s motion subspace defined by the spatial joint axis ŝ m . The joint actuation force associated with the strand deformation model Q m  is accounted for, which also acts in the joint&#39;s motion subspace (the deformation model is explained below): 
         [0000]        rs   m   =ŝ   m   S   {circumflex over (r)}ŝ   m   −Q   m    (8) 
         [0032]    Here rs m  has a number of components equal to the dimension of the joint&#39;s motion subspace, i.e. to the number of degrees of freedom of the joint connecting link m  to the rest of the multi-body system. 
         [0033]    The force balance equation for link i , which is the parent of link j , is: 
         [0000]    
       
         
           
             
               
                 
                   
                     
                       rs 
                       ^ 
                     
                     i 
                   
                   = 
                   
                     
                       
                         
                           I 
                           ^ 
                         
                         i 
                       
                        
                       
                         
                           a 
                           ^ 
                         
                         i 
                       
                     
                     + 
                     
                       
                         
                           v 
                           ^ 
                         
                         i 
                       
                        
                       
                         × 
                         ^ 
                       
                        
                       
                         
                           I 
                           ^ 
                         
                         i 
                       
                        
                       
                         
                           v 
                           ^ 
                         
                         i 
                       
                     
                     - 
                     
                       
                         f 
                         ^ 
                       
                       i 
                     
                     + 
                     
                       
                         ∑ 
                         j 
                       
                        
                       
                         
                           
                             X 
                             ^ 
                           
                           j 
                           
                             - 
                             1 
                           
                         
                          
                         
                           
                             rs 
                             ^ 
                           
                           j 
                         
                       
                     
                   
                 
               
               
                 
                   ( 
                   9 
                   ) 
                 
               
             
           
         
       
     
         [0000]    where {circumflex over (F)} i  is the external spatial force acting on link i . Note that the balance equation accounts for the effect of the spatial force residual of link j . To bring this force in link i &#39;s frame the inverse of the transform {circumflex over (X)} j  is used. The sum is over the children of link i . As above, the DAE residual is the projection of the spatial force residual in the joint&#39;s motion space: 
         [0000]        rs   i   =ŝ   i   S   {circumflex over (r)}ŝ   i   −Q   i    (10) 
         [0034]    This process is repeated for each link, progressing towards the base link of the multi-body system. 
         [0035]    The residue associated with joint velocities is simply given by consistency conditions in the representation of the dynamic state of the strand. There is {dot over (q)} appearing in both y and {dot over (y)}, which leads to: 
         [0000]        rs   i+n   ={dot over (q)}   i *−{dot over (q)} i    (11) 
         [0000]    where {dot over (q)} i  and {dot over (q)} i  are the joint velocities in y and {dot over (y)}, respectively. 
         [0036]    The base joint in the multi-body system in  FIG. 2  is typically constrained to some prescribed motion path. In practice, it is often necessary to impose some additional constraints on the system, e.g. constraining a point on some other link to a motion path, or constraining a joint to an animated value in time, see  FIG. 3 . These constraints (“constraint path” in  FIG. 3 ) are transient in nature and often introduce cyclic dependency in the system. The effects of such additional constraints are computed by augmenting the DAE based multi-body formulation. A method for computing the effects of such constraints is known, see E. Kokkevis, “Practical Physics for Articulated Characters,”  Proceedings of Game Developers Conference  2004 (2004) and belongs to a class of techniques called analytical methods. The main problem with this known methodology is that it is not compatible with a DAE solver. A DAE solver integrates the system without explicitly calculating the state derivatives, whereas the known method directly alters the explicitly calculated derivatives to satisfy analytical constraints. 
         [0037]    The present approach computes analytical constraints as external forces. This allows imposing analytical constraints in the framework of DAE. Consider a point constraint p j  as depicted in  FIG. 3 . The trajectory (path) of p j , starting with the initial conditions, can be uniquely defined by the time varying acceleration a j   d n j . As discussed above, here one does not directly evaluate the state derivative vectors {dot over (y)} t+i  in order to integrate the system. Therefore, one cannot simply enforce the acceleration constraint, by directly altering the state derivatives {dot over (y)} t+i  as known. The present method enforces the constraint by applying an external force instead. 
         [0038]    Before every DAE integration step, the precise nature of the force f j n j , is determined analytically using the above known method. The unit constraint direction n j  is treated as constant and is updated for every integration step. There is a linear relationship between the magnitude of the applied force f j  and the resulting desired acceleration a j   d : 
         [0000]    
       
         
           
             
               
                 
                   
                     a 
                     j 
                     d 
                   
                   = 
                   
                     
                       
                         
                           ∂ 
                           
                             a 
                             j 
                           
                         
                         
                           ∂ 
                           
                             f 
                             j 
                           
                         
                       
                        
                       
                         f 
                         j 
                       
                     
                     + 
                     
                       a 
                       j 
                       0 
                     
                   
                 
               
               
                 
                   ( 
                   12 
                   ) 
                 
               
             
           
         
       
     
         [0000]    where, a j   0  is the acceleration in the direction n j  before the force is applied. 
         [0039]    If there is another constraint point p i  with force having magnitude f i  in the direction n i , the resulting accelerations a i   d  and a j   d  will be given by the following linear system: 
         [0000]    
       
         
           
             
               
                 
                   
                     [ 
                     
                       
                         
                           
                             a 
                             i 
                             d 
                           
                         
                       
                       
                         
                           
                             a 
                             j 
                             d 
                           
                         
                       
                     
                     ] 
                   
                   = 
                   
                     
                       [ 
                       
                         
                           
                             
                               
                                 ∂ 
                                 
                                   a 
                                   i 
                                 
                               
                               
                                 ∂ 
                                 
                                   f 
                                   i 
                                 
                               
                             
                           
                           
                             
                               
                                 ∂ 
                                 
                                   a 
                                   i 
                                 
                               
                               
                                 ∂ 
                                 
                                   f 
                                   j 
                                 
                               
                             
                           
                         
                         
                           
                             
                               
                                 ∂ 
                                 
                                   a 
                                   j 
                                 
                               
                               
                                 ∂ 
                                 
                                   f 
                                   i 
                                 
                               
                             
                           
                           
                             
                               
                                 ∂ 
                                 
                                   a 
                                   j 
                                 
                               
                               
                                 ∂ 
                                 
                                   f 
                                   j 
                                 
                               
                             
                           
                         
                       
                       ] 
                     
                     + 
                     
                       [ 
                       
                         
                           
                             
                               a 
                               i 
                               0 
                             
                           
                         
                         
                           
                             
                               a 
                               j 
                               0 
                             
                           
                         
                       
                       ] 
                     
                   
                 
               
               
                 
                   ( 
                   13 
                   ) 
                 
               
             
           
         
       
     
         [0040]    Generalizing, for m such constraints one needs to determine the vector of f ε R m  of unknown force magnitudes by solving the following linear system: 
         [0000]        Kf+a   0 =0.   (14) 
         [0041]    The Jacobian matrix K ε R m×m  (well known in the field) can be evaluated by applying unit test forces at each constraint and evaluating the changes in accelerations. An efficient procedure to evaluate the Jacobian using the framework of Featherstone&#39;s Articulated Body Method is in Kokkevis, referenced above. The constraint forces thus determined are applied to the multibody system over the next integration step via the external force {circumflex over (f)}, as disclosed above. 
         [0042]    One replaces the constraint direction n by a spatial vector {circumflex over (n)} to generalize the type of the constraints that can be represented. For example, {circumflex over (n)} can represent also joint acceleration constraints. One further extends the constraint system to include unilateral constraints such as collisions, friction and joint limits by posing it as a Linear Complementarity Problem (LCP) which is well known in the field. An LCP solver takes a linear system written in the form: 
         [0000]        Mz+q=w    (15) 
         [0000]    and solves for z such that for every i: 
         [0000]      z i ≧0, w i ≧0 and z i w i =0.   (16) 
         [0043]    The following abbreviates these conditions as: 
         [0000]      z i ≧0 complementary to w i ≧0.   (17) 
         [0044]    The LCP of the present constraint system can then be written as: 
         [0000]        Kf+a   0   −a   d ≧0 complementary to f≧0   (18) 
         [0045]    The LCP states that the forces f, applied only in positive constraint direction, would strive to make the resulting constraint accelerations a equal to desired acceleration a d . However, force f i  will be zero if the resulting constraint acceleration a i  is greater than desired acceleration a i   d . Iterative methods are used to solve the LCP. 
         [0046]    Using forces for enforcing the constraints has an advantage here. If the forces are indeterminate, they are projected into the multi-body&#39;s motion null-space, thus always giving valid configuration, without any jumps or jittering in the simulation. Further, as the forces are determined analytically, as compared to, e.g., penalty-based formulation, they are small for most types of the constraints. Thus they are well within the stability zone of the integrator when taking 4 to 8 time-steps per animation frame. The only exception to this is the velocity impulse constraint, discussed below. As the constraint may not be satisfied accurately, one augments the constraint acceleration by a proportional-derivative form. To exemplify, for a positional constraint, the desired acceleration and the constraint direction are: 
         [0000]    
       
         
           
             
               
                 
                   
                     a 
                     i 
                     d 
                   
                   = 
                   
                     
                       
                         p 
                         ¨ 
                       
                       i 
                       d 
                     
                     + 
                     
                       
                         K 
                         p 
                       
                        
                       
                         ( 
                         
                           
                             p 
                             i 
                             d 
                           
                           - 
                           
                             p 
                             i 
                           
                         
                         ) 
                       
                     
                     + 
                     
                       
                         K 
                         d 
                       
                        
                       
                         ( 
                         
                           
                             
                               p 
                               . 
                             
                             i 
                             d 
                           
                           - 
                           
                             v 
                             i 
                           
                         
                         ) 
                       
                     
                   
                 
               
               
                 
                   ( 
                   19 
                   ) 
                 
               
             
             
               
                 
                   
                     a 
                     i 
                     d 
                   
                   = 
                   
                      
                     
                       a 
                       i 
                       d 
                     
                      
                   
                 
               
               
                 
                   ( 
                   20 
                   ) 
                 
               
             
             
               
                 
                   
                     n 
                     i 
                     d 
                   
                   = 
                   
                     
                       a 
                       i 
                       d 
                     
                     
                       a 
                       i 
                       d 
                     
                   
                 
               
               
                 
                   ( 
                   21 
                   ) 
                 
               
             
           
         
       
     
         [0000]    where {umlaut over (p)} i   d , {dot over (p)} i   d , p i   d  are respectively the acceleration, velocity and position of the constraint path and p i , v i  are the current position and velocity of the constraint. 
         [0047]    The unilateral position constraints disclosed above are used to develop a collision response model for the multi-body system. Collision detection is well known in computer graphics. For brevity, the following describes only the requirements from the collision detection system as needed for the present method. Between the current configuration given by the state vector y t  and the extrapolated configuration using the derivative vector {dot over (y)} t  and integration time step h, one finds all the points on the multi-body system that would collide with the obstacles.  FIG. 4  shows two such collision positions—point p i  has already penetrated the obstacle and point p j  is about to collide. There may be more than one colliding point for a link. Let n i  be the collision normal, with direction pointing away from the obstacle, and a i  and v i  be the collision accelerations and velocities respectively, relative to the obstacle. 
         [0048]    One applies collision response in two steps—contacts and impacts. First, the unilateral constraints are computed that would prevent collision points from accelerating towards the obstacle. This is followed by computation of velocity impulses that would prevent collision points from moving towards the obstacle. 
         [0049]    One can decompose the collision acceleration and the collision velocity into the normal components a ni , v ni  and tangential components a ti , v ti . To prevent any acceleration towards the obstacle, a unilateral constraint is inserted along the collision normal direction n i . The unilateral constraint corresponding to the friction acts in the tangent plane defined by the collision normal. Instead of sampling the tangent plane into discrete set of tangents to formulate a complex and numerically expensive friction model based of the discrete frictional cone, the present method uses a technique as follows. The unilateral constraint direction corresponding to friction is set as: 
         [0000]    
       
         
           
             
               
                 
                   
                     t 
                     i 
                   
                   = 
                   
                     unit 
                      
                     
                       ( 
                       
                         
                           a 
                           ti 
                         
                         + 
                         
                           
                             v 
                             ti 
                           
                           h 
                         
                       
                       ) 
                     
                   
                 
               
               
                 
                   ( 
                   22 
                   ) 
                 
               
             
           
         
       
     
         [0000]    If both a ti  and v ti  are zero, the previously determined tangent vector is used for the contact. 
         [0050]    Finally, the LCP formulation corresponding to the two unilateral constraints in the direction n i  and t i  at collision position p i  is: 
         [0000]        a   ni   −a   ni   d ≧0 complementary to  f   ni ≧0   (23) 
         [0000]      μ i   f   ni   −f   ti ≧0 complementary to λ i ≧0   (24) 
         [0000]      ( a   ti   −a   ti   d )+λ i ≧0 complementary to f ti ≧0   (25) 
         [0051]    The desired normal acceleration a ni   d  is set to be proportional to the penetration depth d i  if the point is penetrating, see Equations 19-21, or zero if the collision point is outside the obstacle. The desired tangential acceleration a ti   d  is set to (−v ti /h). The LCP formulation will compute the required amount of normal force f ni  to remove the normal acceleration a ni . The tangential force f ti  will be at most μ i f ni , and it will try to remove any tangential non-zero velocity component (the dynamic friction case) or if the tangential velocity is zero it will try to remove any tangential acceleration (the static friction case). 
         [0052]    The present method uses impulses to arrest the collision normal velocity v ni . Only those contacts are considered that have the normal velocity component v n i&lt;0. For the impulse computations, one can use the same acceleration constraints discussed above by setting: 
         [0000]        a   ni   d =−(1 +v )v ni ,   (26) 
         [0000]    where v is the coefficient of restitution associated with the collision. Instead of applying potentially large forces, one alters the joint velocities {dot over (q)} t . Although theoretically this would invalidate the consistent solution set y t , {dot over (y)} t , in practice the solver is tolerant to the error. 
         [0053]    The dynamic strand is composed from a hierarchy of input curves defined in a local frame, that defines the initial rest shape of the corresponding multi-body system. One embodiment of the present method provides the user (animator) with high-level control over the direct numerical simulation by means of relevant physical parameters of the dynamic strand, such as mass per unit length, strand radius, bend stiffness and damping, torsional stiffness and damping, gravity, and air drag. 
         [0054]    The user can animate all the parameters and specify their length-wise anisotropic variation. The collision parameters coefficient of restitution, coefficient of static friction, and coefficient of dynamic friction are defined per obstacle surface. The strand may have additional anisotropic weights over the collision parameters, along with their length-wise variation. 
         [0055]    Three-degrees-of-freedom spherical joints are used in one implementation of the present strand system. The joint variable of the i-th joint is expressed as a quaternion q i  and the joint velocity as a vector w i . One decomposes the quaternion defining the relative transformation between two links into a twist component θ t  around the local y-axis and a pure bend component θ b  around a bend axis b. A nonlinear stiffness model is defined as follows: 
         [0000]    
       
         
           
             
               
                 
                   
                     Q 
                     b 
                   
                   = 
                   
                     
                       K 
                       b 
                     
                      
                     b 
                      
                     
                         
                     
                      
                     
                       tan 
                        
                       
                         ( 
                         
                           
                             ( 
                             
                               
                                 θ 
                                 b 
                               
                               - 
                               
                                 θ 
                                 b 
                                 0 
                               
                             
                             ) 
                           
                           2 
                         
                         ) 
                       
                     
                   
                 
               
               
                 
                   ( 
                   27 
                   ) 
                 
               
             
             
               
                 
                   
                     Q 
                     t 
                   
                   = 
                   
                     
                       K 
                       t 
                     
                      
                     y 
                      
                     
                         
                     
                      
                     
                       tan 
                        
                       
                         ( 
                         
                           
                             ( 
                             
                               
                                 θ 
                                 t 
                               
                               - 
                               
                                 θ 
                                 t 
                                 0 
                               
                             
                             ) 
                           
                           2 
                         
                         ) 
                       
                     
                   
                 
               
               
                 
                   ( 
                   28 
                   ) 
                 
               
             
           
         
       
     
         [0000]    where θ b   0  and θ t   0  describe the rest configuration, K t  is the torsional stiffness coefficient and K b  is the anisotropic bend stiffness coefficient. The response model is almost linear for small deformations. However, the non-linear response prevents excessive deformations. 
         [0056]    General external force fields are supported using the plug-in aspect of the well known Maya animation software and that of other animation software. The user can attach any complex combination of time-varying fields such as wake, turbulence, fluid simulations and general event driven scripted force fields. The user can further specify length-wise anisotropic weights for the external force fields. The user can optionally include these fields in computing the Jacobians numerically, as discussed above. 
         [0057]    Often the dynamical simulation is driven by very extreme and brisk animated motions. Although the present method is a robust solution able to cope with such constraints, animators sometimes want the resulting motion to be selectively less violent (abrupt). A time scale β is provided here to control the amount of energy pumped into the strand system. β is a factor by which the velocity and acceleration of the base joint can be scaled and it is typically between zero and one; however the user can set it to values larger than one to accentuate the motion. 
         [0058]    The local dynamics γ is another similar parameter which blends out velocity and acceleration of some local dynamics reference frame. 
         [0000]        â   0 =β( â   0   −γâ   ref )   (29) 
         [0000]        {circumflex over (v)}   0 =β( {circumflex over (v)}   0   −γ{circumflex over (v)}   ref )   (30) 
         [0059]    One scenario that is frequent in animation is a braid or long hairs of a character that fly away when the character starts running or rides a horse. In this case the local dynamics reference frame can be simply set to the character&#39;s hip joint, and with an appropriate value of the local dynamics parameter one can control the amount of flyaway in the result of the simulation. 
         [0060]    Often the animator wants to reach a certain pose for an animated character at a certain time to make the character expressive. In order to support this posable dynamics requirement the present method inserts “springs” between a dynamic strand and the desired animated pose at the tip of each segment. This results in introducing penalty based constraints. The user will typically animate the stiffness and damping of the springs, to change gradually from a purely dynamics simulation to key-frame animation and back. 
         [0061]    The present method is typically implemented in the form of computer software (code) executed on a computer. Coding the software, in light of this disclosure, would be routine; a suitable computer language is C+.  FIG. 5  shows the structure of such a program  68  in terms of a group of modules each corresponding to activity in one of the steps of the  FIG. 1  method. Thus in  FIG. 5 , memory  70  stores the strand data structures. A motion path component  74  (module) defines the motion path for the base link of each strand. A dynamics force module  78  evaluates the dynamic forces for each strand. A collision detector module  82  performs the collision detection for the strands, while collision response module  86  evaluates the collision response forces DAE solver  88  solves the differential algebraic equations for all strands in the current frame. The software itself would typically be stored on a computer readable medium (memory), such as a tape, disk, hard disk, flash memory, etc., for execution by the computer processor. 
         [0062]    This disclosure is illustrative and not limiting; further modification will be apparent to those skilled in the art in light of this disclosure and are intended to fall within the scope of the appended claims.