Online Remeshing for Soft Tissue Simulation in Surgical Training

Virtual and Augmented Reality Supported Simulators

Online Remeshing for Soft Tissue Simulation in Surgical Training urgical procedures most often entail modifying a patient’s macroscopic tissue structure. Cutting, tearing, and carving are some of the most common tasks that surgeons perform in both open surgery and minimally invasive surgery. A surgical training system must allow the practice of such surgical tasks. While the ability to interactively simulate the accurate deformation of soft tissue is important, a system that does not include the capability to modify the simulated tissue To graphically model and has limited utility. animate the realistic behavior Researchers have done much to improve soft-tissue-deformation of deformable tissue in surgical simulations, the authors’ system modeling. Previous research projects have focused on designing adapts tetrahedra resolution by more powerful and efficient mathedynamically retessellating the matical models and designing mesh in and around the regions frameworks to optimize the model of interest. This technique representation. There have been several attempts to locally adapt a overcomes limitations of model in and around the regions of previous methods that made it interest.1,2 However, these models difficult to modify the mesh’s rely on a preprocessing phase that topology online. depends on the object’s geometry. The major drawback to this is that performing dynamic structural modifications online becomes a challenging issue. In this article, we propose a methodology for locally modifying a deformable model’s topology online.

S

Deformable models for soft tissue simulation Surgical simulation is a graphical application that requires real-time interaction with deformable 3D objects. A surgical training system must allow the trainee to interact with virtual soft tissue. Therefore, the challenging task is to develop efficient models of living tissues so that the realistic simulation of tool–tissue interaction can be performed in real time. Among the potential applications of deformable objects, human tissue modeling is one of the most demanding because of the complexity of tissue mechanics. From a purely mechanistic point of view, soft tissues exhibit complex material properties. They are nonlin-

24

November/December 2006

Céline Paloc VICOMTech Alessandro Faraci and Fernando Bello Imperial College London

ear, anisotropic, viscoelastic, and nonhomogeneous. Soft tissues deform considerably under the application of relatively small loads. In addition, it’s quite difficult to obtain material properties of tissues in vivo. In areas such as surgical planning, the main objective is to reproduce tissue properties as accurately as possible to predict a surgery’s exact outcome; the computation time is a negligible factor. In surgical simulation, however, a tradeoff has to be considered. On the one hand, computational issues have to be examined to allow physical simulation in real time; on the other hand, the model must be realistic enough to allow practitioners to acquire relevant surgical skills. The most frequently used methods to simulate tissue deformation in real time are mass-spring systems (MSSs) and finite-element models (FEMs). MSSs achieve real time but lack accuracy when precise modeling of physiological behavior is needed. An FEM is more accurate for rendering the mechanical characteristics of human organs but is more computationally expensive. Many schemes attempt to overcome these limitations, often balancing execution speed and physical accuracy (see the “Local Adaptation Approaches” sidebar). We will demonstrate how the adaptive framework presented in this article can be applied to these methods to improve their respective performance and efficiency.

Online remeshing The model presented in this article uses remeshing techniques to locally adapt a 3D mesh on the fly. Refining and simplifying a 3D mesh have been widely studied in other research areas such as computational geometry. Our method uses existing algorithms and adapts them to deformable modeling. This section gives details of these algorithms and demonstrates how they can be applied to optimize soft tissue simulation.

Tetrahedral mesh representation The method we describe here assumes a tetrahedral mesh representation of the tissue model to provide the greatest flexibility and be less restrictive than other mesh topologies. The mesh representation offers the following advantages:

Published by the IEEE Computer Society

Authorized licensed use limited to: University of Minnesota. Downloaded on September 6, 2009 at 20:11 from IEEE Xplore. Restrictions apply.

0272-1716/06/$20.00 © 2006 IEEE

Local Adaptation Approaches Both mass-spring systems (MSSs) and finite-element models (FEMs) require a discretization of the deformable object with a 2D or 3D mesh. In mesh-based representation, the accuracy of a mesh in representing a given object is related to the density (size and number) of its cells; a parameter that we call the resolution of the mesh. For MSSs as well as FEMs, a high resolution or high number of cells of small size is required to produce accurate physical behavior. However, the computational load increases with the resolution. The concept of local adaptation uses the fact that the highest possible accuracy is not always required in every part of an object. We can achieve a sufficiently high accuracy for a specific task by locally adapting a mesh’s resolution in different parts, thus reducing processing costs and memory space. Two alternative ways can be used to tackle the problem of local mesh adaptation: ■ On-the-fly refinement. The mesh is locally refined online,

whenever and wherever the application requirements concerning accuracy change. ■ Use of a multiresolution model. A comprehensive structure is built offline in a preprocessing step that organizes a collection of alternative mesh representations of a spatial object. The hierarchy can then be queried efficiently according to parameters that an application task specifies. The concept of local mesh adaptation has only recently been applied to deformable modeling. We can classify the proposed methods into two main groups: ■ surface-based deformable models refined on the fly, and ■ volumetric deformable models based on a multiresolu-

