Real-time, interactive animation of deformable two- and three-dimensional objects

[0002] This application claims the benefit of U.S. Provisional Application No. 60/181,381, filed Feb. 9, 2000, which is hereby fully incorporated by reference herein as though set forth in full.

[0004] 1. Field of the Invention

[0005] This invention generally relates to real-time, interactive animation of objects, and, more specifically, to real-time, interactive animation of two- and three-dimensional deformable objects.

[0006] 2. Background

[0007] A mass-spring system is a useful way to model a deformable two- or three-dimensional body in an environment in which the body is subject to external forces. According to a technique known as explicit Euler integration, the forces on the mass points of the object at the beginning of a time period Δt are used to estimate the velocities and locations of the mass points at the end of the time period Δt. The explicit Euler integration step may be used to estimate the locations of the mass points of the object over successive time periods. Deformation of the object resulting from the external forces can be animated by displaying the object as represented by these updated mass points over successive time periods.

[0008] A problem with the explicit Euler technique, however, is that, in order to achieve stability, the duration of the time period over which the explicit integration is performed should be inversely proportional to the stiffness of the mass-spring system used to model the object. Since most real-world objects will require relatively stiff springs for modeling purposes, the duration of the time period over which the integration must be performed will be relatively small. The smaller the time period over which the integration is performed, the slower the resulting animation of the object. The result is that, using the explicit Euler integration technique, the majority of real-world objects cannot be animated in real-time. Moreover, all known techniques which have been attempted to overcome this problem yield results which lack realism.

[0009] For further background on the subject, the reader is referred to “Interactive animation of structured deformable objects,” Mathieu Desbrun et al., Graphics Interface '99 Proceedings, pp. 1-8, June, 1999, Kingston, Ontario; “Interactive Animation of Cloth-like Objects for Virtual Reality,” Mathieu Desbrun et al., published in “Cloth Modeling and Animation,” A. K Peters, Ltd., Publisher, Donald H. House & David E. Breen, Editors, Kathryn Maier, Associate Editor, Apr. 30, 2000; and “Interactive Animation of Cloth-like Objects in Virtual Reality,” Mark Meyer et al., submitted for publication in Journal of Visualization and Computer Graphics. All of these papers or publications are hereby incorporated by reference herein as through set forth in full.

[0010] Accordingly, there is a need for a system of and method for animating two- and three-dimensional objects which overcomes one or more disadvantages of the prior art.

[0011] There is also a need for a system of and method for animating two- and three-dimensional objects in real-time while maintaining realistic animation quality.

[0012] A further need is for a system of and method for animating two- and three-dimensional objects which fulfills any of the foregoing needs, singly or in combination.

[0013] In accordance with the purpose of the invention as broadly described herein, there is provided a method for updating the locations and velocities of the mass points of a two- or three-dimensional object to reflect passage of a time period having a duration Δt. The object is represented by a mass-spring system comprising a plurality of mass points interconnected with springs. In one embodiment, the method comprises the steps of: performing a modified implicit integration step to determine new locations of the mass points of the object after a time period Δt, optionally correcting the new locations of the mass points of the object after the time period Δt in order to preserve angular momentum, and optionally further correcting the new locations of the mass points of the object after the time period Δt to correct for high stiffness or long time step conditions.

[0014] In one embodiment, the method further comprises displaying the object as represented by the corrected locations of the mass points, and successively repeating the foregoing steps over at least one successive time period.

[0015] In one configuration, the modified implicit integration step comprises the substeps of calculating the internal force on each mass point of the object at the beginning of the time period Δt due to the action of springs; adding an artificial viscosity term to each of these internal forces; filtering out high frequencies from the resultant forces; calculating new velocities for each of the mass points of the object using the filtered forces and any external force applied to the mass points; and calculating new positions for each of the mass points of the object from the updated velocities.

[0016] In one implementation, the modified implicit integration step is preceded by a precomputation step in which an n×n matrix H is precomputed, where n is the number of mass points in the object. The matrix H is defined such that H_ij, for i≠j, is equal to k_ij, where k_ijis the spring constant for the spring between the ith and jth mass points of the object, and H_ijis equal to $-\sum _{j\ne i}{k}_{\mathrm{ij}}.$

[0017] This latter term is the negative of the sum of the spring constants for all mass points connected to mass point i with a spring.

[0018] Once the matrix H has been precomputed, the matrix W is precomputed in accordance with the following expression: $W={\left(I_n-\frac{{\mathrm{dt}}^2}{m}H\right)}^{-1}$

[0019] where I_nis the n×n identity matrix, dt is the duration time period over which the implicit integration is being performed, m is the mass of a single mass point, and H is the hessian matrix previously discussed.

