2D XFEM-based modeling of retraction and successive resections for preoperative image update

Abstract

This paper considers an approach to improving outcomes for neurosurgery patients by enhancing intraoperative navigation and guidance. Currently, intraoperative navigation systems do not accurately account for brain shift or tissue resection. We describe how preoperative images can be incrementally updated to take into account any type of brain tissue deformation that may occur during surgery, and thus to improve the accuracy of image-guided navigation systems. For this purpose, we have developed a non-rigid image registration technique using a biomechanical model, which deforms based on the Finite Element Method (FEM). While the FEM has been used successfully for dealing with deformations such as brain shift, it has difficulty with tissue discontinuities. Here, we describe a novel application of the eXtended Finite Element Method (XFEM) in the field of image-guided surgery in order to model brain deformations that imply tissue discontinuities. In particular, this paper presents a detailed account of the use of XFEM for dealing with retraction and successive resections, and demonstrates the feasibility of the approach by considering 2D examples based on intraoperative MR images. To evaluate our results, we compute the modified Hausdorff distance between Canny edges extracted from images before and after registration. We show that this distance decreases after registration, and thus demonstrate that our approach improves alignment of intraoperative images.

• Intraoperative
• non-rigid
• registration
• FEM
• XFEM
• retraction
• resection

Introduction

The main goal of brain surgery is to remove as much of the lesion tissues as possible, while avoiding contact with eloquent areas and fiber tracts located in white-matter tissues. Surgery is planned on the basis of preoperative images from multiple modalities, such as Computed Tomography (CT), structural and functional Magnetic Resonance Imaging (sMRI and fMRI), Diffusion Tensor Imaging (DTI), Positron Emission Tomography (PET), and Magneto-Encephalography (MEG). Surgery is generally performed using an image-guided navigation system that relates the three-dimensional (3D) preoperative images to 3D patient coordinates, thereby allowing the surgeon to position his instruments in the patient's brain while navigating on the preoperative images. However, throughout surgery, the brain deforms, mostly as a result of leakage of the cerebrospinal fluid (CSF) out of the skull cavity, modifications in cerebral perfusion, pharmacological modulation of the extracellular fluid, and surgical actions such as cuts, retractions and resections [1], [2], [3]. Therefore, as surgery progresses, preoperative images become less representative of the actual brain, and navigation accuracy decreases.

These navigation errors could be reduced if one could acquire, throughout surgery, fresh images of the same modalities and quality as those acquired preoperatively. Doing so presents major challenges, since, except at a handful of surgical facilities, intraoperative MR images are usually acquired using low-field MRI scanners that provide lower resolution and contrast than their preoperative counterparts. Moreover, even in those rare facilities that possess high-field intraoperative MRI scanners, the acquisition of intraoperative DTI images remains rather time-consuming. Finally, several useful imaging modalities, such as PET and MEG, cannot yet be acquired intraoperatively. The solution that is generally favored is to acquire limited-quality intraoperative images at several critical points during surgery, track the tissue deformation from one image to the next, update (i.e., deform) the preoperative images accordingly, and feed the resulting images to the navigation system [4].

The deformation of preoperative images so that they conform to intraoperative images falls into the domain of non-rigid registration. One category of such techniques uses biomechanical models based on the Finite Element Method (FEM) [5]. Prior to surgery, a biomechanical brain model specific to the patient is built from preoperative images: the model consists of a volume mesh of finite elements (FEs) and one or more mechanical behavior laws assigned to them. During surgery, a number of key anatomical landmarks are extracted and tracked through successive intraoperative images. The estimated displacements of these landmarks are applied to the biomechanical model and drive its deformation. The resulting displacement field of the biomechanical model is used to deform the preoperative images. So far, most of the mechanical conditions of the brain, such as the volume of CSF flowing out of the skull cavity, or the forces applied by a retractor tool, cannot be estimated in the operating room. The fact that an intraoperative image can provide knowledge of the current state of the brain after some deformation partly eliminates the need for a complete evaluation of these mechanical conditions. The non-rigid registration techniques replace them with the landmark displacements determined from successive intraoperative images.

Most studies of brain deformation based on biomechanical models have focused on shifts (in which the topology of the brain is not modified), i.e., without explicitly taking into account any cut and subsequent deformation [6], [7], [8], [9], [10], [11]. A good review of parameters that can be included in these biomechanical models can be found in reference [12]. The reported accuracy for deformation prediction is approximately one voxel [13]. The situation becomes more complex when the surgeon performs cuts, retractions or resections, the last two necessarily involving a cut. The main difficulty associated with a cut is the discontinuity it implies in the tissues. Indeed, FEM cannot handle discontinuities directly and requires one to realign the discontinuity with element boundaries based on mesh-adaptation or remeshing techniques that provide from an initial mesh a new mesh that is conformed to the discontinuity.

Extensive work has been performed in the domain of surgical simulation on the problem of cutting through an FE mesh, and the methods can be divided into three main classes. The first and simplest method avoids any mesh adaptation or remeshing, and consists of deleting the FEs touched by the cutting tool [14], [15]. This method has the advantage of avoiding the creation of any new element, but it models the cut with some finite thickness, as if some tissue were removed even though this is not actually the case. The second class of algorithms is based on subdivision methods [16], [17], [18], [19], [20]. All FEs intersected by the cut are divided into sub-elements to create FE boundaries matching the cut, while ensuring that the new FEs have good aspect ratios. The third class comprises the methods of mesh adaptation. The idea here is to relocate existing nodes to cling as well as possible to the cut geometry [21]. The mesh can then possibly be relaxed [22] to avoid large distortions. Finally, hybrid methods have been developed that combine subdivision and mesh adaptation to minimize the creation of new FEs [23].