tion representation. Within the first group, the work of Hutchinson, Preston, and Hewitt presented the first attempt to use adaptive resolution in the simulation of draped cloth by locally refining an MSS in regions of high curvature.1 This idea allows the model to converge toward the static equilibrium faster by limiting the number of masses used during calculation. Later, Kähler, Haber, and Seidel increased the smoothness of irregular triangle meshes in regions of high deformation by splitting triangles using a Loop subdivision algorithm.2 The method was used for the rendering of a physics-based facial animation system. Both approaches refine the 2D mesh on the fly. In the case of volumetric deformable models, researchers have recently proposed some methods based on a multiresolution representation. A model based on a multiresolution triangulation and stored in a preprocessed directed acyclic graph was used to refine a volumetric mass-spring network near user-controlled lines.3 Debunne et al. combined a linear finite-volumebased mechanical model with a preprocessed octree representation making it adaptive in both time and space.4 The same author later developed an adaptive framework for deformable models employing an

unstructured hierarchy of tetrahedral meshes.5 Recently, Jerabkova et al. used an FEM based on cubic elements of uniform size and extended it with an octree-based adaptive multiresolution framework.6 Wu et al. proposed a scheme for mesh adaptation based on an extension of the progressive mesh concept for simulation with a nonlinear FEM.7 Similar to Debunne et al.,4 Capell et al. represented deformations of complex objects using a hierarchical basis constructed using volumetric subdivision.8 Because volumetric parameterization is difficult for complex objects, Capell et al. supported the embedding of objects in domains that are easier to parameterize. All the adaptive approaches proposed for volumetric deformable models rely on a multiresolution presentation, built offline in a preprocessing step. Such methods have several drawbacks. The range and accuracy of the stored resolutions are fixed, which limits the flexibility of the method. Moreover, large memory is required to store the multiresolution data representation. Finally, the data structure strictly depends on the object’s topology, which is supposed to be invariable, and thus prohibits structural modifications. In sum, adaptive approaches for surface-based deformable models rely on online remeshing, whereas volume-based models use a multiresolution representation. The work described in the main article proposes to locally adapt a volume-based deformable model using online remeshing.

References 1. D. Hutchinson, M. Preston, and W. Hewitt, “Adaptive Refinement for Mass-Spring Simulations,” Proc. 7th Eurographics Workshop Animation and Simulation, Springer-Verlag, 1996, pp. 31-45. 2. K. Kähler, J. Haber, and H.-P. Seidel, “Dynamic Refinement of Deformable Triangle Meshes for Rendering,” Proc. Computer Graphics Int’l (CGI), IEEE CS Press, 2001, pp. 285-290. 3. F. Ganovelli, P. Cignoni, and R. Scopigno, “Introducing Multiresolution Representation in Deformable Object Modeling,” Proc. Spring Conf. Computer Graphics, Comenius Univ., 1999, pp. 149-158. 4. G. Debunne et al., “Interactive Multiresolution Animation of Deformable Models,” Proc. Eurographics Workshop on Computer Animation and Simulation, Springer, 1999, pp. 133-144. 5. G. Debunne et al., “Dynamic Real-Time Deformations Using Space and Time Adaptive Sampling,” Proc. Computer Graphics, ACM Press, 2001, pp. 31-36. 6. L. Jerabkova et al., “A Voxel Based Multiresolution Technique for Soft Tissue Deformation,” Proc. ACM Symp. Virtual Reality Software and Technology (VRST), ACM Press, 2004, pp. 158-161. 7. X. Wu et al., “Adaptive Nonlinear Finite Elements for Deformable Body Simulation Using Dynamic Progressive Meshes,” Computer Graphics Forum, Blackwell Publishing, 2001, pp. 349-358. 8. S. Capell et al., “A Multiresolution Framework for Dynamic Deformations,” Proc. ACM Siggraph Symp. Computer Animation, ACM Press, 2002, pp. 41-47.

IEEE Computer Graphics and Applications Authorized licensed use limited to: University of Minnesota. Downloaded on September 6, 2009 at 20:11 from IEEE Xplore. Restrictions apply.

25


■ Cartesian, rectilinear, and curvilinear meshes can all ■ ■

■

■ ■

be converted into a tetrahedral mesh. A tetrahedral mesh produces a more accurate representation of a domain with relatively few elements. A tetrahedral mesh gives a triangular surface representation, which is the standard for visualization hardware. In an MSS, a tetrahedral cell has only two zero energy configurations—that is, with all their springs at their rest length. In an MSS, a tetrahedral mesh is more appropriate to preserve volume. Tetrahedral meshes are intensively used in FEMs.

Local refinement based on incremental insertion Starting with a coarse mesh, we can perform arbitrary refinement by inserting new mesh points at any desired location. Such location is not limited to, for example, the midpoints of element edges. However, the real difficulty is that the elements should be relatively round in shape, because elements with large or small angles can degrade the quality of the numerical solution. An important refinement method for angle-bounded tetrahedralizations by arbitrary point insertion is based on the Delaunay criterion.3

1 Incremental insertion of a new vertex in a tetrahedral mesh. (a) In a tetrahedron, (b) on an inner face, and (c) on an inner edge.

(b)

(c)

(a)

(b)

(c)

26

■ if v is found inside a tetrahedron, the tetrahedron

