Advanced Nonrigid Registration Algorithms for Image Fusion

300 downloads 62122 Views 2MB Size Report
Nonrigid brain registration techniques have numerous applications. ... In the fifth section, we present a novel method for producing area preserving surface de- .... that best describes T in terms of S. In a maximum likelihood context, this function ...... domain so that the total potential energy can be written as a sum of potential ...
Advanced Nonrigid Registration Algorithms for Image Fusion Simon K. Warfield Aditya Bharatha Juan Ruiz-Alzola Sigurd Angenent

Alexandre Guimond Alida Tei

Alexis Roche

Florin Talos

Carl-Fredrik Westin Allen Tannenbaum

Jan Rexilius Steven Haker

Ferenc A. Jolesz

Ron Kikinis

1

Introduction Medical images are brought into spatial correspondence, or aligned, by the use of registration algorithms. Nonrigid registration refers to the set of techniques that allow the alignment of datasets that are mismatched in a nonrigid, or nonuniform manner. Such misalignments can result from physical deformation processes, or can be a result of morphological variability. For example, physical deformation in the brain can occur during neurosurgery as a result of such factors as swelling, cerebrospinal fluid (CSF) loss, hemorrhage and the intervention itself. Nonrigid deformation is also characteristic of the organs and soft tissues of the abdomen and pelvis. 1 This chapter appeared in Brain Mapping: The Methods, Second Edition, as chapter 24 on pages 661–690, published by Academic Press of San Diego, USA in 2002.

In addition, nonrigid morphological differences can arise when comparisons are made among image datasets acquired from different individuals. These changes can be a result of normal anatomical variability or the product of pathological processes. Because the gross structure of the brain is essentially similar among humans (and even among related species), the factors described above tend to produce local nonrigid shape differences. Nonrigid brain registration techniques have numerous applications. They have been used to align scans of different brains, permitting the characterization of normal and pathological morphological variation (brain mapping). They have also been used to align anatomical templates with specific datasets, thus facilitating segmentation (i.e. segmentation by registration). More recently, these techniques have been used to capture changes which occur during neurosurgery. With the ongoing development of robust algorithms and advanced hardware platforms, further applications in surgical visualization and enhanced functional image analysis are inevitable. One exciting application of nonrigid registration algorithms is in the automatic registration of multimodal image data. Rigid registration of multimodal data has been greatly facilitated by the framework provided by mutual information (MI). However, MI-based strategies to effectively capture large nonrigid shape differences are still being explored. An alternate approach is to normalize multimodality images and thus reduce the problem to a monomodality match. In the first section, we present a nonrigid registration method which uses an intensity transform which allows a single intensity in one modality to be mapped onto (up to) two intensities. The method is iterative, combining in each iteration an intensity correction and a geometric transform using intensity-similarity criteria. The method is applied in two cases with promising results. In the next section, we turn our attention to the issue of image registration and fusion during neurosurgery. It is common to desire to align preoperative data with images of the patient acquired during neurosurgery. It is now widely acknowledged that during neurosurgical opera-

2

tions, nonrigid changes in the shape of the brain occur as a result of the intervention itself and due to reactive physiological changes. These deformations (“brain shift”) make it difficult to relate preoperative image data to the intraoperative anatomy of the patient. Since preoperative imaging is not subject to the same time constraints and limitations in tissue contrast selection methods as intraoperative imaging, a major goal has been to develop robust nonrigid registration algorithms for matching of preoperative image data onto intraoperative image data. We present our biomechanical modeling algorithm which can capture nonrigid deformations based on surface changes and infer volumetric deformation using a finite element discretization. We also describe our early prospective experience using the method during neurosurgical cases, and provide examples of the enhanced visualizations which are produced. In the third section, we build upon the theme of physics-based models by presenting a novel inhomogeneous elasticity model which uses a local similarity measure to obtain an initial sparse estimate of the deformation field. The method includes automatic feature point extraction using a nonlinear diffusion filter. Correspondence detection is achieved by maximizing a local normalized cross-correlation. The sparse estimates of the deformation field calculated at the feature points are then introduced as external forces, restricting the registration process so that the deformation field is fixed at those points. An advantage of the method is that feature points and correspondences are established automatically. Thus neither segmentation nor the manual identification of correspondences is required. In the fourth section we discuss registration of Diffusion Tensor MRI data and introduce a framework for nonrigid registration of tensor data (including the special case of vector data). The approach is based on a multiresolution template matching scheme followed by interpolation of the sparse displacement field using a Kriging interpolator. After warping the data, the tensors are locally realigned based on information from the deformation gradient of the displacement.

3

In the fifth section, we present a novel method for producing area preserving surface deformations, and more general mass preserving area and volume deformations, based on the minimization of a functional of Monge–Kantorovich type. The theory is based around the optimal mass transport problem of minimizing the cost of redistributing a certain amount of mass between two distributions given a priori. Here the cost is a function of the distance each bit of material is moved, weighted by its mass. The problem of optimal transport is classical and has appeared in econometrics, fluid dynamics, automatic control, transportation, statistical physics, shape optimization, expert systems, and meteorology. We show how the resulting low-order differential equations may be used for image registration. The challenge of nonrigid registration remains one of the outstanding open problems in medical image analysis. New algorithm developments, often targeted toward specific clinical applications, have helped to identify further unsolved issues. This chapter provides an overview of the nonrigid registration algorithms being pursued today at the Surgical Planning Laboratory.

4

1 Inter-Modality and Multi-Contrast Images 1.1

Introduction

Automatic registration techniques of brain images have been developed following two main trends: 1) registration of multimodal images using low to intermediate degree transformations (less than a few hundred parameters), and 2) registration of monomodal images using high-dimensional volumetric maps (elastic or fluid deformations with hundreds of thousands of parameters). The first category mainly addresses the fusion of complementary information obtained from different imaging modalities. The second category’s predominant purpose is the evaluation of either the anatomical evolution process present in a particular subject or of anatomical variations between different subjects. Despite promising early work such as (Hata, 1998), dense transformation field multimodal registration has, so far, remained relatively unexplored. Research on multimodal registration culminated with the concept of mutual information (MI) (Viola and Wells, 1995; Wells et al., 1996b; Hata et al., 1996; Viola and Wells, 1997; Maes et al., 1997), leading to a new class of rigid/affine registration algorithms. In this framework, the registration of two images is performed by maximizing their MI with respect to the transformation space. A significant reason for the success of MI as a similarity measure resides in its generality, as it does not use any prior information about the relationship between the intensities of the images. For instance, MI does not assume a linear relationship as is typically the case in standard optical flow techniques. Also, unlike some earlier approaches, MI does not require the identification of corresponding features in the images to be registered. Significant work has been done in establishing the applicability of MI for nonrigid registration (Gaens et al., 1998; Maintz et al., 1998; Meyer et al., 1999; Likar and Pernus, 2000; Hellier and Barillot, 2000; Rueckert et al., 2000; Hermosillo et al., 2001). Some authors have

