Virtual Reality Simulation of Liver Biopsy with a Respiratory Component

1. Introduction

The use of simulation to train skills has grown exponentially in the last few decades, especially in medicine, where simulator models can provide a platform where trainees can practice procedures without risk or harm to patients. Potential advantages of computer based simulator models over other forms of medical simulators include:

1. anatomical fidelity;

2. the ability to train in an environment wherein physiological behavior is observed, something that is not possible using in-vitro phantoms;

3. variability of anatomy and pathology which is important to the acquisition of experience;

4. quantification of metrics relating to task performance that can be used to monitor trainee performance throughout the learning curve; and

5. cost effectiveness.

The work described herein was conducted as part of the CRaIVE collaboration (Gould et al., 2004), the aim of which is to develop simulator models specific to interventional radiology (IR). Our use of the term ‘interventional radiology’ is intended to encompass all micro-surgical, minimally invasive procedures that use percutaneous needle access for biopsy or for guidewire and catheter introduction and manipulation, and where medical imaging (i.e. ultrasound, fluoroscopy, computed tomography) is used for intra-procedural guidance. A frequently performed, IR technique is liver biopsy, where a needle is inserted into the liver to acquire a tissue sample for histological diagnosis of conditions such as cancer, inflammation, cirrhosis etc.

In this chapter we present a detailed description of a simulator framework for liver biopsy using patient-specific data. Key features of this simulator include the incorporation of modifiable respiratory movements and real-rime haptic interaction. We first present an overview of the current state of the art of medical simulators. Subsequent sections are intended to provide a complete description of parameters required to design a medical simulator for liver biopsy, detailing methods used to produce a computer-based model of both patient respiration and haptic interaction, as well as their integration into the final simulator.

2. State of the art

Various simulators exist to train IR practitioners (interventional radiologists). In the case of liver biopsy, which has the added challenge of respiratory motion, only non computer-based simulators currently exist (Müller et al., 2007). As outlined in the introduction, there are, however particular benefits of using computer-based medical simulators. We focus below on different aspects in the literature that are considered to be integral to the accurate simulation of liver biopsy.

Although liver displacement due to respiratory motion has not been incorporated into existing liver biopsy simulators, previous authors have sought to model respiration for a range of applications, including animated entertainment (Zordan et al., 2004) and medical applications. These latter models described in the literature focus on medical augmented reality (Santhanam et al., 2007) and pre-procedural planning of radiological interventions (Didier et al., 2007; Hostettler et al., 2008). For future simulators to find application within clinical practice, they must be able to offer a realistic real time depiction of anatomical structures and their dynamic interaction. Thus, our aim was a simulator that is fast enough to support real-time, or near real-time, computer graphics techniques and interaction, but realistic enough for the specific medical application.

To date, needle insertion simulators have been developed for a range of applications (Abolhassani et al., 2007). In one example, complex mechanical models are used to simulate needle insertion and soft-tissue behavior in the setting of pre-procedural training and planning (Chentanez et al., 2009). Experimental data has also been used by others to model the force evolution during needle manipulation (Maurin et al., 2004). The requirement for a liver biopsy simulator to incorporate respiratory motion does, however, limit us to a relatively simplistic needle model.

The incorporation of force feedback in simulators has been recently reviewed (Coles et al., 2010). Application of this technology is intended to enhance the user’s immersion in a virtual reality environment. Commercial devices are already available for both basic laparoscopic tasks and interventional radiology procedures. Figure 1 illustrates two examples of these simulators: the LAP Mentor by Simbionix for laparoscopic procedures (www.simbionix.com) and the Mentice VIST for endovascular IR procedures (www.mentice.com). Increased access to commercial haptic devices has helped support the application of this technology in new simulators.

Figure 1.

Existing micro-surgical simulators: (a) LAP Mentor by Simbionix, (b) VIST by Mentice.

3. Construction of the virtual patients

We now focus on the construction of the proposed virtual environment, defining the anatomy used in our simulation as well as the software framework, prior to detailing the simulation of organ behavior in a breathing patient during needle insertion.

3.1. Strategy and input data

The initial objective was to construct virtual patients for use in the proposed liver biopsy simulator. To achieve this, it was necessary to first construct a patient database for visualization of organ position and motion (graphic rendering), and subsequently to adapt this environment to support tool-tissue interactions within a collision detection loop (haptics rendering). This use of patient imaging data was passed by the local research ethics committee.