splits into four (see Figure 1a); ■ if v lies on a face, the face splits into three, thereby

splitting one tetrahedron into three or two tetrahedra into six (the face lies on the outer or inner boundary— see Figure 1b); or ■ if v lies on an edge, the edge splits in two, thereby splitting all tetrahedra around this edge into two (see Figure 1 c). The Delaunay property is restored by using face or edge flips as we describe in the next section.

Flip algorithm The well-known flip algorithm restores the Delaunay property by flipping edges and faces that are found not to be locally Delaunay. Figure 2 shows the two basic flip operations in three dimensions. In Figure 2a, a 2-3 flip transforms the two-tetrahedron configuration on the left into the three-tetrahedron configuration on the right. A 3-2 flip is the reverse transformation. A 3D flip is not possible if the containing hexahedron of the edge or face being considered for elimination is not strictly convex (see Figure 2b). A 3D, flip-based, incremental Delaunay tetrahedralization algorithm must also be able to perform the 4-4 flip demonstrated in Figure 2c. A special case, called a 2-2 flip, appears when the bottom two tetrahedra lie on an exterior boundary, and the top two tetrahedra are not present.

Quality and size control of the generated elements

(a)

2 Tetrahedral flip. (a) 2-3 or 32 flip. (b) Unflippable case. (c) 4-4 flip.

The incremental insertion algorithm is commonly used for constructing Delaunay tetrahedralizations. It operates by maintaining a Delaunay tetrahedralization into which vertices can be inserted one at a time. For each vertex v to be inserted into the existing Delaunay tetrahedralization, v is located by searching the containing tetrahedron. Depending on the location of v,

Despite the good performance of the algorithm described previously, we know that the generation of poor-quality elements (such as flat tetrahedra or slivers) can’t be avoided. However, it’s particularly important to preserve the elements’ quality, because elements with large or small angles can degrade the numerical solution. We can solve this problem by inserting an additional vertex at the circumcenter of a triangle or tetrahedron of poor quality.3 Although in practice this method has always produced a result, there is no guarantee that all the generated new elements are larger than a specific value. This factor can become critical in realtime simulation where the integration time step varies with the size of the mesh elements. Decreasing the time step can drastically slow down the simulation. A tradeoff between size and quality of the mesh elements is therefore required.

November/December 2006 Authorized licensed use limited to: University of Minnesota. Downloaded on September 6, 2009 at 20:11 from IEEE Xplore. Restrictions apply.

During the development of the work described in this article, we observed that controlling the minimum time step is possible by conditioning the smallest edge length to be set above a certain threshold. Moreover, letting Δtmin be the smallest acceptable time step if the critical time step Δtc (≥ Δtmin) for a specific mesh whose smallest edge length is lmin is known, the time step for any mesh in the same physical conditions and with no edges smaller than lmin Δtmin/Δtc will always be larger than Δtmin. We used the previous condition to guarantee real-time simulation by controlling the mesh elements’ size. It’s particularly important when the mesh size and/or topology vary over the simulation, which is the case for any adaptive models, but also when structural modifications are performed on the model. For this reason, our method rejects the insertion of vertices whose local feature size falls below lmin.

Numerical solution with time step control To simulate the physical system’s motion, the acceleration, velocity, and position of each mesh element are evolved according to the classical laws of dynamics at each iteration. The numerical solution is an important part of any real-time simulation (using an MSS or FEM), especially when the elements’ sizes might vary over the simulation. In this context, the choice of the integration method is again a tradeoff between speed and accuracy. Explicit versus implicit methods. Unlike explicit methods, the implicit methods proposed to date rely on an expensive preprocessing step. This step mainly consists of computing matrices that strictly depend on the model discretization’s topology, which is supposed to be invariable, and thus prohibits structural modifications. For this reason, we adopted an explicit integration scheme for our model. During our development, we compared various explicit methods by evaluating the tradeoff between stability and computational efficiency. Even though the higher-order methods are computationally more expensive, the fact that they might take larger time steps makes them far more efficient. They are also clearly more accurate with less risk of divergence. Adaptive time step. Typically, the actual time step used in a simulation is defined empirically in a preprocessing stage through testing. The testing method is straightforward. Given a particular model, the time step increases while the user deforms the model until the system becomes unstable. The numerical solver then slightly reduces the value to achieve robustness. This method works well in practice. However, there is no guarantee that the model will not become unstable under forces different from those used during testing. This fact becomes critical in the case where the mesh topology might vary from the initial configuration. A better solution to determine the time step is to control its value over the course of solving the system of differential equations. The idea is to increase the time step whenever possible without incurring too much error, and to reduce it when needed to avoid excessive error. The implementation of adaptive time step control

requires that the algorithm signals information about its performance and, more importantly, an estimate of its truncation error. This method’s main advantages are that it optimizes the numerical solution, ensures the system’s stability, and automatically determines the time step without the need of any preprocessing stage and/or user input. Minimum time step. Our adaptive time stepping determines the maximum time step for a given mesh to ensure stability and optimization. However, there is no limit on the minimum time step to guarantee real-time simulation. When an excessive error is encountered, the simulation can become extremely slow if the time step is reduced below some threshold. Previous work has rarely mentioned this problem, because the system is empirically set up to ensure that the constant time step is long enough. However, this problem received a great deal of interest with the introduction of topological modifications, as the system starts becoming unstable after an element had been divided into much smaller elements.4 Theoretically, smaller element sizes give rise to smaller time steps. In a simple MSS, the Shannon theorem provides a theoretical justification: Δ tc = π