The field of fracture mechanics has been confronted with an equivalent problem in the study of the growth and propagation of cracks in mechanical parts. Three main methods have been investigated to avoid the drawbacks of using FEM in conjunction with mesh adaptation or remeshing [24]. The first method is the boundary element method (BEM), which is based on the discretization of only the object surface(s) [25], and which has also been also used to simulate minimally invasive surgery involving cutting [26]. The second class of methods is the meshless methods [27], [28], [29]. In FEM, the object is represented by a volume mesh. The nodes interact because they are connected via the elements. In meshless methods, the object is represented by a set of non-connected nodes that interact because their domains of influence overlap. Meshless methods have been used to develop surgical simulation tools [30], [31]. The third method, called the eXtended Finite Element Method (XFEM or X-FEM), appeared in 1999 [32] and has since been the object of considerable research [33]. XFEM works by allowing the displacement field to be discontinuous within some FEs of the mesh. First, a mesh of FEs is built without taking the crack into account. Then, based upon the precise geometry of the crack, new degrees of freedom with associated discontinuous functions are added to some of the existing nodes, allowing a discontinuity in displacement at all points along the crack, without any remeshing. The mesh does not have to be conformed to the discontinuities, which can then be located arbitrarily with respect to the underlying FE mesh.

Until recently, the explicit modeling of retraction and resection for preoperative image update has received much less attention than (brain) shift. Retraction is usually performed when the target, e.g., the tumor, lies deep inside the brain. The surgeon cuts through brain tissues and inserts the blades of a retractor to spread tissues out from the incision and create a free path towards the tumor. Retraction has already been modeled by first duplicating the nodes along non-conformal element boundaries, representing the cut as well as possible by a serrated discontinuity, and then moving the duplicated nodes apart to model the opening of the path [34], [35]. Retraction has also been modeled by deleting the FEs falling into the retraction path, as visualized on the intraoperative image [13]. This way of modeling retraction implies the removal of some tissues, which is not physically correct and does not model explicitly the movement of tissues as they are spread out. Resection, i.e., the removal of tissues, has been modeled by deleting FEs that were tagged before surgery [35], or that were falling within the resection cavity visible on the intraoperative image [13], [36].

Because XFEM allows an accurate representation of the discontinuities while avoiding mesh adaption or remeshing, and because of the similarity between cracks in mechanical parts and cuts in tissues, we have proposed, beginning with reference [37], the use of XFEM for handling cuts, retraction and resection in the updating of preoperative images. The present paper gives a complete and unified account of most of our work to date in this 2D context [37–40]. While we are actively working on the extension to three dimensions [41], the 2D case has the advantage of allowing us to address the main points of the application of XFEM to image-guided surgery without getting into the increased complexities and computation times of the 3D case.

The structure of the paper is as follows. In the next section, we introduce the basic principles of FEM and XFEM. We then describe our general system for updating preoperative images. After this, we consider three cases that illustrate our approach to handling brain shift [40], retraction [38, 39, 37] and successive resections [40]. Finally, we conclude and discuss future work.

Basic principles of FEM and XFEM

FEM discretizes the object into a mesh, i.e., into a set of FEs interconnected by nodes, and approximates the displacement field u(x) by the FEM displacement field u^FEM(x) defined as where I is the set of nodes, the φ_i(x)'s are the nodal shape functions (NSFs), and the u_i's are the nodal degrees of freedom (DOFs). Each φ_i(x) is defined as being continuous on its compact support ω_i, which corresponds to the union of the domains of the FEs connected to node i, with the property where x_j is the position of node j and δ_ij the usual Kronecker signal. In our approach, we use linear NSFs.

FEM requires its displacement field u^FEM(x) to be continuous over each FE. In contrast, XFEM handles a crack, or discontinuity, by allowing the displacement field to be discontinuous within FEs [29], [32], [42], [43]. Arbitrarily shaped cracks can then be modeled without any remeshing. The XFEM displacement field generalizes the FEM displacement field (Equation 1) with

The first term corresponds to the FEM displacement field (Equation 1), where I is the set of nodes, the φ_i(x)'s are the FEM NSFs, and the u_i's are the nodal FEM DOFs. The heart of XFEM is the “enrichment” that adds a number, n^E_i, of DOFs a_ji to each node i of the set J, which is the subset of nodes of I whose support is intersected by the crack of interest. These DOFs are multiplied by the NSFs φ_i(x) and the discontinuous functions g_j(x).

An important consequence of the XFEM enrichment is that, in contrast to FEM, the XFEM displacement field (Equation 3) does not interpolate the nodal displacements for the nodes that are enriched, i.e., for k ∈ J,

Equation 4 implies that nodal displacements cannot be imposed by eliminating the DOFs u_k of nodes that are enriched, but must be applied as a linear combination between several DOFs with, for instance, Lagrange multipliers.

As previously indicated, any node whose support is intersected by the crack is enriched. However, the choice of the function g_j(x) for a node i ∈ J varies according to whether the support of the node is fully or partially intersected by the crack. When the crack fully intersects the support of the node, a simple choice is a piecewise-constant function that changes sign at the boundary across the crack, i.e., the Heaviside function where x is the position of a point of the solid, x* is the position of the point on the crack that is the closest to x, and e_n is the outward normal to the crack at x*, as shown in Figure 1a [32].

Figure 1. (a) Parameters related to the Heaviside function corresponding to a crack (2D case). (b) Parameters related to the crack-tip functions (2D case). (c) 2D XFEM enrichment for a single pair of one crack and one crack tip. Each node enriched with Heaviside DOFs is represented by an open circle (set J), while each node enriched with crack-tip DOFs is represented by a filled square (set K).

When the crack ends within a node support, using H(x) is not sufficient. A possible enrichment relies on functions that incorporate the radial and angular behavior of the asymptotic linear-elastic crack-tip displacement field, which is two-dimensional by nature, i.e.,

