Xu2001.pdf

Medical Image Analysis 5 (2001) 271–279 www.elsevier.com / locate / media

Talairach–Tournoux brain atlas registration using a metalforming principle-based finite element method Meihe Xu*, Wieslaw L. Nowinski Laboratory of Neuro Imaging, Room 4238, UCLA, 710 Westwood, Los Angeles, CA 90095, USA Received 24 August 2000; received in revised form 18 January 2001; accepted 28 June 2001

Abstract In this paper, a novel non-rigid registration method is proposed for registration of the Talairach–Tournoux brain atlas with MRI images and the Schaltenbrand–Wahren brain atlas. A metalforming principle-based finite element method with the large deformation problem is used to find the local deformation, in which finite element equations are governed by constraints in the form of displacements derived from the correspondence relationship between extracted feature points. Some detectable substructures, such as the cortical surface, ventricles and corpus callosum, are first extracted from MRI, forming feature points which are classified into different groups. The softassign method is used to establish the correspondence relationship between feature points within each group and to obtain the global transformation concurrently. The displacement constraints are then derived from the correspondence relationship. A metalforming principle-based finite element method with the large deformation problem is used in which finite element equations are reorganized and simplified by integrating the displacement constraints into the system equations. Our method not only matches the model to the data efficiently, but also decreases the degrees of freedom of the system and consequently reduces the computational cost. The method is illustrated by matching the Talairach–Tournoux brain atlas to MRI normal and pathological data and to the Schaltenbrand–Wahren brain atlas. We compare the results quantitatively between the force assignment-based method and the proposed method. The results show that the proposed method yields more accurate results in a fraction of the time taken by the previous method.  2001 Elsevier Science B.V. All rights reserved. Keywords: Brain atlas; Non-rigid registration; Finite element method; Large deformation; Metalforming principle-based

1. Introduction The non-rigid registration for the automatic localization of neuroanatomy in MR images through a 3D deformable atlas is gaining a great deal of attention in medical imaging, which is significant in several applications, including computational anatomy, functional image analysis, image-guided neurosurgery, and model-enhanced neuroradiology. While rigid and affine transformations can just describe global geometric differences between images, elastic schemes can additionally cope with local differences. Inter-subject non-rigid registration is a challenging *Corresponding author. Tel.: 11-310-206-2101; fax: 11-310-2065518. E-mail address: [email protected] (M. Xu).

problem due to the complexity and variability of brain structures. The numerous methods that have been proposed in the literature can be broadly classified into two types, namely similarity-based matching and feature-based matching. Similarity-based matching approaches seek for a transformation between two images that maximizes some measure of overlap or similarity between two volumes (Xu et al., 1999). Feature-based approaches attempt to find the correspondence and transformation using distinct anatomical features that are first extracted from images. These features include points (Bookstein, 1989; Amit, 1997; Davis et al., 1997), curves (Ge et al., 1995), or a surface model (Thompson and Toga, 1996) of anatomical structures. The simplest set of anatomic features is a set of landmarks. Bookstein (1989) used landmark points based on the thin-plate spline (TPS) to generate an elastic

1361-8415 / 01 / $ – see front matter  2001 Elsevier Science B.V. All rights reserved. PII: S1361-8415( 01 )00045-7

272

M. Xu, W.L. Nowinski / Medical Image Analysis 5 (2001) 271 – 279