A computer-based simulator was constructed based on the database, and depicting the anatomy of a virtual patient. In the current project we chose to extract this information from medical imaging of patients. Three-dimensional computed tomography (CT) scans of patients attending the Centre Léon Bérard (Lyon, France) and the Marie Curie Institute (Paris, France) were acquired. Scan parameters were as follows: size 512x512xN_z with N_z varying from 140 to 200, and voxel dimensions 0.976mm x 0,976mm x 2mm.

In addition to the liver, it is necessary to also include the ribcage and thoracic diaphragm, as they are essential to modeling respiratory motion. The strategy adopted for data acquisition was as follows:

1. retrieval of anonymized scan data for recruited patients after obtaining informed written consent;

2. extraction of organ contours from scans; and

3. computation of a more simplistic representation of organ geometry via a meshing process.

Each of these steps is illustrated in Figure 2 using the ribcage as an example.

Figure 2.

Construction of the virtual patient using the ribcage as an example: (a) 3D CT scan showing the ribcage, (b) Voxels resulting from segmentation of the ribcage, (c) Polygonal mesh used in the simulation.

3.2. Segmentation

Contour data can be extracted from medical imaging either by automatic or manual segmentation. The current project utilized automatic segmentation for bony anatomy (spine and rib cage), and manual segmentation for the liver and diaphragm, where interpretation by a medical expert is required.

3.2.1. Bones, lungs and skin segmentation

All segmentation during the current project was performed using ITK-snap, an open-source program developed by the University of Pennsylvania (Yushkevich et al., 2005). Using this program, automatic segmentation is achieved by a process based on the level set technique. This method is divided into three steps:

1. bounded region is defined in order to drive the contour evolution. This boundary is given by the double-threshold image binarization method.

2. Spherical seeds of different diameters are placed at different locations inside the previously defined boundaries, and

3. an evolution function is applied to these seeds. This time-depending function, C_t, is defined by equation (1).

$$C_t=(\delta g_I-\xi K)\overrightarrow{n}$$ E1

where δ and ξ are tunable parameters that have been respectively set to 1.0 and 0.2; g_I is the force intensity based on the image gradient; K is the force intensity based on the contour curvature; and n is the contour normal.

Due to the high density of bone, segmentation of this tissue was performed using the automatic method with a double threshold (lower and upper) based on high Hounsfield density. The same method is used to segment the lungs. As the lungs are composed primarily of air, a double threshold based on low Hounsfield density is employed. The air-tissue interface at the skin can likewise be segmented by this method based on low Hounsfield density.

3.2.2. Liver segmentation

Segmentation of the liver by automated methods poses a far greater challenge, because of the difficulty in performing the initial step of the level set algorithm due to the similarity between the density of the liver and its surrounding organs and tissues. Whilst methods exist to overcome this problem (Chen et al., 2009; Foruzana et al., 2009; Rusko et al., 2009), they are seldom satisfactory. Consequently, manual segmentation of organs such as the liver, whilst time consuming, is preferred. Manual segmentation can be more reliably and easily performed by using graphic tablets, such as the Cintiq 21UX Tablet PC (Wacom Co Ltd., Saitama, Japan) that was used in the current project. Such tablets, used in combination with a pressure sensitive stylus, permit a trained expert performing the segmentation to trace, freehand, the border of the liver in consecutive CT scan images in the sagittal, coronal and axial view. As with automatic segmentation, manual segmentation was performed using the ITK-snap software.

3.2.3. Diaphragm segmentation

Due to its thickness, typically 3mm (1-2 pixels), it is often difficult to visualize the thoracic diaphragm in standard CT scan images. Moreover, the density of the fibro-muscular tissue which makes up the diaphragm is similar to that of other surrounding organs and tissues. Segmenting of this structure, therefore, presented a major challenge. The most effective way of manually segmenting the diaphragm was to isolate its boundaries with surrounding structures, namely the liver, lungs and thoracic wall. For the diaphragm, manual segmentation in the sagittal, coronal and axial planes is particularly desirable as it helps to ensure a smooth contour. Important anatomical relationships of the diaphragm include:

1. the spine in the region of the twelfth thoracic vertebra; diaphragmatic crura are demonstrable in the axial plane (see Figure 3.a);

2. the inner aspect of the ribs and the costal margin;

3. the overlying visceral and parietal pleura, including the potential space of the pleural reflection;

4. the liver;

5. the heart / pericardium and lungs,

6. stomach and

7. spleen (See Figure 3.b and 3.c).