5

further improved the robustness of the approach by modifying the original MI measure, either by including some prior information on the joint intensity distribution (Maintz et al., 1998; Likar and Pernus, 2000), or by using higher-order definitions of MI which incorporate spatial information (Rueckert et al., 2000). Our approach described here stems from the observation that a number of multimodal rigid registration problems can be solved in practice using other similarity measures than MI, one of which is the correlation ratio (CR) (Roche et al., 1998). The CR is much more constrained than MI as it assumes a functional, though non-linear, relationship between the image intensities. In other words, it assumes that one image could be made similar to the other by a simple intensity remapping. Thus, the CR method amounts to an adaptive estimation strategy where one image is alternately corrected in intensity and in geometry to progressively match the other. For most combinations of medical images, the functional dependence assumption is generally valid for a majority of anatomical structures, but not for all of them. Although this problem does not turn out to be critical in a rigid/affine registration context, we observe that it may seriously hamper the estimation of a high-dimensional transformation. We propose here an extension of the functional dependence model, which we call the bifunctional model, to achieve better intensity corrections. While the bifunctional model is more realistic than the functional one, it remains strongly constrained and thus enables a good conditioning of the multimodal non-rigid registration problem.

1.2

Method

The registration algorithm described here is iterative and each iteration consists of two parts. The first part transforms the intensities of anatomical structures of a source image S so that they match the intensities of the corresponding structures of a target image T . The second part is concerned with the registration of S (after intensity transformation) with T using an optical

6

flow algorithm.

1.2.1

Intensity Transformation

The intensity correction process starts by defining the set C of intensity pairs from corresponding voxels of T and S. Hence, the set C is defined as

 



S x  T x  ; 1  x  N 

C



(1)



where N is the number of voxels in the images. S x  and T x  correspond to the intensity value of the xth voxel of S and T , respectively, when adopting the customary convention of considering images as one-dimensional arrays. We shall now show how to perform intensity correction if we can assume that a single intensity value in S has either 1) exactly one corresponding intensity value in T (monofunctional dependence) or 2) at least one and at most two corresponding intensity values in T (bifunctional dependence).

Monofunctional Dependence Assumption

Our goal is to characterize the mapping from

voxel intensities in S to those in T , knowing that some elements of C are erroneous, i.e. that would not be present in C if S and T were perfectly matched. Let us assume here that the intensity in T is a function of the intensity in S corrupted with an additive stationary Gaussian white noise η,



T x

 



f S x   η x 

(2)

where f is an unknown function to be estimated. This is exactly the model employed by Roche et al. (Roche et al., 2000) which leads to the correlation ratio as the measure to be maximized for registration. In that approach, for a given transformation, one seeks the function that best describes T in terms of S. In a maximum likelihood context, this function is actually a least squares (LS) fit of T in terms of S. 7

The major difference between our respective problems is that we seek a high-dimensional geometrical transformation. As opposed to affine registration where the transformation is governed by the majority of good matches, elastic deformations may be computed using mainly local information (i.e. gradients, local averages, etc.). Hence, we cannot expect good displacements in one structure to correct for bad ones in another; we have to make certain each voxel is moved properly at each iteration. For this, since the geometrical transformation is found using intensity similarity, the most precise intensity transformation is required. Consequently, instead of performing a standard LS regression, we have opted for a robust linear regression estimator which will remove outlying elements of C during the estimation of the intensity transformation. To estimate f we use the least trimmed squares (LTS) method followed by a binary reweighted least squares (RLS) estimation (Rousseeuw and Leroy, 1987). The combination of these two methods provides a very robust regression technique with outlier detection, while ensuring that a maximum of pertinent voxel pairs are taken into account. Different types of functions can be used to model f . In (Guimond et al., 2001) we made use of polynomial functions. The intensity correspondences between T and S is then defined as:



T x where θ

 θ0    θ p 





θ0 θ1 S x  θ2 S x 

2





θ pS x

p



(3)

needs to be estimated and p is the degree of the polynomial function.

This model is adequate to register images that have a vast range of intensities; the restricted polynomial degree imposes intensity space constraints on the correspondences, mapping similar intensities in S to similar intensities in T . In the case that S is a labeled image, neighboring intensities in S will usually correspond to different structures. Hence the intensity space constraint is no longer required. f is then modeled as a piecewise constant function, such that each label of S is mapped to the LTS/RLS estimate of intensities corresponding to that label in T .

8

Bifunctional Dependence Assumption

Functional dependence as expressed in (2) and in (3)

implies that two structures having similar intensity ranges in S should also have similar intensity ranges in T . With some combinations of images, this is a crude approximation. For example, CSF and bones generally give similar intensity values in T1-weighted images, while they appear with very distinct values in PD-weighted scans. Conversely, CSF and gray matter are well contrasted in T1-weighted images, while they correspond to similar intensities in PD weighted scans. To circumvent this difficulty, we have developed a strategy that enables the mapping of an intensity value in S to not only one, but two possible intensity values in T . This method is a natural extension of the previous section. Instead of computing a single function that maps the intensities of S to those of T , two functions are estimated and the mapping becomes a weighted sum of these two functions. We start with the assumption that if a point has an intensity s in S, the corresponding point in T has an intensity t that is normally distributed around two possible values depending on s,





f1 s  and f2 s  . In statistical terms, this means that given s, t is drawn from a mixture of Gaussian distribution,



P t  s





where π1 s  and π2 s 



 

π1 s  N f1 s  σ

2





 

π2 s  N f2 s  σ2 

(4)



1  π1 s  are mixing proportions that depend on the intensity in the

source image, and σ2 represents the variance of the noise in the target image. Consistent with the functional relationship, we will restrict ourselves to polynomial intensity functions, i.e.



f1 s 



θ0 θ1 s θ2 s2   θ p s p , and f2 s 

ψ0 ψ1 s ψ2 s2  ψ p s p .

An intuitive way to interpret this modeling is to state that for any voxel, there is a binary “selector” variable ε  1  2  that would tell us, if it was observed, which of the two functions f1 or f2 actually serves to map s to t. Without knowledge of ε, the best intensity correction to

9

apply to S (in the minimum variance sense) is a weighted sum of the two functions,



f s t 









1  s  t  f1 s  P ε



2  s  t  f2 s 

(5)

in which the weights correspond to the probability that the point be mapped according to either the first or the second function. To estimate the functions, we employ a sequential strategy that performs two successive LTS/RLS regressions as in the monofunctional case. Details on how the other parameters are determined can be found in (Guimond et al., 2001).

1.2.2

Geometrical Transformation