where r and θ are the local polar coordinates of the position of a point x measured in the (e_s(x_tip), e_n(x_tip))-axes, as shown in Figure 1b [29]. Note that the first function in Equation 6 is discontinuous across the crack. If we used H(x) instead of the crack-tip functions F_l(x) = F_l(r, θ), l = 1, …, 4, the crack would effectively be modeled as if it extended to the node support boundary. The functions F_l(r, θ) ensure that the crack terminates precisely at the location of the crack tip.

As explained above, the discontinuous functions g_j(x) allow accurate modeling of the location of the crack, both along and at the crack tip. Moreover, they can be used to take into account the a priori knowledge of the exact solution. The crack-tip functions (Equation 6), for instance, specifically allow representation of the high gradient of stress near the crack tip by enriching the approximation with the known asymptotic behavior of the solution that can be expressed with the basis functions in Equation 6. Depending on the problem to be solved, one can choose other discontinuous functions for enrichment [44], as well as the method of selecting the nodes to be enriched [45].

Based on Equations 3, 5 and 6, the XFEM displacement field for a single pair of one crack and one crack tip is given by

where J now denotes the subset of nodes whose support is fully intersected by the crack and K the subset of nodes whose support is partially intersected by the crack. The a_j's and 's are the new notations for the corresponding DOFs, and H(x) and F_l(x) are the corresponding discontinuous functions. J and K can be formally defined by

where x_c is the position of the crack tip, ω_k is the compact support of node k, and ϖ_k its closure. Equation 7 can easily be generalized to several pairs of cracks and crack tips. Figure 1c illustrates the sets of enriched nodes for a 2D crack with one tip.

The support of any node in J is fully cut into two disjointed pieces by the crack. If, for a given node, one of the two pieces is very small compared to the other, then the Heaviside function is almost constant over the nodal support, which causes the XFEM problem to be ill-conditioned. Consequently, the node has to be removed from the set J, which is equivalent to moving the crack location to the boundary of the support. A criterion should be used to remove from the set J a node for which the ratio of the area of the smallest piece to that of the full support is less than a set tolerance value, e.g., ε = 10⁻⁴ [42].

When minimizing the total deformation energy, the resulting XFEM equations remain sparse and symmetric as for FEM. Whereas FEM requires a remeshing and the duplication of the nodes along the crack to take into account any discontinuity, XFEM requires the identification of the nodes belonging to the sets J and K and the addition of DOFs: (1) any node whose support is not intersected by the discontinuity remains unaffected and thus possesses two DOFs; (2) any node whose support is fully intersected by the discontinuity is enriched with two Heaviside DOFs and thus possesses four DOFs; and (3) any node whose support contains the tip of the discontinuity is enriched with eight crack-tip DOFs and thus possesses ten DOFs. (These numbers of DOFs are for the 2D case.) Because of the additional DOFs in XFEM, the number of equations to be solved is larger than for FEM.

Figures 2 and 3 illustrate two examples of 2D XFEM calculation (carried out in our MATLAB implementation of 2D XFEM). In each case, the first graph shows the meshed object and the position of the crack, and the second shows the result of the XFEM calculation. In Figure 2, the crack ends inside the object. Displacements corresponding to Mode I, i.e., the tensile mode, of fracture [46], are applied to nodes located along the boundary of the mesh. The result of the XFEM calculation indicates that the crack can open even though there was no remeshing; indeed, the FE edges along the crack have become discontinuous. In Figure 3, the crack fully crosses the object. Distinct translations are applied to each part (more specifically to two nodes of each part) of the mesh. The result of the XFEM calculation indicates that the two parts of the mesh can move independently, even though there was no remeshing [47].

Figure 2. First illustration of 2D XFEM calculation: the case of a rectangular object with a crack ending inside the object. (a) Original mesh and position of the crack. (b) Final mesh resulting from application of displacements corresponding to Mode I of fracture to nodes located along the boundary of the mesh.

Figure 3. Second illustration of 2D XFEM calculation: the case of a rectangular object with a crack fully crossing the object. (a) Original mesh and position of the crack. (b) Final mesh resulting from the application of distinct translations to each part (more specifically, to two nodes of each part) of the mesh.

General approach to the updating of preoperative images

The block diagram of Figure 4 shows our general approach to updating preoperative images using successive intraoperative MR (iMR) images acquired at different critical points during surgery. In our present strategy, the preoperative images are updated incrementally. At the end of each update, the preoperative images should be in the best possible alignment with the last iMR image acquired.

Prior to surgery, a patient-specific biomechanical model is constructed from the set of preoperative images. Because the patient does not necessarily lie in the same position during the acquisition of each preoperative image, one may need to perform a rigid registration (involving translations, rotations and scalings) to bring them into correspondence, assuming there is no local brain deformation between preoperative images. As indicated in Figure 4, we assume here that the preoperative images have been registered beforehand. The registered preoperative images are segmented into various regions, e.g., corresponding to cortex, ventricles, tumor, vessels, etc. These regions are then transformed into a mesh of FEs, and possibly distinct behavior laws are assigned to these regions.

Figure 4. Block diagram of our preoperative image-update system based on biomechanical modeling of the brain.

Once the first iMR image has been acquired, we can commence the update of the preoperative images. Ideally, the first iMR image will be acquired prior to the opening of the skull. In the present situation, it is assumed that the brain has the same shape in the first iMR image as in the preoperative images. The set of registered preoperative images and the biomechanical model are registered to the first iMR image via a rigid transformation. The set of registered preoperative images are then transformed accordingly.