Figure 3.

Manual segmentation of the diaphragm: red line depicts (a) diaphragmatic crural attachments to the spine, (b) relationship to the lungs, pleura and ribs within the thorax, and (c) inferior aspect of the heart.

3.3. Meshing

For this project we have segmented scans from five patients with liver cancer, each with variation in both normal anatomy and tumour location. Data acquired from segmentation leads to the creation of ‘stacks’ of voxels that can be labeled according to the organ or structure they depict. These data are subsequently transformed into 3D models via a process of meshing, itself divided into three clear stages:

1. marching cubes algorithm;

2. decimation; and

3. smoothing.

Marching cubes is a computer graphics algorithm for extracting a polygonal mesh of an isosurface from voxels. This isosurface comes directly from the label boundaries that are the result of segmentation. Decimation consists of replacing a large group of polygonal elements by a smaller group, given some constrains of curvature and position accuracy. The last technique, smoothing, is a methods used to diminish sharp edges within the model by “blurring” the angles between triangular elements.

3.4. Implementation

Once fully acquired and preprocessed, this 3D organ framework is implemented in object-oriented C++. We use the application programming interface H3D (SenseGraphics, 2008) as this environment combines both visual and haptic rendering. The H3D visual rendering is a superior layer to OpenGL, permitting the fast depiction of sets of triangles as organ meshes. The haptic rendering is also a superior layer to OpenHaptics, the software provided with the haptic device used in this project, the Phantom Omni (SensAble, Wilmington, USA). H3D uses the XML format called X3D to represent a scene or, in this particular instance, a set of organs, with different programmable behaviors. Figure 4 shows two examples of virtual patients that were developed within this framework.

The steps described up to this point have led to the creation of static models of the organs and tissues required for construction of the proposed liver biopsy simulator. Using a suitable deformation model, this static 3D representation of the organs and tissues can be manipulated without difficulty to allow deformation, either during breathing, or when there is interaction with a surgical instrument. In addition, it is required that force feedback is rendered to the user while inserting a needle through the body of the virtual patient. The next two sections of this chapter detail the methods used to achieve these goals.

Figure 4.

Example of two virtual patients with transparent skin (light brown), lungs (pink), bones (white), liver (dark brown) and thoracic diaphragm (red)

4. Respiration simulation

In order to create a simulator which convincingly recreates conditions of liver biopsy, it is necessary to model liver displacement by respiration. In the current project this was achieved by applying dynamic behavior to previously segmented anatomy. This motion is defined by a sinusoidal curve that mimics normal tidal breathing in a patient. The curve M(t) modeled in equation (2), is used to defined the motion of respiratory muscles.

$$M(t)=A\sin (2\pi ft)$$ E2

Where A is the muscle amplitude and f the respiration frequency. As is common during real liver biopsy, the clinician can instruct the virtual patient to perform an end tidal breath hold, whereupon motion due to respiration stops. In the following section we first introduce the mechanism of normal human respiration, and subsequently discuss details of the soft-tissue deformation technique that was used, as well as the handling of deformation for each organ.

4.1. Description of the natural respiration process

Respiration is a process consisting of the cyclical inflation and deflation of the lungs. The lungs are in effect passive organs whose motion is driven by muscles within the chest wall and diaphragm. The lungs are surrounded by the pleural cavity, bound by the visceral (lung surface) and parietal (chest wall and diaphragm) pleural membrane and containing a small volume of pleural fluid. Surface tension within the pleural fluid ensures close apposition of the lung surfaces with the chest wall and diaphragm, allowing for optional lung inflation. During inhalation, the active downwards contraction of the diaphragm and outwards motion of the chest wall (mainly due to contraction of the intercostal muscles) results in a negative pressure being generated within the thorax, which is in turn responsible for lung expansion. In comparison, exhalation, for the most part, is a passive process facilitated by relaxation of diaphragmatic and chest wall muscles.

The upper portion of the liver is attached to the diaphragm via a collection of ligaments (see Figure 5). The falciform ligament extends upwards on the anterior surface of the liver from a notch in the lower margin, as far back as the posterior surface. On the superior aspect of the liver, the falciform ligament extends medially as the left triangular ligament and laterally as the coronary ligament. The coronary ligament is formed of both an upper and lower layer which meet at the apex of the liver, where they form the right triangular ligament.

Figure 5.

Positions of liver ligaments; (a) anterior and (b) posterior views of the liver