[0020] The precomputation step in this implementation is next followed by a series of steps in which, for each mass point i, the internal force F_ion the mass point due to the action of springs i is first determined, and then an artificial viscosity term is computed and added to this internal force. This series of steps can be represented by the following pseudo-code:

[0021] For each mass point i

[0022] F_i=0

[0023] For each mass point j linked to mass point i with a spring $F_i=F_i+k_{\mathrm{ij}}\left(\left|x_i-x_j\right|-l_{\mathrm{ij}}^0\right)\frac{x_i-x_j}{\left|x_i-x_j\right|}$  F_i=F_i+k_ijdt(v_j−v_i)

[0024] In the foregoing, F_iis the force on mass point i due to springs and viscosity, k_ijis the spring constant for the spring between mass points i and j, x_iis the vector representing the current location of mass point i, x_jis the vector representing the current location of mass point j, l_ij⁰is the rest length of the spring between mass points i and j, dt is the duration of the time step over which the implied integration step is being performed, v_iis the current velocity of mass point i, and v_jis the current velocity of mass point j.

[0025] This series of steps is then followed in this implementation by a series of steps to filter out high frequencies from the force field determined in the previous series of steps, and to determine updated velocities and locations for each of the object mass points. This series of steps can be represented by the following pseudo-code:

[0026] For each mass point i $F_i^{\mathrm{filtered}}=\sum _j{F}_j{W}_{\mathrm{ij}}$ $v_i^{n+1}=v_i^n+\left[F_i^{\mathrm{filtered}}+F_i^{\mathrm{ext}}\right]\frac{\mathrm{dt}}{m}$  x_i^new=x_i+v_iⁿ⁺¹dt

[0027] In the foregoing, F_jis the force on mass point j computed in the previous loop of pseudo-code, W_ijis the element of matrix W located at the intersection of the ith row and jth column, F_i^extis the external force exerted on mass point i, m is the mass of a single mass point, v_iⁿis the velocity of mass point i at the inception of the time period of duration dt, x_iis the location of mass point i at the inception of the time period of duration dt, v_iⁿ⁺¹is the updated velocity of a mass point, and x_i^newis the updated location of the mass point.

[0028] The addition of an artificial viscosity term to the internal forces due to springs, filtering of the resultant forces to remove one or more high frequencies, and use of a constant matrix W to perform this high frequency filtering are all independent aspects of the invention which may be present separately or collectively in various embodiments of the invention.

[0029] Once the modified implicit integration step has been performed, the next step comprises optionally correcting the updated locations of the mass points in order to preserve angular momentum. In one implementation, this step is performed first by calculating the center of gravity of the object x_Gin accordance with the following equation: $x_G=\frac{\left(\sum _i{x}_i\right)}{n},$

[0030] in which n is the number of mass points in the object. Then, the global torque δT created by the filtered force is computed using the following equation: $\delta T=\sum _i{F}_i^{\mathrm{filtered}}\otimes x_i.$

[0031] Next, a corrective displacement for each mass point i, D_i^corr, is computed in accordance with the following equation: D_i^coor=(x_G−x_i){circle around (×)}δT. Finally, the updated location of each mass point, x_i^newis corrected in accordance with the following equation: $x_i^{\mathrm{new}}=x_i^{\mathrm{new}}+D_i^{\mathrm{corr}}\frac{{\mathrm{dt}}^2}{m}.$

[0032] Finally, an optional correction step is performed to accommodate cases involving high stiffness springs or large time periods. In one implementation, this step is performed by iterating through all or substantially all of the pairs of mass points joined by a spring, and checking to see if the spring length l_ijbetween mass points i and j exceeds the rest length, l_ij⁰, by a predetermined maximum amount such as, for example, 10%. If this condition is present, the two mass points are moved closer together in equal measures until the distance between the two equals the rest length plus the predetermined maximum amount. The only exception to this remedy is if movement of one of the mass points is constrained in some way, such as if it is being held for example. In this case, the unconstrained mass point is moved closer to the constrained mass point until the length between the two is at the predetermined maximum.

[0033] In this implementation, this process continues to iterate until one of the following three conditions are present. First, the distance between none of the mass point pairs exceeds the predetermined maximum. Second, the number of iterations equals or exceeds a predetermined maximum number of iterations. Third, a predetermined maximum time limit is equaled or exceeded.

[0034] In addition to the foregoing methods, related processor readable media, network accessible servers, and systems in accordance with the subject invention are also provided.

[0035] FIG. 1 illustrates one embodiment of a method in accordance with the subject invention.

[0036] FIG. 2 illustrates an implementation of one embodiment of a method in accordance with the subject invention.

[0037] FIGS. 3A-3E illustrate an implementation example of one embodiment of a method in accordance with the subject invention.