As each iMR image is acquired, this new image and the preceding iMR image are used to estimate the deformation of the brain. The update of the preoperative images is done incrementally with each new pair of successive iMR images. For each pair, we proceed as follows. A set of common anatomical landmarks are segmented in each image of the pair, and the corresponding landmarks are then matched. In our approach, we match the surfaces of key structures, such as cortex and tumor. Using surface structures rather than volume structures [10] seems more reliable, given the limited quality of typical intraoperative images, and would be more easily adapted to other intraoperative modalities such as ultrasound, for instance. The landmark surface displacements resulting from the matching are applied to the biomechanical model, which is deformed using FEM or XFEM, depending on the type of deformation occurring between the pair of iMR images under consideration, e.g., brain shift, retraction or resection. The resulting volume displacement field of the biomechanical model is used to warp the set of preoperative images in their current state of updating. This process is repeated with each new acquisition of an iMR image.

The initial stress state of the brain is set to zero for each FEM or XFEM computation because it is physically unknown. A plane state of strain is assumed for each FEM or XFEM computation, and an homogeneous linear elastic behavior law is assigned to the mesh. Because displacements, rather than forces, are applied to the model, the FEM or XFEM solution is independent of Young's modulus E [48]. The Poisson ratio ν is set to 0.45 [49].

Below, we illustrate the general strategy just described by considering three particular cases that cover the main sources of deformations: brain shift, retraction and successive resections. In these cases, we make several simplifying assumptions. First, we deal with the problem of updating the first iMR image, which thus plays the role of a substitute preoperative image. This allows us to validate the method easily by comparing the first iMR image, successively deformed, and the subsequent iMR images. These images are of the same modality and quality, and thus show equivalent anatomical features. Except for the rigid registration between the preoperative images, the biomechanical model and the first iMR image, all key aspects of the system are discussed and illustrated. The biomechanical model is thus constructed on the basis of structures segmented in the first iMR image, rather than in the preoperative images, and the structures used in the model are limited to those visible in the intraoperative image. Second, we treat the whole problem in two dimensions. This means that we select a particular 2D slice in the iMR images. The surfaces and volumes of the ultimate 3D formulation respectively become contours and surfaces in the present 2D formulation.

To evaluate our results, we use the modified Hausdorff distance [50] that computes the distance between two sets of points and represents the accuracy of the alignment. For each update, we compare the modified Hausdorff distance between edges extracted with a Canny edge detector from (1) the pair of successive iMR images that are not registered and (2) the first iMR image in its current state of update and the last iMR image acquired. The difference between the two values of the modified Hausdorff distance represents the improvement in the alignment computed with our approach.

Results

For the three cases discussed below, iMR images were acquired with a 0.5-Tesla intraoperative GE Signa scanner. All computation was done offline. We segmented 2D images manually into regions using 3D SLICER (www.slicer.org), we created triangular meshes of these regions using DISTMESH (www-math.mit.edu/~persson/mesh/), and we deformed biomechanical models using our own MATLAB implementation of 2D FEM or 2D XFEM.

Case 1: Deformation due to brain shift

Case 1 is concerned with the brain shift resulting from the opening of the skull and dura. The 2D iMR images acquired before and after the opening are shown in Figures 5(1) and 5(2). iMR image size is 256 × 256 pixels and pixel size is 0.9375 × 0.9375 mm.

Figure 5. Two-dimensional iMR images for Case 1 (brain shift) and segmentation of the first iMR image. (1) The first iMR image, acquired before the opening of the skull and dura, and also used as a substitute preoperative image to be updated. (2) The second iMR image, acquired after the opening of the skull and dura, but prior to further surgical actions. The deformation observed is thus due to brain shift alone. (3) Segmentation of (1) into two regions: the healthy brain and the tumor.

In the present case, no tissue discontinuity is involved in the deformation, so the biomechanical model is deformed using FEM. From the first iMR image (which serves as a substitute preoperative image), we segment out the healthy brain and tumor regions, as illustrated in Figure 5(3). We mesh the whole brain, then choose two contour landmarks (described below) and compute their displacements, which are later used to drive the deformation of the biomechanical model. The first contour landmark is the external boundary of the brain, i.e., the cortex, defined by the set of boundary nodes of the mesh. The second is the internal tumor boundary defined by computing the intersections of the mesh with the boundary between the healthy brain and the tumor, where this boundary is derived from the segmentation. We compute the displacement field of each of these two contours as the closest point of the whole-brain and tumor regions segmented out in the second iMR image. Both contour displacement fields are smoothed to prevent numerical problems.

The resulting surface displacement field (recall that we are working in two dimensions, where objects are planar) of the biomechanical model is used to warp the part of the first iMR image (serving as a substitute preoperative image) [Figure 6(1a)] corresponding to the brain, i.e., with the skull and external cerebrospinal fluid masked out [Figure 6(1b)]. The result of the warping of Figure 6(1b) is shown in Figure 6(2b). To evaluate the quality of the registration, the Canny edge detector is applied to Figures 6(1b) and 6(2b), and to the brain part of Figure 6(2a). Two modified Hausdorff distances are computed between the edges of (1) Figure 6(1b) and the brain part of Figure 6(2a), shown jointly in Figure 6(1c); and (2) Figure 6(2b) and the brain part of Figure 6(2a), shown jointly in Figure 6(2c). The value of the modified Hausdorff distance goes from 1.62 mm to 1.18 mm between Figures 6(1c) and 6(2c), which confirms the improvement in the alignment.

Figure 6. Two-dimensional experimental results for Case 1 (brain shift). (1a) The first iMR image. (2a) The second iMR image. (1b) Brain extracted from (1a). (2b) Deformation of (1b) using the surface displacement field of the biomechanical model, computed via FEM. (We say “surface” as we are working in two dimensions.) (1c) Juxtaposition of Canny edges of (1b) and the brain part of (2a). (2c) Juxtaposition of Canny edges of (2b) and the brain part of (2a).

