Modelling of soft tissue deformation for laparoscopic surgery simulation

Medical Image Analysis 4 (2000) 57–66 www.elsevier.com / locate / media

Modelling of soft tissue deformation for laparoscopic surgery simulation *, Ch. Brechbuhler, ´ ¨ G. Szekely R. Hutter, A. Rhomberg, N. Ironmonger, P. Schmid ¨ , Switzerland ETH Zentrum, Swiss Federal Institute of Technology, Gloriastr. 35, CH-8092 Zurich Received 22 September 1998; received in revised form 12 July 1999; accepted 22 July 1999

Abstract Virtual reality based surgical simulator systems offer a very elegant solution to the development of endoscopic surgical trainers. While the graphical performance of commercial systems already makes PC-based simulators viable, the real-time simulation of soft tissue deformation is still the major obstacle in developing simulators for soft-tissue surgery. The goal of the present work is to develop a framework for the full-scale, real-time, finite element simulation of elastic tissue deformation in complex systems such as the human abdomen. The key for such a development is the proper formulation of the model, the development of scalable parallel solution algorithms, and special-purpose parallel hardware. The developed techniques will be used for the implementation of a gynecological laparoscopic VR-trainer system.  2000 Elsevier Science B.V. All rights reserved. Keywords: Soft tissue deformation; Surgery simulation; Minimal invasive surgery

1. Introduction Endoscopic operations have recently become a very popular technique for the diagnosis as well as treatment of many kinds of human diseases and injuries. The basic idea of endoscopic surgery is to minimize damage to the surrounding healthy tissue, normally caused in reaching the point of surgical intervention for the more inaccessible internal organs. The relatively large cuts in open surgery can be replaced by small perforation holes, serving as entry points for optical and surgical instruments. The small spatial extent of the tissue injury and the careful selection of the entry points result in a major gain in patient recovery after operation. The price for these advantages is paid by the surgeon, who loses direct contact with the operation site. The necessary visual information is mediated by a specialized camera (the endoscope) and is presented on a screen. While *Corresponding author. Tel.: 141-1-632-52-88; fax: 141-1-632-1199. ´ E-mail address: [email protected] (G. Szekely)

preliminary systems experimenting with stereo optics are already available, today’s surgery is usually performed under monoscopic conditions. Due to geometrical constraints posed by the external control of the surgical instruments through the trocar hull, the surgeon loses much of the manipulative freedom usually available in open surgery. Performing operations under these conditions demands very special skills of the surgeon, which can only be gained with extensive training. The basic optical and manipulative skills can be learned today by using inexpensive, traditional training devices. These training units allow one to learn navigating under monoscopic visual feedback, as well as to acquire basic manipulative skills. In this way the surgeon becomes accustomed to completing a particular task, but because the real-life effect is lost, gets only a limited training range in dexterity and problem solving. Additionally, the organs used by these units are generally made of foam or other plastic material, hence there are strong limits on the level of realism which can be achieved. And while experiments on animals are sometimes used for testing new surgical techniques, practical as

1361-8415 / 00 / $ – see front matter  2000 Elsevier Science B.V. All rights reserved. PII: S1361-8415( 00 )00002-5

58

´ et al. / Medical Image Analysis 4 (2000) 57 – 66 G. Szekely

