Reconstruction Algorithms
THE APPLICATION OF HYBRID BEM/FEM METHODS TO SOLVE ELECTRICAL IMPEDANCE TOMOGRAPHY'S FORWARD PROBLEM FOR THE HUMAN HEAD J Sikora1a, S.R. Arridgeb, and R.H. Bayfordc, L. Horeshb a
Institute of Theory of Electrical Engineering, Measurement and Information Systems, Warsaw Univeristy of Technology, Koszykowa 75, 00-661 Warsaw, Poland b UCL, London c Middlesex University, London
ABSTRACT: The forward problem in Electrical Impedance Tomography (EIT) requires an accurate estimation of its solution. Although analytical solutions exist for a limited number of cases, there is no suitable analytical solution for the geometry of the human head. A commonly adopted method is the use of Finite Element Methods (FEM). However, the human head consists of layers, which include the scalp and CSF, which are approximately 1 mm thick. Meshing these layers is technically difficult and introduces significant errors in to the solution. An interesting alternative is the use of hybrid methods combining both BEM and FEM. The Boundary Element Method (BEM) is an alterative to FEM. This method essentially uses Green's theorem to map boundary distributions to boundary distributions, using the explicit form of the Green's functions. Although BEM has the disadvantage of dense matrices, and assumes the region is homogenous, it would be of great advantage for the thin layer of the human such as the scalp and CSF, which are homogenous. The larger regions are better suited to FEM, which are composed of non-uniform conductivity distribution. A logical approach for this forward problem is to combine the BEM and FEM to produce hybrid code, which assumes some region, have contact values. Keywords: BEM, FEM, Hybrid methods, EIT
1. INTRODUCTION The use of FEM to model the forward problem for geometries of the human head has limitation due to the thin shells [1]. However, this problem is ideally suited to a hybrid solution. The advantages of FEM and BEM coupling has been investigated extensively in several engineering fields, such as geomechanics[2], and electromagnetics [3] and there are several different methods of coupling BEM and FEM [3,4,5]. This paper presents the results of simple 2D domain and shows the potential of this method. 1.1. FEM-BEM coupling This problem is closely related to the multi-region problem of the BE method such as presented in Figure 1. The multi-region analysis has to fulfil continuity and equilibrium conditions along the interface line Γ1 between Ω1 and Ω 2 regions. This results in the following two relationships:
Φ (Γ1i ) = Φ (Γ2i ) and
∂Φ (1) ∂n
=− Γi
∂Φ ( 2 ) ∂n
(1) Γi
Figure 1: The multi-region boundary element analysis (left) and hybrid discretization (right).
1
E-mail:
[email protected];
503
XII ICEBI – V EIT 2004 Gdańsk
Let the sub-region Ω1 be discretized by the Finite Elements and the Ω 2 by the Boundary Elements. Along the common interface has to be satisfied Eqn.(1). Continuity of the state function Φ can be maintained by using the same order of basis functions in both FE and BE formulations. Thus, if a three-node isoperimetric quadratic boundary element is used an equivalent finite element such as for example eight node quadrilateral quadratic element or six nodes isoperimetric triangle has to be used for the finite element approximation. The essentiality of the problem lies in the fact that the interpolation for the derivatives of the potential for the FEM lies one order lower than the order of the potential itself, whereas for the BEM formulation developed here, the interpolation functions has the same order not only for the potential but also for its derivatives. Such unequal interpolation of the normal derivatives on the interface implants an error to the resulting system of equations. Because along the interface the continuity and equilibrium conditions have to be fulfilled for the FE approach we have assume that on the interface we have additional unknown flux expressed by Neumann boundary conditions. Normally to solve FE systems the boundary conditions have to be imposed allow us to solve the system of equations. So now the FE system of equations in its matrix form could be expressed as follows:
A ( FE ) Φ ( FE ) = B ( FE ) where Φ
( FE )
∂Φ ( FE ) +F ∂n
in
Ω FE
(2)
∂Φ ( FE ) and are column matrices containing the nodal values for the potential (electric ∂n
potential or photon density) and its normal derivatives (electric current or current photon). The corresponding boundary integral equation for the BEM sub-domain is given by:
A
where Φ ( BE ) and
( BE )
Φ
( BE )
=B
( BE )
∂Φ ( BE ) +q ∂n
on
ΓBE
(3)
∂Φ ( BE ) are column matrices containing the nodal values for the potential and its ∂n
normal derivatives for the BEM sub-region.
1.2. Governing equations The boundary-value problem for the FE sub-region is defined by the second-order differential equation [7]: − ∇ ⋅ D ( r )∇Φ ( r ) + k 2 Φ ( r ) = 0 (4) in conjunction with the boundary conditions:
Φ = Φ0 ∂Φ =Ψ ∂n
on
Γi
on
Γi
(5)
where D(r ) is the conductivity (EIT) or diffusion coefficient (OT), k tends to zero in the case of Laplace’s equation (EIT). The equivalent variational problem for the boundary-value defined above is given by: δF (Φ) = 0 (6) in conjunction with the boundary conditions Eqn.(5), form where:
F (Φ ) =
1 2 [ D ∇Φ + k 2 Φ 2 ]dΩ + ∫ ΦΨ dΓi ∫ 2 Ω ( FE ) Γi
(7)
To discretize the functional Eq.(7), the FE sub region is divided into elements ME with M nodes and interface boundary Γi (see Fig. 1) is broken into MSi segments with Mi number of nodes. Usually, M is
504
Reconstruction Algorithms
much larger than Mi. Within each area element, for example six-node triangle, the field is expressed and as:
{ } {Φ } = {Φ } {N }
6
Φ ( x, y ) = ∑ N ie ( x, y )Φ ie = N e
T
i T
e
e
(8)
i =1
and on each line segment on the interface the field is expressed as:
{ } {Φ } = {Φ } {N }
3
Φ i ( x, y ) = ∑ N ie ( x, y )Φ ie = N e
T
i T
e
i
(9)
i =1
Assuming that the interface boundary Γi is a smooth contour, the normal derivative, which is Ψ , is well defined at each interface node and therefore can also be expressed as:
{ } {Ψ } = {Ψ } {N }
Ψ i ( x, y ) = ∑ N ii ( x, y )Ψii = N i
T
i T
i
i
(10)
This is a weak point of this approach because Φ and its normal derivatives Ψ are approximated by the same shape functions. Substituting Eqn. (8-10) into Eqn. (7), we obtain: T
F=
{ }
{ }
{ }
{ }
MSi T 1 M e e e E Φ [ A ] Φ + Φ i [B i ] Ψ i ∑ ∑ 2 e =1 i =1
(11)
[ ]
where in the 2D space matrix A e will take the form:
⎡⎧ ∂N e ⎫⎧ ∂N e ⎫T ⎧ ∂N e ⎫⎧ ∂N e ⎫T ⎤ e e 2 A = ∫ D ⎢⎨ ⎬⎨ ⎬ +⎨ ⎬⎨ ⎬ ⎥dxdy − ∫ k 0 N N ⎢⎣⎩ ∂x ⎭⎩ ∂x ⎭ ⎩ ∂y ⎭⎩ ∂y ⎭ ⎥⎦ Ωe Ωe
[ ]
{ }{ } dxdy
e
[B ] = ∫ Γ {N }{N } dΓ i
i T
i
i
T
(12)
(13)
i
Provided that the element length of the interface is small, the Jacobian of transformation to local coordinate system may be assumed constant and taken out of the integral sign in Eqn.(13) without causing significant errors. Therefore, by substituting the explicit expressions for the shape functions, it is easy to perform the indicated integrations analytically. So the entries of matrix B i in case of the quadrilateral three nodes isoparametric elements of interface in local coordinate system are defined by:
[ ]
⎡ N 1i N ii ⎢ B i = ⎢ N 2i N 1i ⎢ N 3i N 1i ⎣
[ ]
N 1i N 2i N 2i N 2i N 3i N 2i
N 1i N 3i ⎤ ⎡ 4 2 − 1⎤ 1 ⎢ i i⎥ N 2 N 3 ⎥ J (ξ ) = ⎢ 2 16 − 2⎥⎥ J (ξ ) 15 ⎢⎣− 1 2 − 4⎥⎦ N 3i N 3i ⎥⎦
(14)
Than, performing the assembly, Eqn. (11) can be written as:
F=
[
]
[
]
1 {Φ}T A ( FE ) {Φ} + {Φ}T B ( FE ) {Ψ} 2
(15)
where A ( FE ) is an MxM square matrix, B ( FE ) is an MxM i rectangular matrix, Φ is a column vector representing the nodal value of electric potential or photon intensity and Ψ is a column vector representing the nodal values of
∂Φ on the M i node of the interface. ∂n
2. RESULTS A numerical example is present in figure 2 with the FEM region discretized with 200 elements and 440 nodes (bandwidth equal to 14) and BEM region discretised by 104 elements and 207 nodes. Figure also shows the solution on the perimeter of the region and along the interface. Figure 3 shows internal solutions and the relative error.
505
XII ICEBI – V EIT 2004 Gdańsk
Figure 2: Discretization for the BEM-FEM (left) Solution for the perimeter and interface (right).
Figure 3:Amplitude of the internal solution for BEM-FEM (left) Relative error (right).
3. CONCUCLSIONS The paper presents a simple 2-D coupled domain with only one discontinuous region, which demonstrates the potential for this approach. Further work is in progress in applying this technique to the 3-D problem and dealing with multi shell domain, which forms a more realistic model.
4. REFERENCES 1.
2. 3.
4.
5.
506
Bagshaw AP, Liston AD, Bayford RH, Tizzard A, Gibson AP, Tidswell AT, Sparkes MK, Dehghani H, Binnie CD and Holder DS. Electrical Impedance Tomography of human brain function using reconstruction algorithms based on the finite element method. NeuroImage 2003 G. Beer. Programming the boundary Element Method. An Introduction for Engineers. John Wiley and Sons, 2001. CP Bradley, GM Harris and AJ Pullan. The computational performance of a high-order coupled fem-bem produce in elekropotential problems. IEEE Transactions on Biomedical Engineering, 48(11):1238-1250, November 2001. S Kurz and S Russenschuck. Accurate calculation of magnetic fields in the end regions of superconducting accelerator magnets using the bem-fem coupling method. Proceeding of the 1999 Particle Accelerator Conference. New York.1999. OK Panagouli and PD Panagiotopulos. The fem and bem for fractal boundaries and interfaces, applications to unilateral problems. Computer and Structures, 64(-14):329-339, 1997.