transformation for registration. This method is inaccurate when the number of control points is insufficient. On the other hand, an increased number of control points results in expensive computation due to interpolation. Later works include the use of physical deformation models to constrain the deformation field using optical flow (Bauchemin and Barron, 1995), elastic (Bajesy and Kovacic, 1989; Davatzikos, 1997), or even viscous fluid deformation models (Bro-Nielsen and Gramkow, 1996). The drawback of such methods is that they either require user intervention or other means to compute the forces exerted on the model. The computation of a deformation field using FEM has become more and more popular for various applications such as surgical simulation and non-rigid registration (Ferrant et al., 1999; McInerney and Kikinis, 1998; Xu and Tang, 1996; Koch et al., 1996; Kyriacou and Davatzikos, 1998). The external force field in these methods is derived by minimizing a local similarity measure. It easily leads to divergence of an FEM solution, if either there is a large variation between the reference image and the target image, or for different modalities. Moreover, these methods are carried out under the assumptions of linear elastic, small deflection behavior. However, due to genetic and life-style factors, there are innate variations among individuals in the appearance and location of brain anatomical structures. The inter-subject registration can be categorized as a large deformation problem. Metaxas et al. (1997) directly utilized the distance between matched points to derive the external forces, in which the force strength is defined by the user. The solution encompassing the expected deformation or displacements is usually unreachable using this method, especially in cases in which there are a lot of matched point pairs in a small region, resulting in instability during fitting. Inter-subject registration has an intuitive resemblance to the metalforming operation in the manufacturing industry. In the metalforming operation, the deformation or displacement in each part of the workpiece is known beforehand. A novel finite element approach, based on the metalforming principle, is proposed here for non-rigid registration. Our method has the following original properties. • The method avoids directly deriving the external forces imposed on the model and integrates the deformationdependent constraint condition to finite element equations. • Rearrangement of the FEM equations decreases the degrees of freedom of the system and consequently reduces the computational cost. • The practical non-linearity in shape has been taken into account. The geometrically non-linear finite element method with large deformation problem is used to solve the non-rigid registration. • Some detectable features in MRI data are first extracted and classified by different substructures. The corre-

spondence relationship of feature points between the atlas and MRI data or other atlas can be built up separately for each extracted substructure, according to the softassign method (Chui et al., 1999). The global affine transformation is acquired from the point set of the cortical surface. • The method works well for cases when the variation between the atlas and the data is large. • The proposed method has been applied to match the Talairach–Tournoux brain atlas (Talairach and Tournoux, 1988) to MRI data and to the Schaltenbrand– Wahren brain atlas (Nowinski et al., 2000; Schaltenbrand and Wahren, 1977). The remaining sections are organized as follows. Section 2 gives the non-linear finite element formula with the large deformation problem. Section 3 describes the metalforming principle-based finite element method. Section 4 discusses implementation issues and results. Section 5 contains the conclusion.

2. Finite element method with the large deformation problem The conventional linear finite element method does not work well for large deformation behavior, which occurs when there is a large difference between the atlas and the data. The non-rigid registration addresses cases in which there is a large variation between the atlas model and MRI data. We first review the notion of the finite element technique with the large displacement problem (Zienkiewicz, 1977). According to the virtual work principle, we can obtain the equilibrium equations between internal and external work: dhd jhc (hd j)j 5

E dh´j hs j dV 2 dhd jhR¯ j 5 0, T

(1)

V

where hc j represents the sum of external and internal generalized forces, hd j is the vector of nodal displacements, hR¯ j represents all the external forces due to imposed loads, and h´ j and hs j denote, respectively, the strain and stress vector. We can write for the variation of strains: dh´ j 5 [B¯ ]dhd j,

(2)

where [B¯ ] is the strain-displacement matrix (see Appendix A). If displacements are large, the strains depend nonlinearly on displacement, and the matrix [B¯ ] is now dependent on hd j. By eliminating dhd j, we have

E [B¯ ] hs j dV 2 hR¯ j 5 0. T

hc (hd j)j 5

(3)

V

If the Newton process is to be adopted, we have to find

M. Xu, W.L. Nowinski / Medical Image Analysis 5 (2001) 271 – 279

the relation between increment dhc j and dhd j. Thus taking the appropriate variation of Eq. (3) with respect to dhd j we have dhc j 5

E d[B¯ ] hs j dV 1E [B¯ ] dhs j dV. T

T

V

[K0 ]hd j 5 hR¯ j.

We assume a constitutive equation, suitable for a small strain, large displacement problem, in the form hs j 5 [D]h´ j,

(5)

where [D] is the usual set of elastic constants, which are related to Young’s modulus of elasticity and Poisson’s ratio. Using (2) and (5), we have

5 [D]([B0 ] 1 [BL ])dhd j,