The diaphragm is a sheet-like structure formed of both a tendinous (central) and a muscular (peripheral) part. During contraction, the diaphragm flattens, enlarging the volume of the thoracic cavity and causing the lungs to expand and inflate. When relaxed, the diaphragm is pulled upwards by the elastic recoil of the lung tissue. Due to its ligamentous attachments and close anatomical relationship, during normal respiration the liver follows the downwards and upwards motions of the diaphragm (Moore &Dalley, 2009).

The liver is also surrounded by the lower portion of the chest wall, which comprises ribs (bony elements) and costal cartilage (cartilaginous component of the rib cage), and associated muscle and other soft tissues. All twelve pairs of ribs are attached posteriorly to their respective thoracic vertebrae. Anteriorly, ribs 1 to 10 attach to the sternum via costal cartilage. Ribs 11 and 12 are considered ‘floating ribs’ as they have no anterior connection to the sternum. The rib cage, whilst rigid, is able to move in what is classically described as a combination of ‘bucket handle’ and ‘pump handle’ motions. Contraction of chest wall muscles generates an outwards and upwards motion of the anterior chest wall, causing expansion of the thoracic cavity. The attachments of the ribs to the thoracic vertebrae act as a pivot for this movement. Motions and deformations that are described from the diaphragm, ribcage and liver are shown in Figure 6.

Figure 6.

Physiological behavior and motion of the: a) diaphragm, b) ribcage and c) liver. The black arrows represent the organ deformations during breathing.

Both the diaphragm and the ribcage are therefore independently involved in the process of normal respiration and, as such, their motions influence the position of the liver. Accordingly, a model of respiratory motion is essential to the development of a realistic simulation of liver biopsy.

4.2. Soft-tissue behavior simulation

Modeling soft-tissue behavior for computer-based simulators has been widely studied over the last few decades. Efforts have primarily focused on improving mechanical fidelity so as to best mimic reality. Some algorithms use, directly, the continuous mechanics equations that are adopted in the field of mechanical engineering, whilst others use the time-consuming finite element formulation (Bro-Nielsen, 1996). Algorithms that use fewer equations have also been developed, including the boundary element method (James & Pai, 1999) and the mass-tensor method (Cotin et al., 1998), but their computing time still remains high. Another approach that attempts to minimize computational time, requires that the object being modeled is no longer considered as a continuous body, but rather as one that is discretized. The most recognized example of such a model is the mass-spring method, which consists in linking discrete elements by springs and then using the proportional spring interaction law (Zordan et al., 2004). A similar method is the particle-system, where the elements are no longer connected, but are restricted by a repulsion/attraction law that helps to ensure the systems stability (Amrani et al., 2000). We have chosen to use the Chainmail (Meier et al., 2005) method as it is one of the fastest available for simulation of soft tissue behavior.

A full description of the Chainmail method can be found in (Gibson, 1997). The primary aim of this method is to simulate soft tissues as if they were a sheet of medieval chain mail, composed of small, interconnected metal rings. The expected behavior is as follows: when a chain element is moved, it may influence its neighbors; if the initial motion is too small, only the first element is moved but, if this motion is large enough, neighboring elements may also moved (see Figure 7).

Figure 7.

Chainmail behavior: the blue element is moved and the red elements may move depending on the first element displacement amplitude and direction. Element displacement that is (a) small enough to not influence its neighbors and (b) large enough to influence its neighbors.

Mathematically, this model can be described by the mechanical parameters of compression, γ; stretching, λ; and shearing, μ. If we focus on an axis, i, and on the motion of two elements along this axis, the second element will not move if the length between the elements remains below a quantity i_max and above a quantity i_min defined by γ, λ and μ. Consequently, compression is observed when the length becomes small enough i.e. i_min is equal to γ times the initial length with γ<1. Stretching is observed when the length becomes high enough i.e. i_max is equal to λ times the initial length with λ>1. Finally, shearing is observed when the combined displacement in the other axes becomes high enough i.e. i_max is equal to μ times the initial combined displacement with μ>1. The behavior defined by the Chainmail algorithm is given by equations (3-8).

$$x_{\min }=x_2+(\gamma *\left|x_2-x_1\right|-\mu *(\left|y_2-y_1\right|+\left|z_2-z_1\right|))$$ E3

$$x_{\max }=x_2+(\lambda *\left|x_2-x_1\right|-\mu *(\left|y_2-y_1\right|+\left|z_2-z_1\right|))$$ E4