Having completed the intensity transformation stage, we end up with an intensity corrected



version of the source image, which will be denoted S . In the monofunctional case S x 

 



f S x  and in the bifunctional case S x 

 



f S x  T x  . We may assume that S is roughly

of the same modality as T in the sense that corresponding anatomical structures have similar intensities in S and T . The geometrical transformation problem may then be treated in a monomodal registration context. Many algorithms have been developed that deform one brain so its shape matches that of another (Maintz and Viergever, 1998; Toga, 1999). The procedure used here was influenced by a variety of optical flow methods, primarily the demons algorithm (Thirion, 1995; Thirion,



1998). At a given iteration n, each voxel x of T is displaced according to a vector vn x  so as to match its corresponding anatomical location in Sn . We use the following scheme: vn

 1



x

 Sn  ! hn x   Gσ  vn ! ∇Sn  hn x  2  Sn 



T x   hn x " T x 

where Gσ is a 3D Gaussian filter with isotropic variance σ2 ,





2

∇Sn  hn x  #$

(6)

denotes the convolution,





denotes the composition, ∇ is the gradient operator and the transformation hn x  is related to

10





the displacement by hn x 

x vn x  . As is common with registration methods, we also make

use of multilevel techniques to accelerate convergence. Details about the number of levels and iterations as well as filter implementation issues are addressed in Section 1.3. We show here how our method can be related to three other registration methods: the minimization of the sum of squared difference (SSD) criterion; optical flow; and, the demons algorithm.

1.2.3

Relation to SSD Minimization

In the SSD minimization framework, one searches for the transformation h that minimizes the sum of squared differences between the transformed source image and the target image. The SSD is then defined as:



SSD h 

1 N  S 2 x∑ % 1





h x





T x

 2

(7)

The minimization of (7) may be performed using a gradient descent algorithm. By differ-



entiating the above equation, we get for a given x: ∇SSD h 

  S 





h x & T x   ∇S





h x .

Thus, the gradient descent consists of an iterative scheme of the form: vn 

1

vn α  Sn 







hn x " T x   ∇Sn  hn x 

(8)

where α is the step length. If we set α to a constant value, this method corresponds to a first order gradient descent algorithm. Comparing (8) to (6), we see that our method sets α

!  ∇

 ! 1  S  hn  x  2  T x "



S  hn x 

(9)

2

and applies a Gaussian filter to provide a smooth displacement field. Cachier et al. (Cachier et al., 1999; Pennec et al., 1999) have shown that using (9) closely relates (6) with a second order gradient descent of the SSD criterion, in which each iteration n sets hn

1

to the minimum

of the SSD quadratic approximation at hn . We refer the reader to these articles for a more 11

technical discussion on this subject.

1.2.4

Relation to Optical Flow



T and S are considered as successive time samples of an image sequence represented by I x  t  ,



where x

x1  x2  x3  is a voxel position in the image and t is time. The displacements are

computed by constraining the brightness of brain structures to be constant in time so that



dI x  t  dt

0

(10)

It is well known that (10) is not sufficient to provide a unique displacement for each voxel. In fact, this constraint leads to



f x





 ! ∂I x  t ' ∂t! 2 ∇xI x  t  ∇x I x  t 

(11)

which is the component of the displacement in the direction of the brightness gradient (Horn and Schunck, 1981). Other constraints need to be added to (10) to obtain the displacement components in other directions. Many methods have been proposed to fulfill this purpose and thus regularize the resulting vector field (Barron et al., 1994). One that can be computed very efficiently was proposed by Thirion (Thirion, 1998) in his description of the demons registration method, using a complete grid of demons. It consists of smoothing each dimension of the vector field with a Gaussian filter Gσ . He also proposed to add



  ∂I x  t '

∂t  2 to the denominator of (11) for

numerical stability when ∇x I x  t  is close to zero, a term which serves the same purpose as α2 in the original optical flow formulation of Horn and Schunck (Horn and Schunck, 1981). As is presented by Bro-Nielsen and Gramkow (Bro-Nielsen and Gramkow, 1996), this kind of regularization approximates a linear elasticity transformation model.

12

With this in mind, the displacement that maps a voxel position in T to its position in S is found using an iterative method, vn 

 1

x

Gσ 



vn 

!



 ∂I x  t ' ∂t ! ∇x I x  t  ∇x I x  t  2  ∂I x  t ' ∂t  2 

# 

(12)

Spatial derivatives may be computed in several ways (Horn and Schunck, 1981; Brandt, 1997; Simoncelli, 1994). We have observed from practical experience that our method performs best when they are computed from the resampled source image of the current iteration. As shown in Section 1.2.3, this is in agreement with the SSD minimization. Temporal derivatives are obtained by subtracting the target images from the resampled source image of the current iteration. These considerations relate (12) to (6). The reader should note that the major difference between this method and other optical flow strategies is that regularization is performed after the calculation of the displacements in the gradient direction instead of using an explicit regularization term in a minimization framework.

1.2.5

Relation to the Demons Algorithm

Our algorithm is actually a small variation of the demons method (Thirion, 1995; Thirion, 1998) using a complete grid of demons, itself closely related to optical flow as described in the previous section. The demons algorithm finds the displacements using the following formula: vn

 1



x



 S!  hn x   T x  Gσ   vn ! ∇T x  2  S  hn x "



T x



2

∇T x 

# 

(13)

In comparing (13) and (6), it is apparent that the only difference between our formulation and the demon method is that derivatives are computed on the resampled source image of the current iteration. This modification was performed following the observations on the minimization of the SSD criterion.

13

1.3

Results and Discussion