(6)

where [B0 ] is the same matrix as in linear infinitesimal strain analysis and only [BL ] depends on the displacement (see Appendix A). Finding the derivative of [B¯ ] (i.e., the derivative of BL (d )) by combining the first term of Eq. (4) with Eqs. (A.9), (A.11) and (A.14) in Appendix A, we have

E d[B¯ ] hs j dV 5SE [G] [sˆ ][G] dVDdhd j V

5 [Ks ]dhd j,

(7)

where [G] is associated with the displacement gradient, the nodal displacements and the shape function of the element (see Appendix A), and

3

T

4

,

(8)

where sij are the components of the stress vector and I is the unit matrix. Combining Eq. (8) with Eqs. (4), (6) and (7), we obtain dhc j 5 [KT ]dhd j,

(9)

with [KT ] being the total tangential stiffness matrix, and [KT ] 5 [K0 ] 1 [KL ] 1 [Ks ],

(10)

where [K0 ] represents the usual, small displacements stiffness matrix, i.e., [K0 ] 5

E [B ] [D][B ] dV. T

0

(11)

0

V

The matrix [KL ] is known as the large displacement matrix and is given by

E [B ] [D][B ] 1 [B ] [D][B ] 1 [B ] [D][B ] dV. T

[KL ] 5

0

T

L

L

T

L

(14)

(v) If hci j becomes sufficiently small, the process stops; otherwise go to (ii).

3. Metalforming principle-based finite element technique

3.1. Basis of metalforming analysis

T

s11 I s12 I s13 I [sˆ ] 5 s21 I s22 I s23 I s31 I s32 I s33 I

(13)

(ii) The residual hci j is found using Eq. (3) with the definition of [B¯ ] and strain in Appendix A. (iii) The tangential stiffness matrix [KT ] is estimated and a solution dhd j can be obtained from Eq. (9). (iv) The displacement is updated, hdi 11 j 5 hdi j 1 dhd j.

dhs j 5 [D]dh´ j 5 [D][B¯ ]dhd j

V

In order to find the displacement, an incremental iteration in the Newton method can be applied as follows. (i) The elastic linear solution is obtained as a first approximation hd j using

(4)

V

T

273

L

0

V

(12)

In metalforming analysis, an important feature is that the displacement or deformation on the surface of the workpiece is known while its surface pressures, tractions and overall loads are unknown, prior to processing. In the FEM analysis of the structure problem, it is usually required to find a set of nodal deflections that correspond to a specified set of nodal or distributed loads given that the displacement of supported nodes is zero in particular directions. This is not so in metalforming analysis (Rowe et al., 1991) where deforming loads during the registration are rarely known beforehand and the displacements of certain nodes are directly determined by the boundary conditions or constraints of the problem under consideration. For example, if a node is in contact with a die during a forging operation, the displacement of that node perpendicular to the die surface will be the same as that of the die itself. Although the applied nodal forces are not generally known, the resultant force will be known at many nodes, especially those inside the body and those on the free surface of the workpiece. At these nodes, the resultant force is zero and it is the displacement that must be determined. At the start of the analysis, most of the components of the force vector are known to be zero, but some components are unknown. The components of the displacement vector corresponding to the latter will be specified by the boundary conditions of the problem. The solution of the stiffness equations must therefore take these known components of displacement into account and calculate the corresponding reactions, as well as all the other components of the displacement vector.

M. Xu, W.L. Nowinski / Medical Image Analysis 5 (2001) 271 – 279

274

3.2. Finite element model The Talairach bounding box (Talairach and Tournoux, 1988), encompassing the entire atlas, undergoes an initial global transformation. The transformed bounding box can be assumed to be made of elastic material and partitioned into a 3D grid which constructs a cubical element mesh. The deformable atlas is attached to the FEM mesh. For convenience, we distinguish the boundary condition and constraint condition. The boundary condition denotes the displacements at the frame of the bounding box that have to be set to zero. We refer to the feature points on the atlas model as model points and the feature points in the MRI volume as data points. The constraint condition is specified by the displacements we gather from the model points to data points. In the next subsection, the local displacements are designed to match the corresponding feature points.

