VIEWS: 0 PAGES: 16 CATEGORY: Education POSTED ON: 3/23/2010 Public Domain
International Journal of Bioelectromagnetism http://www.ijbem.org Vol. 10, No. 3, pp. 149 – 164, 2008 Finite Element Modeling of a Realistic Head Based on Medical Images Jun Liua, Shanan Zhub, Yingchun Zhangc, and Bin Hec a College of Automation, Nanchang Hangkong University, P.R. China b College of Electrical Engineering, Zhejiang University, P.R. China c Department of Biomedical Engineering, University of Minnesota, USA Correspondence: Jun Liu, College of Automation, Nanchang Hangkong University, P.R. China. E-mail: liujun1980@sohu.com Abstract. Accurate modeling of head volume conductor plays an important role in EEG/MEG forward and inverse problems. The finite element (FE) modeling of realistic head is a key issue for the finite element analysis of brain electromagnetic field. Various FE modeling methods have been investigated and developed recently to model the realistic head volume conductor, but these methods either use approximation for the tissues, which causes large geometry distortions, or need substantial manual interventions, which makes the modeling process time consuming. In the present study, we have developed a new method to generate subject specific FE head models based on their magnetic resonance (MR) and computer tomographic (CT) image data. The present approach consists of three parts: segmentation of MR and CT images, co-registration of MR and CT images, and mesh generation for all tissue volumes. Numerical experiments using MRI/CT data in two human subjects demonstrate the ability and feasibility of the present FE modeling method to automatically construct realistic head FE models including the scalp, skull, cerebrospinal fluid, gray matter and white matter tissues. Keywords: finite element modeling, image segmentation; image co-registration; finite element generation; EEG, MEG I. Introduction Head modeling plays an important role in electromagnetic forward and inverse problems [B. He, 2005]. While spherical models have been used in finding out the relationship between brain electric sources and electroencephalograms (EEG) and magnetoencephalograms (MEG) [S. Rush et al, 1968; J. P. Ary et al, 1981], studies have shown that more and more details of geometry of the head-brain can provide enhanced performance for EEG/MEG forward and inverse solutions[B. He et al,1987; M. Hamalainen et al, 1989; Y. Yan et al,1991; A. Gevins et al, 1994; K. A. Awada et al ,1997; J. Haueisen et al, 1997; G. Marin et al,1998; J. O. Ollikainen et al, 1999; S. Kim et al,2003; Y. C. Zhang et al, 2004; Y. C. Zhang et al, 2006a; Y. C. Zhang et al, 2006b].. Among the numerical methods, the finite element method (FEM) has the merits of being able to handle both realistic geometry and conductivity inhomogeneity and is numerically efficient. A challenge is the need to construct a subject specific realistic geometry finite element head model, which is needed for FEM computation. Often such procedures are labor and time consuming, preventing FEM head modeling from wide applications. Bonovas et al. segmented head images as the scalp, skull, cerebrospinal fluid （CSF）, gray matter (GM) and white matter (WM), according to the averaged thickness of each tissue. This modeling method is easy to 149 implement, but causes large modeling error due to its approximation of the structure [P. M. Bonovas et al, 2001]. Burguet et al. developed a method for realistic head modeling, containing multiple tissues based on MR images via image segmentation and mesh generation [J. Burguet et al, 2004], which provides good modeling results for the five tissues, i.e. the scalp, skull, CSF, GM and WM, but the model has to be meshed into a large number of elements to improve the element quality. On the other hand, as the segmented skull from MR images has low quality, the accuracy of the model is also influenced. Tizzard et al. developed a modeling method based on a surface model [A. Tizzard et al, 2005], in which the regular bounds of the tissues are first segmented manually, and then meshed by the advancing front method which is based on the surface models as in I-DEAS software. This method has the merits that the shapes of the tissues are smooth and the mesh quality is good, but is limited in that the automatic level is low, so it is inconvenient to repeat for different patients. Several commercial finite element analysis software such as ANSYS（http://www.ansys.com） can be used for finite element model building, but most of them are targeted for industrial engineering design and analysis applications with regular-shape geometry, not providing the ability to model complex realistic geometry head model. In the present study, we introduce an automatic method to generate the finite element (FE) mesh for a realistic geometry inhomogeneous head model. Since MRI has high definition for CSF, GM and WM while CT has high definition for the skull, the present approach uses both MRI and CT for constructing high quality realistic geometry head models. The feasibility of the present approach is evaluated by computer simulations to solve the EEG forward problem. II. Methods In order to get a high quality FE model, both MR and CT images for the same subject were used in the present FE modeling procedure. There are three main functional components illustrated by Fig.1 in our FE modeling procedure: the segmentation component, co-registration component, and mesh generation component. Firstly, the MR images were used to segment the CSF, GM and WM, and the CT images were used to segment the scalp and skull. Then the MR and CT images were co-registered to put all the five head tissues in the same coordinate system. Finally, the mesh generation was implemented to mesh the whole head model with the tetrahedrons. CT MRI Brain Segmentation Brain Scalp, Skull CSF, GM, WM Scalp, Skull, CSF, GM, WM Co-registration Finite element model Mesh generation Fig. 1 Schematic diagram of the finite element model building procedures. 1. Segmentation 150 The brain tissue was firstly extracted automatically from the MR images by using the Brain Extraction Tool (BET) software package [M. S. Smith，2002] and from the CT image by using an adaptive region growing algorithm [H. B. Zhang et al, 2005]. The main work in this segmentation functional component can be divided into two parts, as shown in Fig. 1. In the first part, it is convenient to segment the scalp and skull by using the threshold based segment method, because there are only the scalp and skull tissues left in the CT image after extracting the brain tissue and their gray scale difference is large in CT image. In the second part, the accurate boundaries of the CSF, GM, and WM are difficult to be identified since the gray scales of these three tissues are close and their boundaries are complicated and discontinuous. In this paper, we combined the advantage of K-means (KM) cluster algorithm [J. A. Hartigan et al, 1979]and fuzzy-c means (FCM) algorithm [D. L. Pham et al, 1997] and developed an improved fuzzy-c means (IFCM) algorithm to segment the CSF, GM and WM from extracted MR brain images. The IFCM algorithm is summarized as follows: Step (1) Initialize the parameters, i.e. c, n, m and ε , where c is the number of the clusters and in this paper c＝4 （as there are 4 clusters in the image: GM, CSF, WM and background）, n is the total number of the pixels in the image, m is a fuzzy weighted exponential. A large m may result in a lower convergence speed; however a small m may reduce the accuracy of the segmentation. Normally, m ranged in [1.5, 2.5]. The convergence threshold ε is also an empirical parameter, and we find when it ranged in [0.05, 0.001] we can obtain a reasonable result. In our program, m=2 and ε=0.01. Step (2) Randomly choose the initial cluster centers vi , i = 1Lc , commonly vi is a gray level varies from 0 to 255. Step (3) Class each of all pixels to one of the c cluster centers that has the shortest distance with it, and for each cluster, update the cluster center by M 1 vi = M ∑x k =1 k (1) where x k denotes the gray scale of the pixel belonging to the vi and M the number of the pixels belonging to vi . Step (4) Repeat the operation in Step (3) until all the cluster centers are not changed; Step (5) Set the initial cluster centers of FCM as the results obtained from Step (4); （z） （z +1） Step (6) Update U to U by ⎧ ⎡ c 2 ⎤ −1 ⎪ ⎢ ⎛ d ik ⎞ m −1 ⎥ ⎪u ik = ⎢∑ ⎜ ⎜d ⎟ ⎟ ⎥ , if d ij ≠ 0, j = 1,2,L, c ⎪ ⎢ j =1 ⎝ jk ⎠ ⎥ ⎪ ⎣ ⎦ ⎨ (2) ⎪u ik = 1, if d ij = 0 ⎪ ⎪ ⎪u ik = 0, ⎩ if d rj = 0, r ≠ i; where uik denotes the membership degree of pixel k to vi , d ik = viz − p k , pk is the gray 151 is a c × n matrix composing by （z） degree of pixel k, U uik in cycle z. Step (7) Get the mean vector vi z +1 , i = 1,2,L, c through U ( z +1) as n ∑ (u ik ) m pk viz +1 = k =1 n (3) ∑ (u k =1 ik ) m Step (8) if U ( z) − U ( z +1) ≤ ε , where U ( z ) − U ( z +1) is a Frobenius matrix norm, quit the iterative process and return c sets of fuzzy membership functions, otherwise set z = z + 1 and go back to Step (6). Step (9) A maximum membership segmentation can be obtained by assigning each pixel to the class which has the highest membership value at that location. 2. Co-registration The CSF, GM and WM from the MR image and the scalp and skull from the CT image may not be in the same coordinates because the MR and CT images are acquired with different imaging systems, resulting in different coordinate systems. It is thus necessary to co-register the MR and CT images to put all the five tissues into the same coordinate system. An efficient co-registration algorithm was developed in the present study based on the contour of the images. 2.1 Alignment of the corresponding slices in CT and MRI While the CT and MRI data are 3D volumes, we have taken the strategy to develop a 2D co-registration algorithm, and then perform the alignment of the corresponding slices in CT and MRI volume. Firstly, the CT and MRI volumes are transferred into axial and the top slices in CT and MRI are regarded as the first aligned slices. All the other slices are subsequently determined as all of the slices are parallel and the slice thickness is known. For every slice in MRI, we can find out its corresponding slice in CT by calculating the distance from the current slice to the first slice. As in most of the cases, the slice thickness of MRI and CT may not be the same, so the interpolation of new slices in MRI/CT is also needed. And what should be noted here is as the rotation of the head around the axis through the ears can not be completely identical in the MRI scanner and CT scanner, some corresponding slices we aligned with would not be identical in the incline angle. In order to compensate the error produced by this phenomenon, the affine transformation modal is adopted in the registration as it has the ability to simulate the different incline angle of the slices in MRI and CT. 2.2 Image description by random signal sequences When a pair of corresponding slices was found out, we firstly identified the tissue edges of each compartment from the MR and CT slices by using the morphological operation and Canny edge detect algorithm [J. Canny, 1986]. Then a random pixel was chosen as the starting pixel P ( x1 , y1 ) , and all the other 1 pixels were labeled as P (x1 ,y1 ),P2 (x2 ,y2 ),L, PN (xN ,yN ) clockwise or anti-clockwise. Suppose the centers 1 N N 1 1 of these pixels are PC ( xc , yc ) ( xc = N ∑ xi ， y c = i =1 N ∑y i =1 i ), the distances from each pixel to the center can form a sequence 152 d (i ) = ( xi − xc ) 2 + ( y i − y c ) 2 , i ∈ [1 N ] (4) Notice that the two images to be registered may not have the same number of the edge pixels. In that case, the interpolation or equal distance discretization was performed. Then, we can describe the two corresponding ∗ slices as two random signal sequences d (1), d ∗ (2),L, d ∗ ( N ) and d (1), d (2),L, d ( N ) . 2.3 Objective function construction based on random signal sequences Once the corresponding MR and CT slices are transformed into two random signal sequences, the co-registration between the MR and CT images was achieved by maximizing the objective function (Eq. 5), which is used to evaluate the similarity between the different image methodologies. N ∑ d ( n) d * (n + τ ) N Rd *d (τ ) = n =1 1 1 (5) ⎡N ⎤ ⎡ N ⎤ ⎢∑ d ( n) ⎥ ⎢∑ d ( n + τ ) N ⎥ 2 2 2 *2 ⎣ n=1 ⎦ ⎣ n=1 ⎦ where Rd *d (τ ) is the function that represents the degree of similarity between d (n) and d * (n) , τ is the ⎧n + τ if n + τ <= N spatial shift, and ( n + τ ) N = ⎨ . ⎩ n + τ − N if n + τ > N 2.4 Space transformation matrix construction For any random value of τ ∈ [1, N ] , all the elements in the two random signal sequences will constitute a corresponding aggregate. By iterating those samples in the aggregate, we may obtain a space transformation matrix that satisfies the aggregate. Then we can transform the float image (the image to be transformed in image co-registration) by the space transformation matrix and estimate the co-registration of the mask image (the image to be kept fixed in image co-registration) and the transformed image by Rd *d (τ ) . For any space transformation matrix，we should search all of the values of τ ∈ [1, N ] to get the largest value of Rd *d (τ ) and the matrix that maximizes Rd *d (τ ) is the optimum matrix we are seeking. As the rotation of the head around the axis through the ears can not be completely identical in the MRI scanner and CT scanner, the affine transformation model is adopted in the registration as it has the ability to simulate the different incline angle of the slices in MRI and CT. Let P( x, y) denote the mask image and P ∗ ( x ∗ , y ∗ ) the float image, the affine transformation between the two images can be obtained as ⎡ x ⎤ ⎡ s1 s 2 ⎤ ⎡cos(θ ) − sin(θ )⎤ ⎡ x ∗ ⎤ ⎡Δx ⎤ ⎢ y ⎥ = ⎢s ⎢ ⎥+ s 4 ⎥ ⎢sin(θ ) cos(θ ) ⎥ ⎢ y ∗ ⎥ ⎢Δy ⎥ (6) ⎣ ⎦ ⎣ 3 ⎦⎣ ⎦⎣ ⎦ ⎣ ⎦ Where s1－s 4 denote the scaling in X and Y axes respectively, θ the rotation angle, and Δx and Δy the shift in X and Y axes respectively. Ignoring the physical meaning of Δx , Δy and θ , we can simplify Eq.(6) 153 as: ⎡ x∗ ⎤ ⎡ x ⎤ ⎡r1 r2 r3 ⎤ ⎢ ∗ ⎥ ⎢ y ⎥ = ⎢r r ⎢y ⎥ r6 ⎥ ⎢ ⎥ (7) ⎣ ⎦ ⎣4 5 ⎦1 ⎣ ⎦ In order to use the least squares method, we modify Eq. (7) into another form ART = b ⎡ x1 y1 * * 1 0 0 ⎤0 ⎢ ⎥ ⎢0 0 0 x1 y1 1 ⎥ * * ⎢ x* y * 1 0 0 0 ⎥ ⎢ 2 2 ⎥ R = [r r r r r r ]T where A = ⎢0 0 0 x2 y 2 1 ⎥ * * T 1 2 3 4 5 6 (8) ⎢ ⎥ b = [x1 y1 x2 y 2 L x N y N ]T ⎢M M M M M M ⎥ ⎢ x* y * 1 0 0 0 ⎥ ⎢ N N ⎥ ⎢0 0 0 x* y * 1 ⎥ ⎣ N N ⎦ The least squares solution of Eq. (8) can be given by RT = ( AT A) −1 AT b (9) 3. Finite element mesh generation A mesh generation algorithm was developed, which can automatically generate 3-dimentional high quality finite element meshes from the co-registered MR and CT images, for realistic geometry head model which includes the scalp, skull, CSF, WM, and GM. The tetrahedron was chosen as the basic element type due to its rigorous adaptability to the geometric complexity. The present mesh generation algorithm consists of four steps. (a) (b) （c） Fig. 2 Initialization of the finite element mesh. (a) BCC meshed space; (b) Inner structure; (c) the shape of a typical element. The first step is to generate the initial mesh. A cube enclosing the whole head volume was constructed and meshed by using the Body-Centered Cubic (BCC) lattice [D. Mootz et al, 1981]. The cube was equally divided into small cubes with a user-specified parameter to determine the size of the largest tetrahedron in the FE mesh, then a node is inserted into the center of each cube. After that, two neighboring inserted nodes and two original nodes in a cube edge will compose a tetrahedron so the volume is meshed into many tetrahedrons. Fig.2 shows a 154 BCC meshed volume and its inner configuration. As every tetrahedron is equal in BCC lattice, it has desirable properties in the later tetrahedron refinement. The second step is to re-mesh the tetrahedrons on the tissue boundaries by using the Longest Edge Propagation Path (LEPP) algorithm [M. C. Rivara et al, 1997], to improve the quality of the boundaries in the FE mesh. The LEPP algorithm can be summarized as follows: Step (1): For any tetrahedron t to be re-meshed, find out the last tetrahedron tterminal in its longest edge propagation path. Step (2): Bisect the tterminal by inserting a new node in the middle of its longest edge. Step (3): Find out a next tterminal and repeat Step (2) until all the t terminal tetrahedron are bisected, Step (4): Go back to Step (1) until the tetrahedron t is bisected. The tetrahedron will be bisected by using LEPP in our procedure if it straddles the border of two different tissues. All the tetrahedrons were defined as first degree elements after the initial mesh generation with the BCC lattices. For any element, the degree of the child element will be set to a degree one level higher than the mother element after one re-meshing. The final FE mesh must follow the rule that the degree difference of the neighboring elements should not be more than one level to avoid the severe change of the element size. The third step is to adjust the tetrahedron edges which straddle the border of two different tissues in order to keep the high quality tissue boundaries obtained from the MR and CT images through the image segmentation procedure described in Section 1. Two operations, the node moving and edge splitting, were selectively used to achieve the adjustment. The basic idea of the node moving operation is to move one of the two nodes on the edge which straddles the tissue boundaries to the surface point, as shown in Fig.2(c). The edge N1N2 in the tetrahedron N1N2N3N4 is supposed to straddle tissue boundary and N0 is the cross point of the boundary surface and the edge N1N2. N1 or N2 will be moved to N0 according to the criterion R(t ) f (t ) = 24 ( ) （10） Lmax (t ) where f (t ) , called stretch value, is a target function describing the shape of a tetrahedron t, R(t ) is the radius of the inscribed sphere of t, and Lmax (t ) the length of the longest edge of t. For nodes N1 and N2, the one that leads to larger f (t ) value, which means better tetrahedron quality, for the new tetrahedron is chosen to be moved. Note that the node moving operation should not destroy the topological character of the FE mesh, i.e. the node should not be moved to the other side of its opposite triangle in a tetrahedron element. If this happens in our mesh generation procedure, the edge splitting operation will be used instead of the node moving algorithm. The edge splitting operation is to insert a new node in the location of N0 so that all the tetrahedrons that form the edge N1N2 are bisected. As shown in Fig.2(c), the tetrahedron N1N2N3N4 is bisected into tetrahedrons N1N0N3N4 and N0N2N3N4 when the edge splitting algorithm is implemented. In the 4th step, the edge collapse [J. F. Molinari et al, 2002] and Laplacian smoothing [M. Dorica et al, 2004] operations were used to optimize the quality of the tetrahedral in the FE mesh. The edge collapse operation can 155 be explained as: For any tetrahedron t, if f (t ) < f 0 ( f 0 is a user-specified threshold), then one of the six edges is chosen and its two nodes are moved to one point in the edge according to the criterion of Eq. (10) and the tetrahedron degenerates into a triangle surface of the other tetrahedra. The Laplacian smoothing operation is to move nodes according to Eq. (11), in which N i denotes the coordinate of node i before the moving, N i denotes the coordinate of node i after the moving, N ij denotes ' the neighboring nodes of i, S is the number of the neighboring nodes of i, α is a scale parameter ranging from 0 to 1 and denotes the adjust extent. 1 S N i' = N i + α * ( ∑ N ij −N i ) S j =1 (11) 4. Evaluation In order to test the performance of the present FE modeling procedure, a 3-layer spherical head model was built and meshed by our FE modeling procedure. The EEG forward problem was solved based on this model by using our FEM codes which were verified previously [Y. C. Zhang et al, 2004]. The relative error (RE) and the correlation coefficient (CC) between the FEM results and the analytic results were calculated by using Eq. (12) and (13) to evaluate the accuracy of FE mesh indirectly. T T RE = ∑ (V i =1 i S − Vi A ) 2 ∑ (V i =1 i A 2 ) (12) T ∑ (V S A i S − V i )(Vi A − V i ) CC = i =1 (13) T T ∑ (V ∑ (V S A 2 i S −V ) i 2 i A −V ) i i =1 i =1 where Vi S and Vi A refer to numerical and analytic potentials and T refers to the number of the scalp electrodes, and there are 64 scalp electrodes in the present set-up. We further applied the FE modeling procedure on the MR and CT image data. The finite element model, which includes the scalp, skull, CSF, GM and WM, was built from the MR and CT images by using the present FE modeling procedure and the EEG forward problem was solved based on this FE model by using our verified FEM codes. The FE meshes of different tissue structures in the realistic head model were displayed in different colors and the scalp potential maps from the FEM results with different dipole source configurations were displayed to evaluate the segmentation, co-registration and mesh generation results. III. Results The normalized radius of the scalp, skull and brain was set to 1, 0.92, 0.87 respectively, and the normalized conductivity of the scalp, skull and brain was set to 1.0, 0.05, 1.0 respectively [Y. C. Zhang et al, 2006b]. The 3-layer spherical head model was built and meshed into 199,462 elements and 42,353 nodes automatically by using our FE modeling procedure, as shown in Fig.3(a-b). From the inner to the edge, the size of the tetrahedron changes gradually from the largest to the smallest. This kind of configuration makes the topology of the mesh more reasonable and on the other hand, the quality of the tetrahedron is also guaranteed. A radial or tangential dipole source with different eccentricities was placed in the brain region, and the scalp potential distributions 156 produced by the dipole source were calculated by using the FEM codes, as shown in Fig.3(c-d). The scalp potential distributions were also calculated by analytical method [H. Zhou et al, 1992], and the RE and CC between FEM calculated and analytical scalp potentials were calculated and shown in Fig.3(e-f). We can see that the RE value increases and the CC value decreases when the dipole source is moving from the center to the cortical surface, which is consistent with the FEM calculation literature [Y. C. Zhang et al, 2006a]. We can also see that the FEM results get the maximal RE values of 12.5% and 9%, and minimal CC values of 99.38% and 99.65% for the radial and tangential dipoles, respectively, when the dipole eccentricity reaches 0.8. Considering the fact that the FE model, which consists of only 199,462 elements and 42,353 nodes, is a quite coarse FE mesh, this level of calculation error is reasonable [Y. C. Zhang et al, 2004], which indicates indirectly that the mesh quality of our FE model generated by our FE modeling is acceptable. (a) (b) (c) (d) 14 99.9 12 99.8 10 99.7 CC(%) RE(%) 8 radial 99.6 6 tangential radial 99.5 4 tangential 2 99.4 0 0.2 0.4 0.6 0.8 0 0.2 0.4 0.6 0.8 Eccentricity Eccentricity (e) (f) Fig. 3 A 3-layer spherical model with its application in FEM based potential calculation. (a) The surface of the model; (b) the inner configuration of the model; (c) the potentials produced by a radial dipole; (d) the potentials produced by a tangential dipole; (e) the relative errors between the finite element numerical solution and the analytical solution; (f) the correlation coefficients between the finite element numerical solution and analytical solution. the Then the present FE modeling procedure was used to build the realistic geometry head model from the MR and CT images. The T1-weighted MR images were obtained with a 1.5 T scanner, which composed of 124 157 continuous axial slices each with 1.5 mm thickness. Each slice contained 256×256 pixels, and the field of vision (FOV) was 220×220 (mm). The CT images were obtained from the same subject, which composed of 116 continuous axial slices with slice thickness of 1.25 mm. Each slice contained 512×512 pixels, and the FOV was 250×250 (mm). The segmentation results, shown in Fig.4(a-i)，which include the scalp and skull segmented from the CT images and the CSF, GM and WM segmented from the MR images were obtained automatically by using the present FE modeling procedure. (a) (b) (c) (d) (e) (f) (g) (h) (i) (j) (k) (l) (m) (n) (o) Fig. 4 The segmentation, co-registration and mesh generation procedures from MRI and CT images. (a) Original MRI slice; (b) the brain extracted from MRI; (c) the segmented CSF from MRI; (d) the segmented white matter from MRI; (e) the segmented gray matter from MRI; (f) original CT slice; (g) the brain extracted from CT; (h) the segmented scalp from CT; (i) the segmented skull from CT.; (j) the co-registered image; (k) The meshed scalp with 70,809 elements; (l) the meshed skull with 31,276 elements; (m) the meshed CSF with 19,753 elements; (n) the meshed gray matter with 53,236 elements; (o) the meshed white matter with 38,142 elements. The results of co-registration are shown in Fig.4(j) by a pair of typical corresponding slices (Fig.4(a) and Fig. 4(f)). This typical co-registration result shows that our co-registration algorithm can co-register the MR and CT images for the same patient with about a 2 mm average error and about a 3 mm maximal error, which is calculated according to the maximal distance of the corresponding pixels in the contour of co-registered MRI and CT images. Fig. 4(k-o) shows the five meshed tissues within a FE head model, using the present FE modeling approach, as derived from MR and CT images of a human subject. The head model consists of 213,216 tetrahedrons in total. Compared with the segmentation results in Fig. 4(a-i), we can see that the high quality tissue boundaries 158 obtained from MR and CT images as obtained from segmentation were saved in the final FE mesh. Note that the boundaries of the most of the head tissues, such as the skull, CSF, WM and GM are irregular, complex and sometimes discontinues (for the CSF), but they are all kept in the final FE mesh by using our FE modeling procedure. This result suggests that the present FE modeling procedure has a robust capability of handling the volume model with irregular and arbitrary geometry. We also solved the EEG forward problem using the realistic geometry head model shown in Fig.4(k-o). The original location of the dipole is [0, 0, 60] (mm) which is about in the center of the brain, and then it moves towards the top of the head to [0, 0, 30] (mm) which is about in the boundary of the brain. The conductivity values for the scalp, skull, CSF, GM and WM were set as: 1:0.05:4.62:1:0.43 [R. Van Uitert et al, 2004]. Fig.5 shows the scalp potential distributions calculated in this model when the dipole is located at various depths within the brain. (a) (b) (c) (d) (e) (f) Fig. 5 Potential distributions in the realistic geometry finite element head model (a)-(c) show the scalp potential distributions generated by a radial dipole located at a depth of 60mm, 45mm, and 30mm from the scalp, respectively; (d)-(f) show the scalp potential distributions generated by a tangential dipole located at a depth of 60mm, 45mm, and 30mm from the scalp, respectively. The color bar shows the positive and negative potentials with arbitrary physical units. The present FE modeling procedure was used on another set of MR and CT image data to build another realistic geometry finite element model with the element number of 153,944 and the node number of 33,436, to further test its performance and robustness. The finite difference (FD) model with the node number of 156,233 for the same set of MR and CT image data was also built. Due to the fact that the FD model did not need to include the element information, the element quality of FD model did not have influence on the finite difference method (FDM) calculation results. Therefore, the FDM based EEG forward solution [L. Lemieux et al, 1996; B. Van Rumste et al, 2001] can be used as the relative reference to evaluate the FEM based EEG forward solution, and then to evaluate indirectly the present FE modeling procedure. 159 (a) (b) Fig. 6 A finite element model with 153,944 elements and 33,436 nodes. (a) The surface shape of the model; (b) the inner configuration of the model in which the gray denotes the scalp, the blue denotes the skull, the red denotes CSF, the yellow denotes the gray matter and the green denotes white matter. 12 radial tangential 10 RE(%) 8 6 4 0 10 20 30 40 50 Depth of the dipole (a) 100 99.5 99 CC(%) 98.5 radial 98 tangential 97.5 0 10 20 30 40 50 Depth of the dipole (b) Fig. 7 Relative errors and correlation coefficients between the numerical results using the finite element method and the finite difference method in a realistic geometry head model. (a) The relative errors; (b) the correlation coefficients. 160 Fig. 6 shows the geometry shape of the FE model built from this set of image data by using our FE modeling procedure. The FDM based EEG forward solution and the FEM based EEG forward solution were compared in Fig.7 which shows that the CC and RE between the FDM and FEM solutions change as the dipole source placed in the brain tissue moves from the center to the boundary, and all the RE values are no more than 10.3% (Maximal RE value=10.3% occurs at the location near the boundary of the brain tissue). The statistic analysis of the element quality of the two realistic head FE models built in the present study by using our FE modeling procedure is shown in Table 1, in which the average stretch value is calculated by Eq.(10). Table 1 shows that the average stretch value is around 0.64 for both FE models. Most elements (about 83%) in the two FE models have the average stretch values ranging from 0.5 to 0.7, and only very few numbers of elements (about 0.05%) have the average stretch values below 0.05. Commonly, quantities of elements with stretch values below 0.05 will introduce significant errors in the finite element computation [A. Tizzard et al, 2005]. Table 1 Analysis of models with different numbers of elements Number of Number of Maximal Minimal mean Number of stretch below nodes elements stretch stretch 0.05(%) 33,436 153,944 0.99 0.03 0.637 82(0.053%) 41,368 213,216 0.99 0.04 0.645 96(0.045%) IV. Discussion In the present study, we have developed a new FE modeling approach to generate the 3-dimentional high quality FE head model for solving the EEG forward problem. Compared with the previous FE model building methods, the main contributions of this study are as follows: 1) The present FE modeling procedure is an automatic procedure which is achieved by its three functional parts: image segmentation, co-registration and mesh generation, and can be used to build the realistic geometry head model from the MR and CT images without user interventions. 2) The present FE modeling procedure combines the advantages of the CT and MR imaging modalities to obtain high quality boundaries for different head tissues. 3) The tetrahedron mesh generation algorithm embedded in our FE modeling procedure can automatically generate 3-dimensional finite element mesh for irregular-shaped multiple tissue domains, and is suitable for applications to construct complicated biomedical models. The mesh quality analysis results show that our new FE model building method can generate a realistic FE model which includes multiple irregular-shaped tissue structures, such as the scalp, skull, CSF, GM and WM, with high quality meshes. The simulations were conducted on a 3-layer spherical FE model and a realistic geometry head model built by our FE model method, and the results show that our FE models provide consistent simulations results as assessed by means of the finite difference method. Further more, due to the fact that all of the elements in FDM model are cuboids, so after the segmentation procedure, the accuracy of the FDM model is decided by the number of the elements and has not much to do with the shape of the elements. Therefore, based on the same segmented slices, we have constructed a FDM model with 1615688 elements, and used the FDM based EEG forward solution as the relative reference to evaluate the FEM based EEG forward solution, then to indirectly evaluate the present meshing technique. The CC and RE between the FDM and FEM solutions show that the two results are comparable. One difficulty in biomedical FE modeling research is how to accurately identify different tissue structures from medical images. In our FE modeling procedure, both the MR and CT image data for the same subject are 161 utilized due to the fact that the MR images have high resolution for the CSF, WM and GM, while the CT images have high resolution for the scalp and skull. The MR images are only used to segment the CSF, GM and WM head tissue structures and the CT images are used to segment the scalp and skull. We can see from the segmentation results in Fig.4(a-i) that all the five head tissue structures were easily segmented out with high quality tissue boundaries by using this segmentation protocol. One necessary step in our modeling procedure is the co-registration between the MR and CT images. The MR images and the CT images are fused into one set of images by registration. This step notably improves the accuracy of all the segmented tissues. A limitation of the present study is that we used 2D co-registration instead of 3D co-registration of MR and CT images. However, considering the nature of multiple slices of MR and CT images, with proper transformation, the 2D co-registration procedure is computationally efficient. While further validation studies in a number of clinical data should be conducted in the future, the present evaluation results in two human subjects are encouraging and suggest the applicability of the present FE modeling procedure to construct FE head model for solving the EEG forward problem. Our mesh generation protocol starts from the initial mesh which is uniform mesh with high mesh quality, and all the following element subdivision, which may generate low quality tetrahedrons, strictly follows the rule that the refinement degree difference of the neighbor elements should not be more than one level to avoid the severe change of element size in the finite element mesh. Another way to guarantee the element quality is that the elements have to be re-meshed according to the criterion Eq(10). The model accuracy is improved and the quality of the elements is also guaranteed through the mesh optimization procedure. The total number of the elements can be controlled through the initialization of the lattice and the refinement of the boundary according to user’s requirements. In addition, we optimize the numbering of the nodes iteratively after the model is built, which re-arranges the number of nodes more reasonably, and reduces the bandwidth of the global matrix of the final FE. In summary, we have developed a finite element modeling approach for constructing realistic geometry FE head models from subjects’ MR and CT images. The methods have been tested in two human subjects and satisfactory results obtained, which suggest the applicability of the present method to aid in building finite element head models for solving the EEG forward problem. The method shall also be applicable to aid in MEG forward problem. Acknowledgement We are grateful to Kun Wang for useful discussions and technical assistance. This work was supported in part by NIH RO1EB007920, RO1EB00178, R21EB006070, NSF BES-0411898, NSF-BES-0411480, NSF BES-0602957, and NSFC-50577055. 162 Reference A. Gevins, J. Le, N. Martin and B. Reutter, “High resolution EEG: 124-channel recording, spatial enhancement and MRI integration methods,” Electroenceph. Clin. Neurophysiol, vol.90, pp.337-358, 1994. A. Tizzard, L. Horesh, R. J. Yerworth, D. S. Holder and R. H. Bayford, “Generating accurate finite element meshes for the forward model of the human head in EIT,” Physiol. Meas, vol.26, pp.251-261, 2005. B. He (Ed): “Neural Engineering, ”. Kluwer Academic/Plenum Publishers, 2005. B. He, T. Musha, Y. Okamoto, S. Homma, Y. Nakajima and T. Sato, “Electrical dipole tracing in the brain by means of the boundary element method and its accuracy,” IEEE Trans. Biomed. Eng, vol.34, pp.406-414, 1987. B. Van Rumste, G. Van Hoey, R. Van de Walle, M. Dhave, I. Lemahieu and P. Boon, ”The validation of he finite difference method and reciprocity for solving the inverse problem in EEG dipole source analysis,” Brain Topography, vol.14, pp.83-92, 2001. D. L. Pham, J. L. Prince, A. P. Dagher, and C. Xu, “An Automated Technique for Statistical Characterization of Brain Tissues in Magnetic Resonance Imaging,” International Journal of Pattern Recognition and Artificial Intelligence, vol.11, pp.1189-1211, 1997. D. Mootz and H. G. Wussow, “Crystal structures of pyridine and pyridine trihydrate,” Journal of Chemical Physics, vol.75, pp. 1517-1522, 1981. G. Marin, C. Guerin, S. Baillet, L. Garnero and G. Meunier, “Influence of skull anisotropy for the forward and inverse problem in EEG: simulation studies using FEM on realistic head models,” Human. Brain Mapping, vol.6, pp.250-269, 1998. H. B. Zhang, H. Shen and H. H. Duan, “An automatic and robust algorithm for segmentation of three-dimensional medical images,” Proceedings of the Sixth International Conference on Parallel and Distributed Computing, Applications and Technologies, 2005, pp. 1044-1048. H. Zhou and A. V. Oosterom, “Computation of the potential distribution in a four-layer anisotropic concentric spherical volume conductor,” IEEE Trans. Biomed. Eng. vol.39, pp.154-158, 1992. J. A. Hartigan and M. A.Wong, “A K-Means Clustering Algorithm,” Applied Statistics, vol .28, pp.100-108, 1979. J. Burguet, N. Gadi and I. Bloch, “Realistic models of children heads from 3D-MRI segmentation and tetrahedral mesh construction,” Proceedings of the 2nd International Symposium on 3D Data Processing, Visualization, and Transmission, pp. 631-639, 2004. J. Canny, “A computational approach to edge detection,”. IEEE Transactions on Pattern Analysis and Machine Intelligence. vol.8, pp.679-698, 1986. J. F. Molinari and M. Ortiz, “Three-dimensional adaptive meshing by subdivision and edge-collapse in finite deformation dynamic-plasticity problems with application to adiabatic shear banding,” Int. J. Numer. Meth. Engng, vol.53, pp.1101-1126, 2002. J. Haueisen, C. Ramon, M. Eiselt, H. Brauer and H. Nowak, “Influence of tissue resistivities on neuromagnetic fields and electric potentials studied with a finite element model of the head,” IEEE Trans. Biomed. Eng, vol.44, pp.727-735, 1997. J. O. Ollikainen, M. Vauhkonen, P. A. Karjalainen and J. P. Kaipio, “Effects of local skull inhomogeneities on EEG source estimation,” Med. Eng. Phys, vol.21, pp.143-154, 1999. J. P. Ary, S. A. Klein and D. H. Fender.“Location of sources of evoked scalp potentials: corrections for skull and scalp thicknesses,” IEEE Trans. Biomed. Eng, vol.28, pp.447-452, 1981. K. A. Awada, D. E. Jackson, J. T. Williams, D. R. Wilton and S. B. Baumann, A. C. Papanicolaou, “Computational aspects of finite element modeling in EEG source localization,” IEEE Trans. Biomed. Eng, vol.44, pp.736-752, 1997. L. Lemieux, A. McBride and J. W. Hand, “Calculation of electrical potentials on the surface of a realistic head model by finite difference,” Phys.Med.Biol, vol.41, pp.1079-1091, 1996. M. C. Rivara, “New longest-edge algorithms for the refinement and/or improvement of unstructured triangulations,” Int. J. Numer. Meth. Engng, vol.40, pp.3313-3324. 1997. M. Dorica and D. D. Giannacopoulos, “Toward optimal mesh quality improvements for adaptive finite element electromagnetics with tetrahedral,” IEEE Trans.Magnetics, vol.40, pp.989-992, 2004. M. Hamalainen and J. Sarvas, “Realistic conductor geometry model of the human head for interpretation of neuromagnetic data,” IEEE Trans. Biomed. Eng, vol. 36, pp.165-171, 1989. 163 M. S. Smith, “Fast Robust Automated Brain Extraction,” Human Brain Mapping, vol.17, pp.143-155, 2002. P. M. Bonovas, G. A. Kyriacou and J. N. Sahalos, “A realistic three dimensional FEM of the human head,” Physiol. Meas, vol.22, pp.65-76, 2001. R. Van Uitert, C. Johnson and L. Zhukov, “Influence of head tissue conductivity in forwardand inverse magnetoencephalographic simulations using realistic head models,” IEEE Trans. Biomed. Eng, vol.51, pp.2129-2137, 2004. S. Kim, T. S. Kim,Y. Zhou and M. Singh, “Influence of conductivity tensors on the scalp electrical potential: study with 2-D finite element models,” IEEE Trans. Nuclear Sci, vol.50, pp.133-139, 2003. S. Rush and D. Driscoll, “Current distribution in the brain from surface electrodes,” Anesth. Analg, vol.47, pp.717-723, 1968. Y. C. Zhang, L. Ding, W. van Drongelen, K. Hecox, D. M. Frim and B. He, “A cortical potential imaging study from simultaneous extra- and intracranial electrical recordings by means of the finite element method,” NeuroImage, vol.31, pp. 1513-1524, 2006. Y. C. Zhang, S. A. Zhu and B. He, “A second-order finite element algorithm for solving the three-dimensional EEG forward problem,” Phys. Med. Biol, vol.49, pp.2975-2987, 2004. Y. C. Zhang, W. van Drongelen and B. He, “Estimation of in vivo brain-to-skull conductivity ratio in humans,” Applied Physics Letters, vol.89, pp.2239031-3, 2006. Y. Yan, P. L. Nunez and R. T. Hart, “Finite-element model of the human head: scalp potentials due to dipole sources,” Med. Biol. Eng. Comput, vol.29, pp.475-481, 1991. 164