[0038] FIG. 4 illustrates one embodiment of a processor-based system in accordance with the subject invention.

[0039] FIG. 5 illustrates one embodiment of a client-server architecture in accordance with the subject invention.

[0040] An embodiment of a method in accordance with the subject invention is illustrated in FIG. 1. The method is for updating the locations and velocities of mass points of an object represented by a mass-spring system in which the mass points are interconnected with springs. The updating is to take account of the passage of a time period having the duration dt. One application of this method is for supporting animation of a two- or three-dimensional object in real-time. In this application, the method is advantageously performed at the inception of each of the time periods of duration Δt over which the animation is performed.

[0041] In step 1, a modified implicit integration step is performed to determine new locations of the mass points of the object at the conclusion of the time period of duration Δt. To explain the principles behind this step, the following terminology is adopted:

[0042] x is the geometric state of the system, consisting of all of the positions x_iof the mass points: x=(x₁, x₂, . . . , x_n)^T

[0043] v is the vector containing all of the velocities: v=dx/dt

[0044] F_idenotes the internal forces (due to springs and viscosity) acting on a mass point, i, whereas F_i^extwill denote external forces such as gravity

[0045] The superscript index n or lack thereof indicates the time equal to an arbitrary time t₀plus ndt. For instance, x_iⁿor x_iboth mean x_i(t₀+ndt).

[0046] The superscript index n+1 indicates a time equal to t₀+(n+1)dt

[0047] The backward difference operator Δⁿ⁺¹x refers to xⁿ⁺¹−xⁿ

[0048] dt and Δt are used synonomously

[0049] The purpose of the modified implicit integration step is to estimate, using a variant of implicit Euler integration, new locations of the mass points of the object, x_iⁿ⁺¹, from the previous locations and velocities of these mass points, x_iⁿand v_iⁿ, respectively, and also the forces, F_iⁿ⁺¹, exerted on these mass points at the end of the time period Δt. The following equations describe the implicit Euler integration step: $\begin{array}{cc}{v}_i^{n+1}=v_i^n+F_i^{n+1}\frac{\mathrm{dt}}{m}& \left(1\right)\end{array}$  x_iⁿ⁺¹=x_iⁿ+v_iⁿ⁺¹dt  (2)

[0050] There are three modifications which can be made to the implicit integration process in order to promote computational efficiency. The first two stem from using the following first-order approximation to estimate Fⁿ⁺¹: $\begin{array}{cc}{F}^{n+1}=F^n+\frac{\partial F}{\partial x}{\Delta }^{n+1}x& \left(3\right)\end{array}$

[0051] This approximation avoids the problem, present in the original implicit integration scheme, of having to compute Fⁿ⁺¹without knowing the positions xⁿ⁺¹of the mass points at time t+dt. Letting H refer to $\frac{\partial F}{\partial x},$

[0052] and plugging equation (3) into equations (1) and (2) yields: $\begin{array}{cc}{\Delta }^{n+1}v={\left(I-\frac{{\mathrm{dt}}^2}{m}H\right)}^{-1}\left(F^n+{\mathrm{dtHv}}^n\right)\frac{\mathrm{dt}}{m}& \left(4\right)\end{array}$

[0053] Comparing equation (4) with equation (1), it can be seen that the term dtHvⁿis added to Fⁿ, the force on the mass points due to the springs. For a given mass point i, this additional force can be written as $\sum _i{k}_{\mathrm{ij}}\mathrm{dt}\left(v_j-v_i\right),$

[0054] ∀ j linked to mass point i with a spring. This term represents an artificial viscosity term which is added to the internal force acting on the mass point through springs. It is proportional to the stiffness of the springs to which the mass point is affixed to neighboring mass points, and to the duration dt of the time step over which the implicit integration is to be performed. This is the first modification referred to above.

[0055] The second modification involves multiplying the combined force due to springs and artificial viscosity by the matrix: $W={\left(I-\frac{{\mathrm{dt}}^2}{m}H\right)}^{-1}.$

[0056] This represents a filtering of the force field to remove high frequency components with the type and degree of filtering depending on the stiffness of the springs and the duration of the time step. For low stiffness springs or small time steps, the transfer function of the filter approaches that of a Dirac delta function, which means that a mass point translates independently of its neighbors. For high stiffness springs or large duration time steps, the transfer function of the filter approaches that of a constant value, which means that the mass points all translate together as in rigid body motion.

[0057] The third modification involves approximating the matrix H as a constant matrix, computed assuming that the only internal force acting on the object is the linear portion of the force due to springs. Under this assumption, the matrix is defined as follows: H_ij, where i≠j, equals k_ijwhere k_ijis the spring constant for the spring between mass points i and j, and H_ijequals $-\sum _{j\ne i}{k}_{\mathrm{ij}}.$