3.3. Displacement-dependent constraints In metalforming problems, the effective force acting at all nodes is zero, except for those nodes presumed to be in contact with the dies. In general, the value of the force acting at such a boundary node is not known beforehand, though the value of at least some of the components of its displacement will be. In the same manner, when we use the FEM method to achieve non-rigid registration, the external forces of most nodes are known to be zero at the start of the registration, but that of other nodes (i.e., those that have a common element with feature points) are unknown. The components of displacements corresponding to the latter can be specified by the boundary and constraint conditions of the problem, e.g. the distance of a pair of corresponding feature points can be regarded as the displacement constraints of the FEM system. The solution of the stiffness equations must therefore take these known components of displacement into account and calculate the corresponding reactions as well as all the other components of the displacement. For this purpose, all nodes are classified into constrained nodes, i.e. those having a common element with model points and boundary nodes, and free (unrestricted) nodes. Eq. (12) can be reorganized as

F

K11 K12 K21 K22

GH J H J

d1 R¯ 1 5 ¯ , d2 R2

(15)

[K11 ]hd1 j 5 2 [K12 ]hd2 j.

(17)

The right side of equations can be regarded as external loads. During the process of solution by the Newton method, the tangent stiffness matrix must be rearranged in the same way as Eq. (15). From Eq. (17), we can see that the degrees of freedom of the FEM model is greatly decreased. Consequently, the computational cost of the FEM solution is reduced.

3.4. Displacement assignment The FEM mesh used is generated by cubic elements. Therefore, the nodes of elements cannot correspond to the feature points extracted and the displacement of feature points must be assigned to nodes. We need to distribute the displacement of the feature points to elemental nodes. Computing the displacements from model points to data points and distributing them to the nodes in the volumetric finite elements results in an accurate deformation. The displacements of the deformed substructures are automatically passed on to the neighbors based on the FEM solution. Several neighboring substructures can be deformed concurrently. As described before, the correspondence relationship can be determined using the softassign method separately in each substructure. The displacements are computed for each of the matched point pairs and then distributed to→the nodes. After the global transformation, →we assume P to be a model point (i.e., on the atlas) and Q to be a data point (i.e., in the MRI volume), → → which corresponds to P. The displacement of P is →

→

dP 5Q 2P.

(18) →

We may distribute the displacement of the point P to nodes of the element it locates. →

Nk (P)hdP j hdk j 5 ]]]]] → → , j Nj (P)Nj (P)

O

(19)

→

where d2 is assumed to be the vector of the specified (known) nodal degrees of freedom, and d1 is the vector of the unrestricted nodal degrees of freedom. Then R¯ 2 is the vector of unknown nodal actions and R¯ 1 is the vector of known nodal actions. Eq. (15) can be written as [K11 ]hd1 j 1 [K12 ]hd2 j 5 hR¯ 1 j.

not generally known, the resultant forces are known at many nodes, specifically the unrestricted nodes. At these nodes, the resultant force is zero and it is their displacements that must be determined. Let hR¯ 1 j 5 0; thus

(16)

Although the external forces of the constrained nodes are

where Nk (P) is the shape function that corresponds to the kth node in this element. The subscript j is the number of an arbitrary node in the element. In the implementation, if there are several pairs of feature points in each element, we select one pair, whose distance is the largest. When the displacement of a node, which is the common node of adjacent elements, is distributed by that of these elements, the geometric mean of these displacement vectors is adopted.

M. Xu, W.L. Nowinski / Medical Image Analysis 5 (2001) 271 – 279

4. Implementation and results