In the following section we present registration results involving images obtained from multiple modalities. First, we show a typical example where monofunctional dependence can be assumed: the registration of an atlas (Collins et al., 1998b) with a T1-weighted MR image. We next present an example where bifunctional dependence may be assumed: the registration of a PD-weighted image with the same T1-weighted image. All of the images used in this section have a resolution of 1 ( 1

(

1mm3 and respect the

neurological convention, i.e. on coronal slices, the patient’s left is on the left side of the image. Before registration, images are affinely registered using the correlation ratio method (Roche et al., 1998). The multilevel process was performed at three resolution levels, namely 4mm, 2mm and 1mm per voxel. Displacement fields at one level are initialized from the result of the previous level. The initial displacement field v0 is set to zero. 128 iterations are performed at 4mm/voxel, 32 at 2mm/voxel and 8 at 1mm/voxel. These are twice the number of iterations used for registration of monomodal images using the conventional demons algorithm. We believe that making use of a better stopping criterion, such as the difference of the SSD values between iterations, would probably improve the results shown below. This aspect is presently under investigation. The Gaussian filter Gσ used to smooth the displacement field has a standard deviation of 1 voxel regardless of the resolution. This models stronger constraints on the deformation field at the beginning of the registration process to correct for gross displacements, and weaker constraints near the end when fine displacements are sought. The resampling process makes use of trilinear interpolation, except in the case of the atlas where nearest-neighbor interpolation is used. Computation time to obtain the following results is around 60 minutes on a 450MHz PC with 500MB of RAM (10 minutes at 4mm, 20 minutes at 2mm and 30 minutes at 1mm). Most

14

of this time ( ) 85%) is devoted to the intensity correction part, which has not been optimized in this first version of our program. The other 15% is taken by the standard registration code which is stable and well optimized.

1.3.1

Monofunctional Dependence

We present here the result of registering the atlas with a T1-weighted image. This is a typical example of monofunctional dependence between the intensities of the images to register: since the atlas can be used to generate realistic MR images, it is safe to assume a functional dependence between the intensity of the atlas and those of the T1-weighted image. Also, since the source image S is a labeled image, the function f is modeled as a piecewise constant function. In this case, each intensity level (10 in all) corresponds to a region from which to estimate the constant function. The result of registration is presented in Figure 1. The first image (Figure 1a) shows one slice of the atlas. The second one (Figure 1b) is the corresponding slice of the T1-weighted image. The third and fourth images (Figures 1c and 1d) present the result of registering the atlas with the T1-weighted image using our algorithm. Figure 1c shows the result without the intensity transformation; we have simply applied to the atlas the geometrical transformation resulting from the registration procedure. Figure 1d shows the image resulting from the registration process. It has the same shape as the T1-weighted image (Figure 1b) and intensities have been transformed using the intensity correction. To facilitate the visual assessment of registration accuracy, contours obtained using a Canny-Deriche edge detector (on the T1-weighted image) have been overlaid over each image in Figure 1. Figure 1e shows the joint histogram of intensities after registration. Values have been compressed logarithmically and normalized as is depicted in the color scale. The histogram is color-coded and ranges from dark red representing high point densities to light yellow depict-

15

ing low point densities. The values of the piecewise constant function f are overlaid as white dots.

1.3.2

Bifunctional Dependence

When registering images from different modalities, monofunctional dependence may not necessarily be assumed. We presented in Section 1.2.1 such an example: the registration of PD and T1-weighted images. The main problem in this case is that the CSF/GM intensity of the PD image needs to be mapped to two different intensities in the T1-weighted scan. To solve this problem, we applied the method described in Section 1.2.1 to register PD and T1-weighted images where two polynomial functions of degree 12 are estimated. This polynomial degree was set arbitrarily to a relatively high value to allow significant intensity transformations. As shown in Figure 1f-j, the CSF, which is white in the PD-weighted image (Figure 1f) and black in the T1-weighted image (Figure 1g), is well registered. Also, the intensity transformation is adequate (Figure 1i). The same comments apply to the GM, which is white in the PD-weighted image (Figure 1f) and gray in the T1-weighted image (Figure 1g). Figure 1j presents the joint histogram of the two images after registration. The functions f1 and f2 found during the registration process are superimposed, the blue line corresponds to f1 and the green one to f2 . The line width for a given intensity s is proportional to the value of the



corresponding πε s  . As can be observed in Figure 1j, the polynomial functions f1 and f2 fit well with the high density clusters of the joint histogram. In particular, we see that the CSF/GM intensity values from the PD-weighted image (with values around 220) get mapped to two different intensity values in the T1-weighted scan: 75 and 45. The mapping to 75 represents the GM (red polynomial) while the mapping to 45 represents CSF (blue polynomial).

16

Note that in the registration of the atlas with the T1-weighted image and the PD- with the T1- weighted image, we selected as the source image the one which has the best constrast between structures. This is simply because our algorithm permits many structures of the source image to be mapped to a single intensity. But a single intensity in the source image can be mapped to at most two intensities in the target image. Hence, it is always better to use the image with the greater number of visible structures as the source image.

1.4

Conclusion

We have presented an original method to perform non-rigid registration of multimodal images. This iterative algorithm is composed of two steps: the intensity transformation and the geometrical transformation. Two intensity transformation models were described which assume either monofunctional or bifunctional dependence between the intensity values in the images being matched. Both of these models are built using robust estimators to enable precise and accurate transformation solutions.

17

1

250

150

0.5

100 50 gro un Gr d a C y SF W M hit at e M ter att Mu e scl F r e / at Sk i Sk n in Gl ial Sku Co Mat ll nn ter ect ive

d

200

0

Ba

ck

c

b

T1 Intensities

a

Atlas

e

h

g

i

1

250

T1 Intensities

f

200 150

0.5

100 50 0

75

150

225

0

PD Intensities

j Figure 1: Results of 3D registration. Registration of Atlas with T1-weighted image: (a) Atlas (Source). (b) T1weighted image (Target). (c) Atlas without intensity correction, after registration with T1. (d) Atlas with intensity correction, after registration with T1. (e) The joint histogram of the atlas and T1-weighted image after registration; values range from dark red representing high point densities to light yellow depicting low point densities; the white dots correspond to the intensity transformation found by registering the atlas with the T1-weighted image and assuming monofunctional dependence (piecewise constant function). Registration of PD-weighted with T1-weighted image: (f) PD-weighted image (Source). (g) T1-weighted image (Target). (h) PD-weighted image without intensity correction, after registration with T1-weighted image. (i) PD-weighted image with intensity correction after registration with T1-weighted image. (j) The joint histogram of PD-weighted image and T1weighted imagem, after registration; the blue line corresponds to f1 and the green one to f2 ; the line width for a given intensity value s in T1 corresponds to the value of the corresponding πε * s + . Contours were obtained using a Canny-Deriche edge detector on the targets (b and g) and overlaid on the other images to better assess the quality of registration. The joint histograms values have been compressed logarithmically and normalized as is depicted 18 in the color scale.

2 Image Fusion During Neurosurgery with a Biomechanical Model of Brain Deformation Introduction A critical goal of neurosurgery is to accurately locate, access and remove intracranial lesions without damaging healthy brain tissue. The overriding goal is to preserve neurological function. This requires the precise delineation of the functional anatomy and morphology of the patient’s brain, as well as lesion margins. The similar visual appearance of healthy and diseased brain tissue (e.g. as with infiltrating tumors) and the inability of the surgeon to see critical structures beneath the brain surface can pose difficulties during the operation. Some critical structures, such as white matter fiber tracts, may not be visible at all. Moreover, the difficulty in perceiving lesion (e.g. tumor) boundaries makes complete resection extremely difficult (Jolesz, 1997). Over the last decade, advances in image-guided neurosurgery (IGNS) techniques have contributed to the growth of minimally-invasive neurosurgery. These procedures must be carried out in operating rooms which are specially-equipped with imaging systems. These systems are used to acquire images intraoperatively, as necessitated by the procedure. The improved visualization of deep structures and the improved contrast between the lesion and healthy tissue (depending on the modality) allow the surgeon to plan and execute the procedure with greater precision. IGNS has largely been a visualization-driven task. In the past, it had not been possible to make quantitative assessments of intraoperative imaging data, and instead physicians relied on qualitative judgments. In order to create a rich visualization environment which maximizes the information available to the surgeon, previous work has been concerned with image acquisition, registration and display. Algorithm development for computer-aided IGNS has focussed on improving imaging quality and speed. Another major focus has been to develop sophisticated

19

multimodality image registration and fusion techniques, to enable fusion of preoperative and intraoperative images. However, clinical experience with IGNS involving deep brain structures has revealed the limitations of existing rigid registration approaches. This motivates the search for nonrigid techniques that can rapidly and faithfully capture the morphological changes that occur during surgery. In the future, the use of computer-aided surgical planning will include not only three dimensional (3D), models but also information from multiple imaging modalities, registered into the patient’s reference frame. Intraoperative imaging and navigation will thus be fully integrated. Various imaging modalities have been used for image guidance. These include, among others, digital subtraction angiography (DSA), computed tomography (CT), ultrasound (US), and intraoperative magnetic resonance imaging (IMRI). IMRI represents a significant advantage over other modalities because of its high spatial resolution and superior soft tissue contrast. However, even the most advanced intraoperative imaging systems cannot provide the same image resolution or tissue contrast selection features as preoperative imaging systems. Moreover, intraoperative imaging systems are by necessity limited in the amount of time available for imaging. Multimodality registration can allow preoperative data that cannot be acquired intraoperatively [such as nuclear medicine scans (SPECT/PET), functional MRI (fMRI), MR angiography (MRA), etc.] to be visualized alongside intraoperative data.

2.1

Nonrigid Registration for IGNS

During neurosurgical operations, changes occur in the anatomical position of brain structures and adjacent lesions. The leakage of cerebrospinal fluid (CSF) after opening the dura, hyperventilation, the administration of anesthetic and osmotic agents, and retraction and resection of tissue all contribute to shifting of the brain parenchyma. This makes information based upon preoperatively acquired images unreliable. The loss of correspondence between pre- and in-

20

traoperative images increases substantially as the procedure continues. These changes in the shape of the brain have been widely recognized as nonrigid deformations called “brain shift” (see (Nabavi et al., 2001)). Suitable approaches to capture these shape changes and to create integrated visualizations of preoperative data in the configuration of the deformed brain are currently in active development. Previous work aimed at capturing brain deformations for neurosurgery can be grouped into two categories. In the first categorey are those approaches that use some form of biomechanical model (recent examples include (Hagemann et al., 1999; Skrinjar and Duncan, 1999; Miga et al., 1999; Skrinjar et al., 2001; Ferrant et al., 2000b)). In the second category are those approaches that use phenomenological methods, relying upon image related criteria (recent examples include (Hill et al., 1998; Hata, 1998; Ferrant et al., 1999b; Hata et al., 2000)). Purely image-based matching may be able to achieve a visually pleasing alignment, once issues of noise and intensity artifact are accounted for. However, in our work on intraoperative matching we favor physics-based models which ultimately may be expanded to incorporate important material properties (such as inhomogeneity, anisotropy) of the brain, once these are determined. However, even among physics-based models, there exist a spectrum of approaches, usually involving a trade-off between physical plausibility and speed. A fast surgery simulation method has been described in (Bro-Nielsen, 1996). Here, high computational speeds were obtained by converting a volumetric finite element model into a model with only surface nodes. The goal of this work was to achieve very high graphics speeds consistent with interactive computation and display. This is achieved at the cost of simulation accuracy. This type of model is best suited to computer graphics-oriented display, where high frame rates are needed. A sophisticated finite element based biomechanical model for two-dimensional brain deformation simulation has been proposed in (Hagemann et al., 1999). In this work, correspon-

21

dences were established by manual interaction. The elements of the finite element model were the pixels of the two dimensional image. The manual determination of correspondences can be time consuming, and is subject to human error. Moreover, when methods are generalized to three dimensions, the number of points which must be identified can be very large. Due to the realities of clinical practice, two-dimensional results are not practical. A (three dimensional) voxelwise discretization approach, while theoretically possible, is extremely expensive from a computational standpoint (even considering a parallel implementation) because the number of voxels in a typical intraoperative MRI dataset leading to a large number of equations to solve (256x256x60 = 3,932,160 voxels, which corresponds to 11,796,480 displacements to determine). Downsampling can lead to fewer voxels, but leads to a loss of detail. Edwards et al. (Edwards et al., 1997) presented a two dimensional three component model for tracking intraoperative deformation. This work used a simplified material model. However, the initial 2D multigrid implementation required 120–180 minutes when run on a Sun Microsystems Sparc 20, which may limit its feasibility for routine use. Skrinjar et al. (Skrinjar and Duncan, 1999) have presented a very interesting system for capturing real-time intraoperative brain shift during epilepsy surgery. In this context, brain shift occurs slowly. A very simplified homogeneous brain tissue material model was adopted. Following the description of surface based tracking from intraoperative MRI driving a linear elastic biomechanical model in (Ferrant et al., 2000b), Skrinjar et al. presented a new implementation (Skrinjar et al., 2001) of their system using a linear elastic model and surface based tracking from IMRI with the goal of eventually using stereoscopic cameras to obtain intraoperative surface data and hence to capture intraoperative brain deformation. Paulsen et al. (Paulsen et al., 1999) and Miga et al. (Miga et al., 1999; Miga et al., 2001) have developed a sophisticated finite element model to simulate brain deformation. Their model is interesting because it incorporates simulations of forces associated with tumor tissue, as well

22

as those resulting from retraction and resection. A limitation of the existing approach is that the preoperative segmentation and tetrahedral finite element mesh generation currently require around five hours of operator time. Nevertheless, this approach holds promise in actually predicting brain deformation. The real-time needs of surgery dictate that any algorithm used for prospective image matching must rapidly, reliably and accurately capture nonrigid shape changes in the brain which occur during surgery. Our approach is to construct an unstructured grid representing the geometry of the key structures in the image dataset. This technique allows us to use a finite element discretization that faithfully models key characteristics in important regions while reducing the number of equations to solve by using mesh elements that span multiple voxels in other regions. The algorithm allows the projection of preoperative images onto intraoperative images, thereby allowing the fusion of images from multiple modalities and spanning different contrast mechanisms. We have used parallel hardware, parallel algorithm design and efficient implementations to achieve rapid execution times compatible with neurosurgery.

2.2

Method

Figure 2 is an overview, illustrating the image analysis steps used during intraoperative image registration. The critical image processing tasks include segmentation (the identification of anatomical structures) and registration. Segmentation data are used both for preoperative planning, and to create intraoperative segmentations. The segmentation data are used to calculate an initial affine transformation (rotation, translation, scaling) which rigidly registers the images, thus initializing the data for nonrigid matching using our biomechanical simulation. Using the biomechanical model, the volumetric deformation is inferred through a mechanical simulation with boundary conditions established via surface matching. This sophisticated deformation model can be solved during neurosurgery, providing enhanced intraoperative visualization.

23

Im age A naly sis During Im age Guided Neurosurgery :

Preoperativ e data

I nitial Registration

T issue Segmentation

I ntraoperativ e MRI Surface Matching

Biomechanical Simulation

Visualization Neurosurgery progresses

Figure 2: Schematic illustrating image analysis tasks carried out during neurosurgery. 2.2.1

Preoperative Data Acquisition and Processing

The time available for image processing during surgery is extremely limited compared to that available preoperatively. Consequently, preoperative data acquisition can be more comprehensive, and more extensive image analysis (for example segmentation) can be performed. A variety of manual (Gering et al., 1999), semi-automated (Kikinis et al., 1992; Yezzi et al., 2000) and automated (Warfield et al., 2000a; Kaus et al., 1999; Warfield et al., 2000b) segmentation approaches are available. We select the most accurate, robust approach based upon the preoperative data and the particular critical structures. For the matching experiments which will be described below, we have used an anatomical atlas, although other preoperative data such as magnetic resonance angiography or diffusion tensor images could ultimately also be used. The atlas was constructed from a high resolution scan of a single patient, in which over 200 structures were segmented (Kikinis et al., 1996) using a combination of automated and interactive techniques. During surgery, we are especially interested in the corticospinal tract, a region of white matter which can be difficult or impossible to directly observe with

24

conventional MRI, and which must be preserved. We have previously shown that we can project the corticospinal tract from the atlas onto patient scans for preoperative surgical planning (Kaus et al., 2000).

2.2.2

Intraoperative Image Processing

Intraoperative image processing consists of: 1.) acquiring one or more intraoperative volumetric data sets; 2.) constructing a segmentation of the intraoperative acquisition; 3.) computing an affine registration of the preoperative data onto the new acquisition; 4.) identifying the correspondences between key surfaces of the preoperative and intraoperative data; 5.) solving a biomechanical model to infer a volumetric deformation field; 6.) applying the deformation to the preoperative data and constructing a new visualization merging critical structures from the preoperative data with the intraoperative data.