$$y_{\min }=y_2+(\gamma *\left|y_2-y_1\right|-\mu *(\left|x_2-x_1\right|+\left|z_2-z_1\right|))$$ E5

$$y_{\max }=y_2+(\lambda *\left|y_2-y_1\right|-\mu *(\left|x_2-x_1\right|+\left|z_2-z_1\right|))$$ E6

$$z_{\min }=z_2+(\gamma *\left|z_2-z_1\right|-\mu *(\left|x_2-x_1\right|+\left|y_2-y_1\right|))$$ E7

$$z_{\max }=z_2+(\lambda *\left|z_2-z_1\right|-\mu *(\left|x_2-x_1\right|+\left|y_2-y_1\right|))$$ E8

These equations give the minimum and maximum location of elements linked at a moving element, whose position goes from (x₁, y₁, z₁) before displacement, to (x₂, y₂, z₂) after displacement.

4.3. Rib cage simulation

Due to their rigid nature, it is not necessary to use soft-tissue deformation to model movement of the ribs. Their motion can instead be based on kinematics laws. Furthermore, we assume each rib to have only rotation and no translation, unlike the finite helical axis model used by previous authors (Didier et al., 2007). As this model is only defined by rotational motion, it is comparably simpler (Wilson et al., 2001). In this project, rotational parameters are independently defined for each rib; and are composed of a rotation center, two rotation axes, a minimal rotation angle and a maximal rotation angle. The first axis defines a ribs ‘bucket handle’ rotation (angle α) and the second axis defines the ‘pump handle’ rotation (angle β), as seen in Figure 6. These angles have been previously defined by other authors (Wilson et al., 2001) using patient data. The statistical amplitude angles, α_max and β_max, are presented on Table 1.

Rib number	αmax	βmax
2	14,3	13,7
3	11,4	13,3
4	10,7	10,1
5	9,6	8,9
6	9,4	6,9
7	7,9	6,6
8	7,9	6,2
9	6	6,3

Table 1.

List of angle amplitude for each rib and in the two rotation plans

We have applied such methodology during the current project on datasets acquired from patient scans. To achieve this, we initially computed the rotation parameter during a preprocessing stage and chose to use a mechanical approach by computing the inertia matrix. This matrix, M, is defined by equation (9). It corresponds to integrations over the volume of each rib; x, y and z being the position of a within a rib.

$$M=\left(\begin{array}{ccc}{\displaystyle \int y^2+z^{2}dv}& {\displaystyle \int xydv}& {\displaystyle \int xzdv}\\ {\displaystyle \int xydv}& {\displaystyle \int x^2+z^{2}dv}& {\displaystyle \int zydv}\\ {\displaystyle \int xzdv}& {\displaystyle \int zydv}& {\displaystyle \int y^2+x^{2}dv}\end{array}\right)$$ E9

These moments of inertia represent the rigidity of a 3D object and are widely used in solid mechanics. The Eigen vectors of this matrix directly give the principal axes of the objects. The two highest Eigen values give the axis that contains the largest proportion of matter and is therefore of most interest to us. We define the ‘rib plane’ as that which is given by these two Eigen vectors (see Figure 8).

Figure 8.

Rib rotation algorithm. Inputs are: (i) Eigen vectors; (ii) rotation angles α and β; (iii) rotation center; and (iv) initial rib location. The output is the rib at its final position.

The rotation centers are found by a heuristic method. First, we compute the higher Eigen value of the spine inertia matrix, the associated axis gives the axial line of the spine. The rotation center is then given by the intersection between this line and the rib plane. From this pre-computation stage, we can determine for each rib the rotation center and the two axes. The next step is to apply a dynamic behavior to each rib based on these patient-customized parameters. At each time step, the new coordinates of the rib points are computed with Rodrigues' rotation formula (Koks, 2006). These two rotations are time-dependant and their angle values are derived from equation (2) taking α and β as parameters. The results are given by equations (10) and (11). Acquired behavioral data was stored and accessed in X3D format.

$$\alpha (t)={\alpha }_{\max }\sin (2\pi ft)$$ E10

$$\beta (t)={\beta }_{\max }\sin (2\pi ft)$$ E11

4.4. Diaphragm simulation

Diaphragmatic behavior has been modeled using various techniques, including geometrical methods based on external sensors (Hostettler et al., 2008); Mass-spring systems with incompressibility constraints (Zordan et al., 2004); and nature-inspired tensegrity (Villard et al., 2008). These methods are time-consuming and their accuracy is beyond that which is required by the current application. Accordingly, we have developed a faster method that still reflects the physiological parameters discussed in section 4.1.