m k

(1)

Equation 1 implies that increasing the stiffness and decreasing the mass decreases the critical time step Δtc. If the number of points increases in a mesh—in other words, edge length decreases—the mass of every point must be decreased to conserve the model’s total mass. The time step thus decreases as a consequence of the Shannon theorem. Moreover, in a model where the spring stiffness is assigned so that it varies inversely to the edge length, a reduction of the length directly implies a larger stiffness and, by consequence, a smaller time step. In a simple FEM, the stability criteria of CourantFriedrich-Lewy gives a theoretical limit of the time step that is proportional to the ratio of the minimum edge length lmin to maximum wave velocity in the material νmax:4 l Δ tc = min ν max We implemented the size control (described in the previous section) to overcome this issue. Any topological operations performed on the model are constrained to generate only elements whose size is larger than a minimum threshold to ensure that the time step does not fall below an unacceptable value. Our high-order integration method with time step control permits the optimization of the numerical solution while preserving its stability and guaranteeing the highest possible frame rate within predefined conditions.

Local mesh decimation A mesh region that was refined at some stage of the simulation by the insertion of new vertices might not


27


Volumetric MSS Our initial MSS consists of a model subdivided into tetrahedra. The vertices of each tetrahedron are assigned mass points and the mass points are connected by springs and dampers. Newton’s law of dynamics governs the system’s evolution as follows:

3

Local mesh decimation based on vertex removal.

Fi = miai

T (αi)

T i(αi)

4

Interpolation of vertex position. The position of a new vertex is computed by linear interpolation of its neighbors’ positions.

require the additional elements at a latter stage of the simulation. The model implementation must thus resimplify the mesh region that has previously been refined. Because the refinement method described previously inserts new vertices in the tetrahedralization, we focused on algorithms that simplify the tetrahedralization by performing the reverse operation—that is, removing these same vertices from the mesh. We used the vertex decimation algorithm proposed by Renze and Oliver.5 A general tetrahedralization algorithm attempts to fill the hole created by the removal of a vertex or a set of adjacent vertices (see Figure 3). If the algorithm is unable to tetrahedralize the hole, the candidate decimation vertex cannot be removed at the current decimation step. Instead, the vertex is added to the end of the current decimation list for reconsideration later.

Regions of adaptation Here we examine where the mesh needs to be adapted. Clearly, the mesh should be refined where the representation is too coarse to exhibit a realistic behavior. On the other hand, the mesh should be simplified where a high accuracy is no longer required. During a simulation, different regions might require a higher accuracy and should be refined to exhibit a more realistic behavior. We have identified such regions as those of high deformation, high curvature, and contact. In the context of surgical simulation, tissue tearing is an example of high deformation. The regions subject to the highest deformation should be highly refined to exhibit a realistic modification of the tissue structure. The area of an organ being modified by the environment or by its own deformation should be refined to model a plausible deformation at the flexion point. Finally, contact can occur between soft tissues, soft tissue and a surgical instrument, and soft tissue and other anatomical structures such as bones. The soft tissue model should be refined in the area of contact to model the collision accurately.

Application to mass-spring systems We applied our methodology to a traditional MSS to demonstrate how local remeshing can contribute to increase its performance.

28

(2)

where mi is the mass of each point Pi, and ai is its acceleration caused by the force Fi. Fi is the sum of all internal and external forces. The internal forces are the resultant of the tensions of the springs and dampers linking Pi to its neighbors Pj: ⎡

( ) ∑ ⎢⎢ k (l − r ) −c ⎢

Fint Pi =

ij

j

ij

ij

ij

⎣

⎤ vij⋅ rij ⎥⎥⎥⎥ rij ⎥ ⎥ ⎥⎥⎥ ⎦

⎥⎥ rij ⎥⎥⎥⎥ r ⎥ ij

(3)

where rij and vij are the relative position and relative velocity of the point Pi with respect to its neighbor Pj. The stiffness, damping, and natural length of the spring between Pi and Pj are given by kij, cij, and lij, respectively. The nature of the external force varies according to the kind of load to which the model is exposed. Omnipresent loads represent forces such as gravity and a viscous damping. Occasional loads model the interactions with the surrounding bodies: Fext (Pi) = mig + Finterest

(4)

where mi is the mass of the point Pi and g is the acceleration of gravity. To evolve the system over time, the internal and external forces that act on each mass point Pi are computed according to Equations 3 and 4 during each iteration. The new acceleration ai due to the acting forces is deduced from Equation 2. The acceleration and the velocity are then integrated through time to compute the new velocity vi and the new position ri.

Physical consistency Elsewhere, we presented a quantitative and qualitative study of the parameter assignments of a volumetric MSS to guarantee its most consistent behavior despite resolution changes.6 Based on our previous work, we assign the different parameters (mass, stiffness, and damping, respectively) of the MSS. For point mass, the mass mi of each vertex i according to the volume Vj of its adjacent tetrahedron j, with D the material density, is mi