Segmentation of Intraoperative Volumetric Images In the experiments conducted below, a rapid segmentation of the brain and ventricles was obtained using a binary curvature driven evolution algorithm (Yezzi et al., 2000). The region identified as brain or ventricle was then interactively corrected to remove misclassified tissue using the software described by Gering et al. (Gering et al., 2001). This approach allows the surgeon to inspect and interactively edit the segmentation data, increasing its accuracy. Alternatively, we have experimented with automated intraoperative segmentation (Warfield et al., 1998b; Warfield et al., 2000a) utilizing tissue classification in a multi-channel feature space using a model of expected anatomy as an initial template. Automated approaches will likely be preferable once they have been validated.

25

Unstructured Mesh Generation and Surface Representation We have implemented a mesh generator which is optimized for use with biomedical structures, building upon previously described techniques (Schroeder et al., 1996; Geiger, 1993). During mesh generation, we extract an explicit representation of the surface of the brain and ventricles based on the preoperative segmentation. We also create a volumetric unstructured mesh using a multiresolution version of the marching tetrahedra algorithm. The mesher begins by subdividing the image into cubes, which are then divided into 5 tetrahedra using an alternating split pattern which prevents diagonal crossings on the shared faces. The mesh is iteratively refined in the region of complex boundaries, and then a marching tetrahedra-like approach is applied to this multiresolution mesh. For each cell of the final mesh, the label value of each vertex is checked, and if different, the tetrahedron is divided along the edge having different node labels. A detailed description can be found in (Ferrant et al., 1999b; Ferrant et al., 2000a). The meshing process is extremely robust, allowing us to generate tetrahedral meshes of the brain and ventricles from rapid segmentations of each volumetric intraoperative acquisition carried out during the surgery. This facilitates intraoperative matching from one volumetric acquisition to the next.

