Direct calculation of the electrode movement Jacobian for 3D EIT.pdf

Loading ...
Global Do...
News & Politics
Direct calculation of the electrode movement Jacobian for 3D EIT
Camille Gómez-Laberge1 and Andy Adler1
1
Systems and Computer Engineering, Carleton University, Ottawa, Canada
Abstract—Electrical Impedance Tomography (EIT) of me-
dia with deformable boundaries is very sensitive to electrode
movement. This is especially important for EIT images of the
thorax, which become distorted with breathing and posture
change. Previously, we proposed a reconstruction method for
imaging conductivity change and electrode movement based
on an indirect perturbation Jacobian calculation, involving the
re-computation of the forward solution. Although suitable for
2D and small 3D imaging, the reconstruction accuracy of this
method gradually decreases, while the computation time grows
rapidly for large 3D problems. We propose an efficient, novel
method of calculating the Jacobian matrix directly from the
Finite Element Method (FEM) system equations, without the
re-calculation of the forward solution.
The implemented algorithm is based on asymmetric rank-
one perturbations of the admittance matrix. We show that the
measurement sensitivity calculations, due to the displacement
of each electrode, reduce to operations on small submatrices
of the FEM system matrix. The proposed algorithm is applied
to simulated data using 3D FEM reconstruction models of
various element densities. The computation speed and the
reconstruction fidelity are compared between the proposed and
previous methods.
Keywords—electrical
impedance tomography, electrode
movement, inverse problems.
I. INTRODUCTION
In clinical EIT, boundary deformation occurs mostly in
thoracic measurement due to breathing and posture change
during measurement [1]. The former causes the expansion
and the contraction of the rib cage; the latter causes a
displacement of the skin and the electrodes with respect to
the lungs. The difficulty arises since EIT measurements are
imposed onto a geometric reconstruction model that approx-
imates the shape of the body being observed. The movement
of electrodes reduces the accuracy of the reconstruction
model and introduces spurious difference measurements
producing broad artefacts in the images [2]. Related are the
errors from inadequate geometric models and inconsistent
electrode placement. These problems are well known in the
literature.
Blott et al. studied the effects of electrode movement
in 2D EIT and compensated for reconstruction error by
re-adjusting measurement data using a spatial smoothness
constraint [3]. Lionheart showed that in 3D EIT, the mea-
surement data is sufficient to determine both the conductivity
distribution and the boundary shape [4]. In order to attempt
accounting for the errors due to electrode movement, EIT
imaging algorithms have been developed to reconstruct
both the conductivity change and the electrode movement
[5], [6]. These approaches are promising and show dra-
matic decreases in image artefacts associated with electrode
movement. One limitation of these techniques it that the
Jacobian (or sensitivity) matrix needs to be calculated for
both conductivity change and electrode movement. While
good algorithms for the conductivity change Jacobian are
available [7], the electrode movement portion has been
calculated using numerical perturbation techniques. This
computation is slow and can be numerically inaccurate as
the FEM models get larger. This article proposes an efficient
method of directly calculating the Jacobian matrix, i.e.,
without using numerical perturbation and re-calculation of
the forward problem.
II. METHODS
A. The system model
The conductivity distribution is represented within the 3D
FEM volume. The FEM used in this paper is constructed by
extruding a 2D circular FEM into a 3D cylinder with tetra-
hedral elements. The conductivity distribution is constant on
each element, and is defined as the vector σ ∈ RNk , where
Nk is the number of elements in the FEM.
Sixteen electrodes are placed around the volume’s cir-
cumference in correspondence with the experiment’s elec-
trode configuration. Point electrodes are considered, which
are directly connected to the nodes on the 3D model
boundary. Figure 1 illustrates a 7,680-element FEM model
with the electrode configuration used in this paper. The
electrodes are numbered in a “zig-zag” pattern [8], where
electrode 1 (light), electrode 2 (medium), and the remaining
electrodes (dark) are shown by green stars in figure 1. The
potential difference between electrode i and the reference is
defined as φi. The two blue elements shown in figure 1 are
contrast elements used to compare image reconstructions.
While two electrodes excite the medium, the remaining
−1
−0.5
0
0.5
1
−1
−0.5
0
0.5
1
−1
−0.8
−0.6
−0.4
−0.2
0
0.2
0.4
0.6
0.8
1
x
y
zFig. 1 A forward model FEM with 7,680 tetrahedral elements.
electrodes measure the voltage between adjacent pairs. The
kth measurement vk = φi − φj , is the voltage between two
electrodes i and j. This excitation process is repeated for
each electrode pair; therefore, for Ne electrodes, we obtain
Nv = Ne(Ne − 3) voltage measurements v ∈ RNv .
In the forward problem, we also consider each electrode’s
position vector r ∈ R3, since we expect displacements to
occur during the measurements. Hence, σ and r are used to
construct the admittance matrix Y ∈ RNp×Np , where Np
is the number of nodes in the FEM. This matrix associates
each FEM element i = 1, . . . , Nk to its conductivity σi and
its constituent nodes. The excitation pattern used to inject
current into the medium is represented by the matrix Q ∈
RNp×Ne . Each column of Q indicates which electrodes are
injecting current into the medium. The matrix product
V = Y(σ, r)−1Q
(1)
yields the nodal potential difference matrix V ∈ RNp×Ne ,
such that for excitation j = 1, . . . , Ne, we haveVij = φi for
nodes i = 1, . . . , Np. Hence, v can be extracted fromV by a
series of element-wise subtractions Vij−Vkl corresponding
to the voltage measurement pattern. Although an extraction
operator T is required to obtain v = T [V], we omit this to
simplify the formulation.
B. Electrode movement Jacobian
The electrode movement Jacobian Jm is a submatrix
which augments the standard, conductivity Jacobian Jc.
Thus, the system matrix is formed as
J = [Jc, Jm].
The conductivity Jacobian is calculated exactly as in [11]
and has dimension RNv×Nk . This algorithm is available in
the EIDORS software base under the GNU Public Licence
[9]. By definition, Jm is the sensitivity of the measurement
data V to the movement of an electrode with position vector
r = (rx, ry, rz). That is, for electrodes 1, . . . , Ne, we have
Jm =
[
∂V
∂r1
· · · ∂V
∂rNe
]
∈ RNv×3Ne .
(2)
Hence, we re-write the admittance matrix from equation (1)
as a separable product
Y(σ, r) = C>D(σ)S(r)C,
(3)
where C ∈ RNk×Np
is the connectivity matrix, which
associates each element to its vertices. D(σ) ∈ R4Nk×4Nk
represents the conductivity distribution and is formed from
D(σ) = diag(σ)⊗ I4
(4)
where ⊗ is the Kronecker product, and I4 is the 4×4 identity
matrix. We refer to the matrix S(r) ∈ R4Nk×4Nk as the FEM
system matrix, since it encapsulates the geometric properties
of the FEM model. It is block-diagonal and defined as
S =

