- Email: [email protected]
A non-iterative linear inverse solution for the block approach in EIT
Journal of Computational Science 1 (2010) 190–196
Contents lists available at ScienceDirect
Journal of Computational Science journal homepage: www.elsevier.com/locate/jocs
A non-iterative linear inverse solution for the block approach in EIT A. Abbasi ∗ , B. Vosoughi Vahdat Department of Electrical Engineering, Sharif University of Technology, Tehran, Iran
a r t i c l e
i n f o
Article history: Received 30 May 2010 Received in revised form 23 September 2010 Accepted 23 September 2010
Keywords: Computer modelling Computational method Inverse problem Numerical solution Electrical impedance tomography
a b s t r a c t Electrical impedance tomography (EIT) is a simple, economic and healthy technique to capture images from the internal area of the body. Although EIT is cheaper and smaller than other imaging systems and requires no ionizing radiation, the resolution associated with this technique is intrinsically limited and the image reconstruction algorithms proposed up to now are not efficient enough. In addition to low resolution EIT is an ill-posed inverse problem. Block method in EIT is based on electrical properties of materials and used to enhance image resolution and also to improve the reconstruction algorithm. Recently an inverse solution for EIT based on block method has been developed, however, this method uses non-linear algorithm. The present article provides a non-iterative linear inverse solution for the block approach on EIT. Using linear equations in this new approach provides a fast algorithm and the ability to solve complicated block problems. We have assumed that the subject has a 2D rectangular shape and is made up of identical fixed size blocks and all of the particles of each block have the same electrical conductivities. It is shown by computer simulations that this linear reconstruction algorithm employing the block method results in an accurate identification. Crown Copyright © 2010 Published by Elsevier B.V. All rights reserved.
1. Introduction There are a variety of medical applications for which it would be useful to know the distribution of electrical properties inside the body. By the term of “electrical properties”, we mean both the electrical conductivity and permittivity, which are of interest in the medical applications [1]. Different tissues have different conductivities and permittivities, on the other hand, the knowledge of the map of the internal electrical properties has a number of advantages in the many of medical diagnosis. EIT is a useful method for medical imaging of pulmonary embolism and blood clots in the lungs [1,2], breast [3], neural system studies [4], breath system studies [5], and other medical issues. Existing reconstruction algorithms for EIT can be categorized as following: (1) (2) (3) (4)
Iterative non-linear algorithms. Iterative algorithms. Layer-stripping algorithms. Direct algorithms.
In iterative non-linear algorithms it is assumed that the conductivity differs smoothly. Examples of iterative non-linear algorithms
∗ Corresponding author. E-mail addresses: [email protected] (A. Abbasi), [email protected] (B.V. Vahdat).
include back-projection methods [6], Calderon’s approach [7], onestep Newton methods [8], and moment methods [9]. However, in some physiological applications the conductivity variations are considerable and iterative non-linear algorithms are less useful in these cases, such as the detection of breast tumours, which are known to be two to four times more conductive than healthy breast tissues [10]. Iterative algorithms may solve the nonlinear problem. Examples of iterative algorithms include methods based on output least-squares [11,12], the equation-error formulation [13,14], or statistical inversion [15]. These algorithms are promising for obtaining accurate reconstructed values; however they often are slow on convergence. In layer-stripping algorithms it is assumed that the higher frequencies in injected current have low penetration depth. This method is introduced and generated from inverse scattering view [16]. Convergence of this algorithm for radially symmetric problems is shown [17], however, there is no proof for its convergence in all problems and also this method is sensitive to measurement noise. Direct algorithms represent the other class of reconstruction algorithms for EIT. They solve the nonlinear problem, so they have the potential of finding the conductivity values with high accuracy. Examples of direct algorithms based on numerical techniques may be found at Mueller’s work [18] or in block method for a rectangular shape subject [19]. Block method is a new approach and in this method reconstructed conductivity values are accurate; however, equations of
1877-7503/$ – see front matter. Crown Copyright © 2010 Published by Elsevier B.V. All rights reserved. doi:10.1016/j.jocs.2010.09.001
A. Abbasi, B.V. Vahdat / Journal of Computational Science 1 (2010) 190–196
191
and (i,j) is specific conductivity in the whole parts of B(i,j) as shown in Fig. 1B. Jx (i,j) and Jy (i,j) are the current densities entering the B(i,j) block from X and Y directions respectively. Similarly, ex (i,j) and ey (i,j) are the voltages of the edge-centres for the B(i,j) block. For B(i,j) block following equations are true [19]:
∀i ∈ N,
1≤i≤m
and
∀j ∈ N,
1≤j≤n
ex (i, j + 1) = ex (i, j) −
1 Jx (i, j) + Jx (i, j + 1) 2 (i, j)
(1)
ey (i + 1, j) = ey (i, j) −
1 Jy (i, j) + Jy (i + 1, j) 2 (i, j)
(2)
Jx (i, j) + Jy (i, j) = Jx (i, j + 1) + Jy (i + 1, j) ex (i, j) − ey (i, j) =
(3)
1 (3Jx (i, j) + Jx (i, j + 1) − 3Jy (i, j) 8 (i, j) −Jy (i + 1, j))
(4)
where is the lengths of the block in both X and Y directions (the transverse and vertical sizes of the block). In forward problem of EIT, (i,j) are known. In this problem if the current densities on the boundaries are known the voltages on the same boundaries could be found. On the other hand if the voltages on the boundaries are known the current densities on the same boundaries could be found. The forward problem is solved as follows:
Fig. 1. (A) A schematic of a rectangular shape subject divided to m × n similar blocks. (B) Block B(i,j) with (i,j) specific conductivity and its current and voltage components.
reconstructed algorithm are nonlinear. The number of equations is very much and solving many nonlinear equations simultaneously encounters difficulties. In this paper, we first present the equations governing the physical problem and then describe the steps of the reconstruction algorithm in Section II. In Section III, we propose more measuring test to have a new set of data and we present linear equations for EIT based on electrical properties of material blocks. Section 4 contains a numerical example and reconstructed conductivity values from both nonlinear and linear methods. Simulation results show the validity and speed of this algorithm. 2. Block model and nonlinear reconstruction solution To generate an EIT image, a series of electrodes are attached to a subject. Various currents can be injected through these electrodes and the produced voltages can be measured. By current injection, voltage measuring and using reconstruction algorithm, conductivity distribution inside the subject would be calculated [20,21]. EIT forward problem involves constructing a block model and calculating the voltages (or currents) produced on the boundary when currents are injected (or voltages are applied) on the same boundary [21]. In block method, a rectangular shape subject is divided into m × n similar size blocks with same electrical impedances in each block [19] as shown in Fig. 1A. In Fig. 1A the rectangular subject has been aligned in Cartesian system and each block has been named as B(i,j) (block in the ith row and jth column of Cartesian system). For a single block, Jx (i,j), Jx (i,j + 1), Jy (i,j) and Jy (i + 1,j) are current density components and ex (i,j), ex (i,j + 1), ey (i,j) and ey (i + 1,j) are voltage components for B(i,j)
1. 2. 3. 4.
Getting row-number and column-number. Generating (i,j) for all blocks. Generating Jx (i,1), Jy (1,j), ex (i,1) and ey (1,j). Calculating Jx (i,j), Jy (i,j), ex (i,j) and ey (i,j) by the steps 2 and 3, using Eqs. (1)–(4). 5. Repeating steps 3 and 4 to generate enough tests. In inverse problem of EIT, currents and voltages of boundaries are known and (i,j) of blocks must be calculated. In nonlinear solution for block method for a subject with m × n blocks (Fig. 1A), in the first row from B(1,1) to B(1,n) we have the following equation: (1, j)ex (1, j) − 2(1, j)ey (1, j) + (1, j)ex (1, j + 1) + Jy (1, j) =0
1≤j≤n
(5)
where ex (1,1), Jx (1,1), ex (1,n + 1), Jx (1,n + 1), ey (1,j) and Jy (1,j) are known. (1,j) and ex (1,j) are the unknown parameters. Eq. (5) shows n equations with n + (n − 1) = 2n − 1 unknown parameters for the first row. If the measurement test is repeated by new values in the first row, different boundary currents and voltages would be resulted. Therefore new n equations with 2n − 1 unknown parameters would be obtained. It should be known that (1,j) are the same in all tests, while ex (1,j) are different in each test. Therefore, each test generates n equations with 2n − 1 unknown parameters where n unknowns are common in all tests. If we try the test for n times, n2 equations and n + n(n − 1) = n2 unknown parameters are obtained. n2 unknown parameters can be solved by numerical methods and (1,j) for 1 ≤ j ≤ n of the first row can be found. For the second row from B(2,1) to B(2,n), we need the boundary values of this row of the n previously achieved tests. Boundary values ex (2,1), Jx (2,1), ex (2,n + 1) and Jx (2,n + 1) are known from the measurements. ey (2,j) and Jy (2,j) are calculated by Eqs. (1)–(4) and (2,j) from n2 equations in the second row can be calculated. To find the parameters of all blocks, this procedure can be repeated for all rows. In nonlinear solution of block method with n tests, n2 equations are available for each row and by solving these equations the conductivities would be obtained. In next section we propose to
192
A. Abbasi, B.V. Vahdat / Journal of Computational Science 1 (2010) 190–196
have more measurement test and solve inverse problem by linear equations.
Writing this equation for j = n and rewriting X(i,j) from previous equations results:
3. Proposed linear inverse solution
X(i, n + 1) = A(i, n)A(i, n−1)· · ·A(i, 6)A(i, 5)A(i, 4)A(i, 3)A(i, 2)A(i, 1) X(i, 1) + A(i, n)A(i, n − 1)· · ·A(i, 6)A(i, 5)A(i, 4)A(i, 3)
In this section we describe the procedure to arrange a set of linear equations and show how these equations can be solved. In inverse solution, currents and voltages of boundaries are known and (i,j) of blocks must be calculated. Currents and voltages of boundaries are obtained from measuring tests. Each test has a new set of currents and voltages of boundaries, but the same (i,j) of blocks. In nonlinear inverse solution, n tests generate n2 equations and solving n2 nonlinear equations acquire (i,j) of blocks in each row. We propose to have 2 × n + 2 measurement tests instead of n tests. By 2 × n + 2 tests and categorising currents and voltages of boundaries in different equations, linear inverse solution would be resulted. The number of measurement tests does not have considerable cost for EIT systems. Therefore, increasing test numbers from n tests to 2 × n + 2 tests is acceptable. Following equation show the equations of a single block (Fig. 1B) as a matrix form. Eq. (6) can be driven from Eqs. (1)–(4):
Jx (i, j + 1) ex (i, j + 1)
−1
=
0
⎡ +
⎣
4(i, j) −1 2
− (i, j)
Jx (i, j) ex (i, j)
−4(i, j) 2
⎤
⎦ Jy (i, j) ey (i, j)
X(i, j) =
U(i, j) =
Jx (i, j) ex (i, j) Jy (i, j) ey (i, j)
−1
A(i, j) =
0
⎡ 2
B(i, j) = ⎣ − (i, j)
(8)
(9)
⎤ ⎦
(10)
where A(i,j) and B(i,j) are unknown parameters containing (i,j) and . is a known parameter and (i,j) is main unknown parameter and the aim of the inverse solution is finding all of (i,j). (i,j) is constant in all measurement tests. Therefore, A(i,j) and B(i,j) are unknowns and constant parameters in all measurement tests. X(i,j) and U(i,j) show currents and voltages of the blocks. U(1,j) and X(i,1) are known parameters, because they are currents and voltages on boundary; however U(i,j) for 2 ≤ i ≤ m and X(i,j) for 2 ≤ j ≤ n are unknown parameters, because they are currents and voltages on inner blocks. Using definitions in Eqs. (7)–(10) onto Eq. (6) yields: X(i, j + 1) = A(i, j)X(i, j) + B(i, j)U(i, j)
B(i, n − 2)U(i, n − 2) + A(i, n)B(i, n − 1)U(i, n − 1) +B(i, n)U(i, n)
(12)
We define P(i,j) as a secondary equation: P(i, 1) = A(i, n)A(i, n − 1)· · ·A(i, 6)A(i, 5)A(i, 4)A(i, 3)A(i, 2)A(i, 1) P(i, 2) = A(i, n)A(i, n − 1)· · ·A(i, 6)A(i, 5)A(i, 4)A(i, 3)A(i, 2)B(i, 1) P(i, 3) = A(i, n)A(i, n − 1)· · ·A(i, 6)A(i, 5)A(i, 4)A(i, 3)B(i, 2) P(i, 4) = A(i, n)A(i, n − 1)· · ·A(i, 6)A(i, 5)A(i, 4)B(i, 3) (13) ··· P(i, n − 1) = A(i, n)A(i, n − 1)B(i, n − 2) P(i, n) = A(i, n)B(i, n−1) P(i, n + 1) = B(i, n)
X(i, n + 1) = P(i, 1)X(i, 1) + P(i, 2)U(i, 1) + P(i, 3)U(i, 2)
(7)
−4(i, j) 2
A(i, 5)A(i, 4)B(i, 3)U(i, 3) + · · · + A(i, n)A(i, n − 1)
By appending of Eqs. (12) and (13) to Eq. (11) we have:
4(i, j) −1
A(i, 4)A(i, 3)B(i, 2)U(i, 2) + A(i, n)A(i, n − 1)· · ·A(i, 6)
(6)
For next step, we define X(i,j), U(i,j), A(i,j) and B(i,j) as auxiliary equations according to the following equations:
A(i, 2)B(i, 1)U(i, 1) + A(i, n)A(i, n − 1)· · ·A(i, 6)A(i, 5)
(11)
+ P(i, 4)U(i, 3) + · · · + P(i, n − 1)U(i, n − 2) + P(i, n)U(i, n − 1) + P(i, n + 1)U(i, n)
(14)
X(i,1), X(i,n + 1) and U(1, j), 1 ≤ j ≤ n are boundary conditions and known parameters from measurement test. Therefore, X(i, n + 1) only depends on P(i, j), 1 ≤ j ≤ n + 1 and P(i,j) includes A(i,j) and B(i,j). On the other hand, Eq. (14) illustrates X(i, n + 1) only depends on A(i,j) and B(i,j)which are unknown and also constant parameters during tests containing (i,j). Nonlinear Eq. (5) is rearranged as a linear Eq. (14). In Eq. (14), P(i, j), 1 ≤ j ≤ n + 1 is unknown parameter in matrix format. The number of unknown parameters in single format is 2 × n + 2. Solution of Eq. (14) would be available by 2 × n + 2 measurement tests. In each measurement test P(i, j), 1 ≤ j ≤ n + 1 is unknown parameter and X(i,1), X(i, n + 1) and U(1, j), 1 ≤ j ≤ n are equation coefficients. If measurement test parameter is added to Eq. (14) following equation appears: X(i, n + 1, tn) = P(i, 1)X(i, 1, tn) + P(i, 2)U(i, 1, tn) + P(i, 3)U(i, 2, tn) + P(i, 4)U(i, 3, tn) + · · · + P(i, n − 1)U(i, n − 2, tn) + P(i, n)U(i, n − 1, tn) + P(i, n + 1)U(i, n, tn)
(15)
tn = test number; 1, 2, 3, . . . , 2n + 2 where tn describes measurement test number. Eq. (15) is the representation of linear relationship between unknown parameters and known coefficients which are obtained by measurement tests. By replacing variables in Eq. (15), it can be rewritten as following:
A. Abbasi, B.V. Vahdat / Journal of Computational Science 1 (2010) 190–196
193
∀1 ≤ i ≤ m and m = number of rows
Jx (i, n + 1, 1) ex (i, n + 1, 1)
Jx (i, n + 1, 2) ex (i, n + 1, 2)
Jx (i, n + 1, 3) ex (i, n + 1, 3)
Jx (i, n + 1, 4) ex (i, n + 1, 4)
... Jx (i, n + 1, 2n + 2) ... ex (i, n + 1, 2n + 2)
⎡
= P(i, 1)2×2
P(i, 2)2×2
Jx (i, 1, 1) ⎢ ex (i, 1, 1) ⎢ e (i, 1, 1) ⎢ y ⎢ Jy (i, 2, 1) ⎢ ⎢ × ⎢ Jy (i, 1, 1) · · · P(i, n + 1)2×2 2×(2n+2) ⎢ ey (i, 2, 1) ⎢ .. ⎢ . ⎢ ⎣ J (i, n, 1) y ey (i, n, 1)
2×(2n+2)
⎤
Jx (i, 1, 2) ex (i, 1, 2) Jy (i, 1, 2) ey (i, 1, 2) Jy (i, 2, 2) ey (i, 2, 2)
Jx (i, 1, 3) ex (i, 1, 3) Jy (i, 1, 3) ey (i, 1, 3) Jy (i, 2, 3) ey (i, 2, 3)
Jx (i, 1, 4) ex (i, 1, 4) Jy (i, 1, 4) ey (i, 1, 4) Jy (i, 2, 4) ey (i, 2, 4)
... Jx (i, 1, 2n + 2) ... ex (i, 1, 2n + 2) ⎥ ⎥ ... Jy (i, 1, 2n + 2) ⎥ ... ey (i, 1, 2n + 2) ⎥ ⎥ ... Jy (i, 2, 2n + 2) ⎥ ⎥ ... ey (i, 2, 2n + 2) ⎥
ey (i, n, 2)
ey (i, n, 3)
ey (i, n, 4)
...
⎥ ⎥ ⎥ Jy (i, n, 2) Jy (i, n, 3) Jy (i, n, 4) ... Jy (i, n, 2n + 2) ⎦ ey (i, n, 2n + 2)
(2n+2)×(2n+2)
(16) where the number of rows and J(i, j, tn) and e(i, j, tn) are the current and the voltage for the block B(i,j). Also P(i, j)2×2 are 2 × 2 matrix as following:
⎡ 2
P(i, n + 1) = ⎣ − (i, n)
−1
P(i, n) =
0
4(i, n) −1 −1
P(i, n − 1) =
0
P(i, n − 2) = .. .
−1
P(i, 2) =
P(i, 1) =
0 −1 0
−1 0
−4(i, n) 2
⎡ ⎣
4(i, n) −1 4(i, n) −1
4(i, n) −1 4(i, n) −1
⎦
2 − (i, n − 1)
−1 0
−1 0
−1
⎤
0 −1 0
−4(i, n − 1) 2
4(i, n − 1) −1 4(i, n − 1) −1
4(i, n − 1) −1 4(i, n − 1) −1
⎡ ⎣
⎦ 2
− (i, n − 2) −1 0
−1
⎤
0 −1 0
−4(i, n − 2) 2
4(i, n − 2) −1
4(i, n − 2) −1 4(i, n − 2) −1
In inverse problem of EIT, currents and voltages of boundaries are known and (i,j) of blocks must be calculated. In Eq. (16) currents and voltages of boundaries are Jx (i, 1, tn), Jx (i, n + 1, tn), ex (i, 1, tn), ex (1, j, tn), Jy (1, j, tn), ey (1, j, tn) for 1 ≤ i ≤ m, m = number of rows and 1 ≤ j ≤ n, n = number of columns and (i,j) of blocks must be calculated by P(i, k)2×2 matrix. Following procedure is proposed for EIT inverse problem solution: (1) Doing 2 × n + 2 independent measuring test and saving boundary conditions. (2) Assuming i = 1 (first row). (3) Using Eq. (16) to find P(i, j)2×2 matrix. (4) Using Eq. (17) to calculate (i,n), (i,n − 1), . . ., (i,1). (5) Calculating Jy (i + 1, j, tn), ey (i + 1, j, tn) for 1 ≤ j ≤ n by Eqs. (1)–(4) and results of step 4. (6) Assuming i = i + 1 up to i = m (going to the next row). (7) Doing steps 3–6 and calculating all (i,j). In next section simulation results of linear and non-linear solution for inverse problem in EIT has been compared.
⎡ ⎣
− (i, n − 3)
−1
···
⎦
2
···
⎤
0 −1 0
4(i, 2) −1
−4(i, n − 3) 2
⎡ ⎣
⎤
(17)
⎦
2
− (i, n − 3) 4(i, 2) 4(i, 1) −1 −1 0 −1
−4(i, 1) 2
⎤ ⎦
4. Numerical results and discussion In this section four different numerical examples has been generated and examined. The first example is a 3 × 3 block with similar conductivity distribution except in the central block. The second example describes a 4 × 4 block with similar conductivity distribution except in two of the central blocks. The third example conditions are the same as the second example with a 5 × 5 block structure. The fourth example represents a 5 × 5 block with nonuniform conductivity distribution. In EIT system, there are electrodes on the surface of the subject. These electrodes can inject the current to the subject and also measure the surface voltages of the subject. For a 3 × 3 block there are 9 blocks, so 12 electrodes should cover the whole boundary blocks’ surface (see Fig. 1A). Therefore a 4 × 4 block is 16 electrode EIT system and a 5 × 5 block is a 20 electrode EIT system. MATLAB version 7.0.0 is used for simulation. For each example, the algorithm includes of three following parts. In the first part, using the EIT forward solution, a phantom model has been constructed with known boundary values. In the second part, using
194
A. Abbasi, B.V. Vahdat / Journal of Computational Science 1 (2010) 190–196
Table 1 Difference between linear and non-linear solution results for four numerical example. Measurement parameter A simple 3 × 3 block
RMSE Algorithm time RMSE Algorithm time RMSE Algorithm time RMSE Algorithm time
A simple 4 × 4 block A simple 5 × 5 block A complicated 5 × 5 block
Non-linear solution results (inverse problem) −8
3.2232776 × 10 1.75 s 2.6140932 × 10−5 5.313 s 0.02404753 14.469 s 2.381684132 15.563 s
the data in the first part, inverse problem is solved by non-linear solution method [19]. In the third part using the data in the first part inverse problem is solved by proposed method on this article (linear method). Then the results of two methods have been compared and discussed in accuracy and time consuming difference. In the first part (forward problem) we have the following steps: 1. 2. 3. 4.
Getting row-number and column-number. Generating random (pre-defined) (i,j) for all blocks. Generating random values for Jx (i,1), Jy (1,j), ex (i,1) and ey (1,j). Calculating Jx (i,j), Jy (i,j), ex (i,j) and ex (i,j) using the values in steps 2 and 3 and also Eqs. (1)–(4). 5. Repeating steps 3 and 4 to generate enough tests (2 × n + 2 independent measuring test).
In the second part (non-linear inverse solution) we have the following steps [19]: 1. Using the boundary values and “MATLAB solve function” for the first row ((1,j) in the first row would be obtained). 2. Using Eqs. (1)–(4) and the results in step 1 the boundary conditions for the second row would be calculated. 3. Using step 1 for the second row and calculating (2,j). 4. Repeating steps 2 and 3 for all of the other rows. In the third part (proposed linear solution) we have the following steps: 1. 2. 3. 4.
Assuming i = 1 (first row). Using Eq. (16) to find P(i, j)2×2 matrix. Using Eq. (17) to calculate (i,n), (i,n − 1), . . ., (i,1). Calculating Jy (i + 1, j, tn), ey (i + 1, j, tn) for 1 ≤ j ≤ n by Eqs. (1)–(4) and results of step 4. 5. Assuming i = i + 1 up to i = m (going to the next row). 6. Doing steps 3–6 and calculating all (i,j). 7. Comparing accuracy and time consuming between two methods. 4.1. Numerical example 1: 3 × 3 block
Linear solution results (inverse problem) 6.4023258 × 10−11 0.04 s 2.8884456 × 10−8 0.156 s 2.0819278 × 10−7 0.438 s 0.016674409 0.437 s
4.2. Numerical example 2: 4 × 4 block In this example a 4 × 4 block with similar conductivity distribution ( = 20) has been considered, however, in two of the central blocks conductivity value is = 200. Table 1 describes both solutions have acceptable RMSE however linear solution is about 33 times faster than non-linear method. 4.3. Numerical example 3: 5 × 5 block with two different block conductivity In this example a 5 × 5 block with similar conductivity distribution ( = 20) has been considered, however, in two of the central blocks conductivity value is = 200. Table 1 shows that both methods have acceptable RMSE however linear solution is about 33 times faster than non-linear method. 4.4. Numerical example 4: 5 × 5 block with non-uniform conductivity distribution Table 1 shows a complicated problem (a 5 × 5 block with nonuniform conductivity distribution) and results of inverse linear and non-linear solutions. As the numerical example gets bigger and complicated, results of inverse solution get worse and RMSE value gets smaller. Fig. 2 shows results of RMSE for the previous examples. For a better view we have used RMSE logarithm to base 10. Although linear and non-linear solutions both get worse as the numerical examples get bigger and complicated, linear solution has considerable better RMSE values than non-linear solution in all examples. EIT inverse solution is a time consuming problem and fast algorithms are more favourite. As it has been shown in Fig. 3 proposed linear solution for EIT is faster than non-linear solution. Also linear solution has a very fast and almost constant algorithm time for different examples, however non-linear solution gets very slow on complicated problems. In most algorithms for EIT inverse problem presented in literature, there is no numerical or quantitative comparison between real and calculated conductivities. However in some papers RMSE has been adopted as a comparison factor [19,22–26]. For example
In this example a 3 × 3 block with similar conductivity distribution ( = 20) has been considered, however, the central block conductivity is = 200. Table 1 shows differences of non-linear and linear inverse solution results. Root mean square error (RMSE) has been used to compare real and calculated values of conductivities.
rmse =
m n i=1
j=1
(r (i, j) − c (i, j)) m×n
2
(18)
where r is the real value of the conductivity and c is the calculated one. Table 1 describes both solutions are accurate and have acceptable RMSE however linear solution is very fast. For this example proposed linear solution is about 44 times faster than non-linear method.
Fig. 2. Difference of RMSE logarithm to base 10 for the previous examples.
A. Abbasi, B.V. Vahdat / Journal of Computational Science 1 (2010) 190–196
195
this method is theoretical, it may have better results than other common ways in literature. Development of this method from 2D to 3D can be suggested as new field for future investigations.
References
Fig. 3. Algorithm time in linear and non-linear solutions for the previous examples. Table 2 RMSE comparison between methods. Algorithm for EIT
RMSE
Shape estimation and state estimation for dynamic EIT [26] Shape estimation by EKF for dynamic EIT [24] Employing IMM, EKF and CCEKF for dynamic EIT [23] Neural networks and front point approach [25] Interpolation of front points method [22] Non-linear block method [19] Proposed linear block method (this paper)
10−1 10−2 10−2 10−2 10−5 10−5 10−8
to 10 to 10−1 to 10−1 to 10−1 to 10−3 to 10−4 to 10−7
in a work conducted by Kim for 2D EIT inverse problem, applying interpolation of front points method for approximation of regions of inner object results in RMS on order of 10−5 to 10−3 for reconstructed image [22]. Also, using neural networks and front point approach for estimation of 2D EIT image results 10−2 to 10−1 in RMSE [25]. Tossavainen et al. applied shape estimation and state estimation formulation for a dynamic EIT problem and reported RMSE from 10−1 to 10 in conductivity [26]. In other work conducted by Ijaz, a dynamic reconstruction algorithm has been presented to monitor the concentration distribution inside the fluid vessel based on EIT by employing interacting multiple model (IMM), extended Kalman filtering (EKF) and covariance compensation extended Kalman filtering (CCEKF). For this study RMSE of conductivity has been reported between 10−2 to 10−1 [23]. Using shape estimation of regions of known resistivities based on extended Kalman filtering for dynamic EIT results in RMSE from 10−2 to 10−1 [24]. Table 2 represents these results. Therefore RMSE from 10−8 to 10−7 for this method is in a proper order and indicates that real and calculated conductivities are almost the same. Although the computer simulations proved that this linear reconstruction algorithm employing the block method results in an accurate identification, the block method has some limitations too. The assumption that: “object has a 2D rectangular shape and is made up of identical fixed size blocks”, causes some difficulties in the practical examples. The block method reconstruction algorithm can be applied for special cases with rectangular shapes. Also it can be used for the other subjects by increasing the number of block (m × n → ∞) which makes more complicated equations. In block method, measurement noise can be propagated easily and influence the results [27]. From the application point of view, error propagation on EIT block method is another limitation of this method. 5. Conclusion The proposed method is an accurate solution for 2D electrical impedance tomography. Low error and high resolution of this method is clear. Fast solution and almost the same algorithm time for complicated problems is the other advantage of linear solution. This method can lead a great step in EIT problem solution. Although
[1] M. Cheney, D. Isaacson, J.C. Newell, Electrical impedance tomography, SIAM Review 41 (1) (1999) 85–101. [2] I. Frerichs, Electrical impedance tomography (EIT) in applications related to lung and ventilation: a review of experimental and clinical activities, Physiological Measurement 21 (2000) 1–21. [3] K.S. Osterman, T.E. Kerner, D.B. Williams, Multifrequency electrical impedance imaging: preliminary in vivo experience in breast, Physiological Measurement 21 (2000) 99–110. [4] N. Polydorides, W.R. Lionheard, H. McCann, Krylov subspace iterative techniques: on the detection of brain activity with electrical impedance tomography, IEEE Transactions on Medical Imaging 21 (6) (2002) 596–603. [5] X. Wu, Y. Sun, J. Cen, Circuit for measuring breath waveform with impedance method and method and device for resisting interference of electrical fast transient, Shenzhen Mindray Bio-Medical Electronics Co., Ltd., United States patent and trademark office pre-grant publication, 2007. [6] F. Santosa, M. Vogelius, A backprojection algorithm for electrical impedance imaging, SIAM Journal on Applied Mathematics 50 (1) (1990) 216–243. [7] M. Cheney, D. Isaacson, An overview of inversion algorithms for impedance imaging, Inverse Scattering and Applications Contemporary Mathematics 122 (1991) 29–39. [8] M. Cheney, D. Isaacson, J.C. Newell, J. Goble, S. Simske, Noser: an algorithm for solving the inverse conductivity problem, International Journal of Imaging Systems and Technology 2 (1990) 66–75. [9] A. Allers, F. Santosa, Stability and resolution analysis of a linearized problem in electrical impedance tomography, Inverse Problems 7 (4) (1991) 515–533. [10] A.J. Suroviec, S.S. Stuchly, J.R. Barr, A. Swarup, Dielectric properties of breast carcinoma and the surrounding tissues, IEEE Transactions on Bimedical Engineering 35 (4) (1988) 257–263. [11] T.J. Yorkey, J.G. Webster, W.J. Tompkins, Comparing reconstruction algorithms for electrical impedance tomography, IEEE Transactions on Bimedical Engineering 34 (11) (1987) 843–852. [12] J.S. Kallman, J.G. Berryman, Weighted least-squares criteria for electrical impedance tomography, IEEE Transactions on Medical Imaging 11 (2) (1992) 284–292. [13] M. Soleimani, Computational aspects of low frequency electrical and electromagnetic tomography: a review study, International Journal of Numerical Analysis and Modeling 5 (3) (2008) 407–440. [14] R.V. Kohn, A. McKenney, Numerical implementation of a variational method for electrical impedance imaging, Inverse Problems 9 (1990) 389–414. [15] J. Kaipio, V. Kolehmainen, E. Somersalo, M. Vauhkonen, Statistical inversion and Monte Carlo sampling methods in electrical impedance tomography, Inverse Problems 16 (5) (2000) 1487–1522. [16] M. Cheney, G. Krislensson, Three dimensional inverse scattering: layerstripping formulae and ill-posedness results, Inverse Problems 4 (1988) 625–642. [17] J. Sylvester, A convergent layer stripping algorithm for radially symmetric impedance tomography problem, Communications in Partial Differential Equations 17 (1992) 1955–1994. [18] J.L. Mueller, S. Siltanen, D. Isaacson, A direct reconstruction algorithm for electrical impedance tomography, IEEE Transactions on Medical Imaging 21 (6) (2002) 555–559. [19] A. Abbasi, B. Vosoughi Vahdat, G. Ebrahimi Fakhim, An inverse solution for 2D electrical impedance tomography based on electrical properties of material blocks, Applied Sciences 9 (10) (2009) 1962–1967. [20] W.R. Lionheart, EIT reconstruction algorithms: pitfalls, challenges and recent developments, Physiological Measurement 25 (2004) 125–142. [21] S Babaeizadeh, D.H. Brooks, D. Isaacson, 3-D electrical impedance tomography for piecewise constant domains with known internal boundaries, IEEE Transactions on Bimedical Engineering 54 (1) (2007) 2–10. [22] M.C. Kim, S. Kim, K.Y. Kim, K.H. Seo, H.J. Jeon, J.H. Kim, B.Y. Choi, Estimation of phase boundary by front points method in electrical impedance tomography, Inverse Problems in Science and Engineering 14 (5) (2006) 455–466. [23] U.Z. Ijaz, J.H. Kim, A.K. Khambampati, M.C. Kim, S. Kim, K.Y. Kim, Concentration distribution estimation of fluid through electrical impedance tomography based on interacting multiple model scheme, Flow Measurement and Instrumentation 18 (1) (2006) 47–56. [24] U.Z. Ijaz, K.Y. Kim, Kinematic models for non-stationary elliptic region boundary in electrical impedance tomography, in: Proceedings of the 2006 International Conference on Hybrid Information Technology (ICHIT 2006), Jeju, Korea, 2006, pp. 9–11. [25] J.H. Kim, B.C. Kang, S.H. Lee, B.I. Choi, M.C. Kim, B.S. Kim, U.Z. Ijaz, K.Y. Kim, S. Kim, Phase boundary estimation in electrical resistance tomography with weighted multi-layered neural networks and front point approach, Measurement Science and Technology 17 (10) (2006) 2731–2739. [26] O.P. Tossavainen, M. Vauhkonen, V. Kolehmainen, K.Y. Kim, Tracking of moving interfaces in sedimentation processes using electrical impedance tomography, Chemical Engineering Science 61 (23) (2006) 7717–7729.
196
A. Abbasi, B.V. Vahdat / Journal of Computational Science 1 (2010) 190–196
[27] A. Abbasi, F. Pashakhanlou, B.V. Vahdat, Error propagation in non-Iterative EIT block method, in: IEEE International Symposium Signal Processing and Information Technology, 2007, pp. 678–681. Ata Abbasi is an assistant professor of Sahand University of Technology, Tabriz, Iran. He received the B.Sc. degree in Biomedical Engineering from Sahand University of Technology and M.Sc. degree from Sharif University of Technology, Iran, in 2003 and 2005 respectively. From 2003 to 2005, he was a research assistant in Sharif University of Technology. He graduated Ph.D. in Biomedical Engineering Sharif University of Technology, 2006–2010.
Bijan Vosughi Vahdat is an assistant professor of Sharif University of Technology, Tehran Iran. He graduated B.Sc. in Electronics, Mashad University, 1983–1987 and M.Sc. in Biomedical Engineering, Sharif University of Technology, 1987–1991 and Ph.D. in Biomedical Engineering, University of New South Wales, 1993–1996. His areas of research interest include Medical Imaging Systems, Advanced Programming and Modelling.
Contents lists available at ScienceDirect
Journal of Computational Science journal homepage: www.elsevier.com/locate/jocs
A non-iterative linear inverse solution for the block approach in EIT A. Abbasi ∗ , B. Vosoughi Vahdat Department of Electrical Engineering, Sharif University of Technology, Tehran, Iran
a r t i c l e
i n f o
Article history: Received 30 May 2010 Received in revised form 23 September 2010 Accepted 23 September 2010
Keywords: Computer modelling Computational method Inverse problem Numerical solution Electrical impedance tomography
a b s t r a c t Electrical impedance tomography (EIT) is a simple, economic and healthy technique to capture images from the internal area of the body. Although EIT is cheaper and smaller than other imaging systems and requires no ionizing radiation, the resolution associated with this technique is intrinsically limited and the image reconstruction algorithms proposed up to now are not efficient enough. In addition to low resolution EIT is an ill-posed inverse problem. Block method in EIT is based on electrical properties of materials and used to enhance image resolution and also to improve the reconstruction algorithm. Recently an inverse solution for EIT based on block method has been developed, however, this method uses non-linear algorithm. The present article provides a non-iterative linear inverse solution for the block approach on EIT. Using linear equations in this new approach provides a fast algorithm and the ability to solve complicated block problems. We have assumed that the subject has a 2D rectangular shape and is made up of identical fixed size blocks and all of the particles of each block have the same electrical conductivities. It is shown by computer simulations that this linear reconstruction algorithm employing the block method results in an accurate identification. Crown Copyright © 2010 Published by Elsevier B.V. All rights reserved.
1. Introduction There are a variety of medical applications for which it would be useful to know the distribution of electrical properties inside the body. By the term of “electrical properties”, we mean both the electrical conductivity and permittivity, which are of interest in the medical applications [1]. Different tissues have different conductivities and permittivities, on the other hand, the knowledge of the map of the internal electrical properties has a number of advantages in the many of medical diagnosis. EIT is a useful method for medical imaging of pulmonary embolism and blood clots in the lungs [1,2], breast [3], neural system studies [4], breath system studies [5], and other medical issues. Existing reconstruction algorithms for EIT can be categorized as following: (1) (2) (3) (4)
Iterative non-linear algorithms. Iterative algorithms. Layer-stripping algorithms. Direct algorithms.
In iterative non-linear algorithms it is assumed that the conductivity differs smoothly. Examples of iterative non-linear algorithms
∗ Corresponding author. E-mail addresses: [email protected] (A. Abbasi), [email protected] (B.V. Vahdat).
include back-projection methods [6], Calderon’s approach [7], onestep Newton methods [8], and moment methods [9]. However, in some physiological applications the conductivity variations are considerable and iterative non-linear algorithms are less useful in these cases, such as the detection of breast tumours, which are known to be two to four times more conductive than healthy breast tissues [10]. Iterative algorithms may solve the nonlinear problem. Examples of iterative algorithms include methods based on output least-squares [11,12], the equation-error formulation [13,14], or statistical inversion [15]. These algorithms are promising for obtaining accurate reconstructed values; however they often are slow on convergence. In layer-stripping algorithms it is assumed that the higher frequencies in injected current have low penetration depth. This method is introduced and generated from inverse scattering view [16]. Convergence of this algorithm for radially symmetric problems is shown [17], however, there is no proof for its convergence in all problems and also this method is sensitive to measurement noise. Direct algorithms represent the other class of reconstruction algorithms for EIT. They solve the nonlinear problem, so they have the potential of finding the conductivity values with high accuracy. Examples of direct algorithms based on numerical techniques may be found at Mueller’s work [18] or in block method for a rectangular shape subject [19]. Block method is a new approach and in this method reconstructed conductivity values are accurate; however, equations of
1877-7503/$ – see front matter. Crown Copyright © 2010 Published by Elsevier B.V. All rights reserved. doi:10.1016/j.jocs.2010.09.001
A. Abbasi, B.V. Vahdat / Journal of Computational Science 1 (2010) 190–196
191
and (i,j) is specific conductivity in the whole parts of B(i,j) as shown in Fig. 1B. Jx (i,j) and Jy (i,j) are the current densities entering the B(i,j) block from X and Y directions respectively. Similarly, ex (i,j) and ey (i,j) are the voltages of the edge-centres for the B(i,j) block. For B(i,j) block following equations are true [19]:
∀i ∈ N,
1≤i≤m
and
∀j ∈ N,
1≤j≤n
ex (i, j + 1) = ex (i, j) −
1 Jx (i, j) + Jx (i, j + 1) 2 (i, j)
(1)
ey (i + 1, j) = ey (i, j) −
1 Jy (i, j) + Jy (i + 1, j) 2 (i, j)
(2)
Jx (i, j) + Jy (i, j) = Jx (i, j + 1) + Jy (i + 1, j) ex (i, j) − ey (i, j) =
(3)
1 (3Jx (i, j) + Jx (i, j + 1) − 3Jy (i, j) 8 (i, j) −Jy (i + 1, j))
(4)
where is the lengths of the block in both X and Y directions (the transverse and vertical sizes of the block). In forward problem of EIT, (i,j) are known. In this problem if the current densities on the boundaries are known the voltages on the same boundaries could be found. On the other hand if the voltages on the boundaries are known the current densities on the same boundaries could be found. The forward problem is solved as follows:
Fig. 1. (A) A schematic of a rectangular shape subject divided to m × n similar blocks. (B) Block B(i,j) with (i,j) specific conductivity and its current and voltage components.
reconstructed algorithm are nonlinear. The number of equations is very much and solving many nonlinear equations simultaneously encounters difficulties. In this paper, we first present the equations governing the physical problem and then describe the steps of the reconstruction algorithm in Section II. In Section III, we propose more measuring test to have a new set of data and we present linear equations for EIT based on electrical properties of material blocks. Section 4 contains a numerical example and reconstructed conductivity values from both nonlinear and linear methods. Simulation results show the validity and speed of this algorithm. 2. Block model and nonlinear reconstruction solution To generate an EIT image, a series of electrodes are attached to a subject. Various currents can be injected through these electrodes and the produced voltages can be measured. By current injection, voltage measuring and using reconstruction algorithm, conductivity distribution inside the subject would be calculated [20,21]. EIT forward problem involves constructing a block model and calculating the voltages (or currents) produced on the boundary when currents are injected (or voltages are applied) on the same boundary [21]. In block method, a rectangular shape subject is divided into m × n similar size blocks with same electrical impedances in each block [19] as shown in Fig. 1A. In Fig. 1A the rectangular subject has been aligned in Cartesian system and each block has been named as B(i,j) (block in the ith row and jth column of Cartesian system). For a single block, Jx (i,j), Jx (i,j + 1), Jy (i,j) and Jy (i + 1,j) are current density components and ex (i,j), ex (i,j + 1), ey (i,j) and ey (i + 1,j) are voltage components for B(i,j)
1. 2. 3. 4.
Getting row-number and column-number. Generating (i,j) for all blocks. Generating Jx (i,1), Jy (1,j), ex (i,1) and ey (1,j). Calculating Jx (i,j), Jy (i,j), ex (i,j) and ey (i,j) by the steps 2 and 3, using Eqs. (1)–(4). 5. Repeating steps 3 and 4 to generate enough tests. In inverse problem of EIT, currents and voltages of boundaries are known and (i,j) of blocks must be calculated. In nonlinear solution for block method for a subject with m × n blocks (Fig. 1A), in the first row from B(1,1) to B(1,n) we have the following equation: (1, j)ex (1, j) − 2(1, j)ey (1, j) + (1, j)ex (1, j + 1) + Jy (1, j) =0
1≤j≤n
(5)
where ex (1,1), Jx (1,1), ex (1,n + 1), Jx (1,n + 1), ey (1,j) and Jy (1,j) are known. (1,j) and ex (1,j) are the unknown parameters. Eq. (5) shows n equations with n + (n − 1) = 2n − 1 unknown parameters for the first row. If the measurement test is repeated by new values in the first row, different boundary currents and voltages would be resulted. Therefore new n equations with 2n − 1 unknown parameters would be obtained. It should be known that (1,j) are the same in all tests, while ex (1,j) are different in each test. Therefore, each test generates n equations with 2n − 1 unknown parameters where n unknowns are common in all tests. If we try the test for n times, n2 equations and n + n(n − 1) = n2 unknown parameters are obtained. n2 unknown parameters can be solved by numerical methods and (1,j) for 1 ≤ j ≤ n of the first row can be found. For the second row from B(2,1) to B(2,n), we need the boundary values of this row of the n previously achieved tests. Boundary values ex (2,1), Jx (2,1), ex (2,n + 1) and Jx (2,n + 1) are known from the measurements. ey (2,j) and Jy (2,j) are calculated by Eqs. (1)–(4) and (2,j) from n2 equations in the second row can be calculated. To find the parameters of all blocks, this procedure can be repeated for all rows. In nonlinear solution of block method with n tests, n2 equations are available for each row and by solving these equations the conductivities would be obtained. In next section we propose to
192
A. Abbasi, B.V. Vahdat / Journal of Computational Science 1 (2010) 190–196
have more measurement test and solve inverse problem by linear equations.
Writing this equation for j = n and rewriting X(i,j) from previous equations results:
3. Proposed linear inverse solution
X(i, n + 1) = A(i, n)A(i, n−1)· · ·A(i, 6)A(i, 5)A(i, 4)A(i, 3)A(i, 2)A(i, 1) X(i, 1) + A(i, n)A(i, n − 1)· · ·A(i, 6)A(i, 5)A(i, 4)A(i, 3)
In this section we describe the procedure to arrange a set of linear equations and show how these equations can be solved. In inverse solution, currents and voltages of boundaries are known and (i,j) of blocks must be calculated. Currents and voltages of boundaries are obtained from measuring tests. Each test has a new set of currents and voltages of boundaries, but the same (i,j) of blocks. In nonlinear inverse solution, n tests generate n2 equations and solving n2 nonlinear equations acquire (i,j) of blocks in each row. We propose to have 2 × n + 2 measurement tests instead of n tests. By 2 × n + 2 tests and categorising currents and voltages of boundaries in different equations, linear inverse solution would be resulted. The number of measurement tests does not have considerable cost for EIT systems. Therefore, increasing test numbers from n tests to 2 × n + 2 tests is acceptable. Following equation show the equations of a single block (Fig. 1B) as a matrix form. Eq. (6) can be driven from Eqs. (1)–(4):
Jx (i, j + 1) ex (i, j + 1)
−1
=
0
⎡ +
⎣
4(i, j) −1 2
− (i, j)
Jx (i, j) ex (i, j)
−4(i, j) 2
⎤
⎦ Jy (i, j) ey (i, j)
X(i, j) =
U(i, j) =
Jx (i, j) ex (i, j) Jy (i, j) ey (i, j)
−1
A(i, j) =
0
⎡ 2
B(i, j) = ⎣ − (i, j)
(8)
(9)
⎤ ⎦
(10)
where A(i,j) and B(i,j) are unknown parameters containing (i,j) and . is a known parameter and (i,j) is main unknown parameter and the aim of the inverse solution is finding all of (i,j). (i,j) is constant in all measurement tests. Therefore, A(i,j) and B(i,j) are unknowns and constant parameters in all measurement tests. X(i,j) and U(i,j) show currents and voltages of the blocks. U(1,j) and X(i,1) are known parameters, because they are currents and voltages on boundary; however U(i,j) for 2 ≤ i ≤ m and X(i,j) for 2 ≤ j ≤ n are unknown parameters, because they are currents and voltages on inner blocks. Using definitions in Eqs. (7)–(10) onto Eq. (6) yields: X(i, j + 1) = A(i, j)X(i, j) + B(i, j)U(i, j)
B(i, n − 2)U(i, n − 2) + A(i, n)B(i, n − 1)U(i, n − 1) +B(i, n)U(i, n)
(12)
We define P(i,j) as a secondary equation: P(i, 1) = A(i, n)A(i, n − 1)· · ·A(i, 6)A(i, 5)A(i, 4)A(i, 3)A(i, 2)A(i, 1) P(i, 2) = A(i, n)A(i, n − 1)· · ·A(i, 6)A(i, 5)A(i, 4)A(i, 3)A(i, 2)B(i, 1) P(i, 3) = A(i, n)A(i, n − 1)· · ·A(i, 6)A(i, 5)A(i, 4)A(i, 3)B(i, 2) P(i, 4) = A(i, n)A(i, n − 1)· · ·A(i, 6)A(i, 5)A(i, 4)B(i, 3) (13) ··· P(i, n − 1) = A(i, n)A(i, n − 1)B(i, n − 2) P(i, n) = A(i, n)B(i, n−1) P(i, n + 1) = B(i, n)
X(i, n + 1) = P(i, 1)X(i, 1) + P(i, 2)U(i, 1) + P(i, 3)U(i, 2)
(7)
−4(i, j) 2
A(i, 5)A(i, 4)B(i, 3)U(i, 3) + · · · + A(i, n)A(i, n − 1)
By appending of Eqs. (12) and (13) to Eq. (11) we have:
4(i, j) −1
A(i, 4)A(i, 3)B(i, 2)U(i, 2) + A(i, n)A(i, n − 1)· · ·A(i, 6)
(6)
For next step, we define X(i,j), U(i,j), A(i,j) and B(i,j) as auxiliary equations according to the following equations:
A(i, 2)B(i, 1)U(i, 1) + A(i, n)A(i, n − 1)· · ·A(i, 6)A(i, 5)
(11)
+ P(i, 4)U(i, 3) + · · · + P(i, n − 1)U(i, n − 2) + P(i, n)U(i, n − 1) + P(i, n + 1)U(i, n)
(14)
X(i,1), X(i,n + 1) and U(1, j), 1 ≤ j ≤ n are boundary conditions and known parameters from measurement test. Therefore, X(i, n + 1) only depends on P(i, j), 1 ≤ j ≤ n + 1 and P(i,j) includes A(i,j) and B(i,j). On the other hand, Eq. (14) illustrates X(i, n + 1) only depends on A(i,j) and B(i,j)which are unknown and also constant parameters during tests containing (i,j). Nonlinear Eq. (5) is rearranged as a linear Eq. (14). In Eq. (14), P(i, j), 1 ≤ j ≤ n + 1 is unknown parameter in matrix format. The number of unknown parameters in single format is 2 × n + 2. Solution of Eq. (14) would be available by 2 × n + 2 measurement tests. In each measurement test P(i, j), 1 ≤ j ≤ n + 1 is unknown parameter and X(i,1), X(i, n + 1) and U(1, j), 1 ≤ j ≤ n are equation coefficients. If measurement test parameter is added to Eq. (14) following equation appears: X(i, n + 1, tn) = P(i, 1)X(i, 1, tn) + P(i, 2)U(i, 1, tn) + P(i, 3)U(i, 2, tn) + P(i, 4)U(i, 3, tn) + · · · + P(i, n − 1)U(i, n − 2, tn) + P(i, n)U(i, n − 1, tn) + P(i, n + 1)U(i, n, tn)
(15)
tn = test number; 1, 2, 3, . . . , 2n + 2 where tn describes measurement test number. Eq. (15) is the representation of linear relationship between unknown parameters and known coefficients which are obtained by measurement tests. By replacing variables in Eq. (15), it can be rewritten as following:
A. Abbasi, B.V. Vahdat / Journal of Computational Science 1 (2010) 190–196
193
∀1 ≤ i ≤ m and m = number of rows
Jx (i, n + 1, 1) ex (i, n + 1, 1)
Jx (i, n + 1, 2) ex (i, n + 1, 2)
Jx (i, n + 1, 3) ex (i, n + 1, 3)
Jx (i, n + 1, 4) ex (i, n + 1, 4)
... Jx (i, n + 1, 2n + 2) ... ex (i, n + 1, 2n + 2)
⎡
= P(i, 1)2×2
P(i, 2)2×2
Jx (i, 1, 1) ⎢ ex (i, 1, 1) ⎢ e (i, 1, 1) ⎢ y ⎢ Jy (i, 2, 1) ⎢ ⎢ × ⎢ Jy (i, 1, 1) · · · P(i, n + 1)2×2 2×(2n+2) ⎢ ey (i, 2, 1) ⎢ .. ⎢ . ⎢ ⎣ J (i, n, 1) y ey (i, n, 1)
2×(2n+2)
⎤
Jx (i, 1, 2) ex (i, 1, 2) Jy (i, 1, 2) ey (i, 1, 2) Jy (i, 2, 2) ey (i, 2, 2)
Jx (i, 1, 3) ex (i, 1, 3) Jy (i, 1, 3) ey (i, 1, 3) Jy (i, 2, 3) ey (i, 2, 3)
Jx (i, 1, 4) ex (i, 1, 4) Jy (i, 1, 4) ey (i, 1, 4) Jy (i, 2, 4) ey (i, 2, 4)
... Jx (i, 1, 2n + 2) ... ex (i, 1, 2n + 2) ⎥ ⎥ ... Jy (i, 1, 2n + 2) ⎥ ... ey (i, 1, 2n + 2) ⎥ ⎥ ... Jy (i, 2, 2n + 2) ⎥ ⎥ ... ey (i, 2, 2n + 2) ⎥
ey (i, n, 2)
ey (i, n, 3)
ey (i, n, 4)
...
⎥ ⎥ ⎥ Jy (i, n, 2) Jy (i, n, 3) Jy (i, n, 4) ... Jy (i, n, 2n + 2) ⎦ ey (i, n, 2n + 2)
(2n+2)×(2n+2)
(16) where the number of rows and J(i, j, tn) and e(i, j, tn) are the current and the voltage for the block B(i,j). Also P(i, j)2×2 are 2 × 2 matrix as following:
⎡ 2
P(i, n + 1) = ⎣ − (i, n)
−1
P(i, n) =
0
4(i, n) −1 −1
P(i, n − 1) =
0
P(i, n − 2) = .. .
−1
P(i, 2) =
P(i, 1) =
0 −1 0
−1 0
−4(i, n) 2
⎡ ⎣
4(i, n) −1 4(i, n) −1
4(i, n) −1 4(i, n) −1
⎦
2 − (i, n − 1)
−1 0
−1 0
−1
⎤
0 −1 0
−4(i, n − 1) 2
4(i, n − 1) −1 4(i, n − 1) −1
4(i, n − 1) −1 4(i, n − 1) −1
⎡ ⎣
⎦ 2
− (i, n − 2) −1 0
−1
⎤
0 −1 0
−4(i, n − 2) 2
4(i, n − 2) −1
4(i, n − 2) −1 4(i, n − 2) −1
In inverse problem of EIT, currents and voltages of boundaries are known and (i,j) of blocks must be calculated. In Eq. (16) currents and voltages of boundaries are Jx (i, 1, tn), Jx (i, n + 1, tn), ex (i, 1, tn), ex (1, j, tn), Jy (1, j, tn), ey (1, j, tn) for 1 ≤ i ≤ m, m = number of rows and 1 ≤ j ≤ n, n = number of columns and (i,j) of blocks must be calculated by P(i, k)2×2 matrix. Following procedure is proposed for EIT inverse problem solution: (1) Doing 2 × n + 2 independent measuring test and saving boundary conditions. (2) Assuming i = 1 (first row). (3) Using Eq. (16) to find P(i, j)2×2 matrix. (4) Using Eq. (17) to calculate (i,n), (i,n − 1), . . ., (i,1). (5) Calculating Jy (i + 1, j, tn), ey (i + 1, j, tn) for 1 ≤ j ≤ n by Eqs. (1)–(4) and results of step 4. (6) Assuming i = i + 1 up to i = m (going to the next row). (7) Doing steps 3–6 and calculating all (i,j). In next section simulation results of linear and non-linear solution for inverse problem in EIT has been compared.
⎡ ⎣
− (i, n − 3)
−1
···
⎦
2
···
⎤
0 −1 0
4(i, 2) −1
−4(i, n − 3) 2
⎡ ⎣
⎤
(17)
⎦
2
− (i, n − 3) 4(i, 2) 4(i, 1) −1 −1 0 −1
−4(i, 1) 2
⎤ ⎦
4. Numerical results and discussion In this section four different numerical examples has been generated and examined. The first example is a 3 × 3 block with similar conductivity distribution except in the central block. The second example describes a 4 × 4 block with similar conductivity distribution except in two of the central blocks. The third example conditions are the same as the second example with a 5 × 5 block structure. The fourth example represents a 5 × 5 block with nonuniform conductivity distribution. In EIT system, there are electrodes on the surface of the subject. These electrodes can inject the current to the subject and also measure the surface voltages of the subject. For a 3 × 3 block there are 9 blocks, so 12 electrodes should cover the whole boundary blocks’ surface (see Fig. 1A). Therefore a 4 × 4 block is 16 electrode EIT system and a 5 × 5 block is a 20 electrode EIT system. MATLAB version 7.0.0 is used for simulation. For each example, the algorithm includes of three following parts. In the first part, using the EIT forward solution, a phantom model has been constructed with known boundary values. In the second part, using
194
A. Abbasi, B.V. Vahdat / Journal of Computational Science 1 (2010) 190–196
Table 1 Difference between linear and non-linear solution results for four numerical example. Measurement parameter A simple 3 × 3 block
RMSE Algorithm time RMSE Algorithm time RMSE Algorithm time RMSE Algorithm time
A simple 4 × 4 block A simple 5 × 5 block A complicated 5 × 5 block
Non-linear solution results (inverse problem) −8
3.2232776 × 10 1.75 s 2.6140932 × 10−5 5.313 s 0.02404753 14.469 s 2.381684132 15.563 s
the data in the first part, inverse problem is solved by non-linear solution method [19]. In the third part using the data in the first part inverse problem is solved by proposed method on this article (linear method). Then the results of two methods have been compared and discussed in accuracy and time consuming difference. In the first part (forward problem) we have the following steps: 1. 2. 3. 4.
Getting row-number and column-number. Generating random (pre-defined) (i,j) for all blocks. Generating random values for Jx (i,1), Jy (1,j), ex (i,1) and ey (1,j). Calculating Jx (i,j), Jy (i,j), ex (i,j) and ex (i,j) using the values in steps 2 and 3 and also Eqs. (1)–(4). 5. Repeating steps 3 and 4 to generate enough tests (2 × n + 2 independent measuring test).
In the second part (non-linear inverse solution) we have the following steps [19]: 1. Using the boundary values and “MATLAB solve function” for the first row ((1,j) in the first row would be obtained). 2. Using Eqs. (1)–(4) and the results in step 1 the boundary conditions for the second row would be calculated. 3. Using step 1 for the second row and calculating (2,j). 4. Repeating steps 2 and 3 for all of the other rows. In the third part (proposed linear solution) we have the following steps: 1. 2. 3. 4.
Assuming i = 1 (first row). Using Eq. (16) to find P(i, j)2×2 matrix. Using Eq. (17) to calculate (i,n), (i,n − 1), . . ., (i,1). Calculating Jy (i + 1, j, tn), ey (i + 1, j, tn) for 1 ≤ j ≤ n by Eqs. (1)–(4) and results of step 4. 5. Assuming i = i + 1 up to i = m (going to the next row). 6. Doing steps 3–6 and calculating all (i,j). 7. Comparing accuracy and time consuming between two methods. 4.1. Numerical example 1: 3 × 3 block
Linear solution results (inverse problem) 6.4023258 × 10−11 0.04 s 2.8884456 × 10−8 0.156 s 2.0819278 × 10−7 0.438 s 0.016674409 0.437 s
4.2. Numerical example 2: 4 × 4 block In this example a 4 × 4 block with similar conductivity distribution ( = 20) has been considered, however, in two of the central blocks conductivity value is = 200. Table 1 describes both solutions have acceptable RMSE however linear solution is about 33 times faster than non-linear method. 4.3. Numerical example 3: 5 × 5 block with two different block conductivity In this example a 5 × 5 block with similar conductivity distribution ( = 20) has been considered, however, in two of the central blocks conductivity value is = 200. Table 1 shows that both methods have acceptable RMSE however linear solution is about 33 times faster than non-linear method. 4.4. Numerical example 4: 5 × 5 block with non-uniform conductivity distribution Table 1 shows a complicated problem (a 5 × 5 block with nonuniform conductivity distribution) and results of inverse linear and non-linear solutions. As the numerical example gets bigger and complicated, results of inverse solution get worse and RMSE value gets smaller. Fig. 2 shows results of RMSE for the previous examples. For a better view we have used RMSE logarithm to base 10. Although linear and non-linear solutions both get worse as the numerical examples get bigger and complicated, linear solution has considerable better RMSE values than non-linear solution in all examples. EIT inverse solution is a time consuming problem and fast algorithms are more favourite. As it has been shown in Fig. 3 proposed linear solution for EIT is faster than non-linear solution. Also linear solution has a very fast and almost constant algorithm time for different examples, however non-linear solution gets very slow on complicated problems. In most algorithms for EIT inverse problem presented in literature, there is no numerical or quantitative comparison between real and calculated conductivities. However in some papers RMSE has been adopted as a comparison factor [19,22–26]. For example
In this example a 3 × 3 block with similar conductivity distribution ( = 20) has been considered, however, the central block conductivity is = 200. Table 1 shows differences of non-linear and linear inverse solution results. Root mean square error (RMSE) has been used to compare real and calculated values of conductivities.
rmse =
m n i=1
j=1
(r (i, j) − c (i, j)) m×n
2
(18)
where r is the real value of the conductivity and c is the calculated one. Table 1 describes both solutions are accurate and have acceptable RMSE however linear solution is very fast. For this example proposed linear solution is about 44 times faster than non-linear method.
Fig. 2. Difference of RMSE logarithm to base 10 for the previous examples.
A. Abbasi, B.V. Vahdat / Journal of Computational Science 1 (2010) 190–196
195
this method is theoretical, it may have better results than other common ways in literature. Development of this method from 2D to 3D can be suggested as new field for future investigations.
References
Fig. 3. Algorithm time in linear and non-linear solutions for the previous examples. Table 2 RMSE comparison between methods. Algorithm for EIT
RMSE
Shape estimation and state estimation for dynamic EIT [26] Shape estimation by EKF for dynamic EIT [24] Employing IMM, EKF and CCEKF for dynamic EIT [23] Neural networks and front point approach [25] Interpolation of front points method [22] Non-linear block method [19] Proposed linear block method (this paper)
10−1 10−2 10−2 10−2 10−5 10−5 10−8
to 10 to 10−1 to 10−1 to 10−1 to 10−3 to 10−4 to 10−7
in a work conducted by Kim for 2D EIT inverse problem, applying interpolation of front points method for approximation of regions of inner object results in RMS on order of 10−5 to 10−3 for reconstructed image [22]. Also, using neural networks and front point approach for estimation of 2D EIT image results 10−2 to 10−1 in RMSE [25]. Tossavainen et al. applied shape estimation and state estimation formulation for a dynamic EIT problem and reported RMSE from 10−1 to 10 in conductivity [26]. In other work conducted by Ijaz, a dynamic reconstruction algorithm has been presented to monitor the concentration distribution inside the fluid vessel based on EIT by employing interacting multiple model (IMM), extended Kalman filtering (EKF) and covariance compensation extended Kalman filtering (CCEKF). For this study RMSE of conductivity has been reported between 10−2 to 10−1 [23]. Using shape estimation of regions of known resistivities based on extended Kalman filtering for dynamic EIT results in RMSE from 10−2 to 10−1 [24]. Table 2 represents these results. Therefore RMSE from 10−8 to 10−7 for this method is in a proper order and indicates that real and calculated conductivities are almost the same. Although the computer simulations proved that this linear reconstruction algorithm employing the block method results in an accurate identification, the block method has some limitations too. The assumption that: “object has a 2D rectangular shape and is made up of identical fixed size blocks”, causes some difficulties in the practical examples. The block method reconstruction algorithm can be applied for special cases with rectangular shapes. Also it can be used for the other subjects by increasing the number of block (m × n → ∞) which makes more complicated equations. In block method, measurement noise can be propagated easily and influence the results [27]. From the application point of view, error propagation on EIT block method is another limitation of this method. 5. Conclusion The proposed method is an accurate solution for 2D electrical impedance tomography. Low error and high resolution of this method is clear. Fast solution and almost the same algorithm time for complicated problems is the other advantage of linear solution. This method can lead a great step in EIT problem solution. Although
[1] M. Cheney, D. Isaacson, J.C. Newell, Electrical impedance tomography, SIAM Review 41 (1) (1999) 85–101. [2] I. Frerichs, Electrical impedance tomography (EIT) in applications related to lung and ventilation: a review of experimental and clinical activities, Physiological Measurement 21 (2000) 1–21. [3] K.S. Osterman, T.E. Kerner, D.B. Williams, Multifrequency electrical impedance imaging: preliminary in vivo experience in breast, Physiological Measurement 21 (2000) 99–110. [4] N. Polydorides, W.R. Lionheard, H. McCann, Krylov subspace iterative techniques: on the detection of brain activity with electrical impedance tomography, IEEE Transactions on Medical Imaging 21 (6) (2002) 596–603. [5] X. Wu, Y. Sun, J. Cen, Circuit for measuring breath waveform with impedance method and method and device for resisting interference of electrical fast transient, Shenzhen Mindray Bio-Medical Electronics Co., Ltd., United States patent and trademark office pre-grant publication, 2007. [6] F. Santosa, M. Vogelius, A backprojection algorithm for electrical impedance imaging, SIAM Journal on Applied Mathematics 50 (1) (1990) 216–243. [7] M. Cheney, D. Isaacson, An overview of inversion algorithms for impedance imaging, Inverse Scattering and Applications Contemporary Mathematics 122 (1991) 29–39. [8] M. Cheney, D. Isaacson, J.C. Newell, J. Goble, S. Simske, Noser: an algorithm for solving the inverse conductivity problem, International Journal of Imaging Systems and Technology 2 (1990) 66–75. [9] A. Allers, F. Santosa, Stability and resolution analysis of a linearized problem in electrical impedance tomography, Inverse Problems 7 (4) (1991) 515–533. [10] A.J. Suroviec, S.S. Stuchly, J.R. Barr, A. Swarup, Dielectric properties of breast carcinoma and the surrounding tissues, IEEE Transactions on Bimedical Engineering 35 (4) (1988) 257–263. [11] T.J. Yorkey, J.G. Webster, W.J. Tompkins, Comparing reconstruction algorithms for electrical impedance tomography, IEEE Transactions on Bimedical Engineering 34 (11) (1987) 843–852. [12] J.S. Kallman, J.G. Berryman, Weighted least-squares criteria for electrical impedance tomography, IEEE Transactions on Medical Imaging 11 (2) (1992) 284–292. [13] M. Soleimani, Computational aspects of low frequency electrical and electromagnetic tomography: a review study, International Journal of Numerical Analysis and Modeling 5 (3) (2008) 407–440. [14] R.V. Kohn, A. McKenney, Numerical implementation of a variational method for electrical impedance imaging, Inverse Problems 9 (1990) 389–414. [15] J. Kaipio, V. Kolehmainen, E. Somersalo, M. Vauhkonen, Statistical inversion and Monte Carlo sampling methods in electrical impedance tomography, Inverse Problems 16 (5) (2000) 1487–1522. [16] M. Cheney, G. Krislensson, Three dimensional inverse scattering: layerstripping formulae and ill-posedness results, Inverse Problems 4 (1988) 625–642. [17] J. Sylvester, A convergent layer stripping algorithm for radially symmetric impedance tomography problem, Communications in Partial Differential Equations 17 (1992) 1955–1994. [18] J.L. Mueller, S. Siltanen, D. Isaacson, A direct reconstruction algorithm for electrical impedance tomography, IEEE Transactions on Medical Imaging 21 (6) (2002) 555–559. [19] A. Abbasi, B. Vosoughi Vahdat, G. Ebrahimi Fakhim, An inverse solution for 2D electrical impedance tomography based on electrical properties of material blocks, Applied Sciences 9 (10) (2009) 1962–1967. [20] W.R. Lionheart, EIT reconstruction algorithms: pitfalls, challenges and recent developments, Physiological Measurement 25 (2004) 125–142. [21] S Babaeizadeh, D.H. Brooks, D. Isaacson, 3-D electrical impedance tomography for piecewise constant domains with known internal boundaries, IEEE Transactions on Bimedical Engineering 54 (1) (2007) 2–10. [22] M.C. Kim, S. Kim, K.Y. Kim, K.H. Seo, H.J. Jeon, J.H. Kim, B.Y. Choi, Estimation of phase boundary by front points method in electrical impedance tomography, Inverse Problems in Science and Engineering 14 (5) (2006) 455–466. [23] U.Z. Ijaz, J.H. Kim, A.K. Khambampati, M.C. Kim, S. Kim, K.Y. Kim, Concentration distribution estimation of fluid through electrical impedance tomography based on interacting multiple model scheme, Flow Measurement and Instrumentation 18 (1) (2006) 47–56. [24] U.Z. Ijaz, K.Y. Kim, Kinematic models for non-stationary elliptic region boundary in electrical impedance tomography, in: Proceedings of the 2006 International Conference on Hybrid Information Technology (ICHIT 2006), Jeju, Korea, 2006, pp. 9–11. [25] J.H. Kim, B.C. Kang, S.H. Lee, B.I. Choi, M.C. Kim, B.S. Kim, U.Z. Ijaz, K.Y. Kim, S. Kim, Phase boundary estimation in electrical resistance tomography with weighted multi-layered neural networks and front point approach, Measurement Science and Technology 17 (10) (2006) 2731–2739. [26] O.P. Tossavainen, M. Vauhkonen, V. Kolehmainen, K.Y. Kim, Tracking of moving interfaces in sedimentation processes using electrical impedance tomography, Chemical Engineering Science 61 (23) (2006) 7717–7729.
196
A. Abbasi, B.V. Vahdat / Journal of Computational Science 1 (2010) 190–196
[27] A. Abbasi, F. Pashakhanlou, B.V. Vahdat, Error propagation in non-Iterative EIT block method, in: IEEE International Symposium Signal Processing and Information Technology, 2007, pp. 678–681. Ata Abbasi is an assistant professor of Sahand University of Technology, Tabriz, Iran. He received the B.Sc. degree in Biomedical Engineering from Sahand University of Technology and M.Sc. degree from Sharif University of Technology, Iran, in 2003 and 2005 respectively. From 2003 to 2005, he was a research assistant in Sharif University of Technology. He graduated Ph.D. in Biomedical Engineering Sharif University of Technology, 2006–2010.
Bijan Vosughi Vahdat is an assistant professor of Sharif University of Technology, Tehran Iran. He graduated B.Sc. in Electronics, Mashad University, 1983–1987 and M.Sc. in Biomedical Engineering, Sharif University of Technology, 1987–1991 and Ph.D. in Biomedical Engineering, University of New South Wales, 1993–1996. His areas of research interest include Medical Imaging Systems, Advanced Programming and Modelling.