well as ethical reasons strongly restrict their use in everyday surgical training. Virtual reality based surgical simulator systems offer a very elegant solution to this training problem. A wide range of VR simulator systems have been proposed and implemented in the past few years. Some of them are restricted to purely diagnostic endoscopic investigations (Satava, 1996; Vining, 1996; Alyassin and Lorensen, 1998), while others, for example, allow the training of surgical procedures for laparoscopic (Cover et al., 1993; ¨ Kuhnapfel et al., 1995; Cotin et al., 1996; Baur et al., 1998; Downes et al., 1998), arthroscopic (Ziegler et al., 1995; Gibson et al., 1997), or radiological (Hahn et al., 1998) interventions. While the graphical performance of commercial systems already renders PC-based simulators possible (Alyassin and Lorensen, 1998), the real-time simulation of soft tissue deformation is still the major obstacle while developing simulator systems for soft tissue surgery. Different methods in use for deformation modelling include: • Free-form deformation techniques of computer graphics (Barr, 1984; Sederberg and Parry, 1986) use parametric interpolative models (such as polynomial models, splines, or superquadrics) for deformation estimation of solid primitives. While analogy to physical deformation processes is not always obvious, such techniques have become very popular in surgical simulators (Basdogan et al., 1998; Baur et al., 1998) due to the resulting fast deformation calculation. • Different more or less justified simple physical analogies have also been used for tissue deformation modelling. Most popular are mass-spring models ¨ (Kuhnapfel et al., 1995; Bro-Nielsen et al., 1998; Downes et al., 1998), but other alternatives like spacefilling spheres have also been implemented (Suzuki et al., 1998). • Elastically deformable surface models introduced in computer graphics and computer vision (Terzopoulos et al., 1987) calculate surface deformations by solving the linear elasticity equations using different numerical techniques. These methods allow simulation of tissue deformation based on physical principles (Cover et al., 1993). Full 3D extensions of these techniques (BroNielsen and Cotin, 1996; Cotin et al., 1996) already represent the first attempts for finite element-based tissue deformation modelling. The finite element method (FEM) is a very common and accurate way to solve continuum-mechanical boundaryvalue problems (Zienkiewicz and Taylor, 1994; Bathe, 1996). In the case of biological tissue, we have to deal with large deformations and also with anisotropic, inhomogeneous, and nonlinear materials. Furthermore, organs and surgical instruments interact, therein leading to numerical contact problems. Nonetheless, provided an adequate formulation is chosen, even in these cases the FEM is a very powerful tool.

Unfortunately, finite element calculations are notoriously slow, making them not very appealing for real-time applications like endoscopic surgery simulations. Accordingly, tissue deformation for surgical training systems is, up to now, only calculated by the FEM for fairly simple physical systems (Sagar et al., 1994; Bro-Nielsen and Cotin, 1996; Cotin et al., 1996). The goal of the presented work is to develop a framework for the full-scale real-time finite element simulation of elastic tissue deformation in complex systems such as the human abdomen. The key for such a development is the proper formulation of the model (see Section 2), the development of scalable parallel solution algorithms (see Section 4) as well as dedicated parallel hardware (see Section 3). A first modular prototype system implementing the basic algorithms of the developed techniques for VRtraining in gynecological laparoscopy is also presented (see Section 5).

2. Elastomechanical modelling To enable the virtual reality simulation of surgical operations, we must achieve the following: • calculation of realistic deformations and contact forces of the organs; • performance of the calculations in real-time; • stable calculations for the time range of the entire simulation.

2.1. Explicit finite element formulation A continuum can be described mathematically with the following set of partial differential equations, satisfied for each material point: div s 1 f 5 r u¨

(momentum equations)

(1)

div( r u~ ) 1 r~ 5 0 (continuity equation)

(2)

s 5 f1 (´) (constitutive law)

(3)

´ 5 f2 (u) (strain formulation)

(4)

where s is the Cauchy stress tensor, ´ is the strain tensor, f is the vector of the volume forces, u stands for the displacements, and r is the density. A superposed dot denotes a derivative with respect to time. Within the FEM a body is subdivided by a finite number of well-defined elements (e.g., hexahedrons, tetrahedrons, quadrilaterals). Displacements and positions in the element are interpolated from discrete nodal values. A trilinearly interpolated hexahedron consists of eight nodes lying in the corners. Fig. 5(left) shows a 2D set of elements. Neighbor elements share some nodes. For every element

´ et al. / Medical Image Analysis 4 (2000) 57 – 66 G. Szekely

the equations (1)–(4) can be formulated resulting in the following discrete system of differential equations: Mu¨ 1 Cu~ 1 Kd u 5 f 2 r ,

(5)