S1
0
· · ·
0
0 S2
...
...
. . .
0
0
· · ·
0 SNk
 .
(5)
Each block is given by
Si =
2
Nd!
1
|detAi|
B>
i Bi,
where the model dimension Nd = 3 for 3D. The matrix Ai
contains the geometry data of element i
Ai =

1 r1x
r1y
r1z
1 r2x
r2y
r2z
1 r3x
r3y
r3z
1 r4x
r4y
r4z

−1
,
(6)
where (rj,i
x , r
j,i
y , r
j,i
z ) is the (x, y, z) Cartesian co-ordinates
of vertex j of element i. Each element i, is defined by its
four vertices 1 ≤ j ≤ 4. The matrix Bi is the submatrix of
Ai with the top row deleted, denoted by Bi = [Ai]\row1.
In general, for an electrode with position vector r, we
have, from equations (1) and (2),
∂V
∂r
=
∂
∂r
(
Y−1Q
)
(7)
= Y−1
∂Y
∂r
Y−>Q,
(8)
where, from equation (3),
∂Y
∂r
= C>D
∂S
∂r
C.
(9)
In order to evaluate equation (9), we need only consider the
blocks Si where r is a vertex node of the perturbed element
i; the other blocks vanish because they are constant. For
each non-zero block, we require the derivative
∂Si
∂r
=
∂
∂r
[
2
Nd!
1
|detAi|
B>
i Bi
]
.
The product rule of this equation yields
∂Si
∂r
=
2
Nd!
[
∂
∂r
(
1
|detAi|
)
B>
i Bi +
1
|detAi|
(
∂B>
i
∂r
Bi +B>
i
∂Bi
∂r
)]
.
(10)
The three partial derivative terms can be calculated using
an asymmetric rank one perturbation technique similar that
shown in lemma 2 of [12]. The proof can be found in [13].
Asymmetric rank one perturbations: Let α ∈ R, a,b ∈
Rd, and X ∈ Rd×d be an invertible matrix. Then if α
6=
−(b>Xa)−1, the rank one perturbed matrix X+ αab> is
invertible and
(X+ αab>)−1 = X−1 − αX
−1ab>X−1
1 + αb>X−1a
.
(11)
In addition, the perturbation determinant is
det(X+ αab>) = (1 + αb>X−1a) detX.
(12)