[0058] This modification avoids the problem of having to compute the matrix H at each step Δt over which the object is animated. The matrix H would otherwise have to be computed at each step Δt over which animation is performed since it varies with the force field, and the locations of the mass points of the object. That can be seen from the formula for $\frac{\partial F_{\left(i,j\right)}}{\partial x_i}$

[0059] (where F_(ij)is the force due to springs between mass points i and j) as follows: $-k_{\mathrm{ij}}\frac{\left|x_i-x_j\right|-l_{\mathrm{ij}}^o}{\left|x_i-x_j\right|}{I}_3+l_{\mathrm{ij}}^o\frac{{\left(x_i-x_j\right)}^T\left(x_i-x_j\right)}{\left|x_i-x_j\right|}$

[0060] The internal force due to springs acting on the object can be split up into a linear force and a non-linear force, where the linear portion of the force, F_(i,j)^linearcan be represented as: −k_ij(x_i−x_j), and the non-linear portion of the force, F_(i,j)^non-linear, can be represented as: $k_{\mathrm{ij}}{l}_{\mathrm{ij}}^0\frac{\left(x_i-x_j\right)}{\left|x_i-x_j\right|}.$

[0061] The matrix H for the linear portion of the force is as defined previously: H_ij, where i≠j, equals k_ij, where k_ijis the spring constant for the spring between mass points i and j, and H_ijequals $-\sum _{j\ne i}{k}_{\mathrm{ij}}.$

[0062] As to the non-linear portion of the force, the absolute value of this force is constant, and only the phase varies over time. Consequently, this force only causes rotation of the object during the time period Δt. If this rotation is ignored, the matrix H for the non-linear force will be the null matrix since the non-linear force is assumed not to change. Given this assumption, the matrix H for the combined linear and non-linear portions of the force can be taken to be that of the linear force taken alone. Since this matrix is a constant, the computational burden of having to recompute H at each time period Δt over which the animation is performed is avoided.

[0063] These three modifications represent independent aspects of the invention. Accordingly, the subject invention includes embodiments in which one or any combination of the foregoing modifications are performed. For example, the subject invention includes a modified implicit integration step in which the artificial viscosity term is considered, the matrix H is assumed to be constant, but high frequency filtering of this combined force is either avoided entirely or implicitly performed in one of the other steps. Similarly, the invention includes an embodiment in which high frequency filtering of the internal forces is performed, but the addition of an artificial viscosity term is avoided.

[0064] Turning back to FIG. 1, in optional step 2, the new locations of the mass points are corrected in order to compensate, at least in part, for one or more of the approximations made in step 1. In one embodiment, this step comprises correcting the new locations of the mass points as computed in step one in order to preserve angular momentum, thereby achieving more realistic animation. In this embodiment, a correction to conserve linear momentum is not necessary since such is automatically preserved in the implicit integration step. However, such is not the case with angular momentum. This factor is not conserved, particularly in cases involving stiff springs, since the rotation resulting from the non-linear force is ignored in computing H. Fortunately, satisfactory results can be achieved by correcting, in a post implicit integration step, the new locations of the mass points as computed in step 1.

[0065] In one embodiment, this step involves first determining the center of gravity of the object, x_G, which can be represented as follows: $x_G=\frac{\sum _i{x}_i}{n}.$

[0066] Next, the global torque δT is computed in accordance with the following: $\delta T=\sum _i{F}_i^{\mathrm{filtered}}\otimes x_i.$

[0067] Thirdly, the new locations of the mass points of the object as determined in step 1 are corrected so that the global torque is 0.

[0068] In one implementation, this last substep involves first computing the approximate corrective displacement for each mass point i in accordance with the following: D_i^corr=(x_G−x_i){circle around (×)}δT. Next, the new locations of the object are updated as follows: $x_i^{\mathrm{new}}=x_i^{\mathrm{new}}+D_i^{\mathrm{corr}}\frac{{\mathrm{dt}}^2}{m}.$

[0069] In step 3, the new locations of the mass points are optionally further corrected to accommodate high stiffness or long time step conditions. The motivation for this step is the recognition that springs are only an imperfect model for deformable objects. With springs, it is assumed that the elongation is proportional to the force applied. That may result in implausibly large deformations for real world objects. In step 3, a form of constraint enforcement may be applied to ensure that the displacements of the mass points of the object, as determined in the previous two steps, are not too great.

[0070] A second embodiment of a method in accordance with the subject invention involves performing the foregoing method steps, but, in addition, further comprises displaying the object as represented by the updated mass points. Advantageously, the method may be repeated over successive time periods of the same or varying duration, thereby achieving animation of the object in real time.