Case 2: Deformation due to brain shift and retraction

Case 2 is concerned with the brain shift resulting from the opening of the skull and dura, and with a subsequent retraction. The 2D iMR images acquired before the opening of the skull and dura, and after some retraction, are shown in Figures 7(1) and 7(2). iMR image size is 256 × 256 pixels, and pixel size is 0.8593 × 0.8593 mm. The beginning of resection that is also visible in the second iMR image will be ignored. Ideally, we would have preferred to handle the brain shift first and then the retraction; however, this is not possible since an iMR image corresponding to the brain shift alone was not available. We thus handle both effects simultaneously. (The possibility of doing so shows the flexibility of our general approach.)

Figure 7. Two-dimensional iMR images for Case 2 (brain shift and retraction). (1) The first iMR image, acquired before the opening of the skull and dura, and also used as a substitute preoperative image to be updated. (2) The second iMR image, acquired after the opening of the skull and dura, the retraction, and the beginning of resection. The deformation observed is principally due to brain shift and retraction.

The main difference between Case 1 and Case 2 is the presence in the latter of a discontinuity due to the cut required to retract tissues. To build the biomechanical model, we proceed as follows. First, as in Case 1, the brain region is segmented out in the first iMR image (which again serves as a substitute preoperative image). We estimate the position of the cut on the first iMR image with the help of a neurosurgeon. For simplicity, the cut is assumed to be linear; its estimated position is shown in Figure 8a. If we were to use FEM, we would mesh the brain region in such a way as to align the tissue discontinuity (due to the cut) with FE edges, and would duplicate the nodes lying on these edges. Here, however, we use XFEM rather than FEM. Thus, we start by meshing the brain region into FEs without taking the cut into account. We then compute the intersections of the mesh with the cut (Figure 8b). The nodes are enriched with XFEM DOFs in accordance with Equation 7: (a) any node whose support is not intersected by the discontinuity remains unaffected; (b) any node whose support is fully intersected by the discontinuity is enriched with two Heaviside DOFs; (c) any node whose support contains the tip of the discontinuity is enriched with eight crack-tip DOFs. The mesh in Figure 8b consists of 828 nodes, which corresponds to 1,656 FEM DOFs. Enrichment adds 30 XFEM Heaviside DOFs and 24 crack-tip DOFs, leading to a total of 1,710 DOFs.

Figure 8. (a) The first iMR image with our estimate of the location of the cut. (b) The mesh created from the brain region of the first iMR image and cut; the cut defines intersections with the mesh.

We choose two contour landmarks (whose displacement fields will drive the deformation of the biomechanical model): the cortex (as in Case 1) and the cut (along which the blades of the retractor are introduced). The cortex contour is defined by the set of external nodes of the mesh, as in Case 1, except for the nodes that are located close to the cut. Indeed, we do not want to use those nodes whose displacements would bring them into the hole created by the retraction. The displacements of the selected cortex nodes are obtained in the same way as in Case 1. The cut is defined by the set of cut intersections defined above. To estimate the displacement of the cut intersections, we assume that the tissues on the left and right sides of the cut were displaced perpendicular to the line of the cut. We estimate the displacement of each cut intersection using the second iMR image; one such displacement is shown in Figure 9. Both contour displacement fields are smoothed to prevent numerical problems. In particular, the displacement applied at the intersection of the element containing the tip should not be too large to avoid element flipping.

Figure 9. The second iMR image with the cut and, for one particular intersection of the cut with the mesh, the corresponding displacement due to retraction.

The displacements and obtained earlier for the left and right lips of each cut intersection n are applied by imposing, based on Equation 7, the conditions

for the left lip, where the Heaviside function H(x) (see Equation 5) is equal to +1, and, similarly, the conditions

for the right lip, where the Heaviside function H(x) is equal to −1. (The assignments of the values +1 and −1 to the left and right lips can be exchanged.) In Equations 9 and 10, is the position of cut intersection n, A_n and B_n are the two nodes of the FE edge on which cut intersection n is located, φ_Am(x) and φ_Bm(x) are the FEM NSFs of nodes A_n and B_n, u_An and u_Bn are the nodal FEM DOFs of nodes A_n and B_n, and a_An and a_Bn are the nodal XFEM Heaviside DOFs of nodes A_n and B_n. Conditions 9 and 10 are applied with Lagrange multipliers. If nodes A_n and/or B_n are enriched with crack-tip DOFs rather than with Heaviside DOFs, the displacements are applied in a similar way.

Figures 10a and 10b show, respectively, the initial mesh and the final mesh resulting from the XFEM calculation. The surface displacement field of the biomechanical model is used to warp the part of the first iMR image (serving as a substitute preoperative image) [Figure 11(1a)] corresponding to the brain, i.e., with the skull and external cerebrospinal fluid masked out [Figure 11(1b)]. The result of the warping of Figure 11(1b) is shown in Figure 11(2b). To evaluate the quality of the registration, the Canny edge detector is applied to Figures 11(1b) and 11(2b) and to the brain part of Figure 11(2a). Two modified Hausdorff distances are computed between edges of (1) Figure 11(1b) and the brain part of Figure 11(2a), shown jointly in Figure 11(1c); and (2) Figure 11(2b) and the brain part of Figure 11(2a), shown jointly in Figure 11(2c). The value of the modified Hausdorff distance goes from 2.17 mm to 1.78 mm between Figures 11(1c) and 11(2c), which confirms the improvement in the alignment.

Figure 10. (a) The initial mesh with the cut. (b) The final mesh. The deformation is the result of applying the displacements due to brain shift and retraction to the biomechanical model using XFEM.