where M is the mass matrix, C is the damping matrix, K is the incremental stiffness matrix, u is the vector of the nodal displacements, f are the external node forces, and r are the internal node forces. All these matrices and vectors may be time dependent. One possibility to solve equations in (5) is a quasi-static manner (Crisfield, 1991; Hinton, 1992). In this case the ¨ ~ dynamic part of the equations is neglected (u5u50) and the solution is reached iteratively. Within every iteration a huge set of linear algebraic equations has to be solved. In problems with large deformation and contact interactions the iteration hardly converges and a stable simulation over several minutes is nearly impossible. Otherwise, the dynamic equations have to be integrated with respect to time (Hinton, 1992). The time integration of the equations can, in principle, be performed using implicit or explicit integration schemes. Using the implicit method, the solution has to be calculated iteratively at every discrete time step, like in the quasi-static problem. Contrary, the explicit time integration can be performed without iteration and without solving a system of linear algebraic equations. This integration scheme is only conditionally stable; that is, only very small time steps lead to a stable solution (Flanagan and Belytschko, 1984). An estimation for the critical time step is made by Dt 5 DL /c, where c is the maximal wave propagation speed in the medium and DL is the smallest element length of the model. In the case of our model of a uterus, this equation leads to 10 000 time steps per second (Dt5100 ms). These short time steps increase the computational effort but lead to a very stable contact formulation. The disadvantage of this method — that many steps have to be taken — is not so significant since each step is much less time consuming than in an implicit algorithm. Because of the considerations mentioned above, and based on several numerical tests, we decided to solve the problems with an explicit finite element formulation. Unfortunately, this integration scheme is not as accurate as the implicit one and time discretization errors will usually accumulate. In the next section we will show a way of avoiding this problem.

2.2. Element formulation /constitutive equations The most time consuming part in the explicit formulation is the computation of the internal element forces (see Section 2.3), including the calculation of stresses and strains. These state variables are related to well-defined reference configurations. In the following we have to distinguish the updated and the total Lagrange formulation. Within the updated Lagrange formulation a state is

59

related to the last successfully calculated configuration. Here, strains are usually of an incremental form. This leads to an accumulation of discretization errors and in consequence to the above-mentioned lack of accuracy. Additionally, an incremental strain formulation (even in the case of elastic materials) usually leads to remaining deformations in a stress-free state after a load cycle (Kojic and Bathe, 1987; Roy et al., 1992). These errors drastically influence the robustness of a simulation and have to be eliminated. When a total Lagrange formulation is chosen every state is related to the initial configuration. In this case absolute strain formulations have to be considered (e.g., GreenLagrange or Hencky). Even with the explicit time integration scheme this leads to an exact static solution, i.e., when the transient terms have been eliminated by damping effects. Contrary to incremental strain formulations, these strains lead to correct results after a load cycle and no error accumulation occurs. These considerations suggested to use a hyperelastic material law. The Mooney–Rivlin law for nearly incompressible materials leads to an elegant finite element formulation and is frequently used in biological tissue modelling. The Mooney–Rivlin material law is a first attempt, and will probably be replaced by a neo-Hookean law, completed with an appropriate damping model.

2.3. Volume integration To obtain the internal forces of an element we have to solve the following integral:

EF B S

f Ii 5

il

I k kl

dV 0 ,

(6)

V0

where repeated subscripts imply summation over the range of that subscript. f Ii are the internal forces in direction i of node I of the element, Fil are the components of the deformation gradient, B Ik are the derivatives of the interpolation functions with respect to coordinate k, Skl are the components of the second Piola–Kirchhoff stress tensor and V 0 is the initial volume of the element. Commonly, this integral is evaluated by a numerical eight-point quadrature. That is, the integrand will be computed with respect to eight different integration points, and the eight values are added with a well-defined weighting factor. This computation is very time consuming. The reduced volume integration is a method to circumvent this overhead (Belytschko and Ong, 1984). Within this scheme the integrand has to be calculated just once for each element. Considering only the mean strains within an element, some deformation modes (hourglass modes) are not influenced by the resulting forces. These modes have to be controlled with additional stabilization forces (hourglass control): f Ii 5 Fil B kI SklV 0 1 stab f iI .

(7)

60

´ et al. / Medical Image Analysis 4 (2000) 57 – 66 G. Szekely

2.4. Simulation experiments The developed methods have been extensively tested by off-line simulations of deformations using a first model of the uterus. Fig. 1 shows the generated FEM mesh, containing about 2000 elements. The abdominal cavity has been modeled as a rigid surface. Fig. 2 illustrates the initial state of the uterus within the abdominal cavity (left) compared to an inter-operative image of the operation site (right). An example for the deformation of the corpus (left) as well as for the manipulation of the fallopian tube (right) of the uterus is shown in Fig. 3.

3. Design of a real-time FEM computation engine

Fig. 1. Finite element model of the uterus containing about 2000 elements.