4.1. Implementation In the previous sections we mainly discussed the local deformation based on a finite element method with the large deformation problem and metalforming principle. Our registration algorithm relies on the assumption that a large set of macroscopic structures is common between individuals. Fig. 1 shows the details of our implementation. • Feature extraction. After performing a fuzzy segmentation (Ahmed et al., 1999) of the MRI volume, we can obtain thresholds for some detectable tissues, such as the cortical surface, ventricles, and corpus callosum region near the midsagittal plane. We first extract contours of the external surface of the skin (i.e., the intensity of the background is the threshold). Based on the distance constraints and the threshold method, we extract their contours for different tissues and resample the 3D feature points representing these tissues. Thus, the feature points of each group represent the respective substructure. Because our electronic Talairach–Tournoux atlas (Nowinski et al., 1997a,b) is fully colorcoded and labelled, the feature points of each substructure can easily be extracted. • Establishment of points correspondence and global affine transformation. A correspondence relationship can be built up for each segmented substructure using the softassign method (Chui et al., 1999). In this approach, the correspondence and the spatial affine transformation are co-determined from hundreds of feature points of each substructure. In addition, the anatomical variability between subjects can create many outliers, i.e. some feature points which do not match, which are given up. Because different substructures can create different affine transformations, the global transformation is obtained by a point set from the cortical feature. • The global affine transformation and FEM model construction. After the global transformation, we determine a bounding box of the atlas and partition the box into elements. As an initial attempt, the model is assumed to be a uniform, linear elastic material. The displacements of the feature points are distributed to the related elemental nodes. • Stiffness matrix computation, rearrangement and solu-

275

tion. According to the formula described in the previous section, both the small displacement and tangent stiffness matrix can be obtained and rearranged. The displacements of the feature points can be found and distributed to the nodal displacements of the relevant nodes. Based on the Newton method, we obtain the final solution. Suppose that the model is made of isotropic linear elastic material and is uniform, we set Young’s modulus of elasticity to be 0.6 and Poisson’s ratio to be 0.1.

4.2. Application and results We carried out various experiments to test our method. These included matching the Talairach–Tournoux (TT) brain atlas (Talairach and Tournoux, 1988) to normal MRI data, to pathological MRI data and to the Schaltenbrand– Wahren brain atlas (SW) (Schaltenbrand and Wahren, 1977). This work is of vivid clinical and research importance. Our brain atlases are integrated with several commercial systems and are being used clinically, particularly for stereotactic functional neurosurgery (Nowinski et al., 2000). Our Talairach–Tournoux brain atlas is also integrated with Mayo’s ANALYZE version 3.0. Our experiments were performed on an SGI / O2 workstation.

4.2.1. Atlas-to-MRI registration In the experiment of matching the TT atlas to MRI data, the deformable bounding box of the TT atlas is subdivided into 2250 rectangular prism elements and there are 8448 degrees of freedom. After the boundary condition and displacement constraints are added and the FEM equations are rearranged, the number of degrees of freedom is reduced to 5216 and the computation cost is greatly decreased. The deformable finite elements with 4362 model points are matched to 5674 data points (both model points and data points are classified into three groups because they are from different substructures or tissues). The TT atlas was organized into atlas volumes (Nowinski et al., 1997a,b), of which three basic atlas volumes correspond to the orthogonal sections. Because slices of the TT atlas are not continuous, resamplings are obtained by linear interpolation in the MRI volume, after the global affine transformation. Fig. 2 shows the TT atlas matched to normal data. The matching process took approximately 16 minutes. Fig. 3 shows that the TT atlas is matched to pathological MRI data. The first and second row in the

Fig. 1. Flowchart of the algorithm.

276

M. Xu, W.L. Nowinski / Medical Image Analysis 5 (2001) 271 – 279

Fig. 2. Matching of the TT atlas to a normal MRI volume: (a) MRI image; (b) TT atlas; (c) TT superimposed on the MRI space.

figure are from respectively the TT coronal volume and the TT sagittal volume.

4.2.2. Atlas-to-atlas registration The TT atlas and SW atlas (Nowinski et al., 1997a,b, 2000) are two different important atlases. Their combina-

tion by registration has several advantages. First, each has its own strengths and complementary information is merged and provided to the user. Second, direct registration of the SW atlas with MRI data may not be feasible because the features of the SW atlas are very fine and difficult to extract automatically from MRI. In this way,