During normal inhalation, the peripheral muscular portion of the diaphragm contracts pulling downwards on the central tendon and increasing intrathoracic volume. In comparison, during exhalation, the central tendon of the diaphragm retracts upwards into the chest as the muscle relaxes. We have, however, proposed a model whereby the central tendon of the diaphragm assumes an active role, with a sinusoidal up and down motion derived from equation (2), whilst the muscular portion remains passive. As the central tendon of the diaphragm is naturally stiff, our model assumes it to be a rigid structure. The motion function, T(t), of the model, where t is the time as given by equation (12). The translation amplitude, d_max, was set at 2cm.

$$T(t)=d_{\max }\sin (2\pi ft)$$ E12

In contrast, the muscle portion of the diaphragm follows an elastic behavior based on the Chainmail algorithm detailed in section 4.2. This also preserved the natural deformation of the muscle during the contraction-relaxation cycle. The portion of the diaphragm, which is attached to the spine, is programmed to remain motionless. Illustrations of diaphragm behavior, as programmed by our model, are presented in Figure 9.

Figure 9.

Programmed diaphragmatic behavior: the central tendon of the diaphragm, presented here in white, is considered to have a stiff translation in the longitudinal plane. In comparison, the muscular portion (red), with the exception of regions directly attached to the spine, passively deform in response to tendon displacement. The initial and final positions of the diaphragm are depicted above by opaque and transparent boundaries respectively.

4.5. Liver simulation

The final organ whose behavior needs to be detailed is the liver. Whilst the liver is entirely passive during respiration, its upper border lies immediately adjacent to the diaphragm and is attached to its inferior surface via a series of ligaments (refer to details in section 4.1). To model this interaction we chose to apply coupling forces as opposed to displacement constrains, allowing continuous smooth deformation that would otherwise have been hard to reproduce. Coupling between the liver and diaphragm was characterized in a pre-computing step using X3D data. A threshold distance of 5mm was established, below which a connection force is applied. Internal deformations of the liver due to respiration were modeled by the Chainmail algorithm.

Our model also allowed deformation of the liver under the influence of an external force applied by the biopsy needle. This interaction primarily occurs between the tip of the needle and what is in real life a thin fibrous layer covering the liver, referred to as the liver capsule. This layer is deformed as the operator attempts to introduce the biopsy needle into the body of the liver. Once perforated, however, this layer no longer shows deformation. Again, we use the Chainmail algorithm to model this interaction by applying the needle tip displacement to the three closest points on the liver surface. In reality, the greatest influence on needle insertion is not the visual displacement of the liver capsule, but haptic interaction. In the next section we focus on the method used to render this important ‘force feedback’ experienced by the operator during needle insertion.

5. Haptic rendering

As in many simulators, we have attempted to enhance the user’s immersion through the addition of haptic rendering. This is achieved by the operator interacting with the virtual environment through a haptic device that resembles tools used during the actual task. In this section we introduce the synchronization between the visual feedback and the haptic device, as well as explain the computation and rendering to the user of the physical forces.

5.1. Needle avatar

In this project a haptic device is used to simulate the liver biopsy needle. This device generates two important parameters used by our framework; namely, the position (P) and the direction (D) of the needle. The needle avatar is then linked to the stylus that the user is manipulating. At each time step, P and D are recorded and the needle is displayed. When the needle is inserted into the patient’s body, a force feedback (F) is computed and rendered back to the user (see Figure 10).

Figure 10.

Commercially available haptic device. The haptic device is driven by the position (P) and direction (D) of the stylus that are determined by the user’s movements, and renders back a force (F) processed by the computer.

5.2. Soft-tissue force feedback

Haptic rendering is often based on collision detection and collision response algorithms. The combination of these methods could not, however, be directly applied during the current project as it would have resulted in a model that was too complex due to the number of organs and tissues potentially traversed by the needle. In addition, it would need haptic sub-layers, for example the organs inside the skin or the tumor inside the liver. Such an option is currently not fully supported by haptics APIs.

We chose instead the method implemented in (Vidal et al., 2008), where force is no longer dependent on collision detection, but is always computed while the needle is inside the patient’s body and never while outside. In the same way, no classical collision response force algorithm is performed. Instead, force is computed using the initial 3D CT scan data, without the mesh environment information. As the biopsy needle is inserted into the virtual patient, it is possible to access the corresponding voxel Hounsfield densities, provided the medical image remains appropriately registered to the polygonal information. Longitudinal force is then determined on line by means of the voxel gradient as it is similar to matter density. The operator can therefore feel the boundaries between tissue planes as they are traversed by the advancing needle

