Finite Element Modeling of a Realistic Head Based on by uyk41809


									              International Journal of Bioelectromagnetism                                   
              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
                             College of Automation, Nanchang Hangkong University, P.R. China
                                College of Electrical Engineering, Zhejiang University, P.R. China
                            Department of Biomedical Engineering, University of Minnesota, USA
      Correspondence: Jun Liu, College of Automation, Nanchang Hangkong University, P.R. China. E-mail:

     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

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( 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

                   CT                                              MRI

                                                                  Brain                           Segmentation

        Brain           Scalp, Skull                          CSF, GM, WM

                                 Scalp, Skull, CSF, GM, WM

                                       Finite element model                                      Mesh generation

  Fig. 1 Schematic diagram of the finite element model building procedures.
1. Segmentation

     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

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
                                                                       vi =
                                                                                  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

                                     ⎪       ⎢ ⎛ 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

                                             is a c × n matrix composing by
           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


                                                                 ∑ (u     ik   ) m pk
                                                      viz +1 =   k =1
                                                                  ∑ (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

pixels were labeled as        P (x1 ,y1 ),P2 (x2 ,y2 ),L, PN (xN ,yN ) clockwise or anti-clockwise. Suppose the centers

                                                        N                                   N
                                                  1                             1
of these pixels are      PC ( xc , yc ) ( xc =
                                                      ∑ xi , y c =
                                                       i =1                     N
                                                                                        ∑y  i =1
                                                                                                   i   ), the distances from each pixel to the center

can form a sequence

                                        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.

                                                                ∑ d ( n) d     *
                                                                                   (n + τ ) N
                                         Rd *d (τ ) =           n =1
                                                                       1                          1
                                                                  ⎤ ⎡      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 ⎥
                                      ⎣ ⎦ ⎣ 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)

                                                                                ⎡ x∗ ⎤
                                                          ⎡ x ⎤ ⎡r1 r2     r3 ⎤ ⎢ ∗ ⎥
                                                          ⎢ y ⎥ = ⎢r r          ⎢y ⎥
                                                                           r6 ⎥ ⎢ ⎥
                                                          ⎣ ⎦ ⎣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
                        ⎢                                ⎥ 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

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
     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

 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

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
                                                                                   − Vi A ) 2             ∑ (V
                                                                                                          i =1
                                                                                                                             A 2
                                                                                                                                 )                               (12)


                                                           ∑ (V
                                                                                               S                             A
                                                                                       − V i )(Vi A − V i )
                                       CC =                i =1                                                                                                  (13)
                                                    T                                              T

                                                   ∑ (V                                            ∑ (V
                                                                                       S                                             A 2
                                                                           −V )        i
                                                                                                                         −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

        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

          8                                                                                                       radial
          6                                                                                                       tangential
                                          radial                                  99.5
          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

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

     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

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)

                                                 (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.

                           (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.




                       0         10        20             30             40                 50
                                          Depth of the dipole




                       0         10        20             30             40                 50
                                          Depth of the dipole

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.

     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,
                         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

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.

      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.

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,
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.
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.

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,
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.


To top