A shortcoming of commonly used formulations for the stabilization forces is their incremental character. Obviously, this leads to accumulated time integration errors and to incorrect results after load cycles. Consequently we formulate these forces with respect to the initial configuration. This method leads to very stable results even in cases of large deformations. Furthermore, no integration errors will be accumulated. If finite elements are not distorted in the initial configuration, absolute formulations have an additional advantage. Many multiplications can be saved, and even the hourglass control needs considerably fewer multiplications.

The only way to provide the necessary computational power for the real-time solution of complex finite element systems is to build a parallel computer which supports fully parallel algorithms for the explicit time integration scheme. Algorithmic design as described in Section 4 allows to scale the computation to the necessary speed with the selection of an appropriate number of processor units or processing elements (PE). The section summarizes the basic requirements and design principles for implementing special-purpose, parallel hardware to perform explicit finite element calculations.

3.1. Performance requirements 3.1.1. Computation Once the necessary FE calculations and the corresponding parallel algorithms are determined, we have to

Fig. 2. Initial state of the model uterus (left) compared to an inter-operative image of the operation site (right).

´ et al. / Medical Image Analysis 4 (2000) 57 – 66 G. Szekely

61

Fig. 3. Simulated deformation of the corpus (left) and the fallopian tube (right) of the uterus.

analyse the performance requirements in detail to find out what kind of computer is needed to satisfy our demands. Analysis of optimized explicit finite element algorithms shows that approximately 700 floating point operations per element are needed in each time step. Additional computation time has to be reserved for collision detection and handling. This leads to 10 MFLOPS (million floating point operations per second) per element for the chosen time step of 100 ms. For 2000 elements, a total of 20 GFLOPS sustained are needed. This amount of computational power is much too high for a state of the art workstation. Implicit FE calculation can be implemented well on vector-supercomputers, since the main task is to solve huge sparse systems of linear equations (Taylor, 1991; Pommerell, 1992). However, the time steps and the vectors of the explicit method are too short, therefore a high-performance parallel machine is needed. The latest processors are able to deliver . 300 MFLOPS with optimized programming, so 64 of these will meet our demands. As finite element calculations have a tendency to become unstable, accurate computation using double precision floating point operations is necessary. During the time steps, data have to be sent and received with minimal processor involvement, requiring that very light weight protocols be used.

3.1.2. Communication There are three different types of communication: exchange of forces, collision detection and handling, and data transfer to the graphics engine. For a FE model with 2000 elements on a 4 3 4 3 4 processor machine, 5600 force vectors have to be exchanged each time step. This results in a communication bandwidth of 1.35 GByte / s.

This is a huge number, but the fact that data only have to be transferred between neighbors eases our requirements. In collision communication, we note that only part of the surface nodes are potentially colliding with other parts of the FE model, and only about 300 vectors have to be exchanged for a 2000-element model (72 MByte / s). Finally, the information for all the surface nodes must be sent to the graphics engine. This takes place once every 40 ms, but to avoid adding too much latency to the system, data should be sent in a shorter time (20 time steps). On this basis, 18 MByte / s have to be sent to the graphics engine interface.

3.2. Parallel processing architecture According to the previous algorithmic considerations, a three-dimensional mesh of processors can be optimally used, where every PE is connected to its six neighbors (Fig. 4). Every PE calculates a cube of the whole FE model, so it has to exchange data from the faces of its cube with PEs that share the same faces. In addition there are also edges and corners, where data have to be shared among processors that are not connected by a direct channel. Such data have to be routed through other PEs. This is actually no drawback, since in most cases the intermediate PEs have also data to transmit to the same target PE, so they just add their values before transmitting the received data block.

3.2.1. Local communication To achieve the necessary low latency, we have to exploit the specialties of our communication pattern. If data have to be exchanged between processors, it must be communicated in every time step. Fortunately it is known at

62

´ et al. / Medical Image Analysis 4 (2000) 57 – 66 G. Szekely

although data for collision detection have to be exchanged between two processors that may be far away within the mesh. Part of these data values has to be broadcast to all other processors. Additionally, several processors can send and receive data at the same time. Routing such data through the local channels is very inefficient. The routing may consume processor time which can be better utilized for FE calculation and routing through several PEs increases latency to an unacceptable level. Therefore, such messages are better sent over a global network, over which any PE can exchange data with all others. Data must be identified, thus packets of data must be sent with information concerning the receiver as well as an identification.