Soft tissues are not only responsible for forces generated along the long axis of the needle (needle axial forces), they also apply a constraining force that inhibits lateral movement of the needle. Consequently, once the needle has been inserted more than a few millimeters into the virtual patient, lateral movement of needle is greatly constrained by the influence of the surrounding tissue(s), such that it is almost impossible to deviate from the initial trajectory. In our model we have added constraining forces in order to reproduce this phenomenon. The direction of the force (F_a) acting upon the needle is computed from the initial position (P₀) and direction (D₀) of insertion and from current position (P_t). Thus, the relationship describing the soft-tissue force feedback is given by equation (13).

$$F_a=(P_t-P_0)*\sin (\frac{D_0}{P_t-P_0})$$ E13

Figure 11.

Needle insertion force in Newtons with the penetration depth in millimeters (equation (14).

5.3. Liver force feedback

Haptic feedback is of particular importance during needle insertion into the substance of the liver. As was described in section 4.5, the liver is covered by a fibrous capsule that offers considerable resistance to continued passage of the needle. This resistance is, however, dramatically reduced at the moment of capsule penetration. Thereafter, continued needle insertion deeper into the liver is subjected to increasing resistance. In this project we have simulated this behavior by adding a needle axial force whose value is extracted from (Maurin et al., 2004). Experimental data were acquired with a force sensor during IR procedures (Karuppasamy et al., 2008) and fitted to a mathematical model. It defines the force, F_n, varying with the penetration depth d. This force is defined piecewise, before and after penetration with two sets of parameters τ, ω, υ and φ that define equation (14). The corresponding curve is shown on Figure 11.

$$F_n(d)=(\nu +\phi )e^{\tau (d-\omega )}+\phi$$ E14

6. Resulting simulator

The final simulator encompasses all of the elements described in the preceding sections: a collection of virtual patients and their organs, segmented from CT scans of actual patients and represented by meshes; cyclical deformation of theses meshes in accordance with a model of respiration; a virtual needle that can be manipulated with six degrees of freedom via the arm of a haptic device. As the needle passes through the different tissues, the operator is able to determine additional information about their nature through haptic feedback. Furthermore, as this needle is inserted into the capsule of the liver, deformation of the organ’s surface can be observed and felt. The biopsy is triggered by the operator pressing one of the buttons in the Omni haptic interface. In the following sections we present a complete description of the final simulator, including patient-customised anatomy and physiology. We also add a module to measure the user performance. Finally, we present a review of feedback received from both expert and trainee interventional radiologists.

Figure 12.

Liver biopsy simulator using the Sensegraphics Immersive Workbench and Phantom Omni haptic device

6.1. Immersive Workbench

In developing this simulator particular emphasis was made on selecting algorithms that are suitably fast. This has, ultimately, allowed the simulator to be real time, offering 500Hz and 20Hz for haptic and graphic rendering respectively. To achieve a superior ergonomic and 3D immersive experience, the simulator framework has been integrated with the Sensegraphics Immersive Workbench

http://www.sensegraphics.com/

(see Figure 12). This system required active 3D glasses with a shuttering system synchronized with a CRT monitor to visualize the virtual environment in 3D. As is shown in Figure 12, the monitor is positioned at a 45° angle above a horizontal semi-transparent mirror. With this configuration, the user is able to view the 3D target object collocated with their hand and the haptic device.

6.2. Patient specific physiology

The final 3D simulator has been designed to take into account variability in the speed and strength of factors contributing to motion. This has allowed the simulator to match diaphragmatic or thoracic breathing patterns and the overlap that exists between these two dynamic processes (see Figure 13).

Figure 13.

Patient customized anatomy with inhalation motion and needle avatar controlled via a haptic device

The respiration is mathematically monitored by time-dependant functions defined by equations (10), (11) and (12). The breathing frequency could be changed by modifying the shared parameter f. The user can change from quiet breathing to hyperventilation. In the same way the influence of the diaphragm is linked to the d_max parameter of equation (12) and the thoracic behavior is tuned with a scalar value linked to the α_max and β_max parameters of respectively equations (10) and (11). Figure 14.a shows the panel displayed to the user in order to set a respiration pattern.

6.3. Performance metrics

As previously stated, simulators such as the one described here, provide a valuable and safe environment in which trainees may practice and be assessed. To support this application, a number of computer-based metrics are provided, including: time taken; number of times contact is made with designated ‘no-go’ areas (representing errors of needle tip location); number of attempted needle insertions; and biopsy accuracy (adequacy of tissue sample obtained). The time taken for the procedure begins when the task is initialized and finishes when the biopsy sample is acquired. The number of times voxels designated as ‘no go’ areas are crossed by the biopsy needle is also determined. ‘No go’ voxels include the lungs, costophrenic recess, colon, gallbladder, neurovascular bundles (located below the ribs), and bone. A further metric counts the number of attempted needle insertions; these are repeated whenever the operator perceives a failure to locate the biopsy target on the needle trajectory, as defined by the insertion point and the direction of the virtual needle. The final metric computed by the simulator was biopsy precision. In the current simulator, the tumor is modeled by a spherical body of 15mm. Four levels of precision are defined depending on the distance between the site of biopsy and the center of the tumor. Biopsy samples taken from within 5mm of the centre of the tumor are considered to contain “Necrotic debris”, that is a feature of tumors as they outgrow their blood supply. Biopsies taken from within 5-10mm of the tumor centre are judged as being a “good specimen” because this is the region of actively growing tumour cells, whilst tissues 10-15mm from the center may contain “some tumor cells”. Any biopsy that is >15mm from the centre of the tumor will contain normal cells and hence is regarded as “not pathologic”. All metrics described are displayed to the operator on completion of the procedure, allowing objective assessment of their performance (see Figure 14.b)

Figure 14.

Simulator user interface. a) Tuning of patient breathing pattern b) Metrics display.