DΣ j V j (5)

4

Then, for spring stiffness (with E as the material elastic modulus) we use kij =

EΣij Vij lij2

(6)

Finally, we use the following equation for spring damping:


cij =

(

2 kij mi + m j

)

lij

ment. We repeated each experiment three times, using (7)

In this context, it’s then possible to change the mesh density while maintaining the same global mechanical properties.

Online remeshing of volumetric MSS We applied our online remeshing to our volumetric MSS discretized at a middle-level resolution. The initial mesh is refined by inserting new vertices, as described previously. One important issue is that we do not know a priori a new vertex’ attributes, such as its velocity and its position in the reference mesh. To keep a consistent mesh, we consider the local position of the vertex relative to its neighbors. The easiest method computes the barycentric coordinates αi of the vertex relative to its neighbors. Once we know the barycentric coordinates, we can infer position p of the vertex by linear interpolation of the neighbor’s position pi. Let n be the number of neighbors:

■ a low-resolution mesh; ■ a high-resolution mesh; and ■ an adaptive mesh; dynamically refined in the vicini-

ty of the instrument and simplified in previously refined distant regions (see Figure 5). We repeated each set of experiments for several tissue models of differing complexity. Figure 6 shows the computational time during one of the experiments for 100 iterations. Other similar measurements can be found elsewhere.7

n

p=

∑α p

i i

i= 0

And n is equal to 4 if the vertex is located in a tetrahedron, 3 if located on a face, and 2 if located on an edge. As Figure 4 illustrates, the interpolation function allows the computation of the new vertices’ position in a mesh under deformation while preserving consistency. Similarly, we assume that the attribute a of a vertex is a linear combination of its neighbor’s attribute ai, then,

5 Local mesh refinement of a liver model based on a mass-spring system. 90

70

i= 0

For some attributes, the assumption of linearity might seem restrictive. In practice, however, the method is highly satisfactory as no visual or physical artifacts are observed when new vertices are inserted into the mesh. According to Equations 5, 6, and 7, the model’s physical parameters, such as the point masses and the spring values, strongly depend on the mesh’s topology. These parameters must then be updated after each refinement or simplification. The computation depends on the number of parameters to be updated and can quickly become expensive if this number is high. However, the remeshing procedures are local operations and only modify the mesh topology locally. Consequently, only the masses or springs whose neighborhood the adaptive procedure has modified need to be updated, including those which have just been inserted. This property significantly reduces the computational load due to updating the parameters.

Evaluation We integrated the resulting model within an in-house VR surgical training system. This minimally invasive system consists of a laparoscopic instrument, an endoscopic camera, a monitor, and the physical simulation module. We carried out comparative experiments consisting of probing a tissue model with a surgical instru-

60 50 40 30 20 10

(a)

0

0

20

1.2

40 60 Iterations

80

100

High resolution Low resolution Adaptive resolution

1.0 0.8 Δt

∑

αi ai

Time (ms)

n

a=

High resolution Low resolution Adaptive resolution

80

0.6 0.4 0.2

(b)

0

0

20

40 60 Iterations

80

100

6 (a) Computational time and (b) time step during probing simulations of a mass-spring system for different resolutions: low (91 tetrahedra), high (415 tetrahedra), and adaptive (91 to 415 tetrahedra).


29


For each tissue model we observed that the highresolution and the adaptive meshes exhibit similar behavior, while the low-resolution meshes suffer from a lack of accuracy in the interaction region. Moreover, due to the small number of points and springs, the lowresolution mesh is more subject to inverted tetrahedra that can lead, in extreme cases, to an inaccurate rest configuration. The high-resolution mesh tends to prohibit interactive simulation frame rates, whereas the adaptive mesh preserves a low computational time despite the fact that the hierarchy needs to be updated in the area of contact. The computational time increases with the number of mesh elements. A mesh of n mass points will require 2n integrations per time step. It also depends on the time step’s size, as for a step α times smaller than Δtmax, the numerical solver performs 2αn integrations per iteration. Moreover, the time step depends on the minimum edge length in the mesh. This property can again be observed by comparing the diagrams of the highresolution and the adaptive mesh in Figure 6, where Δt always stays over the same threshold, as lmin is set to the same value for both meshes. Therefore, the adaptive mesh has the same α value as the high-resolution mesh and the computational time difference only depends on the difference between the number of elements. For each iteration, if the number of mass points in the adaptive and the high-resolution mesh is na and nhr, respectively, the relaxation of the adaptive mesh is nhr/na faster than the relaxation of the high-resolution mesh. The system allows interacting in real time (30 Hz) with an initial mesh of about a hundred mass points by setting the elastic modulus to 10 kilopascals (kPa), the density to 103 kilograms per cubic meter (kg/m3), and the minimum edge length to a quarter of the initial minimum edge length. This performance is similar to the multiresolution model of DeBunne et al.6 However, this model suffered from troublesome popping artifacts during changes of resolution, whereas the adaptation of our model resolution is totally transparent to the user. In addition, the DeBunne et al. model relies on precomputed meshes, and the capability to perform structural modifications has not been demonstrated, unlike our model, which allows such modifications.