Fig. 3. Matching of the TT atlas to pathological data: (a) a coronal slice of the MRI volume; (b) a slice of the TT coronal volume after global alignment; (c) atlas superimposed on the MRI slice; (d) a sagittal slice of the MRI volume; (e) a slice of the TT sagittal volume after global alignment; (f) the atlas superimposed on the MRI slice.

M. Xu, W.L. Nowinski / Medical Image Analysis 5 (2001) 271 – 279

any data registered with, for example, the TT atlas is automatically registered with another atlas, for example the SW brain atlas. In our implementation, the SW brain atlas is regarded to be the match target. Unlike matching the atlas to data, we may set up the correspondence between model points and ‘‘data points’’ manually (as explained in (Nowinski et al., 1997a,b), for the thalamic region and the basal ganglia). This preprocessing step is just performed once because both model points and data points come from the TT and SW atlases. There are 16 pairs of feature points to be picked up. As the first step, we perform a 2D registration. After performing the global affine transformation, the bounding box is subdivided into 400 rectangular elements and there are 882 degrees of freedom. Fig. 4 shows the result of registration of the TT atlas with the SW atlas. For displaying the results, a 2D interpolation is obtained by the nearest corner-point method because the TT atlas is labelled by colors.

4.3. Performance and accuracy We synthesize a sample 2D image by deforming 2D MRI data. The image has a size of 2563256. The reference image is divided into 400 elements. There are

277

Table 1 Comparison of results between the linear FEM method with force assignment and the proposed scheme (SGI / O2) Scheme

No. of elements No. of nodes No. of degrees of freedom No. of iterations Time of each iteration (s) Total time (min) Error

Force-assignment based method

Proposed method

400 441

400 441

882 22

536 6

47.2 17.3 112.7

39.4 3.94 14.8

441 nodes in total. The registration is performed from the synthesized image back to the original image (target image). The registration error is estimated by computing the squared sum of the intensity difference between the deformed image of the synthesized image and target image. The squared sum of the intensity difference (after intensity normalization) in our method is 14.8. Table 1 gives a comparison of the performance and accuracy

Fig. 4. Matching of the TT atlas to the SW atlas (the thalamic region and the basal ganglia only, so other regions do not match well): (a) outlines of the SW atlas; (b) the TT atlas; (c) the relative position of slices of the two atlases before registration; (d) image of the SW atlas; (e) TT atlas is superimposed on the SW atlas.

278

M. Xu, W.L. Nowinski / Medical Image Analysis 5 (2001) 271 – 279

between the proposed method and a linear FEM method with force assignment (Metaxas et al., 1997). The number of iterations is referred to the steps of iterations of the process repeated until convergence. We quantitatively compare the results between the force-assignment based method and the proposed method. The results show that the proposed method yields more accurate results in a fraction of the time taken by the previous method. From the table we can see that the cost of each iteration in our method is less than that of the forceassignment based method, because of the reduction in the degrees of freedom, although our method needs computation of the tangent stiffness matrix at each iteration. The error of the registration is much smaller than that of the force-assignment based method.

5. Conclusion The work presented here introduces metalforming principle-based non-rigid registration. In conjunction with a finite element model with large deformation, we demonstrated that we can match efficiently the Talairach–Tournoux brain atlas to MRI normal and pathological data and to the Schaltenbrand–Wahren brain atlas. In the presented method, the expected displacements of the matched points have been set to be constraints of the FEM model. The easily detectable substructures are first extracted from MRI. Then, the correspondence relationship of these feature points and the global transformation are co-determined. Finally, a metalforming principle-based finite element method with the large deformation problem was used, and the FEM equations were reorganized and simplified by integrating the displacement constraints into the system. The results show that the proposed method cannot only efficiently match the model to data, but can also decrease the degrees of freedom of the system and reduce the computational expense. In the future, we will set up an FEM model based on a brain atlas, in which all model points will be nodes of the FEM model and involve updating point correspondences in incremental iterations. The remeshing technique will also be exploited. We hope this work will extend the scope of use of brain atlases for research and clinical applications.