Figure 11. Two-dimensional experimental results for Case 2 (brain shift and retraction). (1a) The first iMR image. (2a) The second iMR image. (1b) Brain extracted from (1a). (2b) Deformation of (1b) using the surface displacement field of the biomechanical model, computed via XFEM. (We say “surface” as we are working in two dimensions.) (1c) Juxtaposition of Canny edges of (1b) and the brain part of (2a). (2c) Juxtaposition of Canny edges of (2b) and the brain part of (2a).

Case 3: Deformation due to brain shift and successive resections

Case 3 is concerned with the brain shift resulting from the opening of the skull and dura, and with a subsequent series of successive resections. (No retraction is involved.) The 2D iMR images acquired before the opening of the skull and dura, and after successive resections, are shown in Figure 12. iMR image size is 256 × 256 pixels and pixel size is 0.9375 × 0.9375 mm.

Figure 12. Two-dimensional iMR images for Case 3 (brain shift and successive resections). (1) The first iMR image, acquired before the opening of the skull and dura, and also used as a substitute preoperative image to be updated. (2) The second iMR image, acquired after the opening of the skull and dura, and thus after some brain shift has occurred. (3) The third iMR image, acquired after a first resection. (4) The fourth iMR image, acquired after a second resection. (5) The fifth iMR image, acquired after a third resection and the postoperative brain shift.

Brain shift was already handled in the previous section, based on the two iMR images of Figures 12(1) and 12(2). (Cases 1 and 3 correspond to the same surgical case.) This resulted in (1) a biomechanical model deformed in accordance with the brain shift; and (2) the part of the first iMR image (serving as a substitute preoperative image) [Figure 12(1)] corresponding to the brain, warped to be in registration with Figure 12(2)]. Both elements will be used as a starting point to model the successive resections.

Recall that, for brain shift deformation, we estimated two contour displacement fields using the first and second iMR images, one for the cortex, and one for the internal tumor boundary. Matching two contours to get a displacement field makes sense only if the contours correspond to the same physical entity. For deformations involving some resection, we cannot rely on the totality of the cortex, since a part of it is now missing. Consequently, for modeling resection, we evaluate the displacement field for the contour resulting from the combination of the intact cortex and the internal tumor boundary that constitutes the contour of the healthy brain.

The first, second and third resections are modeled using different techniques, as detailed in the following three subsections.

Modeling of first resection.

Based on the second and third iMR images, we cannot determine the amount and location of the tissues removed by the first resection. This is because the third iMR image shows the combined effect of tissue removal and subsequent deformation. In fact, because the second and third iMR images do not show the same amount of brain tissue, the problems of modeling resection and modeling brain shift are intrinsically different. Nevertheless, we decided to model the first resection by still relying on the displacement field of key contour landmarks, in this case just the healthy brain contour, to deform the biomechanical model. Indeed, these fields appear to be the only reliable source of information that we can extract from the second and third iMR images concerning the deformation caused by this first resection. Consequently, we do not explicitly model the removal of tissues, but do model directly the deformation resulting from it, without introducing any tissue discontinuity. The biomechanical model is deformed using FEM, in accordance with the contour displacement field of the healthy brain.

The surface displacement field of the biomechanical model is used to warp the part of the first iMR image (serving as a substitute preoperative image) corresponding to the brain in its current state of update [Figure 13(1b)]. Figure 13(1b) should be found to be in registration with the second iMR image [Figure 13(1a)] as a result of the brain-shift modeling. The result of the warping of Figure 13(1b) is shown in Figure 13(2b), which is now registered to the third iMR image, except outside the healthy brain boundary, i.e., except for the tumor. Finally, we alter the result of warping to reflect the effect of resection. For this, we assign the background color (here black) to the pixels of Figure 13(2b) corresponding to the resected tissues that are “absent” in the third iMR image. The result of the warping with resection, shown in Figure 13(2c), is performed by masking the warped image [Figure 13(2b)] with the brain region segmented out from the third iMR image [Figure 13(2a)]. To evaluate the quality of the registration, the Canny edge detector is applied to Figures 13(1b) and 13(2c), and to the brain part of Figure 13(2a). Two modified Hausdorff distances are computed between edges of (1) Figure 13(1b) and the brain part of Figure 13(2a), shown jointly in Figure 13(1d); and (2) Figure 13(2c) and the brain part of Figure 13(2a), shown jointly in Figure 13(2d). The value of the modified Hausdorff distance goes from 2.07 mm to 1.18 mm between Figures 13(1d) and 13(2d), which confirms the improvement of the alignment.

Figure 13. Two-dimensional experimental results for modeling of the first resection in Case 3. (1a) The second iMR image. (2a) The third iMR image. (1b) The image resulting from brain-shift modeling [identical to Figure 6(2b)]. (2b) Deformation of (1b) using the surface displacement field of the biomechanical model, computed via FEM. (We say “surface” as we are working in two dimensions.) (2c) Deformation of (1b) with resection, performed by masking (2b) with the brain region segmented from the third iMR image (2a). (1d) Juxtaposition of Canny edges of (1b) and the brain part of (2a). (2d) Juxtaposition of Canny edges of (2c) and the brain part of (2a).

Modeling of second resection.

The significant feature of the second resection is that some tissues have already been removed by the first resection, which means that these tissues can have no physical influence on subsequent brain deformations because they no longer “exist”. Consequently, the first resection must be reflected in the biomechanical model. Recall that, to model the first resection, the biomechanical model was deformed to be registered to the third iMR image. The boundary of the first resection, i.e., the tissue discontinuity to be included in the biomechanical model, can then be defined using the third iMR image (Figure 14a). With an FEM-based calculation of deformation, we would remesh the biomechanical model to take into account the discontinuity. Then, we would just remove the part of the mesh corresponding to the resected tissues, and use the other part of the mesh to model the second resection. Instead, with XFEM, the nodes whose support is intersected by the discontinuity are enriched with XFEM Heaviside DOFs (Figure 14b). (Since the resection boundary fully crosses the brain, there is no need to use crack-tip DOFs in the enrichment of these nodes.) Once the discontinuity splitting the mesh in two parts is modeled with XFEM, the two parts of the mesh are free to deform independently. Consequently, when the biomechanical model is deformed using XFEM, the part corresponding to the tissues removed by the first resection has no influence on the deformation of the remaining part of the brain.