Therefore, by equation (11), the second term of equation
(10) can be calculated for a perturbation of r along any
direction in the limit α → 0+. Recall equation (6), and let
Ai = P−1
i
. Then, for any element i,
dBi
dα
=
dAi
dα
∣∣∣∣
\row1
=
d
dα
[(
Pi + αab>
)−1]∣∣∣∣ α=0
\row1
=
[
−P−1
i ab
>P−1
i
]
\row1
=
[
−Aiab>Ai
]
\row1
(13)
and similarly
dB>
i
dα
=
(
dBi
dα
)>
=
[
−A>
i ba
>A>
i
]
\row1.
Table 1 FEM Forward & Inverse Model Pairs
Model
Elements
Nodes
pair
Forward
Inverse
Forward
Inverse
A
7,680
1,536
1,595
369
B
15,360
3,072
3,045
697
The first term of equation (10) is calculated using equation
(12):
d
dα
1
|detAi|
=
=
−sgn(detAi)
detA2i
d
dα
detAi
=
−sgn(detAi)
detA2i
d
dα
det
[
(P−1
i + αab
>)−1
]∣∣∣∣
α=0
=
1
|detAi|
b>Aia
(14)
Substituting equations (13) and (14) into (10) yields the
sensitivity term for each element i with vertex represented
by r
∂Si
∂r
=
2
Nd!
[
1
|detAi|
(
b>AiaB>
i Bi +
∂B>
i
∂r
Bi +B>
i
∂Bi
∂r
)]
.
(15)
Equation (15) models the deformation of an element i
incurred by the perturbation of the adjoining electrode with
position vector r. For each deformed element, the pertur-
bations along x, y, and z are calculated by choosing the
unit vectors a and b in equations (13) and (14) that modify
the corresponding node’s coordinates. These coordinate data
are stored as shown in equation (6). Thus, the FEM system
matrix S is constructed as in equation (5) where only the
element blocks Si touching the displaced electrode are non-
zero. Finally, the matrix ∂Y/∂r can be directly computed,
and the movement Jacobian is obtained by re-arranging the
elements of ∂V/∂r such that Jm ∈ RNv×3Ne .
III. RESULTS
The computations of the conductivity distribution depend
on quantities derived from the FEM geometry. Hence, we
consider two model densities differing by a factor of 2: 1,536
and 3,072 elements, respectively. When simulating EIT data,
we also rely on the FEM to solve the forward problem, i.e.,
calculating v given S and σ. Table 1 summarizes the two
FEM model pairs used, labeled A and B.
Image reconstructions were computed using the proposed
method and were compared with the previously used indirect
Table 2 Movement Jacobian Results
Model
Computation Time (ms)
Relative Norm
pair
Direct
Indirect
||Jindir − Jdir|| / ||Jdir||
A
840
22,420
1.48 × 10−6
B
1,460
41,430
1.57 × 10−6
Direct J
Indirect J
F ig. 2 Reconstructions of the FEM model in figure 1 using the direct
and indirect methods.
perturbation method in [6]. Table 2 shows the average
computation time required in milliseconds to compute the
Jacobian matrix for each FEM model and the relative norm-
error between both methods. A significantly longer computa-
tion time is required for the indirect method. The electrodes
were slightly displaced (less than 1% of the model diameter)
in simulations using the FEM model shown in figure 1.
Figure 2 shows the reconstructions obtained with the direct
and indirect methods. Three slices at z = −0.4, 0,+0.4 are
shown where the electrode positions are indicated by the
small green circles. The position of both contrast elements
in blue has been recovered; however, each method shows a
slightly different artefact.
IV. DISCUSSION
In this study, we formulate a more efficient and accurate
computation of the electrode movement Jacobian by elimi-
nating the need of re-calculation of the forward problem and
the numerical error associated with numerical perturbations.
This was accomplished by implementing an asymmetric
rank one perturbation technique similar to that described
in [12] directly on the FEM system matrix. The model
developed here applies to point electrodes; however the
underlying approach can be easily modified for use with
more complex electrode models. Results show significantly
faster computations that can complete a 15,000-element 3D
FEM model with 16 electrodes in under 20 seconds using
a standard desktop PC. Imaging results are similar to the
previous, indirect method described in [6] and show no sig-
nificant change in reconstruction fidelity. Such a calculation
method may be useful in reconstruction algorithm design
for large FEM inverse problems.
REFERENCES
1.
Lozano A, Rosell J, Pallás-Areny R (1995) Errors in prolonged
electrical impedance measurements due to electrode repositioning and
postural changes. Physiol Meas 16:121–130.
2.
Adler A, Guardo R, Berthiaume Y (1996) Impedance imaging of lung
ventilation: do we need to account for chest expansion? IEEE Trans
Biomed Eng 43:414–420.
3.
Blott BH, Daniell GJ, Meeson S (1998) Electrical impedance tomog-
raphy with compensation for electrode positioning variations. Phys
Med Biol 43:1731–1739.
4.
Lionheart WRB (1998) Boundary shape and electrical impedance
tomography. Inv Prob 14:139–147.
5.
Kolehmainen V, Lassas M, Ola P (2005) The inverse conductivity
problem with an imperfectly known boundary. SIAM J Appl Math
66:365–383.
6.
Soleimani M, Gómez-Laberge C, Adler A (2006) Imaging of con-
ductivity changes and electrode movement in EIT. Physiol Meas
27:S103–S113.
7.
Polydorides N (2002) Image reconstruction algorithms for soft-field
tomography. Ph.D. Thesis. University of Manchester Institute of
Science and Technology pp.262.
8.
Graham BM, Adler A (2006) Electrode placement strategies for 3D
EIT. 7th Conf Biomed EIT.
9.
Adler A, Lionheart WRB (2006) Uses and abuses of EIDORS: and
extensible software base for EIT. Physiol Meas 27:S25–S42.
10.
Cheng KS, Isaacson D, Newell JC et al. (1989) Electrode models
for electric current computed tomography. IEEE Trans Biomed Eng
36:918–924.
11.
Adler A, Guardo R (1996) Electrical impedance tomography: reg-
ularized imaging and contrast detection. IEEE Trans Med Imag
15:170–179.
12.
Olsen PA, Gopinath AR (2004) Modeling inverse covariance matri-
ces by basis expansion. IEEE Trans Speech Audio Processing 12:37–
46.
13.
Golub GH, Loan CFV (1996) Matrix computations. Johns Hopkins
Univ. Press, Baltimore.
Address of the corresponding author:
Author: Camille Gómez-Laberge
Institute: Carleton University
Street: 1125 Colonel By Drive
City: Ottawa
Country: Canada
Email: cgomez@sce.carleton.ca