described in terms of the initial particle coordinates X at time t 5 0, in which the particle motion can be written as x 5 x(X,t).

(A.1)

The displacement u may be expressed in terms of X or x, giving x 5 X 1 u(X,t).

(A.2)

The deformation gradient tensor F is defined as

F

G

≠x ≠u F 5]5 I 1] , ≠X ≠X

(A.3)

where I is the unit tensor. The Green–Lagrange strain tensor is defined as 1 ´ 5 ][F T ? F 2 I] 2 1 ≠u ≠u ] 1 ] 5] 2 ≠X ≠X

≠u ≠u FS D S D 1S] D ?S] DG. ≠X ≠X T

T

(A.4)

So

´ 5 ´0 1 ]12 A uu,

(A.5)

where

F

S DG

DS

D

≠u 1 ≠u 2 ≠u 3 ≠u 3 ≠u 2 ≠u 3 ≠u 1 ´0 5 ], ], ], ] 1 ] , ] 1 ] , ≠X1 ≠X2 ≠X3 ≠X2 ≠X3 ≠X1 ≠X3 T ≠u 2 ≠u 1 ]1] (A.6) ≠X1 ≠X2

S

and

u 5 [u T1 , u T2 , u T3 ],

F

≠u 1 ≠u 2 ≠u 3 ui 5 ], ], ] ≠Xi ≠Xi ≠Xi

(A.7)

G

T

,

u1 0 0 0 u3 u2 A u 5 0 u2 0 u3 0 u1 0 0 u3 u2 u1 0

3

(A.8) T

4

.

(A.9)

When the shape function N(X) of an element is given, we have u 5 N(X)d,

(A.10)

where d is the nodal displacement. On the derivative of u, we can express

Acknowledgements

u 5 Gd,

(A.11)

This work was supported by Kent Ridge Digital Labs, 21 Heng Mui Keng Terrace, Singapore 119613.

≠N(X) G 5 ]]. ≠X

(A.12)

Appendix A For convenience, a tensorial representation is adopted. In a material or Lagrangian description, all behavior is

Using Eq. (A.10), Green’s strain can be expressed in terms of nodal displacements d, as

´ 5 [B0 1 ]12 BL (d )]d,

(A.13)

B(d ) L 5 A u G,

(A.14)

M. Xu, W.L. Nowinski / Medical Image Analysis 5 (2001) 271 – 279

where B0 can be determined from (A.6) and (A.8), respectively. Combining with (A.4), (A.11) and du 5 Ndd, we obtain d´ 5 [B0 1 BL (d )]dd 5 [B¯ ]dd,

(A.15)

where B0 is the same matrix as in linear infinitesimal strain analysis and only BL (d ) depends on the displacement. In general, BL will be found to be a linear function of such displacements.

References Ahmed, M.N., Yamany, S.M., Mohamed, N.A., Farag, A.A., 1999. A modified fuzzy C-means algorithm for MRI bias filed estimation and adaptive segmentation. In: Taylor, C. (Ed.), Proceedings of the Second International Conference on Medical Image Computing and ComputerAssisted Intervention, MICCAI’99, Cambridge, UK, pp. 72–81. Amit, Y., 1997. Graphical shape templates for automatic anatomy detection with applications to MRI brain scans. IEEE Trans. Med. Imaging 16 (1), 28–40. Bajesy, R., Kovacic, S., 1989. Multi resolution elastic matching. Comput. Vis., Graph., Image Process. 46, 1–21. Bauchemin, S., Barron, J.L., 1995. The computation of optical flow. ACM Comput. Surv. 27 (3). Bookstein, F.L., 1989. Principal warps: Thin-plate splines and the decomposition of deformations. IEEE Trans. Pan. Anal. Mach. Intell. 11 (6), 567–585. Bro-Nielsen, M., Gramkow, C., 1996. Fast fluid registration of medical images. In: Proceedings Visualization in Biomedical Computing, pp. 267–276. Chui, H., Rambo, J., Duncan, J., Schultz, R., Rangarajan, A., 1999. Registration of cortical anatomical structures via robust 3D point matching. In: Kuba, A. (Ed.), Proceedings of the 16th International Conference on Information Processing in Medical Imaging, IPMI’99. Hungary, pp. 168–181. Davatzikos, C., 1997. Spatial transformation and registration of brain images using elastically deformable models. Comput. Vis. Image Understand. 16 (6), 864–877. Davis, M.N., Khotanzad, A., Flamig, D.P., Harms, S.E., 1997. A physics based coordinate transformation for 3D image matching. IEEE Trans. Med. Imaging 16 (3), 317–328. Ferrant, M., Warfield, S.K., Guttmann, C.R.G., Mulkern, R.V., Jolesz, R.A., Kikinis, R., 1999. 3D image matching using a finite element based elastic deformation model. In: Taylor, C. (Ed.), Proceedings of the Second International Conference on Medical Image Computing