Fig. 4. Three-dimensional communication architecture.

programming time which data values have to be sent, in what order, and to which destination. Secondly, when the elements are assigned to the processors, only those processors with elements sharing nodes on the common face of the cubes have to exchange data (collision detection is an exception to this). Bearing this in mind, the best way to provide sufficient communication bandwidth with low latency and minimal software overhead is to use point to point communication and to transmit pure data without additional address information. There is no need for routing information when point to point channels are used, and sending only data values minimizes software overhead on the sender side. But the data still have to be recognized by the receiver. The address information of the data, however, is implicitly given by the order within the time step it is sent. It is known at compile time which values have to be exchanged, so the tool that partitions the elements also sets up the order in which data values are sent and received.

3.2.2. Global communication Communication takes place mostly between neighbors,

4. Partitioning the model for parallel computation Explicit FE computation suggests element-wise parallelism, decomposing the spatial domain. Each element is mapped to exactly one processor, and one processor computes the internal forces of several elements (see Fig. 5, left). Whenever elements residing on different processors share common nodes, these nodes must be represented on all involved processors. The resulting force acting on such a distributed node emerges from the communication among the contributing processors, which treat the foreign forces like external forces, adding them to their own.

4.1. Parallel computation Due to the inherently 3D nature of the problem, threedimensional processor networks appear to be the optimal topology for a multiprocessor architecture. To organize the processors and their communication, they are labeled with ‘coordinates’, which generally do not correspond to geometric coordinates (see Fig. 5, middle). In the following, directions like x, y and z always mean communication directions.

Fig. 5. Left: a small 2D model serves for illustrating the concepts on paper. The six different shades of gray illustrate an arbitrary assignment of elements to 2 3 3 processors. Middle: the coordinates of the processors to which parts of the model (groups of elements) are assigned. Dashed lines show x borders, dotted lines y borders. Right: the model is partitioned in sub-models. Grey shades indicate batches for recovery, arrows mark nodes on the border.

´ et al. / Medical Image Analysis 4 (2000) 57 – 66 G. Szekely