Application to a nonlinear finite-element model Here we describe how to apply our local remeshing methodology to a nonlinear implementation of an FEM.

over time. Let M be the mass matrix, C the damping matrix, u the vector of the displacements, f the external force vector, and F the nodal point force vector. The general Lagrangian equation of motion for the object is Mü + Cu ˙=f–F

where the dots indicate temporal derivatives. Equation 8 represents a system of differential equations of second order; the solution to the equations is obtained by standard procedures for solving differential equations. The vectors u, f, and F depend on time while the mass matrix M and the damping matrix C are independent of time. These two matrices are usually sparse and, after applying the mass lumping simplification, become diagonal. We can obtain the solution of the nonlinear dynamic response of a finite-element system employing the most common explicit time integration operator used in nonlinear dynamic analysis, the central difference operator. We consider the finite-element system’s equilibrium at time t to calculate the displacements at time Δt. The main advantage of the method is that with M and C diagonal matrices, the solution of u(t + Δt) does not involve a triangular factorization of a coefficient matrix. We can calculate the nodal force vector F in a nonlinear model in different ways, depending on the type of nonlinearity needed. In the nonlinear generalization of Hooke’s law with second Piola-Kirchkoff stress and Green-Lagrange strain tensors, we calculate the element force vector Fe at each time step as Fe = BeSeVe where Be is the kinematic matrix, Se is the second PiolaKirchkoff stress vector, and Ve is the initial volume of the tetrahedron before the simulation. We calculate the second Piola-Kirchkoff stress tensor using the GreenLagrange strain tensor through the elastic material matrix. One characteristic of the second Piola-Kirchhoff stress and the Green-Lagrange strain tensors is that they don’t change under rigid body rotations. Thus, only the actual straining of material will yield an increase in the components of the stress tensor. We express the Green-Lagrange strain tensor as t

An FEM is the canonic way to solve continuum mechanics equations providing a mathematical framework for discretizing a continuous problem. We employed this approach to approximate a dynamic nonlinear elastic model that can simulate the behavior of nonlinear, elastic, soft tissue. We considered the tetrahedron finite element, with each of its vertices corresponding to an object’s point. Applying the Lagrangian equation of motion to simulate the deformation of human organs allows for the realistic rendering of mechanic and dynamic behavior

30

∈ ij =

(

1 t uij + t u j ,i + t uk ,it uk , j 2

)

where t

Volumetric nonlinear FEM

(8)

ui, j =

∂ t ui ∂ xj

To obtain a simple and widely used elastic material description for large deformation analysis, we can generalize the linear elastic relations as follows: t S 0 ij

= t0 Cijrs0t ∈ t S 0 ij

rs

t ∈ 0 rs

where and are the components of the second Piola-Kirchhoff stress and Green-Lagrange strain tent sors, and 0 Cijrs are the components of the constant elasticity tensor.


Considering 3D stress conditions, we have t C 0 ijrs

(

= λ δ ijδ rs + μ δ ir δ js + δ isδ jr

)

where λ and μ are the Lamé constants and δij is the Kronecker delta, λ=

Ev

( 1+ v) ( 1−2v)

μ =

(

E

2 1+ v

)

The method explained here is the formulation of a geometrical implicit FEM. This is based on first building the element matrices and vectors and then constructing the overall matrices and vectors appearing in Equation 8 to find the solution for the displacement vector u. We define as a geometrical explicit FEM the approach where, after finding the expressions of the element matrices, we solve for the element displacement vector ue and then we calculate the general displacement vector u.

7

Online selection of active tetrahedra that share the longest closest edge.

Selection of active tetrahedra A geometrical explicit FEM only considers the deforming tetrahedra that the propagation of internal forces affect. Therefore, a geometrical explicit FEM is suitable for dynamic tetrahedra selection and deselection. We define an active node as a mesh point whose internal force vector ensures a configuration change of neighboring tetrahedra—that is, when its norm is greater than a certain threshold value. Also, a node becomes active when we apply an external force to it. We define a tetrahedron as active if at least one of its four nodes is active. If at a certain time step all the nodes of an active tetrahedron present small displacements, velocities, and accelerations, the tetrahedron becomes inactive and its nodes’ displacement, velocity, and acceleration values are set to zero. The system can use a smaller number of active tetrahedra by disregarding those tetrahedra that deform very little—that is, whose nodes are moving less than a given threshold and accelerations and velocities are also below certain values. Figure 7 illustrates the online selection of active tetrahedra.

Online refinement We can interactively add nodes online in the region where the surgical tool touches the object. The insertion of such nodes increases the density of the point cloud only on the mesh surface that creates inhomogeneities in distributing the mesh’s nodes. Volumetric refinement is necessary to increase the nodes’ density inside the object. This refinement maintains the mesh’s standard density throughout the object to ensure that the simulation’s behavior is mechanically uniform. Using the refinement procedure described previously, Figure 8 shows (in red) the new internal nodes that the system added in real time inside a liver’s mesh. We specifically added these nodes close to the contact region. We will add more internal nodes further away from the contact region if the deformation spreads further. During deformation, adding new nodes means that mechanical properties of nodes and tetrahedra have to be updated. Every time we insert new nodes in the mesh, the quantities of each tetrahedron—the element