279

and Computer-Assisted Intervention, MICCAI’99, Cambridge, UK, pp. 202–209. Ge, Y., Fitzpatrick, J.M., Kessler, R.M., Jeske-Janicka, M., 1995. Intersubject brain image registration using both cortical and subcortical landmarks. In: Proceedings Image Processing, SPIE 2434, pp. 81–95. Koch, R.M., Gross, M.H., Carls, E.R., von Buren, D.R., Fankhauser, G., Parish, Y.I.H., 1996. Simulating facial surgery using finite element models. In: Proceedings Siggraph 96, New Orleans, LA. Kyriacou, S.K., Davatzikos, C., 1998. A biomedical model of soft tissue deformation with applications to non-rigid registration of brain images with tumour pathology. In: Wells, W.M. (Ed.), Proceedings of the First International Conference on Medical Image Computing and ComputerAssisted Intervention, MICCAI’98, pp. 531–538. McInerney, T., Kikinis, R., 1998. An object-based volumetric deformable atlas for the improved localization of neuroanatomy in MR images. In: Wells, W.M. (Ed.), Proceedings of the First International conference on Medical Image Computing and Computer-Assisted Intervention, MICCAI’98, pp. 861–869. Metaxas, D., Koh, E., Badler, N.I., 1997. Multi-level shape representation using global deformations and locally adaptive finite elements. Int. J. Comput. Vis. 25 (1), 49–61. Nowinski, W.L., Bryan, R.N., Raghavan, R., 1997a. Electronic Clinical Brain Atlas. Thieme, New York. Nowinski, W.L., Fang, A. et al., 1997b. Multiple brain atlas database and atlas-based neuroimaging system. Comput. Aided Surg. 2, 42–66. Nowinski, W.L., Yang, G.L., Yeo, T.T., 2000. Computer-aided stereotactic functional neurosurgery enhanced by the use of the multiple brain atlas database. IEEE Trans. Med. Imaging 19 (1), 1–8. Rowe, G.W., Sturgess, C.E.N., Hartley, P., Pillinger, I., 1991. Finiteelement Plasticity and Metal-forming Analysis. Cambridge University Press, Cambridge. Schaltenbrand, G., Wahren, W., 1977. Atlas for Stereotaxy of the Human Brain. Georg Thieme, Stuttgart. Talairach, J., Tournoux, P., 1988. Co-planar Stereotaxic Atlas of the Human Brain. Thieme Medical, New York. Thompson, P.M., Toga, A.W., 1996. A surface-based technique for warping 3 dimensional images of the brain. IEEE Trans. Med. Imaging 15 (4), 402–417. Xu, M., Tang, Z., 1996. Two computational models for animation of the deformable objects. In: Proceedings of the ACM Symposium on Virtual Reality Software and Technology VRST’96, pp. 151–155. Xu, M.H., Rajagopalan, S., Nowinski, W.L., 1999. A fast mutual information method for multi-modal registration. In: Kuba, A. (Ed.), Proceedings of the 16th International Conference on Information Processing in Medical Imaging, IPMI’99, Visegrad, Hungary, pp. 466–471. Zienkiewicz, O.C., 1977. The Finite Element Method, 3rd Edition. McGraw-Hill, London.