Affine Registration of Preoperative to Intraoperative Image Datasets For affine registration (rotation, translation, scaling), we use a fast parallel implementation of a robust algorithm which is based upon aligning segmented image data, using a rapidlyconverging multiresolution search strategy (Warfield et al., 1998a). Applying the resulting transform, segmentations and greyscale data from the preoperative and intraoperative scans are rigidly registered.

26

Volumetric Biomechanical Simulation of Brain Deformation During the procedure, the brain undergoes nonrigid shape changes for the reasons described above. During IGNS the surgeon is able to acquire a new volumetric MRI when he wishes to review the current configuration of the entire brain. A volumetric deformation field relating earlier acquisitions to this new scan is computed by first matching surfaces from the earlier acquisition to the current acquisition, and then calculating the volumetric displacements by using the surface displacements as boundary conditions. The critical concept is that forces are applied to the volumetric model that will produce the same surface displacements as were obtained by the surface matching. The biomechanical model can then be used to compute a volumetric deformation map.

Establishing Surface Correspondences

The surfaces of the brain and lateral ventricles are

iteratively deformed using a dual active surface algorithm. Image-derived forces are applied iteratively to an elastic membrane surface model of the early scan, thereby deforming it so as to match the boundary of the current acquisition. The derived forces are a decreasing function of the image intensity gradients, so as to be minimized at the edges of objects in the volume. We have included prior knowledge about the expected gray level and gradients of the objects being matched to increase the convergence rate of the process. This algorithm is fully described in (Ferrant et al., 1999a).

Biomechanical Simulation of Volumetric Brain Deformation We treat the brain as a homogeneous linearly elastic material. The deformation energy of an elastic body Ω, under no initial stresses or strains, and subject to externally applied forces, can be described by the following model (Zienkiewicz and Taylor, 1994):



E u

1 2,



σ - ε dΩ

27

 ,



u- F dΩ 

(14)

where the variables are given in terms of the stress vector, σ, the strain vector, ε, the forces



F

F x  y z  applied to the elastic body (forces per unit volume, surface forces or forces con-

 

centrated at the nodes of the mesh) and u





u x  y z  v x  y z  w x  y z   - , the displacement

vector field we wish to compute. Since we are using a linear elasticity framework, we assume small deformations. Hence the strain vector ε is given by ε



∂u ∂v ∂w ∂u ∂v ∂v ∂w ∂w ∂u #      ∂x ∂y ∂z ∂y ∂x ∂z ∂y ∂x ∂z

which can be written as ε



(15)

L u where L is a linear operator. The elastomechanical relation