[0071] An implementation of a method in accordance with the subject invention is illustrated in FIG. 2. The method begins with step 4, comprising the step of calculating a force field representing the sum of the internal force on each mass point of the object at the beginning of the time period of duration Δt, assuming the object is represented by a mass-spring model. In one implementation example, for each mass point i, this step can be represented as follows, where the index j represents all the mass points linked to mass point i with a spring: $F_i=\sum _j\left[k_{\mathrm{ij}}\left(\left|x_i-x_j\right|-l_{\mathrm{ij}}^0\right)\frac{\left(x_i-x_j\right)}{\left|x_i-x_j\right|}\right]$

[0072] Next, in step 5, an artificial viscosity term is added to each term of the internal force field to provide a corrected force field. In one implementation, this involves adding to the internal force applied to each mass point i the following term for each of the mass points j connected to mass point i with a spring: k_ijdt(v_j−v_i). In one implementation example, this step can be represented as follows: $F_i=F_i+\sum _j{k}_{\mathrm{ij}}dt\left(v_j-v_i\right)$

[0073] Then, in optional step 6, high frequencies are filtered out from the corrected force field determined in step 5. In one implementation, this step involves pre-computing the constant matrix $W={\left(I_n-\frac{{\mathrm{dt}}^2}{m}H\right)}^{-1},$

[0074] where H_ij, i≠j, equals k_ij, and H_iiequals $-\sum _{j\ne i}{k}_{\mathrm{ij}}.$

[0075] Next, this implementation involves computing the filtered force F_i^filteredfor each mass point i as follows: $F_i^{\mathrm{filtered}}=\sum _j{F}_i{W}_{\mathrm{ij}}$

[0076] In step 7, this implementation comprises calculating the new velocity for each mass point i. In this implementation, this step involves, for each mass point i, determining the total force on each mass point due to the filtered force field, and also external factors, such as wind, contact with another object, grasping by a human, etc. Then, this total force is used to compute the new velocity of the object in accordance with the standard implicit integration equation (1) presented earlier. This step can be represented as follows: $v_i^{n+1}=v_i^n+\left[F_i^{\mathrm{filtered}}+F_i^{\mathrm{ext}}\right]\frac{\mathrm{dt}}{m}$

[0077] An alternative implementation, applicable in the case in which filtering step 6 is avoided, involves determining the updated velocities from the unfiltered force F_ifor each mass point i as determined in step 5.

[0078] A second alternate implementation, also applicable in the case in which filtering step 6 is avoided, involves determining the updated velocities from the unfiltered forces in a manner which implicitly filters out one or more high frequencies from the unfiltered forces.

[0079] In step 8, the new position for each mass point i is determined from the new velocity of the mass point as determined in step 7. This step may be represented as follows:

[0080] x_i^new=x_iⁿ+v_iⁿ⁺¹dt

[0081] In step 48, the new position of each mass point is corrected to preserve angular momentum. In one implementation, this step includes first determining global torque δT and the center of gravity x_Gas follows: $\delta T=\sum _i\left(F_i^{\mathrm{filtered}}\otimes x_i\right)$ $x_G=\frac{\sum _i{x}_i}{n}$

[0082] Then, it includes determining the corrective displacement which needs to be applied to each of the mass points i in order to compensate for the non-zero global torque. This substep can be represented as follows:

[0083] D_i^corr=(x_G−x_i){circle around (×)}δT

[0084] Finally, this step includes the substep of correcting the new location of each mass point as determined in step 8. This substep can be represented as follows: $x_i^{\mathrm{new}}=x_i^{\mathrm{new}}+D_i^{\mathrm{corr}}\frac{{\mathrm{dt}}^2}{m}$

[0085] Next, optional step 49 is performed. As previously discussed, this optional step involves iteratively correcting the new locations of the mass points as determined in step 48 in order to more accurately model the elongation characteristics of real world objects. In one implementation, this step proceeds until (1) none of the pairs of mass points connected by springs are extended by an excessive amount; (2) a predetermined limit on the number of iterations is equaled or exceeded; or (3) a predetermined time limit for this step is equaled or exceeded.

[0086] In one implementation example, this step is performed by iterating through all or substantially all of the pairs of mass points joined by a spring, and checking to see if the spring length l_ijbetween mass points i and j exceeds the rest length, l_ij⁰, by a predetermined maximum amount such as, for example, 10%. If this condition is present, the two mass points are moved closer together in equal measures until the distance between the two equals the rest length plus the predetermined maximum amount. The only exception to this remedy is if movement of one of the mass points is constrained in some way, such as if it is being held for example. In this case, the unconstrained mass point is moved closer to the constrained mass point until the length between the two is at the predetermined maximum.