8 Initial surface nodes (blue) and new internal nodes (red) inserted locally in real time during simulation. current displacement, element previous displacement, and element mass vectors—are updated. For each new tetrahedron Tsi (i = 1 … 4) we calculate the entries corresponding to the newly inserted node P for the element e displacement at the previous time step u t−Δ t , the elee ment displacement at current time step u t , and the element mass vectors Me as follows: u te−Δ t ( P) = u te( P) =

( )

u t−Δ t P 4

( )

(9)

ut P

( )

Me P =

(10)

4

( )

M P 4

(11)

where ut−Δt, ut, and M are vectors for the overall mesh. Equations 9 through 11 ensure consistency between the element and general quantities since the following equalities have to be satisfied: 4

∑ u ( P) = u ( P) Tsi t −Δ t

t −Δ t

e= 1


31

Virtual and Augmented Reality Supported Simulators 4

∑ u ( P) = u( P) Tsi t

e= 1 4

∑ M ( P) = M( P) Tsi

e= 1

Our system updates the element stiffness matrix and element damping matrix for each new tetrahedron.

Evaluation We have used a geometrical explicit nonlinear model and dynamic active tetrahedra selection of a liver in real time—with interaction through the haptic device—see Figure 9.8 The data structure of the tetrahedral mesh contained 1,019 external nodes; 557 internal nodes (1,576 total nodes in the mesh); and 5,983 tetrahedra. Elastic modulus, Poisson’s ratio, damping, density, and time step were set respectively to 4 kPa, 0.475, 200, 1,050 kg/m3, and 0.0005 second. During the simulation, our system had an average of 150 tetrahedra simultaneously active in the vicinity of the surgical tool. It performed refinement in that same active region. The frame rate per second (fps) was higher than the requested value of 30 for an acceptable visual rendering. The system maintained the haptic frame rate per second (hps) at 1,000 throughout the simulation, ensuring the required rate necessary for effective haptic feedback.

Application to surgical tasks Here we show that, in a manner similar to probing, the adaptive remeshing permits the performance of various surgical tasks with high accuracy while optimizing the computations. We demonstrate the ability to extend our method to tissue grasping and cutting.

Grasping Beyond probing, another important task that a surgeon performs during minimally invasive surgery is grasping a piece of tissue with forceps. For example, a

surgeon grasps the gall bladder prior to its removal. Any system can easily simulate this task by determining the vertex closest to the forceps and by constraining this vertex to the instrument’s extremity after the forceps are closed. When the force exerted on that vertex exceeds a threshold, the system removes the constraint. We implemented this procedure into our system and tested it on different models using an adaptive mesh. While using a low-resolution mesh without refinement, the closest vertex can be located relatively far from the forceps, which generates visual artifacts and unrealistic grasping. Therefore, refining the mesh in the forceps’ vicinity allows for accurate tissue grasping. After removing the constraint, the system resimplifies the mesh, thus accelerating the model’s relaxation to its rest position. Figure 10 shows an example of grasping the liver using our FEM implementation in which the tetrahedra are progressively selected while the tool moves along the liver’s surface—the majority of elements are not involved in the simulation. We have conducted this simulation online using a haptic device.

Cutting The adaptive on-the-fly remeshing of the underlying mesh permits dynamic structural modifications of the model undergoing deformation. This ability directly leads to a new concept for performing tasks involving structural modifications, such as cutting. Building on our adaptive model, we propose a new approach to the cutting problem using online remeshing in the cut region. We can refine and arrange the mesh along the incision by inserting new elements (vertices and external faces). We integrated and tested this approach in our MSS implementation (see Figure 11).7 Our proposed method allows the simulation to start with a coarse mesh in contrast to existing methods requiring a high-resolution mesh to simulate a cut of the same accuracy.4 Furthermore, serious drawbacks encountered in previous methods, such as the generation of small and degenerate elements, are overcome thanks to our quality and size control. We have also shown that the computational cost of our approach is suitable for real-time simulation. Those advantages contribute to a major improvement of the overall efficiency during cutting simulation.

Conclusion

fps = 38.9 hps = 999.9

9 Deformation of a 3D liver model interacting in real time with the haptic device. 32

In this article we have proposed a methodology to build a soft tissue simulation system that applies adaptive techniques to efficiently support sophisticated surgical interactions, including structural modifications such as cutting. The results obtained with this system were demonstrated in the context of the two main modeling methods: MSS and FEM. For a mesh of moderate size, the adaptive approach offers the same accuracy as a fine mesh while providing real-time interaction that would have been impossible otherwise. Moreover, unlike previous models based on multiresolution, the change of the model resolution is totally transparent to the user (without any popping artifacts) and allows structural modifications.