6.4. Qualitative evaluation

In order to evaluate the simulator, we invited both expert and trainee interventional radiologists to use it to perform a number of tasks. A review of the main comments made by clinicians is presented below.

Feedback received for each module of the simulator was generally positive. Most clinicians reported that the inclusion of movement due to respiration was an important feature of the simulator and considered that the realism of the procedure was enhanced when this function was active. Many comments were made about the haptic rendering, which was considered to satisfactorily reflect the forces experienced during liver biopsy in patients.

It was suggested that the simulator was most suitable for first and second year trainees who generally have limited opportunity to practice and maintain their technical skills in this procedure. By using the simulator, clinicians felt that trainees would be able to acquire and maintain both core practical skills and confidence, within a safe environment, remote from patients. The general consensus of all clinicians was that simulators, such as this, are needed and that there is a place for them in training curricula. Indeed the Royal College of Radiologists and British Society of Interventional Radiologists in the UK are actively developing training opportunities using simulation

We also received a number of comments relating to potential improvement of the simulator. Several clinicians noted that, during this procedure, they would tend to rest their hands on the patients body, something that was not possible with a virtual patient. One clinician reported difficulty in pressing the button on the Omni haptic device used to trigger biopsy, whilst another found it difficult to manouvre the same device.

7. Future work

Future work must both address issues raised by clinicians during the evaluation stage, as well as modifying the existing framework to allow simulation of other interventional radiology procedures within the abdomen. We also intend to develop new methods of automatic segmentation that will aid in the isolation of anatomical structures, including the diaphragm. It is intended that this will enhance the retrieval of 3D data from CT scans, widening the application of this technology to allow patient specific pre-procedural training and planning. Finally, refinement of the respiratory model will allow the described framework to be used in other applications including the administration of targeted radiotherapy. Further thorough and systematic validation is also required and to this end a skills transfer study has recently commenced.

8. Conclusion

In conclusion, we have presented here a new framework of a liver biopsy simulator that includes respiratory motion. This framework is composed of 3D virtual organs extracted from patient data. In addition, organ motion due to respiration and haptic feedback during needle insertion is computed on-line and in real time within a virtual reality environment. An immersive workbench permits the operator to perform the procedure in a 3D environment that allows co-location of their hands with the virtual target. The feedback of experienced and trainee clinicians who were invited to use the simulator was obtained to evaluate the simulator’s strengths and weaknesses. As a direct consequence of its efficacy, the simulator framework described here has been applied in two other medical simulators that are based on fluoroscopic (Villard et al., 2009) and ultrasound (Bello et al., 2009) image guidance.

Acknowledgments

The author would like to acknowledge the financial support of the UK Department of Health under the Health Technology Devices (HTD) Program during the course of this study.