[0087] In this implementation example, this process continues to iterate until one of the following three conditions are present. First, the distance between none of the mass point pairs exceeds the predetermined maximum. Second, the number of iterations equals or exceeds a predetermined maximum number of iterations. Third, a predetermined maximum time limit is equaled or exceeded.

[0088] Further information regarding aspects of this step are available in “Deformation constraints in a mass-spring model to describe rigid cloth behavior,” X. Provot, Graphics Interface, pp. 147-154, June 1995, which is hereby fully incorporated by reference herein as though set forth in full.

[0089] An implementation example of a method in accordance with the subject invention is illustrated in FIGS. 3A-3E.

[0090] In step 9, the matrix W is precomputed in accordance with the following: $W={\left(I_n-\frac{{\mathrm{dt}}^2}{m}H\right)}^{-1}$

[0091] where I_nis the n×n identity matrix, Δt is the duration of the time period over which the implicit integration is being performed, m is the mass of a single mass point, and H is the an n×n marix in which n is the number of mass points in the matrix, and H_ij, an element of the matrix where i≠j, is equal to k_ijwhere k_ijis the spring constant for the spring between the ith and jth mass points of the object, and H_ijis equal to $-\sum _{j\ne i}{k}_{\mathrm{ij}}.$

[0092] This latter term is the negative of the sum of the spring constants for all mass points connected to mass point i with a spring. Step 9 is followed by step 10, in which an index i and a variable x_Gare initialized to 0.

[0093] Steps 11-19 are then performed. In step 11, the variable F_iis set to 0. Step 11 is followed by step 12, in which the value x_iis added to x_G, and the result then replaces x_G. Step 12 is followed by step 13 in which the index j is initialized to 0. Step 14 is then performed. In step 14, the variable F_iis updated in accordance with the following: $F_i=F_i+k_{\mathrm{ij}}\left(\left|x_i-x_j\right|-l_{\mathrm{ij}}^0\right)\frac{x_i-x_j}{\left|x_i-x_j\right|}$

[0094] In the foregoing, F_iis the internal force on mass point i due to the action of springs, k_ijis the spring constant for the spring between mass points i and j, x_iis the vector representing the location of mass point i, x_jis the vector representing the location of mass point j, l_ij⁰is the rest length of the spring between mass points i and j, and dt is the duration of the time step over which the implied integration step is being performed.

[0095] Step 14 is followed by step 15. In step 15, the variable F_iis updated with an artificial viscosity term as follows:

[0096] F_i=F_i+k_ijdt(v_j−v_i)

[0097] Step 15 is followed by step 16 in which the index j is incremented by 1.

[0098] Step 16 is followed by step 17 in which a query is made whether the index j exceeds n, the number of mass points in the object. If not, a jump is made to step 14, whereupon steps 14-17 are repeated. If so, step 18 is performed. In step 18, the index i is incremented by 1. Step 18 is followed by step 19, in which a query is made whether the index i exceeds n, the number of mass points in the object. If not, ajump is made to step 11, whereupon steps 11-19 are repeated. If so, step 20 is performed.

[0099] Steps 11-19 can be represented by the following pseudo-code:

[0100] For each mass point i

[0101] F_i=0

[0102] For each mass point j linked to mass point i with a spring $F_i=F_i+k_{\mathrm{ij}}\left(\left|x_i-x_j\right|-l_0^{\mathrm{ij}}\right)\frac{x_i-x_j}{\left|x_i-x_j\right|}$  F_i=F_i+k_ijdt(v_j−v_i)

[0103] In step 20, the variable x_Gis updated as follows:

[0104] x_G=x_G/n

[0105] At this point, x_Grepresents the center of mass of the object.

[0106] Step 20 is followed by step 21, in which the index i is re-initialized to 0, and the variable δT is initialized to 0. Step 21 is followed by step 22, in which the variable F_i^filteredis computed in accordance with the following expression: $F_i^{\mathrm{filtered}}=\sum _j{F}_j{W}_{\mathrm{ij}}$

[0107] In the foregoing, F_jis the force on mass point j computed in the previous loop of pseudo-code, W_ijis the element of matrix W located at the intersection of the ith row and jth column, and F_i^extis the external force exerted on mass point i, m is the mass of a mass point. This step represents high frequency filtering of the force field using the matrix W.

[0108] Step 22 is followed by step 23, in which the variable δT is updated in accordance with the following expression:

[0109] δT=δT+F_i^filtered{circle around (×)}x_i

[0110] In this step, δT represents global torque.