Figure 14. (a) The final mesh from the modeling of the first resection superimposed on the healthy brain and tumor regions segmented out from the third iMR image. This superimposition allows us to define the tissue discontinuity due to the first resection. (b) Illustration of the XFEM enrichment process for modeling the tissue discontinuity due to the first resection. XFEM Heaviside DOFs are added to each node whose support is intersected by the discontinuity.

Except for the fact that the biomechanical model is deformed with XFEM rather than FEM, the modeling of the second resection is identical to that of the first resection. The bottom part of the mesh, representing the tissues remaining after the first resection, is deformed according to the contour displacement field of the healthy brain derived from the third and fourth iMR images, while the top part, representing the piece of tissues removed by the first resection, is subjected to a translation, but only for visualization purposes. (The two parts of the mesh could indeed overlap around the discontinuity after deformation of the bottom part.) Figure 15 shows the final mesh resulting from the XFEM calculation.

Figure 15. The final mesh resulting from the modeling of the second resection using XFEM.

The surface displacement field of the biomechanical model is used to warp the part of the first iMR image (serving as a substitute preoperative image) corresponding to the brain in its current state of update [Figures 16(1b) and 16(1c)]. Figures 16(1b) and 16(1c) should be found to be in registration with the third iMR image [Figure 16(1a)] as a result of the brain shift and first resection modelings. The result of the warping of Figure 16(1b) is shown in Figure 16(2b). The result of the warping with resection, shown in Figure 16(2c), is achieved by masking the warped image [Figure 16(2b)] with the brain region segmented out from the fourth iMR image [Figure 16(2a)]. To evaluate the quality of the registration, the Canny edge detector is applied to Figures 16(1c) and 16(2c), and to the brain part of Figure 16(2a). Two modified Hausdorff distances are computed between edges of (1) Figure 16(1c) and the brain part of Figure 16(2a), shown jointly in Figure 16(1d); and (2) Figure 16(2c) and the brain part of Figure 16(2a), shown jointly in Figure 16(2d). The value of the modified Hausdorff distance goes from 1.72 mm to 1.21 mm between Figures 16(1d) and 16(2d), which confirms the improvement in the alignment.

Figure 16. Two-dimensional experimental results for modeling of the second resection in Case 3. (1a) The third iMR image. (2a) The fourth iMR image. (1b) The image resulting from brain shift and first resection modelings, identical to Figure 13(2b). (2b) Deformation of (1b) using the surface displacement field of the biomechanical model, computed via XFEM. (We say “surface” since we are working in two dimensions.) (1c) The image resulting from brain shift and first resection modelings, identical to Figure 13(2c). (2c) Deformation of (1b) with resection, performed by masking (2b) with the brain region segmented from the fourth iMR image (2a). (1d) Juxtaposition of Canny edges of (1c) and the brain part of (2a). (2d) Juxtaposition of Canny edges of (2c) and the brain part of (2a).

Modeling of third resection

If we had complete freedom, we would model the third resection and postoperative brain shift based on the bottom part of the output mesh resulting from the modeling of the second resection (Figure 15). However, even though the mesh is displayed as two separate parts, it is, in fact, a single entity. Indeed, a main feature of XFEM is its ability to handle the effect of a discontinuity without modifying the underlying mesh, i.e., without remeshing. For modeling the second resection, the edges of FEs straddling the discontinuity have been “cut” and their nodes moved apart. However, it is not possible to start an XFEM calculation with a mesh already deformed because the FEs are not supposed to be already cut. Consequently, XFEM, expressed in a linear formulation, does not permit us to proceed with a new discontinuity from a configuration already deformed by XFEM. Indeed, a nonlinear formulation of XFEM [51], [52] should be used to proceed as desired.

To temporarily avoid going nonlinear, we devised a new method that continues to preserve the spirit of XFEM, which is mainly to avoid remeshing in the presence of a discontinuity. However, this method implies additional processing. The idea of our method is to “reconnect the pieces” of the deformed mesh (Figure 15), but without remeshing, in such a way that it can be further resected. Since the deformed positions of the nodes are correct for the current bottom part of the mesh, we only need to find a way to “restore” a meaningful top part. This current top part cannot simply be lowered and reconnected with the current bottom part, because the two parts may not fit. Indeed, the bottom part has been deformed during the modeling of the second resection, with the consequence that the discontinuity boundary of the bottom part has changed shape. While the following might not be the best from a computational standpoint, our current solution is to use, for the position of the nodes of the top part of the mesh, the position that we would have obtained if we had modeled the second resection based on FEM, rather than on XFEM, i.e., without taking into account the removal of tissues by the first resection. These positions are erroneous, but this is not an issue, since these nodes will again be resected, given that the third resection is necessarily deeper than the second. Figure 17a shows the final mesh resulting from the modeling of the second resection using XFEM (Figure 15), while Figure 17b shows the final mesh resulting from the modeling of the second resection using FEM. The mesh resulting from the “reconnection” is shown in Figure 17c.

Figure 17. (a) The final mesh resulting from the modeling of the second resection using XFEM, as shown in Figure 15. (b) The final mesh resulting from the modeling of the second resection using FEM. (c) The reconnected mesh used for modeling the third resection. The positions of the nodes in (c) corresponding to the bottom part of mesh (a) come from mesh (a), while the positions of the nodes corresponding to the top part of mesh (a) come from mesh (b).