Previous adaptive refinement methods for deformable surfaces were based on simple techniques, such as subdivision. The first advantage of the tetrahedral Delaunay refinement proposed in this article is to allow arbitrary refinement by inserting new mesh points at any desired location not limited to, for example, the midpoints of element edges. Moreover, the main characteristic of Delaunay refinement is to preserve the mesh quality, or in this particular context, the quality of the area of interest. While the mesh quality is a primary concern in other research areas, such as finite-element analysis, it has received little interest in the field of soft tissue simulation. However, it’s well known that degenerated elements degrade the numerical solution. In this work, we paid particular attention to preserving a reasonably good quality mesh in the region of interest. Another factor degrading the numerical solution is the size of the mesh elements, proportional to the time step. While this factor is not as important in applications where computational time is not the main constraint, it can dramatically decrease the performance of interactive applications. However, similarly to mesh quality, it has received little attention. Most previous authors adapt the time step of the numerical integration in a heuristic manner. In our work, we introduce another approach that consists of adapting the time step depending on the local error while controlling the minimum size of the mesh elements so that the time step does not fall below an unacceptable value. This permits the optimization of the integration while preserving its stability in order to guarantee the highest possible frame rate within predefined conditions. Besides the adaptive refinement, the adaptive decimation procedure is based on vertex removal, which might not be the most efficient decimation method, but we chose it because it is the reverse operation of the refinement based on vertex insertion. This ensures that the complete decimation always leads to the initial constrained Delaunay tetrahedralization. Despite a significant reduction of the computational cost thanks to the adaptive framework, the mesh complexity compatible with real-time requirements is still low. Building a complex and realistic surgical scene remains challenging. However, we have demonstrated the promising potential of the adaptive approach. While hardware acceleration will keep on improving, we believe that dynamic adaptivity can form the foundation of real-time deformable modeling entailing complex interactions, making possible the development of high-fidelity surgical training systems. ■

References 1. G. Debunne et al., “Adaptive Simulation of Soft Bodies in Real-Time,” Computer Animation, May 2000, pp. 133-144. 2. L. Jerabkova et al., “A Voxel Based Multiresolution Technique for Soft Tissue Deformation,” Proc. ACM Symp. Virtual Reality Software and Technology (VRST), ACM Press, 2004, pp. 158-161. 3. J. Shewchuk, “Tetrahedral Mesh Generation by Delaunay Refinement,” Proc. Symp. Computational Geometry, ACM Press, 1998, pp. 86-95.

10 Online active tetrahedra selection taking place while the virtual tool grasps the liver using our finite-element model implementation.

11

Grasping and cutting a mass-spring system liver model based on online remeshing. 4. A. Mor, “Progressive Cutting with Minimal New Element Creation of Soft Tissue Models for Interactive Surgical Simulation,” doctoral dissertation, Robotics Inst., Carnegie Mellon Univ., 2001. 5. K. Renze and J. Oliver, “Generalized Unstructured Decimation,” IEEE Computer Graphics and Applications, vol. 16, no. 6, 1996, pp. 24-32. 6. C. Paloc et al., “Online Multiresolution Volumetric Mass Spring System for Real Time Soft Tissue Deformation,” Proc. Medical Imaging Computing and Computer-Assisted Intervention (MICCAI), Springer, 2002, pp. 219-226. 7. C. Paloc, “Adaptive Deformable Model (Allowing Topological Modifications) for Surgical Simulation,” doctoral dissertation, Dept. of Biological and Medical Systems, Imperial College London, 2003.


33


Céline Paloc is the head of the Biomedical Applications department at VICOMTech. Her main research interests include computer vision, image processing, and computer graphics applied to biomedical environments. Paloc has an MSc in biomedical engineering from the University of Technology of Compiègne, France, and a PhD in computer science from Imperial College London. Contact her at [email protected].

Fernando Bello is a lecturer in surgical graphics and computing in the Department of Biosurgery and Surgical Technology at Imperial College London. His research interests include medical visualization and analysis, medical virtual environments, and haptic interaction. Bello has a BSc (Hon) in electronic systems engineering from Monterrey Institute of Technology, Monterrey, Mexico,and a PhD in biomedical systems from Imperial College London. He is a member of the IEEE, the Web3D Medical Working Group, and the Collaboration in Radiological Interventional Virtual Environments Steering Committee. He is also a full member of the Higher Education Academy and a Business Fellow for the London Technology Network. Contact him at [email protected].

Alessandro Faraci is a technology analyst in the physical science team of the London Technology Network. His research interests are in VR, medical imaging, image processing, real-time simulation, and haptics interaction applied to medical environments. Faraci has a BSc in applied mathematics and an MSc in statistics from the University of Genoa, Italy, and a PhD in surgical imaging and computing from Imperial College London. Contact him at [email protected].

For further information on this or any other computing topic, please visit our Digital Library at http://www. computer.org/publications/dlib.

8. A. Faraci, “A Multiresolution Nonlinear Finite Element Approach to Real-Time Simulation of Soft Tissue Deformation with Haptic Feedback,” doctoral dissertation, Dept. of Surgical Oncology and Technology, Imperial College London, 2005.

Here now from the IEEE Computer Society

IEEE ReadyNotes Looking for accessible tutorials on software development, project management, and emerging technologies? Then have a look at ReadyNotes, another new product from the IEEE Computer Society.

34

ReadyNotes are guidebooks that serve as quick-start references for busy computing professionals.

Available as immediately downloadable PDFs (with a credit card purchase), ReadyNotes sell for $19 or less. www.computer.org/ReadyNotes