between stresses and strains can be expressed by the generalized Hooke’s law as σ



σx  σy  σz  τxy  τyz  τzx 

-

Dε 

(16)

Assuming isotropic material properties for each point, we obtain a symmetric elasticity matrix D in the form

.//

D

1

E ν  1  2ν 

33

1 ν

ν

ν

0

0

0

//

ν

1 ν

ν

0

0

0

ν

ν

1 ν

0

0

0

33

0

0

33

1 2ν 2

0

// 

2433

//

//0

//

//

1

0

0

0

1 2ν 2

0

0

0

0

0

0

0

0

1

0

1

1 2ν 2

33 33 33

(17)

5

with physical parameters E (Young’s modulus) and ν (Poisson’s ratio). See (Zienkiewicz and Taylor, 1994) for the full details. For the discretization, we use the finite element method applied over the volumetric image domain so that the total potential energy can be written as a sum of potential energies for each

28



element: E u 



e e ∑e % nodes 1 E u  . The mesh is composed of tetrahedral elements and thus each N

element is defined by four mesh nodes. The continuous displacement field u everywhere within element e of the mesh is defined as a function of the displacement at the element’s nodes uei weighted by the element’s interpolating functions



Nnodes



u x

%

Nie



x ,



I Nie x  uei



(18)

i 1

Linear interpolating functions are used to define the displacement field inside each element. The interpolating function of node i of tetrahedral element e is defined as



1 ae bei x cei y die z 7 6V e 6 i

Nie x 

(19)

The computation of the volume of the element V e and the interpolation coefficients are detailed in (Zienkiewicz and Taylor, 1994, pages 91–92). The volumetric deformation of the brain is found by solving for the displacement field that minimizes the deformation energy described by Equation (14). For our finitite element approach this is described by



M

∑ δE

δE u 

%

e



ue 

0

(20)

e 1

where



δE u e

e



Nnodes

∑ %

i 1

∂ e e E u  δui e e ∂ui



Nnodes

∑ %

i 1

∂ e e E u  δvi e e ∂vi



Since δui e  δvi e and δwi e are independent, defining matrix Be

29

Nnodes

∑ %

i 1



∂ e e E u  δwi e e ∂wi

Bei 

Nnodes i 1

%

with Bei



(21)

LNie for

every node i of each element e, yields in the following equation:

0

,



Be - D Be ue dΩ

with the element stiffness matrix Ke

98

e ΩB

-

 ,



Ne - Fe dΩ 

(22)

D Be dΩ. An assembly of the equations for

all elements finally leads to a global linear system of equations, which can be solved for the displacements resulting from the forces applied to the body: F

Ku

(23)

The displacements at the boundary surface nodes are fixed to match those generated by the active surface model. Let u : be the vector representing the displacement to be imposed at the boundary nodes. The elements of the rows of the stiffness matrix K corresponding to the nodes for which a displacement is to be imposed are set to zero and the diagonal elements of these rows are set to one. The force vector F is set to equal the displacement vector for the boundary nodes: F

:u (Zienkiewicz and Taylor, 1994). In this way solving Equation (23) for

the unknown displacements will produce a deformation field over the entire volumetric mesh that matches the prescribed displacements at the boundary surfaces.

Hardware and Implementation The volumetric deformation of the brain is computed by solving for the displacement field that minimizes the energy described by Equation (14), after fixing the displacements at the surface to match those generated by the active surface model. Three variables, representing the x, y and z displacements, must be determined for each element of the finite element mesh. Each variable gives rise to one row and one column in the global K matrix. The rows of the matrix are divided equally amongst the CPUs available for

30

computation and the global matrix is assembled in parallel. Each CPU assembles the local Ke matrix for each element in its subdomain. Although each CPU has an equal number of rows to process, because the connectivity of the mesh is irregular, some CPUs may do more work than others. Following matrix assembly, the boundary conditions determined by the surface matching are applied. The global K matrix is adjusted such that rows associated with variables that are determined consist of a single non-zero entry of unit magnitude on the diagonal. The volumetric biomechanical brain model system of equations (and the active surface membrane model equations) are solved using the Portable, Extensible Toolkit for Scientific Computation (PETSc) package (Balay et al., 1997; Balay et al., 2000a) using the Generalized Minimal Residual (GMRES) solver with block Jacobi preconditioning. During neurosurgery, the system of equations was solved on a Sun Microsystems SunFire 6800 symmetric multiprocessor machine with 12 750MHz UltraSPARC-III (8MB Ecache) CPUs and 12 GB of RAM. This architecture gives us sufficient compute capacity to execute the intraoperative image processing prospectively during neurosurgery.

Intraoperative Visualization Once the volumetric deformation field has been computed, it can be applied to earlier data to warp it into the current configuration of the patient anatomy. The imaging data can then be displayed by texture mapping onto flat planes to facilitate comparisons with current intraoperative data as well as prior scans. Triangle models of segmented structures (i.e. based on registered volumetric data) can be used to display surface renderings of critical anatomical structures, overlaid on intraoperative image data. This allows ready appreciation of the 3D anatomy of these segmented structures together with the imaging data in the form of planes passing through or over the 3D triangle models (Gering et al., 2001). This augments the sur-

31

Figure 3: Open-configuration magnetic resonance scanner during neurosurgery. geon’s ability to see critical structures which must be preserved (such as the corticospinal tract) and to better appreciate the lesion and its relationship to other structures. Figure 3 shows the open-configuration magnetic resonance scanner optimized for imaging during surgical procedures (Jolesz, 1997; Black et al., 1997). The image we constructed was presented on the LCD and increased the information available to the surgeon as the operation progressed.

2.3

Results and Discussion

The image registration strategy described here has been applied prospectively during several neurosurgical cases. We present here illustrative results which demonstrate the ability of our algorithm to capture intraoperative brain deformations. The enhancement provided by intraoperative nonrigid registration to the surgical visualization environment is shown by our matching the corticospinal tract of a preoperatively prepared anatomical atlas to the initial and subsequent intraoperative scans of a subject. This matching

32

Figure 4: This figure shows the corticospinal tract from our anatomical atlas in blue, projected into the shape of the brain of the subject shown in Figure 3. was carried out prospectively during the neurosurgery, demonstrating the practical value of the approach and its ability to meet the real-time constraints of surgery. We have also conducted parallel scaling experiments which have yielded very encouraging results. The entire image analysis process can be completed in less than 10 minutes, which has been adequate to display the information to the surgeon. Interestingly, the most computationally intensive task (the biomechanical simulation) has also been optimized the most and is now the fastest step. We anticipate that segmentation techniques requiring less user interaction will result in significant improvements in speed.