In an array of n x 3 n y 3 n z processors, one processor is identified by its coordinates i x , i y , i z , where 0 < i j , n j , ; j [ hx,y,zj. This processor has the running number i x 1 n x (i y 1 n y i z ), corresponding to the conventional layout of a 3D array proc[n z ][n y ][n z ] in C. It has up to six local communication channels, namely 2 in x direction to the neighbor processors (i x 2 1,i y ,i z ) and (i x 1 1,i y ,i z ) and similarly 2 each in y and z. Processors on the border of the block lack some neighbors.

4.2. Distributed nodes The mapping of elements to processors should balance the computational load while minimizing the demand for communication bandwidth. This could be solved as a discrete optimization problem by minimizing the overall computation time, using a model for how the mapping influences computation time. Currently we use a heuristic for defining the mapping. Briefly, the Bounding Box of the processors of the incident elements is relevant for the distribution of a certain node. More precisely, the set of processors which know about a node is defined by the extrema of the coordinates of the processors on which any element resides to which this node belongs. Hence a processor may have to know about some nodes which do not belong to any of its elements. These transit nodes are necessary for ‘routing’ during force communication. All other nodes are called real nodes. We call a mapping admissible iff for every node the coordinates of the involved processors differ at most by 1 in any of the three communication directions. In other words, the bounding box must not extend over more than two processors in any direction. This condition allows the forces concerning one node to accumulate using local communication only, because local communication according to the scheme proposed below reaches the direct neighbors (in all three directions) in one time step, but not the next but one. In case of a non-admissible mapping, the offending nodes must be split into several locally admissible parts, and global communication must merge the forces of the parts. A given node on a processor may possibly require communication in direction x, y and / or z, either with the predecessor (minus) or the successor (plus) in that direction. In this case the node belongs to the x, y and / or z border. The border nodes divide all elements into four classes (batches), defined by the minimal border of any node of the element. That is, an element having a node on the x, y or z border belongs to the corresponding batch, x taking precedence over y and y over z. Elements without any border nodes belong to the last, ‘inner’ batch. A force acting on a node must be applied exactly once. This is inherently given for inner forces (from recovery). Inertial and gravity forces are applied as an initialization of the forces at a node only on one PE, where the node is

63

circled black in Fig. 5 (right). Other nodes (circled gray) initialize forces to 0. The requirement must equally be satisfied for external and mass forces. The mapping defines the sequence of computation and communication which have to run on the parallel computer. To make effective use of the parallel architecture, the load should be balanced among the single processors, and waiting for data from other processors should be minimized. The procedure sketched below aims at this goal. The term ‘recovery’ means determining the interior reaction forces within one element (Eq. (6)). The summation over all elements is split into four phases with communication interlaced between them. Nodes on the x border are marked with a horizontal arrow. Before sending in direction x, all elements with such a node must be computed; they are in batch x, and the figure shows them shaded in a darker gray. While the forces are being transferred, elements in the y batch (lighter gray) are recovered. (No z batch exists in the 2D illustration.) The remaining elements are local to the PE; they belong to the ‘inner’ batch (white).

The following structogramg represents recovery of all elementss in one time step ( 5 meaning a gather operation and 2 5 a scatter / decrement).

The for loop is now split in four phases with communication interlaced between them. The grouping on the right indicates how a loop can code this sequence more compactly, almost eliminating code duplication. Numbers 0, 1, 2 and 3 name the batches (and directions) x, y, z and ‘inner’.

When possible, the recovery of some ‘own’ elements is

64

´ et al. / Medical Image Analysis 4 (2000) 57 – 66 G. Szekely

computed between sending partially assembled node forces in a certain direction and then using the foreign forces received from that direction. This alleviates temporarily unbalanced load and gives communication some time for transferring the data.

4.3. Collision detection The next part is collision detection, which uses the global network for data transfer. We adopt an approximation to collision detection that is commonly used with the FEM. This method only determines which nodes of one contact partner, the slave, penetrate any face of the other contact partner, the master. Flipping the roles of master and slave and averaging will balance this obvious asymmetry. To save time on collision detection, the dynamic behavior of the model is taken into account. Elements do not move faster then 10 cm / s, so if two elements are 2 cm apart, the earliest time they can touch is 100 ms or 1000 time steps away. A collision between these elements cannot take place within the next 1000 time steps and therefore needs no attention. Secondly, the system is built in a hierarchical manner. Every PE first checks bounding boxes with every other PE, before checking for collisions on an element level. Thus collision detection takes the following steps: • Every 1000 time steps all bounding boxes are compared. • If close enough bounding boxes are detected, the corresponding PEs establish a communication channel to exchange the displacements of nodes that may be involved in collisions. • During the next 1000 time steps, the displacements of these nodes are sent to and received from the other PEs. Collisions are detected and each PE determines the consequences for its nodes. Since the communication channels are only redefined every 1000 time steps, software overhead is minimized.

5. Prototype simulator for laparoscopic gynecology Although the dedicated hardware is not yet assembled and functional, we have implemented the algorithms both for the preparatory steps of partitioning the mesh and for the actual parallel computation, as well as a complete system, simulating not-yet-existing devices. We can identify the following five units that work together in the simulator (Fig. 6); each is a more or less independent program. The instrument device encodes the position where the user puts the virtual instrument, and it feeds back the resulting contact forces to the user. The optic device encodes the position of the virtual endoscope and hence the viewpoint of the camera. The mechanic engine performs the computationally expensive FEM calculations. The graphics engine renders the scene as seen

Fig. 6. The prototype has five distinct parts, like the final surgery simulator: an optic device, an instrument device, a communication link (ix), a graphics engine, and a mechanic engine.

through the endoscope. These four parts are connected through a link we chose to call ix; it establishes all communications between the above parts. Fig. 7 sketches the flow of data between these five parts of the surgery simulator. The position of the instrument goes both to the mechanic engine, where it may cause contact, and to the graphics engine, which must render the instrument as long as it is visible in the scene. Mechanic reaction forces propagate back to the force feedback device of the instrument. The new positions of the surface nodes of the deformed organs must be communicated to the graphics for each frame, this includes surface normals. And the graphics must learn the position of the optic device for defining the camera position. In actual surgery, the endoscope may well collide with an organ (causing it to deform which would be handled by the mechanic, too), or even with the instrument. But we do not yet model these possibilities computationally, and the corresponding lines are dashed. The modular setup of the simulator makes it simple to reconfigure, allowing easy adjustment to specific user needs. Different versions of the modules have already been developed including dummy versions, which do just the minimal handshaking necessary to comply with the protocol. Fully functional versions allow the integration and support of different hardware components. The instrument device, e.g., has two active versions: one presents a mousebased interface in an X11 window; the other interfaces to a PHANToM force feedback device. Of course we cannot compute a realistic FE model in real time with this prototype. We use a very simplistic model of 21 elements and achieve time steps of about 2 ms. But the prototype demonstrates the interplay of all relevant components of the system and the validity of our parallelization concept.

´ et al. / Medical Image Analysis 4 (2000) 57 – 66 G. Szekely

65

Fig. 7. Flow of data between the parts of the simulator. The fat lines symbolize data flow ∞. The numbers below the boxes are the repetition rates that our design aims for. Numbers next to lines indicate how many real numbers need to be transferred in one cycle of the slower participant.

Partitioning is a preparatory step and can be done even today without the need of any particular hardware. It effectively schedules in all details the steps of the later computation and exchange of partial results in the mechanic engine, which will be implemented as a 3D lattice of PEs. Today, the parallel computation itself is carried out in several UNIX processes, each process simulating one PE. We ran them as several processes on one CPU, as processes on different hosts, or as processes on different CPUs of an eight-processor SGI Origin 2000; the latter being the fastest. The communication channels are modeled using MPI.

been already implemented using commercial hardware, demonstrating the feasibility of the developed algorithms. The highly modular architecture of the simulator system will allow the successive replacement of its building blocks by modules serving more advanced hardware components during ongoing development.

Acknowledgements This work has been supported by grant No. 5480 / 412625.5RH of the Swiss Federal Institut Zurich and the grant No. 309 / 1995 of the EMDO Foundation, Switzerland.

6. Conclusion Explicit finite element analysis leads to very stable simulations even in the case of large deformations and contact problems. The use of a total Lagrange formulation increases the stability and eliminates integration errors in the static solution. Reduced integration decreases the calculation costs by a factor of 5–7 and nevertheless leads to very stable and sufficiently accurate results when the stabilization forces are calculated with respect to the initial configuration. Appropriately designed parallel algorithms enable implementation of the explicit integration method on a large, three-dimensional network of high-performance processor units, allowing the subsequent use of these techniques even in real-time surgery simulator systems. While the presented results on the deformations of the corpus and fallopian tubes of the uterus during diagnostic gynecological laparoscopy still do not provide highly realistic simulation of tissue deformation, they clearly demonstrate the usefulness and future potential of the developed methods. A first fully functional prototype of the simulator has

References Alyassin, A.M., Lorensen, W.E., 1998. Virtual endoscopy software application on a PC. In: Proc. MMVR’98. IOS Press, pp. 84–89. Barr, A.H., 1984. Global and local deformations of solid primitives. Comput. Graphics 18 (3), 21–30. Basdogan, C., Ho, Ch-H., Srinivasan, M.A., Small, S.D., Dawson, S.L., 1998. Force interactions in laparoscopic simulations: Haptic rendering of soft tissues. In: Proc. MMVR’98. IOS Press, pp. 385–391. Bathe, K.J., 1996. In: Finite Element Procedures, Prentice Hall, Englewood Cliffs, NJ. Baur, Ch., Guzzoni, D., Georg, O., 1998. Virgy: A virtual reality and force feedback based endoscopy surgery simulator. In: Proc. MMVR’98. IOS Press, pp. 110–116. Belytschko, T., Ong, S., 1984. Hourglass control in linear and nonlinear problems. Comput. Methods Appl. Mechanics Eng. 10, 251–276. Bro-Nielsen, M., Cotin, S., 1996. Real-time volumetric deformable models for surgery simulation using Finite Elements and condensation. Comput. Graphics Forum 15 (3), 57–66. Bro-Nielsen, M., Helfrick, D., Glass, B., Zeng, X., Connacher, H., 1998. VR simulation of abdominal trauma surgery. In: Proc. MMVR’98. IOS Press, pp. 117–123. Cotin, S., Delingette, H., Clement, J.M., Bro-Nielsen, M., Ayache, N.,

66

´ et al. / Medical Image Analysis 4 (2000) 57 – 66 G. Szekely

Marescaux, J., 1996. Geometrical and physical representations for a simulator of hepatic surgery. In: Proc. MMVR’96. IOS Press, pp. 139–151. Cover, S.A., Ezquerra, N.F., O’Brian, J.F., 1993. Interactively deformable models for surgery simulation. IEEE Comput. Graphics Appl. 13, 68–75. Crisfield, M.A., 1991. Non-linear finite element analysis of solid and structures. In: Essentials, Vol. 1. Wiley, Chichester. Downes, M., Cavusoglu, M.C., Gantert, W., Way, L.W., Tendick, F., 1998. Virtual environments for training critical skills in laparoscopic surgery. In: Proc. MMVR’98. IOS Press, pp. 316–322. Flanagan, D.P., Belytschko, T., 1984. Eigenvalues and stable time steps for the uniform strain hexahedron and quadrilateral. J. Appl. Mechanics 51, 35–40. Gibson, S., Samosky, J., Mor, A., Fyock, C., Grimson, E., Kanade, T., Kikinis, R., Lauer, H., McKenzie, N., Nakajima, S., Ohkami, H., Osborne, R., Sawada, A., 1997. Simulating arthroscopic knee surgery using volumetric object representations, real-time volume rendering and haptic feedback. In: Proc. CVRMed’97. Springer, pp. 369–378. Hahn, J.K., Kaufman, R., Winick, A.B., Carleton, Th., Park, Y., Lindeman, R., Oh, K.-M., Al-Ghreimil, N., Walsh, R.J., Loew, M., Gerber, J., Sankar, S., 1998. Training environment for inferior vena caval filter placement. In: Proc. MMVR’98. IOS Press, pp. 291–297. Hinton, E., 1992. In: NAFEMS – Introduction to Nonlinear Finite Element Analysis. Bell and Bain, Glasgow. Kojic, M., Bathe, K.J., 1987. Studies of finite element procedures-stress solution of a closed elastic strain path with stretching and shearing using the updated lagrangian jaumann formulation. Comput. Structures 26 (1 / 2), 175–179. ¨ ¨ Kuhnapfel, U.G., Krumm, H.G., Kuhn, C., Hubner, M., Neisius, B., 1995. Endosurgery simulations with KISMET: A flexible tool for surgical instrument design, operation room planning and VR technology based

abdominal surgery training. In: Proc. Virtual Reality World’95. Stuttgart, pp. 165–171. Pommerell, C., 1992. Solution of large unsymmetric systems of linear equations. PhD Thesis, Department of Electrical Engineering, ETH ¨ Zurich. Roy, S. et al., 1992. On the use of polar decomposition in the integration of hypoelastic constitutive laws. Int. J. Eng. Sci. 30 (2), 119–133. Sagar, M.A., Bullivant, D., Mallinson, G.D., Hunter, P.J., Hunter, I.W., 1994. A virtual environment and model of the eye for surgical simulation. Comput. Graphics 28, 205–212. Satava, R., 1996. Virtual endoscopy: Diagnosis using 3-D visualisation and virtual representation. Surg. Endosc. 10, 173–174. Sederberg, T.W., Parry, S.R., 1986. Free-form deformation of solid geometric models. Comput. Graphics 20 (4), 151–160. Suzuki, N., Hattori, A., Ezumi, T., Uchiyama, A., Kumano, T., Ikamoto, A., Adachi, Y., Takatsu, A.A., 1998. Simulator for virtual surgery using deformable organ models and force feedback system. In: Proc. MMVR’98. IOS Press, pp. 227–233. Taylor, V.E., 1991. Application-specific architectures for large finiteelement applications. PhD Thesis, Dept. Elect. Eng. and Comput. Sci., Univ. of California, Berkeley, CA. Terzopoulos, D., Platt, J., Barr, A., Fleischer, K., 1987. Elastically deformable models. Comput. Graphics 21 (4), 205–214. Vining, D.J., 1996. Virtual endoscopy: Is it really?. Radiology 200, 30–31. Ziegler, R., Mueller, W., Fischer, G., Goebel, M., 1995. A virtual reality medical training system. In: Proc. 1st Int. Conf. on Comput. Vision, Virtual Reality and Robotics in Medicine, CVRMed’95, Nice. Springer, Lecture Notes in Comp. Sci., Vol. 905, pp. 282–286. Zienkiewicz, O.C., Taylor, R.L., 1994. The Finite Element Method. McGraw-Hill, London.