- Email: [email protected]
VR models for surgical training with realistic biophysical properties
International Congress Series 1230 (2001) 1044 – 1051
VR models for surgical training with realistic biophysical properties J. Gross*, G. Buess Section for Minimally Invasive Surgery, University of Tu¨bingen, Walhoernlestr. 22, 72072 Tu¨bingen, Germany
Abstract The key issue for surgical simulations is fast modeling of deformable objects. We present a generalized particle system approach that can be tuned to match elastic properties of soft tissues based on a tetrahedral mesh. In the small deformation limit, the model behaves like a linear finite element model, but it has the advantage of being invariant under rotations. Also, we use measured viscoelastic properties to precisely model hysteresis, creep and stress relaxation in tissues. D 2001 Elsevier Science B.V. All rights reserved. Keywords: Models; Surgical training; Biophysical
1. Introduction Training of surgeons currently takes place on the job, i.e. through observation in the OR, assisting an experienced surgeon and finally taking charge of interventions. In the future, this method is expected to be complemented and gradually replaced by training in virtual environments. In minimally invasive surgery (MIS) the issue is further emphasized by the need for special abilities on the side of the surgeon, like eye –hand coordination based on an endoscopic video display, extracting 3D structure from a 2D display and manipulation of special instruments with constrained degrees of freedom. While to date, the emphasis of surgical simulators has often been on good visual quality, in the near future realistic force feedback will be equally important. This requires a higher degree of physical fidelity in the models.
*
Corresponding author. Tel.: +49-70-712981; fax: +49-70-1295569. E-mail address: [email protected] (J. Gross).
0531-5131/01/$ – see front matter D 2001 Elsevier Science B.V. All rights reserved. PII: S 0 5 3 1 - 5 1 3 1 ( 0 1 ) 0 0 1 8 5 - 6
J. Gross, G. Buess / International Congress Series 1230 (2001) 1044–1051
1045
The biggest challenge for surgical simulations is and will continue to be the efficient simulation of deformable objects in real time (above 20 frames per second, fps). A mathematically relatively simple method is the use of mass – spring networks to represent soft tissue. In this kind of model the mass of the organ is concentrated in nodes that are coupled by elastic springs. The solution of the equations of motion is relatively straightforward and changes in topology (e.g. by cutting of tissue during the simulation) are easily implemented by removing springs from the model. However, physical material properties (e.g. shear modulus, bulk modulus, Poisson’s ratio) cannot be simulated precisely because assigning the appropriate stiffnesses to the springs is not trivial. In fact, it was shown that on triangular grids mass spring systems cannot represent physical materials precisely [1]. A generalization of mass spring systems are particle systems in which other kinds of forces are introduced. For textile simulations, shear and bending forces were added to model physical properties of clothes and other soft objects [2– 5]. On rectangular structured grids, this approach can be shown to be equivalent to a finite differences discretization of the continuum mechanics equations governing deformation. Thus, physical properties can be modeled and spring stiffnesses depend on the level of discretization in a defined way [6]. However, for arbitrarily shaped objects rectangular grids are inconvenient. In this paper, we present a generalized particle system on a tetrahedral mesh, and a method to adjust the spring stiffnesses to simulate anisotropic materials. An alternative approach to particle systems is the finite element method to discretize the continuum mechanics equations. The completely linearized version using Cauchy’s strain tensor is numerically efficient [7] and can use additional tricks to reduce the computational effort [8], but is not invariant under rotations. For surgical simulations, this is a severe limitation since parts of organs may have to be lifted and rotated out of the way, which in linear FE systems causes unphysical inflation of the model. Using Greens strain tensor is a possible solution to this problem, but no simulator has been implemented using this technique because the resulting fourth-order energy equation is numerically difficult to deal with. We use a linear FE system as a reference to adjust the spring constants in our generalized particle system, as described below. Another important issue in simulation of biological tissues is their viscoelastic properties. Viscoelasticity results in mechanical damping, time-dependent forces (relaxation, creep) and hysteresis. We know of one attempt to account for hysteresis in textile simulation [2], but in general soft materials simulations have failed to treat this issue properly so far. In fact, viscoelasticity is introduced into all of the published work by the use of attenuators, a prerequisite for numerical stability of the systems, but the attenuation coefficients are usually not adjusted as to reflect measured properties of the simulated materials. In the second part of this paper, we present a method to measure viscoelastic properties of tissues, and methods to reflect these properties in the simulated objects.
2. Nonlinear properties and linearization In general, for soft, deformable materials, even relatively small forces can lead to large displacements and deformations. Thus, the usual approximation used in continuum
1046
J. Gross, G. Buess / International Congress Series 1230 (2001) 1044–1051
mechanics that all deformations are small is, strictly speaking, inappropriate. There are two possible ways of linearizing computer models. The more basic one is to use only linear equations in solving the dynamic equations of motion. The linear FE model described above is the most extreme example of this technique and leads to obviously wrong results. A more subtle example is using the complete nonlinear equations in setting up a particle system, but using only one iteration of the Newton method in solving the equations. This amounts to a local linearization with a simplified solution strategy and can equally lead to erroneous results, as was recently shown in a textile simulation experiment [9]. However, avoiding this kind of linearization by correctly solving the nonlinear equations (of either a particle system or an FE system using Greens strain tensor) does not guarantee that the material behaves exactly as its real-world counterpart would. In order for that to be the case, the correct material constants have to be employed and finite deformations have to be considered in the whole process of assembling the equations of motion. This is a difficult task for several reasons. For example, in FE models the integral of the elastic energy over each finite element is minimized to find the equilibrium position of all nodes. This integral is usually extended over the undeformed element, which is correct only for small deformations. More importantly, it is not trivial to measure higherorder elastic constants of living tissues. Since several elastic constants are involved in most deformation modes, it is difficult to separate the contributions of different combinations of these to the measured nonlinear response. Even in an isotropic material, there are two second-order elastic constants in the elastic energy equation, and three third-order constants. In less symmetric materials, the number of constants explodes quickly, increasing both the number of required measurements and the errors incurred in determining the constants. We presently avoid this problem by restricting ourselves to a model that is adjusted to physical correct material properties only for small deformations. However, we do require that the simulated material behaves correctly under rigid transformations (translation and rotation). Introducing higher-order terms only makes sense if all of them are included, that is, including higher-order elastic constants in the constitutive equations.
3. Generalized particle system The result by van Gelder [1], that triangular mass spring systems are unable to precisely reflect physical material properties, also holds for 3D systems and tetrahedral meshes. This is easy to see: a tetrahedron made up of six springs becomes lower when its base is stretched, but Poisson’s ratio in this case only depends on the geometry, not on the spring stiffnesses. Thus, we need more flexibility in the design of the particle system. The deformation in a tetrahedron is determined by 6 degrees of freedom (DOF): The four vertices together have 12 DOF, but there are 3 DOF in rigid translation and an additional 3 for rigid rotation of the whole tetrahedron, so only 6 remain for deformation. Since there are also 6 edges, the deformation can be completely described
J. Gross, G. Buess / International Congress Series 1230 (2001) 1044–1051
1047
by the lengths of the edges. All forces due to deformation can therefore be thought of as dependent on the lengths of the edges. Next, we need a local basis in order to construct force vectors. For each vertex we use the edge vectors to the other three vertices as basis vectors. This way we have 3 base vectors times 6 edge deformations = 18 stiffness constants per vertex. However, symmetry considerations relate constants for different vertices with each other, so ultimately we end up with 36 stiffness constants to adjust. The force on vertex i in the generalized particle system therefore reads as:
6 ! ! ! ! X Fi ¼ Kijk Yj ðjYk j jY0;k jÞ
ð1Þ
j;k¼1
! ! where Y are the edge vectors, Y0 are the undeformed edge vectors and Kijk is a kind of elastic matrix of the form,
0
K1jk
B B B0 B B B B0 B ¼B B B0 B B B B0 B @ 0 0
K2jk
0
k51
k21
k31
0
k22
k32
0
k52
k23
k33
0
k53
k24
k34
0
k54
k25
k35
0
k55
k26
k36
0
k56
k11
B B B k12 B B B B k13 B ¼B B B k14 B B B B k15 B @ k16
0
k31
0
0
0
k32
0
0
0
k33
0
0
0
k34
0
0
0
k35
0
0
0
k36
0
0
0
1
C C 0C C C C 0C C C; C 0C C C C 0C C A 0 k61
1
C C k62 C C C C k63 C C C; C k64 C C C C k65 C C A k66
1048
J. Gross, G. Buess / International Congress Series 1230 (2001) 1044–1051
0
K3jk
B B B k12 B B B B k13 B ¼B B B k14 B B B B k15 B @ k16 0
K4jk
k11
0
B B B0 B B B B0 B ¼B B B0 B B B B0 B @ 0
k41
0 0
0
k42
0
k23
0
k43
0
k24
0
k44
0
k25
0
k45
0
k26
0
k46
0
k21
0
k22
0 0
k41
k51
0 0
k42
k52
0 0
k43
k53
0 0
k44
k54
0 0
k45
k55
0 0
k46
k56
k61
1
C C 0C C C C 0C C C; C 0C C C C 0C C A 0 1
C C k62 C C C C k63 C C C: C k64 C C C C k65 C C A k66
ð2Þ
The columns with zeros reflect the fact that the forces on each vertex are only based on the adjacent edge vectors (but they still depend on the lengths of all edges). Now we need to assign values to the constants in K. We use a linear finite element of the same geometry and with the desired elastic constants to calculate the elastic matrix of the tetrahedron. Six displacement vectors are calculated such that each only changes the length of one of the edges. Then for each of these displacement vectors the forces on all of the vertices are calculated using the FE elastic matrix. The FE forces are projected on the appropriate edge vectors, and from comparison of the force vectors with those obtained from Eq. (1), we get linear equations for the elastic constants kjk, which are easily solved. The linear system of equations is overdetermined by a factor of 2. The fact that nevertheless a solution is found means that the method is consistent. Note that while the parameters of the particle system are adjusted with a linear FE model, the deformation behavior is still correct under rotations, while this is not the case in the FE system. This follows from the fact that we use a local basis in terms of edge vectors to calculate forces, which is by definition independent of rotations. The advantage of our approach over an FE system using Greens strain tensor is the absence of high-order polynomials in the displacements. However, as stated above, we do not take care of correct behavior for large deformations since the higher order elastic constants are unknown.
J. Gross, G. Buess / International Congress Series 1230 (2001) 1044–1051
1049
4. Viscoelasticity During surgery simulation, the most obvious effect of viscoelasticity of tissues is probably creep. For example, when a piece of tissue is pulled with a forceps for a while and then released, it will not return to its original state immediately but will slowly retract. Measurements have shown that many types of tissue in fact exhibit a very broad range of relaxation times [10]. They show the so-called constant Q (CQ) [11] behavior, i.e. the mechanical quality factor Q is roughly independent of frequency. On the other hand, using a damping term in the simulation that is proportional to the spring stiffness leads to a single relaxation time and a strongly frequency-dependent Q. Fig. 1 shows the result of a shear measurement on a bovine liver sample ex vivo. A slice of the sample was glued between two aluminum surfaces and sheared in a mechanical testing machine. The deformation was varied according to a trapezoid scheme as shown in the figure. The resulting force was fitted to the constant Q response. The figure shows that the linear CQ theory describes the behavior of the tissue very well. Unfortunately, the CQ model cannot be directly implemented in a simulation because its response depends on an integral over the deformation history of the tissue. This would be extremely costly to store and calculate at each time step. The CQ behavior can, however, be approximated by a combination of spring – dashpot pairs in a limited frequency range. For a surgical simulation, the important range of frequencies is bounded by the perceptive ability of the user; we estimate that outside a range of 0.01– 10 Hz, the user will not notice unrealistic viscoelastic response. This range can be approximated quite well using three spring – dashpot elements. However, the range could be expanded if it turned out that users are sensitive to attenuation phenomena to a tactile update frequency of up to 1 kHz [12]. Combinations of spring– dashpot pairs can be simulated using a memory variable approach [13]. Fig. 2 compares an exact CQ result with a 1D linear FE implementation using three memory parameters. A block of tissue was sheared with constant force for 1 s, then released. The figure clearly shows that a simulation with only one arbitrarily chosen
Fig. 1. Measurement of viscoelastic properties on bovine liver ex vivo. (°) Measured stress (Pa); (—) CQ fit (Pa), Q = 4.8; ( – – – ) strain ( 103).
1050
J. Gross, G. Buess / International Congress Series 1230 (2001) 1044–1051
Fig. 2. Response of a soft material to a 1-s shearing force: exact CQ result (solid line), FE simulation using three memory parameters (dotted), and simple Kelvin – Voigt solid with one relaxation time (dashed). (Shear modulus 103 Pa at 2 Hz, density 103 kg m 3, thickness 10 cm, Q = 4.8).
relaxation time misses completely both the vibratory response of the sample (due to reflections of shear waves back and forth between the faces of the sample), and the nonexponential creep behavior that is typical for CQ materials. The method was developed for linear FE models, but it can be easily applied to particle systems as well. Note also that the CQ theory is a linear theory by nature. However, using Fung’s [10] quasielastic approach, it is possible to approximately model nonlinearity and viscoelasticity simultaneously.
5. Conclusion We have presented methods to increase the physical realism of tissue models, while preserving a numerically simple approach that should be feasible even in real-time applications. A generalized particle system is used to combine the advantages of rotational invariance, the possibility to use tetrahedral meshes and the ability to reflect physical material properties. The fact that the model is not tuned to reflect nonlinear mechanical behavior currently has no drawback since these nonlinear properties are not yet known with sufficient detail. We have shown elsewhere for a textile simulation, that the high deformation response of particle systems is comparable to an FE model using Green’s strain tensor even without this tuning [6]. Viscoelastic properties of tissues have been measured in pure shear mode. The results confirm that tissues can be described as constant Q materials. Modeling of these materials is possible with high quality using a memory variable approach. Currently, these methods are implemented to be used both for textile and for surgical simulations, together with efficient numerical techniques in order to obtain fast and reliable results.
J. Gross, G. Buess / International Congress Series 1230 (2001) 1044–1051
1051
Acknowledgements This work was supported by the Deutsche Forschungsgemeinschaft (DFG) through the ElastoMedTrain grant. The authors are indebted to Olaf Etzmuß, Michael Hauth, Bernd Eberhardt and Wolfgang Strasser for fruitful cooperation and valuable discussions.
References [1] [2] [3] [4] [5] [6] [7] [8] [9] [10] [11] [12] [13]
A. Van Gelder, J. Graphics Tools 3 (2) (1998) 21 – 41. B. Eberhardt, A. Weber, Comput. Graphics 23 (1999) 599 – 606. D. Terzopoulos, K. Fleischer, Vis. Comp. 4 (1988) 306 – 331. D. Baraff, A. Witkin, SIGGRAPH 98 Conference Proceedings, ACM, New York, NY, USA, 1998, pp. 43 – 54. D.E. Breen, D.H. House, M.J. Wozny, SIGGRAPH 94 Conference Proceedings, ACM, New York, NY, USA, 1994, pp. 365 – 372. O. Etzmuß, J. Gross, Deriving a particle system from continuum mechanics for the animation of deformable objects, subm. to IEEE Trans. Visualization Comp. Graphics. H. Delingette, S. Cotin, N. Ayache, in: N. Thalmann, D. Thalmann (Eds.), Computer Animation, (Computer Animation’99), IEEE Computer Society, 1999, pp. 70 – 81. M. Bro-Nielsen, Proc. IEEE 86 (3) (1998) 490 – 503. M. Hauth, O. Etzmuß, in: A. Chalmers, T.-M. Rhyne (Eds.), A High Performance Solver for the Animation of Deformable Objects using Advanced Numerical Methods, Subm. to Eurographics 2001. Y.C. Fung, Biomechanics: Mechanical Properties of Living Tissues, Springer, New York, 1993. E. Kjartansson, J. Geophys. Res., B 84 (B9) (1979) 4737 – 4748. This value holds for the Phantom force feedback device by Sensable Technologies. I. Kay, E.S. Krebes, Geophysics 64 (1) (1999) 300 – 307.
VR models for surgical training with realistic biophysical properties J. Gross*, G. Buess Section for Minimally Invasive Surgery, University of Tu¨bingen, Walhoernlestr. 22, 72072 Tu¨bingen, Germany
Abstract The key issue for surgical simulations is fast modeling of deformable objects. We present a generalized particle system approach that can be tuned to match elastic properties of soft tissues based on a tetrahedral mesh. In the small deformation limit, the model behaves like a linear finite element model, but it has the advantage of being invariant under rotations. Also, we use measured viscoelastic properties to precisely model hysteresis, creep and stress relaxation in tissues. D 2001 Elsevier Science B.V. All rights reserved. Keywords: Models; Surgical training; Biophysical
1. Introduction Training of surgeons currently takes place on the job, i.e. through observation in the OR, assisting an experienced surgeon and finally taking charge of interventions. In the future, this method is expected to be complemented and gradually replaced by training in virtual environments. In minimally invasive surgery (MIS) the issue is further emphasized by the need for special abilities on the side of the surgeon, like eye –hand coordination based on an endoscopic video display, extracting 3D structure from a 2D display and manipulation of special instruments with constrained degrees of freedom. While to date, the emphasis of surgical simulators has often been on good visual quality, in the near future realistic force feedback will be equally important. This requires a higher degree of physical fidelity in the models.
*
Corresponding author. Tel.: +49-70-712981; fax: +49-70-1295569. E-mail address: [email protected] (J. Gross).
0531-5131/01/$ – see front matter D 2001 Elsevier Science B.V. All rights reserved. PII: S 0 5 3 1 - 5 1 3 1 ( 0 1 ) 0 0 1 8 5 - 6
J. Gross, G. Buess / International Congress Series 1230 (2001) 1044–1051
1045
The biggest challenge for surgical simulations is and will continue to be the efficient simulation of deformable objects in real time (above 20 frames per second, fps). A mathematically relatively simple method is the use of mass – spring networks to represent soft tissue. In this kind of model the mass of the organ is concentrated in nodes that are coupled by elastic springs. The solution of the equations of motion is relatively straightforward and changes in topology (e.g. by cutting of tissue during the simulation) are easily implemented by removing springs from the model. However, physical material properties (e.g. shear modulus, bulk modulus, Poisson’s ratio) cannot be simulated precisely because assigning the appropriate stiffnesses to the springs is not trivial. In fact, it was shown that on triangular grids mass spring systems cannot represent physical materials precisely [1]. A generalization of mass spring systems are particle systems in which other kinds of forces are introduced. For textile simulations, shear and bending forces were added to model physical properties of clothes and other soft objects [2– 5]. On rectangular structured grids, this approach can be shown to be equivalent to a finite differences discretization of the continuum mechanics equations governing deformation. Thus, physical properties can be modeled and spring stiffnesses depend on the level of discretization in a defined way [6]. However, for arbitrarily shaped objects rectangular grids are inconvenient. In this paper, we present a generalized particle system on a tetrahedral mesh, and a method to adjust the spring stiffnesses to simulate anisotropic materials. An alternative approach to particle systems is the finite element method to discretize the continuum mechanics equations. The completely linearized version using Cauchy’s strain tensor is numerically efficient [7] and can use additional tricks to reduce the computational effort [8], but is not invariant under rotations. For surgical simulations, this is a severe limitation since parts of organs may have to be lifted and rotated out of the way, which in linear FE systems causes unphysical inflation of the model. Using Greens strain tensor is a possible solution to this problem, but no simulator has been implemented using this technique because the resulting fourth-order energy equation is numerically difficult to deal with. We use a linear FE system as a reference to adjust the spring constants in our generalized particle system, as described below. Another important issue in simulation of biological tissues is their viscoelastic properties. Viscoelasticity results in mechanical damping, time-dependent forces (relaxation, creep) and hysteresis. We know of one attempt to account for hysteresis in textile simulation [2], but in general soft materials simulations have failed to treat this issue properly so far. In fact, viscoelasticity is introduced into all of the published work by the use of attenuators, a prerequisite for numerical stability of the systems, but the attenuation coefficients are usually not adjusted as to reflect measured properties of the simulated materials. In the second part of this paper, we present a method to measure viscoelastic properties of tissues, and methods to reflect these properties in the simulated objects.
2. Nonlinear properties and linearization In general, for soft, deformable materials, even relatively small forces can lead to large displacements and deformations. Thus, the usual approximation used in continuum
1046
J. Gross, G. Buess / International Congress Series 1230 (2001) 1044–1051
mechanics that all deformations are small is, strictly speaking, inappropriate. There are two possible ways of linearizing computer models. The more basic one is to use only linear equations in solving the dynamic equations of motion. The linear FE model described above is the most extreme example of this technique and leads to obviously wrong results. A more subtle example is using the complete nonlinear equations in setting up a particle system, but using only one iteration of the Newton method in solving the equations. This amounts to a local linearization with a simplified solution strategy and can equally lead to erroneous results, as was recently shown in a textile simulation experiment [9]. However, avoiding this kind of linearization by correctly solving the nonlinear equations (of either a particle system or an FE system using Greens strain tensor) does not guarantee that the material behaves exactly as its real-world counterpart would. In order for that to be the case, the correct material constants have to be employed and finite deformations have to be considered in the whole process of assembling the equations of motion. This is a difficult task for several reasons. For example, in FE models the integral of the elastic energy over each finite element is minimized to find the equilibrium position of all nodes. This integral is usually extended over the undeformed element, which is correct only for small deformations. More importantly, it is not trivial to measure higherorder elastic constants of living tissues. Since several elastic constants are involved in most deformation modes, it is difficult to separate the contributions of different combinations of these to the measured nonlinear response. Even in an isotropic material, there are two second-order elastic constants in the elastic energy equation, and three third-order constants. In less symmetric materials, the number of constants explodes quickly, increasing both the number of required measurements and the errors incurred in determining the constants. We presently avoid this problem by restricting ourselves to a model that is adjusted to physical correct material properties only for small deformations. However, we do require that the simulated material behaves correctly under rigid transformations (translation and rotation). Introducing higher-order terms only makes sense if all of them are included, that is, including higher-order elastic constants in the constitutive equations.
3. Generalized particle system The result by van Gelder [1], that triangular mass spring systems are unable to precisely reflect physical material properties, also holds for 3D systems and tetrahedral meshes. This is easy to see: a tetrahedron made up of six springs becomes lower when its base is stretched, but Poisson’s ratio in this case only depends on the geometry, not on the spring stiffnesses. Thus, we need more flexibility in the design of the particle system. The deformation in a tetrahedron is determined by 6 degrees of freedom (DOF): The four vertices together have 12 DOF, but there are 3 DOF in rigid translation and an additional 3 for rigid rotation of the whole tetrahedron, so only 6 remain for deformation. Since there are also 6 edges, the deformation can be completely described
J. Gross, G. Buess / International Congress Series 1230 (2001) 1044–1051
1047
by the lengths of the edges. All forces due to deformation can therefore be thought of as dependent on the lengths of the edges. Next, we need a local basis in order to construct force vectors. For each vertex we use the edge vectors to the other three vertices as basis vectors. This way we have 3 base vectors times 6 edge deformations = 18 stiffness constants per vertex. However, symmetry considerations relate constants for different vertices with each other, so ultimately we end up with 36 stiffness constants to adjust. The force on vertex i in the generalized particle system therefore reads as:
6 ! ! ! ! X Fi ¼ Kijk Yj ðjYk j jY0;k jÞ
ð1Þ
j;k¼1
! ! where Y are the edge vectors, Y0 are the undeformed edge vectors and Kijk is a kind of elastic matrix of the form,
0
K1jk
B B B0 B B B B0 B ¼B B B0 B B B B0 B @ 0 0
K2jk
0
k51
k21
k31
0
k22
k32
0
k52
k23
k33
0
k53
k24
k34
0
k54
k25
k35
0
k55
k26
k36
0
k56
k11
B B B k12 B B B B k13 B ¼B B B k14 B B B B k15 B @ k16
0
k31
0
0
0
k32
0
0
0
k33
0
0
0
k34
0
0
0
k35
0
0
0
k36
0
0
0
1
C C 0C C C C 0C C C; C 0C C C C 0C C A 0 k61
1
C C k62 C C C C k63 C C C; C k64 C C C C k65 C C A k66
1048
J. Gross, G. Buess / International Congress Series 1230 (2001) 1044–1051
0
K3jk
B B B k12 B B B B k13 B ¼B B B k14 B B B B k15 B @ k16 0
K4jk
k11
0
B B B0 B B B B0 B ¼B B B0 B B B B0 B @ 0
k41
0 0
0
k42
0
k23
0
k43
0
k24
0
k44
0
k25
0
k45
0
k26
0
k46
0
k21
0
k22
0 0
k41
k51
0 0
k42
k52
0 0
k43
k53
0 0
k44
k54
0 0
k45
k55
0 0
k46
k56
k61
1
C C 0C C C C 0C C C; C 0C C C C 0C C A 0 1
C C k62 C C C C k63 C C C: C k64 C C C C k65 C C A k66
ð2Þ
The columns with zeros reflect the fact that the forces on each vertex are only based on the adjacent edge vectors (but they still depend on the lengths of all edges). Now we need to assign values to the constants in K. We use a linear finite element of the same geometry and with the desired elastic constants to calculate the elastic matrix of the tetrahedron. Six displacement vectors are calculated such that each only changes the length of one of the edges. Then for each of these displacement vectors the forces on all of the vertices are calculated using the FE elastic matrix. The FE forces are projected on the appropriate edge vectors, and from comparison of the force vectors with those obtained from Eq. (1), we get linear equations for the elastic constants kjk, which are easily solved. The linear system of equations is overdetermined by a factor of 2. The fact that nevertheless a solution is found means that the method is consistent. Note that while the parameters of the particle system are adjusted with a linear FE model, the deformation behavior is still correct under rotations, while this is not the case in the FE system. This follows from the fact that we use a local basis in terms of edge vectors to calculate forces, which is by definition independent of rotations. The advantage of our approach over an FE system using Greens strain tensor is the absence of high-order polynomials in the displacements. However, as stated above, we do not take care of correct behavior for large deformations since the higher order elastic constants are unknown.
J. Gross, G. Buess / International Congress Series 1230 (2001) 1044–1051
1049
4. Viscoelasticity During surgery simulation, the most obvious effect of viscoelasticity of tissues is probably creep. For example, when a piece of tissue is pulled with a forceps for a while and then released, it will not return to its original state immediately but will slowly retract. Measurements have shown that many types of tissue in fact exhibit a very broad range of relaxation times [10]. They show the so-called constant Q (CQ) [11] behavior, i.e. the mechanical quality factor Q is roughly independent of frequency. On the other hand, using a damping term in the simulation that is proportional to the spring stiffness leads to a single relaxation time and a strongly frequency-dependent Q. Fig. 1 shows the result of a shear measurement on a bovine liver sample ex vivo. A slice of the sample was glued between two aluminum surfaces and sheared in a mechanical testing machine. The deformation was varied according to a trapezoid scheme as shown in the figure. The resulting force was fitted to the constant Q response. The figure shows that the linear CQ theory describes the behavior of the tissue very well. Unfortunately, the CQ model cannot be directly implemented in a simulation because its response depends on an integral over the deformation history of the tissue. This would be extremely costly to store and calculate at each time step. The CQ behavior can, however, be approximated by a combination of spring – dashpot pairs in a limited frequency range. For a surgical simulation, the important range of frequencies is bounded by the perceptive ability of the user; we estimate that outside a range of 0.01– 10 Hz, the user will not notice unrealistic viscoelastic response. This range can be approximated quite well using three spring – dashpot elements. However, the range could be expanded if it turned out that users are sensitive to attenuation phenomena to a tactile update frequency of up to 1 kHz [12]. Combinations of spring– dashpot pairs can be simulated using a memory variable approach [13]. Fig. 2 compares an exact CQ result with a 1D linear FE implementation using three memory parameters. A block of tissue was sheared with constant force for 1 s, then released. The figure clearly shows that a simulation with only one arbitrarily chosen
Fig. 1. Measurement of viscoelastic properties on bovine liver ex vivo. (°) Measured stress (Pa); (—) CQ fit (Pa), Q = 4.8; ( – – – ) strain ( 103).
1050
J. Gross, G. Buess / International Congress Series 1230 (2001) 1044–1051
Fig. 2. Response of a soft material to a 1-s shearing force: exact CQ result (solid line), FE simulation using three memory parameters (dotted), and simple Kelvin – Voigt solid with one relaxation time (dashed). (Shear modulus 103 Pa at 2 Hz, density 103 kg m 3, thickness 10 cm, Q = 4.8).
relaxation time misses completely both the vibratory response of the sample (due to reflections of shear waves back and forth between the faces of the sample), and the nonexponential creep behavior that is typical for CQ materials. The method was developed for linear FE models, but it can be easily applied to particle systems as well. Note also that the CQ theory is a linear theory by nature. However, using Fung’s [10] quasielastic approach, it is possible to approximately model nonlinearity and viscoelasticity simultaneously.
5. Conclusion We have presented methods to increase the physical realism of tissue models, while preserving a numerically simple approach that should be feasible even in real-time applications. A generalized particle system is used to combine the advantages of rotational invariance, the possibility to use tetrahedral meshes and the ability to reflect physical material properties. The fact that the model is not tuned to reflect nonlinear mechanical behavior currently has no drawback since these nonlinear properties are not yet known with sufficient detail. We have shown elsewhere for a textile simulation, that the high deformation response of particle systems is comparable to an FE model using Green’s strain tensor even without this tuning [6]. Viscoelastic properties of tissues have been measured in pure shear mode. The results confirm that tissues can be described as constant Q materials. Modeling of these materials is possible with high quality using a memory variable approach. Currently, these methods are implemented to be used both for textile and for surgical simulations, together with efficient numerical techniques in order to obtain fast and reliable results.
J. Gross, G. Buess / International Congress Series 1230 (2001) 1044–1051
1051
Acknowledgements This work was supported by the Deutsche Forschungsgemeinschaft (DFG) through the ElastoMedTrain grant. The authors are indebted to Olaf Etzmuß, Michael Hauth, Bernd Eberhardt and Wolfgang Strasser for fruitful cooperation and valuable discussions.
References [1] [2] [3] [4] [5] [6] [7] [8] [9] [10] [11] [12] [13]
A. Van Gelder, J. Graphics Tools 3 (2) (1998) 21 – 41. B. Eberhardt, A. Weber, Comput. Graphics 23 (1999) 599 – 606. D. Terzopoulos, K. Fleischer, Vis. Comp. 4 (1988) 306 – 331. D. Baraff, A. Witkin, SIGGRAPH 98 Conference Proceedings, ACM, New York, NY, USA, 1998, pp. 43 – 54. D.E. Breen, D.H. House, M.J. Wozny, SIGGRAPH 94 Conference Proceedings, ACM, New York, NY, USA, 1994, pp. 365 – 372. O. Etzmuß, J. Gross, Deriving a particle system from continuum mechanics for the animation of deformable objects, subm. to IEEE Trans. Visualization Comp. Graphics. H. Delingette, S. Cotin, N. Ayache, in: N. Thalmann, D. Thalmann (Eds.), Computer Animation, (Computer Animation’99), IEEE Computer Society, 1999, pp. 70 – 81. M. Bro-Nielsen, Proc. IEEE 86 (3) (1998) 490 – 503. M. Hauth, O. Etzmuß, in: A. Chalmers, T.-M. Rhyne (Eds.), A High Performance Solver for the Animation of Deformable Objects using Advanced Numerical Methods, Subm. to Eurographics 2001. Y.C. Fung, Biomechanics: Mechanical Properties of Living Tissues, Springer, New York, 1993. E. Kjartansson, J. Geophys. Res., B 84 (B9) (1979) 4737 – 4748. This value holds for the Phantom force feedback device by Sensable Technologies. I. Kay, E.S. Krebes, Geophysics 64 (1) (1999) 300 – 307.