[0111] Step 23 is followed by step 24 in which the external force on a mass point, F_i^extis added to the filtered internal force on the mass point, F_i^filtered. The combined force is then used to compute the new velocity v_iⁿ⁺¹for the mass point in accordance with the following expression: $v_i^{n+1}=v_i^n+\left[F_i^{\mathrm{filtered}}+F_i^{\mathrm{ext}}\right]\frac{\mathrm{dt}}{m}$

[0112] Step 24 is followed by step 25 in which the new location of a mass point, x_i^new, is computed in accordance with the following expression:

[0113] x_i^new=x_i+v_iⁿ⁺¹dt

[0114] Step 26 is then performed in which the index i is incremented by 1. Step 26 is followed by step 27 in which the index i is compared with n, the number of mass points in the object. If the index does not exceed n, a jump is made to step 22, whereupon steps 22-27 are repeated. If the index exceeds n, step 28 is performed.

[0115] In this implementation example, steps 22-27 can be represented by the following pseudo-code:

[0116] For each mass point i $F_i^{\mathrm{filtered}}=\sum _j{F}_j{W}_{\mathrm{ij}}$ $v_i^{n+1}=v_i^n+\left[F_i^{\mathrm{filtered}}+F_i^{\mathrm{ext}}\right]\frac{\mathrm{dt}}{m}$  x_i^new=x_i+v_iⁿ⁺¹dt

[0117] In step 28, the index i is reinitialized to 0. Step 28 is followed by step 29 in which a corrective displacement D_i^corris computed for each of the mass points in accordance with the following expression:

[0118] D_i^corr=(x_G−x_i){circle around (×)}δT

[0119] Step 29 is followed by step 30 in which the new location of a mass point as determined in step 25 is updated in accordance with the corrective displacement D_i^corr. The following expression describes this process: $x_i^{\mathrm{new}}=x_i^{\mathrm{new}}+D_i^{\mathrm{corr}}\frac{{\mathrm{dt}}^2}{m}$

[0120] Step 30 is followed by step 31 in which the index i is incremented by 1. Step 32 is then performed. In step 32, a query is made whether the index i exceeds n, the number of mass points in the object. If not, a jump is made to step 29, whereupon steps 29-32 are repeated. If so, a jump is made to step 33.

[0121] Steps 28-32 can be represented by the following pseudo-code:

[0122] For each index point i

[0123] D_i^corr=(x_G−x_i){circle around (×)}δT

[0124] $x_i^{\mathrm{new}}=x_i^{\mathrm{new}}+D_i^{\mathrm{corr}}\frac{{\mathrm{dt}}^2}{m}$

[0125] Next, in step 33, an optional correction step is performed to accommodate cases involving high stiffness or large time period conditions. This optional step has been adequately described in the previous sections.

[0126] Step 34 is then performed. In step 34, the index is re-initialized to 0. Step 34 is followed by step 35, in which the variable x_iⁿ⁺¹is updated by setting it equal to x_i^new. Step 35 is followed by step 34 in which the variable v_iⁿ⁺¹is updated based on x_iⁿ⁺¹. This step can be represented by the following expression:

[0127] v_iⁿ⁺¹=(x_iⁿ⁺¹−x_iⁿ)/dt

[0128] Step 35 is followed by step 36 in which the index i is incremented by 1. Next, step 38 is performed. In step 38, a query is made whether the index i exceeds n, the number of mass points in the object. If not, a jump is made to step 35, whereupon steps 35-38 are repeated. If so, the process ends.

[0129] Steps 34-38 can be represented by the following pseudo-code:

[0130] For each mass point i

[0131] x_iⁿ⁺¹=x_i^new

[0132] v_iⁿ⁺¹=(x_iⁿ⁺¹−x_iⁿ)/dt

[0133] The entirety of the implementation example of FIGS. 3A-3E can be represented by the following pseudo-code: $\mathrm{Precompute}W={\left(I_n-\frac{{\mathrm{dt}}^2}{m}H\right)}^{-1}$

[0134] At each time step dt

[0135] //Compute internal forces F_i

[0136] //due to springs and artificial viscosity

[0137] x_G=0

[0138] For each mass point i

[0139] F_i=0

[0140] x_G=x_G+x_i

[0141] For each mass point j linked to mass point i with a spring $F_i=F_i+k_{\mathrm{ij}}\left(\left|x_i-x_j\right|-l_{\mathrm{ij}}^0\right)\frac{x_i-x_j}{\left|x_i-x_j\right|}$  F_i=F_i+k_ijdt(v_j−v_i)

[0142] x_G=x_G/n

[0143] δT=0

[0144] //High frequency filtering

[0145] For each mass point i $F_i^{\mathrm{filtered}}=\sum _j{F}_j{W}_{\mathrm{ij}}$  δT=δT+F_i^filtered{circle around (×)}x_i