After the reconnection of the mesh, the modeling of the third resection is identical to that of the second resection. The tissue discontinuity due to the second resection is defined from the fourth iMR image, and is used to appropriately enrich the nodes of the reconnected mesh. The biomechanical model is deformed using XFEM, in accordance with the contour displacement field of the healthy brain computed from the fourth and fifth iMR images. One significant feature of the procedure described for modeling the third resection, i.e., mesh reconnection followed by model deformation, is that it can be applied iteratively for each subsequent resection visible on successive iMR images, no matter how many there are.

In the particular case considered here, a simplification can be made in the modeling of the third resection because by the time the fifth iMR image is acquired the resection is complete (assuming that the surgeon has not removed tissues located further than the internal tumor boundary). This means that only the surface displacement field corresponding to the healthy brain tissues needs to be computed. Since displacements are applied exactly to the contour of the healthy brain, the results would be the same with FEM and XFEM. For the present case, we choose to use FEM. The surface displacement field of the biomechanical model is used to warp the part of the first iMR image (serving as a substitute preoperative image) corresponding to the brain in its current state of update [(Figures 18(1b) and 18(1c)]. Figures 18(1b) and (1c) should be found in registration with the fourth iMR image [Figure 18(1a)] as a result of the brain shift, first resection and second resection modelings. The result of the warping of Figure 18(1b) is shown in Figure 18(2b). The result of the warping with resection, shown in Figure 18(2c), is achieved by masking the warped image [Figure 18(2b)] with the brain region segmented out from the fifth iMR image (Figure 18(2a)]. To evaluate the quality of the registration, the Canny edge detector is applied to Figures 18(1c) and 18(2c), and to the brain part of Figure 18(2a). Two modified Hausdorff distances are computed between edges of (1) Figure 18(1c) and the brain part of Figure 18(2a), shown jointly in Figure 18(1d); and (2) Figure 18(2c) and the brain part of Figure 18(2a), shown jointly in Figure 18(2d). The value of the modified Hausdorff distance goes from 1.78 mm to 1.61 mm between Figures 18(1d) and 18(2d), which confirms the improvement in the alignment.

Figure 18. Two-dimensional experimental results for modeling of the third resection in Case 3. (1a) The fourth iMR image. (2a) The fifth iMR image. (1b) The image resulting from brain shift, first resection, and second resection modelings, identical to Figure 16(2b). (2b) Deformation of (1b) using the surface displacement field of the biomechanical model, computed via FEM (using XFEM is also possible). (We say “surface” since we are working in two dimensions.) (1c) The image resulting from brain shift, first resection, and second resection modelings, identical to Figure 16(2c). (2c) Deformation of (1b) with resection, performed by masking (2b) with the brain region segmented from the fifth iMR image (2a). (1d) Juxtaposition of Canny edges of (1c) and the brain part of (2a). (2d) Juxtaposition of Canny edges of (2c) and the brain part of (2a).

Conclusions

We have developed a flexible approach for updating a set of preoperative images by using a series of limited-quality intraoperative images acquired at different critical points during surgery. The key elements of the approach are the brain model constructed from preoperative images, and the use of solid-mechanics techniques to evaluate the deformation of the brain based on the displacements of a few landmarks tracked through successive iMR images. The landmarks correspond to surface structures, such as the cortex or the tumor boundary, extracted from the intraoperative images. Using surfaces rather than volume information seems more reliable for limited-quality images. The iMR images used in our experiments were acquired with a 0.5-Tesla intraoperative GE Signa scanner. To facilitate the evaluation of our results, we have updated the first iMR image, assumed to be acquired before the opening of the skull and dura, as a substitute preoperative image. We have also dealt with a 2D version of the real 3D problem. All computation was done offline.

We have focused on deformations due to cutting, retraction and resection, all of them implying one or more tissue discontinuities. Not surprisingly, efficient techniques for dealing with such deformations can be found in the field of fracture mechanics, where XFEM began to be developed in 1999 to study cracks. Because of the similarity between cracks in mechanical engineering and cuts in surgery, we have previously proposed [37–41] the use of XFEM to deal with cuts, retractions and resections in image-guided surgery. The present paper describes and illustrates possible strategies for using both FEM and XFEM (as appropriate) to deal with brain shift, retraction and successive resections. One of the important lessons learned is that, while a linear formulation of XFEM can easily deal with a one-time retraction or resection, the same formulation cannot deal with successive discontinuities. This paper nevertheless shows an effective method for going around this problem. For all cases, the results show an improvement in the alignment computed with our approach. Since we handle in two dimensions deformations that actually occur in three dimensions, we expect to obtain better results when modeling the deformations directly in three dimensions. Nevertheless, these results give us confidence that the proposed approach works correctly. A nonlinear formulation of XFEM should provide a more systematic approach for dealing with successive resections and for handling combined retraction and resection.

Future research will thus be directed at the use of nonlinear XFEM, the joint handling of retraction and resection, the generalization of all techniques from two dimensions to three dimensions, the implementation of an end-to-end software system for real-time preoperative image updating in the operating room, and the validation of these techniques.

Acknowledgments

This investigation was supported in part by the Fonds Léon Fredericq, Liège, Belgium; the Fonds National de la Recherche Scientifique (F.N.R.S.), Brussels, Belgium; the Fonds d'Investissements de la Recherche Scientifique (F.I.R.S.) from the CHU of Liège, Belgium; NSF grant ITR 0426558; a research grant from CIMIT; and by NIH grants R03 CA126466, R01 RR021885, R01 GM074068 and R01 EB008015. P.A. Robe is a Research Associate at the F.N.R.S. of Belgium.