Biomechanical Simulation of Volumetric Brain Deformation Figure 4 shows the corticospinal tract from our anatomical atlas projected into the shape of the brain of the subject. This visualization helps the surgeon to better appreciate the 3D relationship of this essential structure to the lesion and other regions of the brain. The corticospinal tract cannot be readily observed in IMRI acquisitions.

33

(a) A single slice from an early 3D IMRI scan.

(b) The corresponding slice in a later 3D IMRI scan, showing significant brain shift has occurred.

(c) The matched slice of the first volume after simulation of the brain deformation.

(d) Visualization of the magnitude of the deformation field computed in matching image (a) to image (b).

Figure 5: Two dimensional slices through three-dimensional data, showing the match of the simulated deformation of the initial brain onto the actual deformed brain. The quality of the match is significantly better than can be obtained through rigid registration alone.

34

Figure 5 is a typical case illustrating the amount of brain deformation which can occur during neurosurgery, as well as the effectiveness of our algorithm in capturing this shift during neurosurgery. As shown, the quality of the match is significantly better than can be obtained through rigid registration alone. Our early experience has shown that our intraoperative biomechanical simulation of brain deformation is a robust and reliable method for capturing the changes in brain shape that occur during neurosurgery. The registration algorithm requires no user interaction and the parallel implementation is sufficiently fast to be used intraoperatively. We intend to incorporate patientspecific preoperative data in place of the anatomical atlas to increase the surgical value of the intraoperative updates. As we refine our approach, we expect to appreciate performance gains based on more automated segmentation methods, and further optimized parallel implementations which address load imbalances. Improvements in the accuracy of the match could result from a more sophisticated model of the material properties of the brain (such as more accurate modeling of the cerebral falx and the lateral ventricles). Sophisticated MR imaging methods such as diffusion tensor MRI now enable the preoperative imaging of inhomogeneous anistropic white matter structure, which could be incorporated into the material model. Ultimately, the prediction of brain deformation, as opposed to the capture of observed deformation described here, will most likely require a nonlinear material model together with extensive monitoring of physiological data. Such prediction could be used to indicate when new intraoperative imaging is necessary to appropriately update both the simulation model and the surgeon’s understanding of the brain shape.

35

2.4

Conclusion

Nonrigid changes in brain morphology occur during neurosurgery and limit the usefulness of preoperative imaging for intra-treatment planning and surgical navigation. Intraoperative nonrigid registration can add significantly to the value of intraoperative imaging. It provides for quantitative monitoring of therapy application, including the ability to make quantitative comparisons with a preoperatively-defined treatment plan and enables preoperative image data to be aligned with the current configuration of the brain of the patient. We have shown that even a relatively complex biomechanical model can be initialized and solved during neurosurgery, providing enhanced surgical visualization. Ultimately, such approaches may provide a truly integrated, multimodality environment for surgical navigation and planning.

36

3 Physics-Based Regularization with an Empirical Model of Anatomical Variability An important issue in nonrigid registration for computer-assisted neurology and neurosurgery is the generation of deformation fields that reflect the transformation of an image in a realistic way with respect to the given anatomy. Due to lack of image structure, noise, intensity artifacts, computational complexity and a restricted time frame (e.g. during surgery), it is not feasible to measure directly the deformation occuring at each voxel. This leads to estimates of the deformation field only at sparse locations which have to be interpolated throughout the image. Recently, physics-based elastic and viscous fluid models for nonrigid registration have become popular (Bajcsy and Kovacic, 1989), since they have the potential to constrain the underlying deformation in a plausible manner. However viscous fluid models (Lester et al., 1999; Wang and Staib, 2000) have to be chosen carefully, since they allow large deformations. This is not always suitable for medical applications concerning the brain. Furthermore, viscous fluid models driven by alignment of similar gray values may allow anatomically incorrect matches of different but adjacent structures through the same mechanism by which large-deformation matches are permitted. For example, one gyrus may flow from the source brain to match two or more different gyri in a target brain, producing an anatomically incorrect match. In terms of physics-based elastic models, recent work has (Davatzikos, 1997; Ferrant et al., 2000b) proposed an active surface algorithm computed at the boundary of a regarded structure as an initial estimate of the deformation field which was then introduced into a volumetric elastic model to infer the deformation inside and outside the surface. A drawback of this method is that although it has been shown to be accurate close to the object’s boundary, away from the boundaries the solution could be less accurate. The work by (Wang and Staib, 2000) represents an improvement in that statistical shape information (based on a set of images with manually-

37

identified boundary points) was included as an additional matching criterion. Even though such methods are promising for specific brain structures, a robust 3D shape representation of the whole brain still remains difficult to achieve. In (Collins, 1994) another nonrigid registration algorithm was proposed, based on an iterative refinement of a local similarity measure using a simplex optimization. In this approach the deformation field was constrained only by smoothing after correspondence estimation, and thus can only be accurate for specific regions of the brain. To achieve better results, the method was improved by introducing various gyri and sulci of the brain as geometric landmarks (Collins et al., 1998a). In order to obtain realistic deformations, we propose here a physics-based elastic model. The method does not require a segmentation and does not have the drawback that initial estimates of the deformation are only generated for the boundary of a considered structure. Instead, these estimates are calculated based on a template matching approach with a local similarity measure. Furthermore we have incorporated a model for inhomogeneous elasticities into our algorithm. The discretization of the underlying equation is done by a finite element technique, which has become a popular method for medical imaging applications (e.g. see (Bro-Nielsen, 1998) and (Ferrant et al., 2000b)).

3.1

Method

The process of registration can be described as an optimization problem that minimizes the deformation energy between a template and a reference image. Assuming that both images represent the same physical object, the deformation that aligns them is therefore related to the theorem of minimum potential energy. The idea of our registration process can now be described as follows: based on a set of points extracted out of an image as described in (3.2), an initial sparse estimate of the deformation field is found by a local normalized cross-correlation

38

(3.3). In a next step nonrigid registration is performed using an elastic model (3.4) which is constrained by the sparse estimates computed in the previous step.

3.2

Feature Point extraction

Let Ω denote the domain of a volume S : Ω ;=< with voxel positions x



x  y z 

-  x>

Ω. In

a first step a set of feature points is extracted out of the reference image. For that purpose we calculate the gradient magnitude out of blurred image intensities. In order to obtain suitable feature points for an initial sparse estimate of the deformation field, only voxel higher than two standard deviations above the mean of the magnitude of the gradient are used for the correspondence detection (3.3). Figure 6 shows this process for one slice of a Magnetic Resonance (MR) scan of the brain. To overcome the poor edge-preserving properties of linear low-pass filters, we use a nonlinear diffusion filter. A filtered version p of volume S can be described as a solution of the partial differential equation (PDE):



div ? g  ∇pσ  2  ∇p @

∂t p

(24)

with Neumann boundary conditions and the original image as initial state (Weickert, 1997). The diffusion function g :