[0146] $v_i^{n+1}=v_i^n+\left[F_i^{\mathrm{filtered}}+F_i^{\mathrm{ext}}\right]\frac{\mathrm{dt}}{m}$

[0147] //Post correction of angular momentum

[0148] For each mass point i

[0149] D_i^corr=(x_G−x_i){circle around (×)}δT

[0150] $x_i^{\mathrm{new}}=x_i^{\mathrm{new}}+D_i^{\mathrm{corr}}\frac{{\mathrm{dt}}^2}{m}$

[0151] //Now, do post correction for high stiffness or large time

[0152] //period cases

[0153] nbIter=0

[0154] do

[0155] Eliminate conditions in which mass points separated by springs are displaced from one another by excessive amount

[0156] until (error<ε) or (nbIter>nbIterMax) or (time is up!)

[0157] //Update new position and velocity

[0158] For each mass point i

[0159] x_iⁿ⁺¹=x_i^new

[0160] v_iⁿ⁺¹=(x_iⁿ⁺¹−x_iⁿ)/dt

[0161] FIG. 4 illustrates one embodiment of a system for animating two- or three-dimensional deformable objects in real-time in accordance with the subject invention. As illustrated, the system includes CPU 42 which can be any processor, including but not limited to a microprocessor, capable of executing software instructions stored in a memory, storage device or medium accessible or readable by the processor. CPU 42 is coupled to one or more memories, storage devices or media, including by way of illustration only and not limitation, Disk 39, RAM 40, ROM 41, and VRAM 43, through Address Bus 46 and Data Bus 45. In addition, VRAM 43 is coupled to Display 44 through Video Bus 47. In this embodiment, Disk 39 functions as archival storage, RAM (Random Access Memory) 40 functions as volatile temporary storage, ROM (Read Only Memory) 41 functions as non-volatile storage, and VRAM (video RAM) 43 functions as pixel-mapping memory for Display 44. By way of illustration only, and not limitation, Disk 39 may be a hard disk, floppy disk, or compact disk.

[0162] In one configuration, software, comprising of series of instructions embodying any of the methods of the subject invention, is stored in Disk 39. When the system is booted up, the software is loaded into RAM 40. A section of RAM 40 is also allocated to the storage of the locations of the mass points of the object to be animated. CPU 42 accesses this software as stored in RAM 40, and executes it to update the locations, also in RAM 40, of the mass points of an object at successive time intervals of duration Δt or of non-uniform duration. CPU 42 also updates the pixels in VRAM 43 at the successive time intervals in accordance with the new locations of the mass points of the object. The updated pixels are displayed on Display 44. Preferably, this process takes place in real time. The result is real-time animation of the object.

[0163] In another configuration, the software embodying a method in accordance with the subject invention is stored in ROM 41. A section of RAM 40 is allocated to the storage of the locations of the mass points of the object to be animated. CPU 42 accesses the software as stored in ROM 41, and executes it to update the locations, stored in RAM 40, of the mass points of an object at successive time intervals of duration Δt or of non-uniform duration. CPU 42 also updates the pixels in VRAM 43 at successive time intervals in accordance with the new locations of the mass points of the object. The updated pixels are displayed on Display 44. Again, this process preferably takes place in real time. Consequently, real-time animation of the object is achieved.

[0164] FIG. 5 illustrates a second embodiment of a system for animating two- or three-dimensional objects in real time in accordance with the subject invention. In this embodiment, a client-server architecture is employed in which a client 53 accesses a server 51 over a network 52, which may be proprietary network or a public network such as the Internet which utilizes a standardized communications protocol such as TCP/IP. In one embodiment, software embodying a method in accordance with the subject invention is maintained on server 51. Client 53 issues a request to animate an object in real-time, and this spawns a request to access the software maintained on server 51. In one configuration, to eliminate excessive latency, the software is then copied to client 53 whereupon it is executed. In order to customize the animation to a particular user, the client 53 may also provide the software with user-specific information about the object to be animated. In another configuration, assuming latency is not a problem, the software is remotely executed on server 51. In this configuration, the server 51 successively provides the client 53 with updated locations of the mass points of the object for display of pixels representative thereof at the client site. In a third configuration, the software is downloaded from server 51 to a second server (not shown) which is local to client 53 or at least accessible by the client 53 without undue latency. The software may then be executed by this second server. As the software executes, the second server may then successively provide the client with updated locations of the mass points of the object for display of pixels representative thereof at the client site. Hence, again, real-time animation of the object is achieved.

[0165] While embodiments, implementations, and implementation examples have been shown and described, it should be apparent that there are many more embodiments, implementations, and implementation examples that are within the scope of the subject invention. Accordingly, the scope of the subject invention is not to be restricted except in light of the appended claims and their equivalents.