NEW HARDWARE AND SOFTWARE DESIGN FOR ELECTRICAL
IMPEDANCE TOMOGRAPHY
NEW HARDWARE AND SOFTWARE DESIGN FOR
ELECTRICAL IMPEDANCE TOMOGRAPHY
By
Mehran Goharian, M.Sc., B.Sc. (Hons.)
A Thesis
Submitted to the School of Graduate Studies
in Partial Fulfilment of the Requirements
for the Degree
Doctor of Philosophy
McMaster University
© Copyright by Goharian, September 2007
DOCTOR OF PHILOSPHY (2007)
(Medical Physics)
McMaster University
Hamilton, Ontario
TITLE:
New Hardware and Software Design for Electrical
Impedance Tomography
AUTHOR:
Mehran Goharian
SUPERVISOR:
Dr. Gerald R.Moran
NUMBER OF PAGES: xii, 132
ii
Abstract
Electrical impedance tomography (EIT) is an imaging technique that reconstructs the internal electrical
properties of an object from boundary voltage measurements. In this technique a series of electrodes is
attached to the surface of an object and alternating current is passed via these electrodes and the resulting
voltages are measured. Reconstruction of internal conductivity images requires the solution of an ill-
conditioned nonlinear inverse problem from the noisy boundary voltage measurements. Such unreliable
boundary measurements make the solutions unstable. To obtain stable and meaningful solutions
regularization is used. This thesis deals with the EIT problem from the perspective of both image
reconstruction and hardware design. This thesis consists of two main parts. The first part covers the
development of 3D image reconstruction algorithms for single and multi-frequency EIT. The second part
relates to the design of novel multi-frequency hardware and performance testing of the hardware using the
designed phantom.
Three different approaches for image reconstruction of EIT are presented:
1) The dogleg algorithm is introduced as an alternative method to Levenberg-Marquardt for solving the EIT
inverse problem. It was found that the dogleg technique requires less computation time to converge to the
same result as the Levenberg-Marquardt.
2) We propose a novel approach to build a subspace for regularization using a spectral and spatial multi-
frequency analysis approach. The approach is based on the construction of a subspace for the expected
conductivity distributions using principal component analysis (peA). The advantage of this technique is
that priori information for regularization matrix is determined from the statistical nature of the multi-
frequency data.
3) We present a quadratic constrained least square approach to the EIT problem. The proposed approach is
based on the trust region subproblem (TRS), which uses L-curve maximum curvature criteria to fmd a
regularization parameter. Our results show that the TRS algorithm has the advantage that it does not require
any knowledge of the norm of the noise for its process.
4) The second part of thesis discuses the designing, implementation, and testing a novel 48-channel multi-
frequency EIT system. The system specifications proved to be comparable with the existing EIT systems
iii
with capability of 3-D measurement over selectable frequencies. The proposed algorithms are [mally tested
under experimental situation using designed EIT hardware. The conductivity and permittivity images for
different targets were reconstructed using four different approaches: dog-leg, principal component analysis
(PCA), Gauss-Newton, and difference imaging.
In the case of the multi-frequency analysis, the PCA-based
approach provided a substantial improvement over the Gauss-Newton technique in terms of systematic
error reduction. Our EIT system recovered a conductivity value of 0.08 Sm-l for the 0.07 Sm-l piece of
cucumber (14% error).
IV
Acknowledgements
I express my deep gratitude to Dr. Gerald Moran for his continued guidance and help through the whole
years of my graduate studies in McMaster University. Without his support this thesis would have never
been completed. I would like to thank Dr. John Macgregor, a member of my Ph.D.committee, for both his
encouragement and fmancial support for the EIT project. Thanks to my other supervisory committee
members, Drs. W. V. Prestwich and Frank. Prato, for their time commitment and suggestions as well.
I would also to thank Dr.Mark-lohn Bruwer for the collaborative work through parts of this project. Thanks
to Kenrick Chin for his efforts fabricating and designing the various parts of hardware for this project. I
have learnt much from him. He is a real problem solver!
I would also like to thank my friends Jovica, Marijia, Mariangela, and Rebecca who also shared an office
with me. I am grateful to my friends Aslam and Aravinthan for their company and generous help on
different occasions.
I wish to express my love to my mother and sisters for their encouragement during my years. I extend my
deepest gratitude to my father and mother in-laws, for their continued support and friendship throughout
the years. You're the best!
I would like to thank the Department of Medical Physics and Applied Radiation Sciences for the funding
that has made my graduate studies possible. Thanks to the administrative staff, Wendy, Fiona, and Vera
and to Dr. Fiona McNeill, the chair of the department, for their assistance and support over my graduate
studies at McMaster University.
I am deeply thankful to my wonderful wife and best friend, Parisa, for her moral support, encouragement
and I am grateful for her continued patience that has made the completion of this work possible. I dedicate
this work to her.
v
Table of Contents
Abstract
iii
Acknowledgment
V
Th~ifro~~
ri
List of Figures
viii
Abbreviations
xii
Chapter I: Introduction
2
3
Background
Forward Model and Data Collection in EIT
2.1
2.2
2.3
2.4
2.5
Physics of the Problem
Electrode Models
Numerical Solution to Forward Model
Algorithms to Solve Forward Problem
Current Patterns and Measurements Procedure
Image Reconstruction Algorithms
3.1
3.2
3.3
Linearization and Jacobian
Regularization
Iterative Method
3
4
5
7
11
12
15
15
18
19
4
References
23
Chapter II: Paper I: Dogleg trust-region application in electrical impedance tomography
28
2.0.1
Introduction to Paper I
28
2.0.2 Contents of Paper I
29
Chapter III: Paper II: A novel approach for EIT regularization via spatial and spectral principal
component analysis
47
3.0.1
Introduction to Paper II
47
3.0.2
Contents of Paper II
48
Chapter IV: Paper III: Regularization of EIT problem using trust region subproblem
64
4.0.1
Introduction to Paper III
64
4.0.2 Contents of Paper III
65
vi
Chapter V: Paper IV: A DSP based multi-frequency electrical impedance tomography system 75
5.0.1 Introduction to Paper IV
75
5.0.2 Contents of Paper IV
76
Chapter VI: Paper V: 3D multi-frequency electrical impedance tomography imaging in phantoms
6.0.1 Introduction to Paper V
100
6.0.2 Contents of Paper V
Chapter VII: Summary and Future Work
8.1
8.2
8.3
8.4
Development of 3D Algorithms for Image Reconstruction of BIT
Design and Building Multi-Frequency BIT System
Future Research for Multi-Frequency BIT
References
Appendix A: Paper VI: Modifying the MR!, elastic stiffness and electrical properties of polyvinyl
101
119
119
121
122
125
alcohol cryogel using irradiation
126
A.O.I Introduction to Paper VI
126
A.0.2 Contents of Paper VI
127
vii
List of Figures
Chapter I: Introduction
Figure 1: 3D meshes for a cylinder with circular electrodes attached on surface. The mesh was ........ .
generated with NETGEN mesh generator. ............................ , '" ............................................. 8
Figure 2: Adjacent current driver. The current is passed through two neighbour electrodes and the .. ..
voltages differences are collected from successive pairs of adjacent electrodes ................................. 13
Chapter II: Paper I
Figure 1: The complete Levenberg-Marquardt algorithm ............................................. '" ......... .32
Figure 2: The complete Powell's dogleg algorithm .. , ...................................................... '" ... .34
Figure 3: Simulated geometry with 48 electrodes .................................................................... 36
Figure 4: Tumor geometry for two abnormalities case ofEIT simulation ..................................... 36
Figure 5: Original conductivity distribution for the case of one abnormality at (-0.5,-0.5,0.5) at different..
heights, z, of cylinder for 30% changes ....................................................................................... .3 7
Figure 6: Reconstructed images with dogleg method with 1 % noise for the case of one abnormality at (-
0.5,-0.5,0.5) at different heights, z ,of cylinder .................................................. , ................... 38
Figure 7: Reconstructed images with LM method with 1 % noise for the case of one abnormality at (-0.5,-
0.5,0.5) at different heights, z ,of cylinder ............................................................................ 38
Figure 8: Original conductivity distribution for the case of one abnormality at (0,0,0.5) at different ...... ..
Heights,z, of cylinder for 30% change ........ , ..................................................................... .39
Figure 9: Reconstructed images with dogleg method with 0.5% noise for the case of one abnormality at .. .
(0,0,0.5) at different heights, z, of cylinder ........................................................................... 39
Figure 10: Reconstructed images with dogleg method with 0.5% noise for the case of one abnormality at (-
0.5,-0.5,0.5) at different heights, z, of cylinder with 10% conductivity change .............. , ........... , ..... .40
Figure
11: Original conductivity distribution for the case of two abnormality both with 30% conductivity
change at (-0.5,-0.5,0.5) and (1,1,0.5) at different heights, z, of cylinder ....................................... .40
Figure 12: Reconstructed conductivity images with dogleg for 1% noise for the case of two abnormality
both with 30% conductivity change at (-0.5,-0.5,0.5) and (1,1,0.5) at different heights, z, of cylinder ..... .41
viii
Figure 13: Original conductivity distribution for the case of two abnormality at (-0.5,-0.5,0.5) and (1,1,0.5)
at different heights, z, of cylinder one with 30% higher and the other one with 50% lower than background
..................... '" '" .......... , ........ , ...................... '" .. , '" ....................... , '" ..... , .......... , .... 41
Figure 14: Reconstructed conductivity images with dogleg for the case of two abnormality at (-0.5,-0.5,0.5)
and (1,1,0.5) at different heights, z, of cylinder at 1 % noise one with 30% higher and the other one with
50% lower than background ....................................................... " .............. , .................... 42
Figure 14: Original conductivity distribution for the case of two abnormality at (-0.5,-0.5,0.5) and
(0.5,0.5,1.5) at different heights, z, of cylinder one with 30% higher and the other one with 50% lower than
background .... , .................. , .......... " ..... , ........ , ............................................................. 42
Figure 16: Reconstructed conductivity images with dogleg for the case of two abnormality at (-0.5,-0.5,0.5)
and (0.5,0.5,1.5) at different heights, z, of cylinder with 1 % noise one with 30% higher and the other one
with 50% lower than background ..................................................... , .............. , ..................... 43
Figure 17: Convergence rate change for objective function verses iteration for the case of two objects at (-
0.5,-0.5,0.5) and (1,1,0.5) with conductivity changes of+30% and -50% at 1% noise leve1. ................................ 45
Chapter III: Paper II
Figure 1: Flowchart for SS-MIA algorithm ................................................... '" ..................... 54
Figure 2: Simulated EIT geometry with electrodes .................................................................. 55
Figure 3: Tumor geometry for EIT simulation .......................................................................... 56
Figure 4: Breast tissue conductivity variation verses frequency (Surowiec
et a/1988) ..... .................. 56
Figure 5: Loading weights for first principal component of MIA PCA model.. .................................... 57
Figure 6: Loading weights for fIrst principal component of SS-MIA PCA model. Model shows correlation
structure between spectral signature and its relationship with its spatial neighbors. The spectral signature
comprises just one variable. Its associated loading weight is 0.78 (the first element - off the scale) ......... 58
Figure 7: Original conductivity distributions for geometry illustrated in Figureure 2 shown as a series of2-
D slices. (at lowest frequency 10 KHz) ................................................................................. 58
Figure 8: Reconstructed conductivity distributions for case of 1 % noise level with both MIA and SS-MIA
shown as a series of 2-D slices. Images reconstructed after 5 iterations. (at lowest frequency 10 KHz) ..... 59
Figure 9: Reconstructed conductivity distributions for case of 3% noise level with both MIA and SS-MIA
as shown series of 2-D slices. Images reconstructed after 5 iterations (at lowest frequency 10 KHz) ........ 59
Figure 10: Reconstructed conductivity distributions for case of3% noise level with only SS-MIA shown as
a series of 2-D slices. In this case the regularization matrix is updated in each iteration based on SS-MIA
PCA model. Images reconstructed after 5 iterations (at lowest frequency 10 KHz) ............................... 60
Figure 11: Reconstructed conductivity distributions for case of 1 % noise level with TV regularization
shown as a series of 2-D slices. Images reconstructed after 5 iterations (at lowest frequency 10 KHz) ...... 60
Figure 12: Reconstructed conductivity distributions for case of 3% noise level with TV regularization
shown as a series of 2-D slices. Images reconstructed after 5 iterations (at lowest frequency 10 KHz) ..... 61
ix
Figure 13: False positive percentage of finite elements that was determined based on an image intensity
threshold that ensures 95% of tumor fmite elements are correctly classified ........................ '" .. , ...... 61
Chapter IV: Paper III
Figure 1: A cylindrical phantom with two test objects ............................................................. 69
Figure 2: Reconstructed images with TRS at different heights of phantom ..................................... 70
Figure 3: Reconstructed images with CGLS at different heights of phantom .................................. 70
Figure 4: L-curve with TRS, CGLS, and best solution with Tikhonov .......................................... 71
Figure 5: Reconstructed images with CGLS at30% standard deviation .......................................... 72
Figure 6: L-curve with TRS and CGLS at 30% standard deviation ................................................ 73
Chapter V: Paper IV
Figure 1: Multi-frequency EIT system ................................................................................. 78
Figure 2: Block Diagram ofDSP based waveform generator ...................................................... 80
Figure 3: Schematic diagram ofDSP waveform generator ......................................................... 81
Figure 4: A modified Howland current driver ..................................... , ............................ '" ... 82
Figure 5: Block diagram of multiplexers ...................................................................... , ....... 82
Figure 6: Schematic diagram ofMCD based multiplexers ............... '" ....................................... 83
Figure 7: Block diagram ofNI-6251 DAQ board (Adapted from M series user manual-371022G-Ol ,
www.ni.com ) ............................................... , ... '" .................................. , .................... 86
Figure 8: Block diagram ofFFT-based phase sensitive detection ................................................. 86
Figure 9: CMRR variation verses frequency for EIT system and DAQ board ............. , ...... , ............. 89
Figure 10: Maximum reciprocity variation verses frequency ............ '" ....................................... 89
Figure 11: Mean SNR variation with frequency for three different cases; three-ring with adjacent-adjacent
Three-ring with opposite-adjacent and single ring with opposite-adjacent configuration ...................... 91
Figure 12: Measurements accuracy variation with frequency relative to full-scale ofDQA measurement. 91
Figure 13: a): Metal and plastic rods inside of saline filled cylindrical tank. b): reconstructed conductivity
images of metal and plastic rods as a series of 2-D slices using subtraction from empty tank at 50 KHz,
metal rod is in red and plastic rod is in blue color c): The original position of two objects relative to
electrodes in lower part of phantom. The metal rod was placed 1.2 em from electrode 1 and the plastic rod
at the same depth near electrode 9 ......................................................................................... 94
Figure 14: a): Cucumber and rod inside of saline filled cylindrical tank, b): reconstructed conductivity
images of cucumber and rod using subtraction from empty tank as a series of 2-D slices at 125 KHz, metal
rod is in red and cucumber is in blue color c): reconstructed permittivity images of cucumber and rod using
x
subtraction from empty tank as a series of 2-D slices at 125 KHz, metal rod is in blue and cucumber is in
red color d): The original position of two objects relative to electrodes in lower part of phantom. The
cucumber was placed close to electrode 3 and the metal rod at the same depth near electrode 11 ............ 96
Chapter VI: Paper V
Figure 1: Hardware block diagram of Multi-frequency EIT ...................................................... 103
Figure 2: a) Difference conductivity images of a metal rod with 1.2 cm diameter at 1 cm far from electrode
1 b) the actual position of metal rod relative to electrodes ............................................. '" ......... 106
Figure 3: a) Difference conductivity images of a metal rod with 1.2 cm diameter at 4 cm far from electrode
1 b) the actual position of metal rod relative to electrodes ......................................................... 107
Figure 4: a): Difference conductivity images of a metal rod with 1.2 cm diameter at 1 cm far from electrode
1 and 1.5 cm diameter plastic rod at same depth from electrode 9, both rods were 6 cm tall. b) Metal and
plastic rods inside of cylindrical phantom, c): The original position of two objects relative to electrodes in
lower part of phantom. The metal rod was placed 1 cm from electrode 1 and the plastic rod at the same
depth near electrode 9 ...... '" .. , . " ... '" ............................... , .................................. '" ......... 1 09
Figure 5: a) A piece of cucumber inside of saline tank, b) conductivity variation verses frequency for
cucumber from 1 kHz up to 1 MHz c) Reconstructed conductivity images of cucumber using Dogleg
approach at 125 KHz. The cucumber was close to electrode 3, d) position of cucumber relative to
electrodes ......... '" ........................... '" ....................................................................... 112
Figure 6: a): Reconstructed conductivity images of cucumber using PCA approach at 125 KHz. The
cucumber was close to electrode 3. Eleven different frequencies were used for multi-frequency analysis b):
Reconstructed conductivity images of cucumber using Gauss-Newton approach at 125 KHz ............... 113
Figure.7: a): Reconstructed conductivity images of cucumber using difference method at 125 KHz. The
cucumber was close to electrode 3. b): Reconstructed permittivity images. The cucumber was close to
electrode number 3 ...................................................................................................... 114
Appendix A: Paper VI
Figure 1: Variation of the elastic stiffness (indentation) of unirradiated and irradiated PV A-3C and PV A-
6C ......................................................................................................................... A-3
Figure 2: Variation of Tiand T2 values for unirradiated PVA-3C and PVA-6C for various field strengths
............................................................................................................................. A-4
Figure 3: Variation of Tland T2 values for irradiated PV A-3C and PVA-6C for various field strengths. A-4
Figure 4: Variation of electrical conductivity for unirradiated and irradiated PVA-3C and PVA-6C versus
frequency ......... '" ... '" ........ , .................. '" ................... , ................ '" ......... '" .............. A-4
Figure 5: Variation of electrical resistance for unirradiated and irradiated PVA-3C and PVA-6C versus
frequency ................................................................................................................. A-4
Figure 6: Variation of permittivity for unirradiated and irradiated PVA-3C and PVA-6C versus frequency
............................................................................................................
................. A-4
Xl
Abbreviations
CGLS
Conjugate Gradient Least Square
CMRR
Common-Mode-Rejection-Ratio
CT
Computed Tomography
DAC
Digital Analogue Converter
DAQ
Data Acquisition
DL
Dogleg
DSP
Digital Signal Processing
EIT
Electrical Impedance Tomography
FEM
Finite Element Method
GCV
Generalized Cross Validation
GN
Gauss Newton
LM
Levenberg-Marquardt
MIA
Multivariate Image Analysis
MRI
Magnetic Resonance Imaging
MCV
Microcontroller Unit
PCA
Principal Component Analysis
PET
Positron Emission Tomography
SNR
Signal-to-Noise
SS-MIA
Spatial-Spectral MIA
SVD
Singular Value Decomposition
TRS
Trust Region Subproblem
VCCS
Voltage Controlled Current Source
xii
CHAPTER I
Introduction
1.
Background
This thesis discusses hardware design and image reconstruction methods in Electrical Impedance
Tomography (EIT). EIT is an imaging modality which reconstructs images of the distribution of the
electrical properties of an object from boundary measurements. In this modality an array of electrodes is
attached to the surface of the object under investigation, and electrical currents are passed to the electrodes.
The electrical current is an alternating current whose amplitude is between I-lOrnA and frequency between
1-100 kHz. The stimulation results in electric potentials that are measured on the surface of the object using
the same or additional electrodes. The voltages measured on the boundary are a function of the object's
electrical properties and the electrical current. An estimation of the spatial distribution of the electrical
properties inside of the volume is reconstructed by implementing varying current injection strategies and
voltage measurement sequences. In a practical situation, the simultaneous measurement of both the
amplitude and the phase of boundary voltages can result in images of electrical conductivity and permittivity
distributions inside the volume.
EIT has numerous applications in both industrial and medical areas. In industry, EIT has been used to image
fluid flows in pipelines and as a non-destructive test tool such as crack detection (Dickin
et al., 1996;
Alessandrini
et al., 1998). A variation of the technique has also been used for landmine detection by the
military (Wexler
et al., 1985). In medical applications, EIT is used for monitoring cardiac function,
monitoring brain function, detection of haemorrhage and possibly stroke, gastric imaging, breast cancer
detection, and functional imaging of the thorax (Cherepenin
et ai., 2002; Eyuboglu
et ai., 1989; Gibson
et
al., 2000; Bagshaw
et al., 2003; Kunst
et al., 1998; Smallwood
et al., 1994). It is the interest of this thesis to
develop EIT for medical imaging. Ultimately our laboratory's focus will be the detection and possibly
diagnosis of breast cancer.
This is motivated by the great potential of EIT to be used in .clinical environments as a diagnostic and
monitoring modality. However currently, due to several limitations of EIT, it has not yet been used as a
routine medical diagnostic tool in everyday clinical practice. Primarily, EIT suffers from a low spatial
resolution, and is susceptible to noise and electrode errors. The spatial resolution of EIT is relatively poor in
comparison to other modalities, such as, magnetic resonance imaging (MRJ) and Computed X-ray
Tomography (CT). In comparison however, EIT is a technique that does not use ionizing radiation, it is safe
and non-invasive, and is very portable. Moreover, it is much more inexpensive compared to other imaging
modalities such as MR!, CT and positron emission tomography (PET).
The mathematical formulation and uniqueness of the EIT problem was first addressed by Calderon
(Calderon, 1980). Since then several different approaches for solving the EIT problem have been proposed.
A problem is called a well-posed problem if it has a unique solution (uniqueness) and if the solution depends
continuously on the data (stability). In other words, small changes in the data must not cause instability in
the solution. A problem is called ill-posed if at least one of the three conditions (existence, uniqueness,
stability) is not satisfied. The EIT problem is categorized as an "ill-posed" problem (Hansen ,1998).
The EIT reconstruction problem is a nonlinear ill-posed inverse problem, and special approaches must be
implemented to recover a stable solution. As EIT is an ill-posed problem, small perturbations in the
measured boundary voltages can cause arbitrarily large errors in the estimated internal electrical
conductivity. One of the most common techniques for EIT image reconstruction is the minimization of the
squared norm of the difference between the measured boundary voltages and the calculated boundary
voltages through a mathematical model (Vauhkonen, 1997). However, because of the ill-posedness of the
EIT problem, this optimization process has to adopt particular techniques in order to obtain a stable solution.
This modification is called
regularization, a process which introduces additional terms into the minimization
(Borsic, 2002). This modification removes the ill-posedness of the problem. When the problem is
regularized by introducing this additional term, a
priori information about the solution is introduced into the
reconstruction process.
The EIT reconstruction algorithms can be classified into several categories; static/absolute imaging,
dynamic imaging, multi-frequency imaging, and difference imaging. Each of these is intended to image a
different aspect of the object's conductivity distribution. The aim of absolute imaging is to reconstruct an
absolute conductivity distribution that is achieved by using only a single data set of boundary voltage
measurements. This kind of image reconstruction needs accurate calculation of simulated voltages for
known conductivity distributions. The calculation of the boundary voltages when the injected currents and
the conductivity distribution are known is referred to as a
forward problem. In dynamic imaging, the time
variation of the conductivity distribution is included into the reconstruction process, and reconstruction is
updated after each current injection (Vauhkonen, 1997). In difference imaging, the difference between
conductivity distributions are created as a result of the difference between two data set measurements
corresponding to two different object conductivity distributions. The purpose of multiple frequency imaging
systems is to reconstruct the frequency variation of conductivity distributions of the target (Griffiths and
Zhang, 1989).
2
The utility of EIT for medical applications comes from the large variation in the tissue conductivity or
resisitivity of different organs. As shown in Table 1, there is a huge dynamic range of resistivity values for
different organs enabling the creation of high contrast images. One such interesting application of EIT for
medical purposes is the detection of significant changes in the resistivities caused by circulating blood or
inspired air. For these situations the difference imaging method would be an ideal choice as the electrode
artefacts are eliminated. Examples where difference imaging has been used for medical applications are:
gastric imaging (Barber 1990, Dijkstra
et at., 1993; Smallwood
et at., 1992) the detection of intrathoracic
fluid volumes (Newell
et al., 1996), the estimation of cardiac and pulmonary parameters (Brown
et aI.,
1992; Eytiboglu
et at., 1989; Harris
et at., 1992; Newell
et ai., 1992) and the detection of haemorrhage
(Sadleir
et ai., 1992; Murphy
et ai., 1987; McArdle
et aI., 1988). In the following sections, the physics of
the EIT difference image reconstruction method is discussed.
Table 1: Resistivity values for mammalian tissues (Barber and Brown, 1984).
Tissue
CSF
Blood
Liver
Human arm (longitudinal)
Human arm (transverse)
Skeletal muscle (longitudinal)
Skeletal muscle (transverse)
Cardiac muscle (longitudinal)
Cardiac muscle (transverse)
Neural tissue
Gray matter
White matter
Lung( out-in)
Fat
Bone
Resistivity (Om)
0.65
1.5
3.5
2.4
6.75
1.25
18.00
1.6
4.24
5.8
2.84
6.82
7.27-23.63
27.2
166
2.
Forward Model and Data Collection in EIT
From the perspective of EIT image reconstruction, a physical model is required to relate the internal
conductivity distribution to the boundary voltage measurement. This implies a relationship between the
measured voltages, the injected currents, and the conductivity distribution. The equation for the EIT
physical model has been derived from Maxwell's equations (Malmivuo and Plonsey 1995; Doerstling,1995).
In order to reconstruct the conductivity distribution through the EIT inverse solution, the forward solution is
required first. The forward problem is to calculate the boundary voltages when the conductivity distribution
and the injected currents are given. Of course initially, this conductivity distribution is the starting value
specified by a particular model. In the following section the different physical models that are used in the
3
forward solution are discussed (Cheng
et at., 1989; Pidcock
et at., 1995a; Somersalo
et at., 1992; Pidcok
et
at., 1992b). The most accurate model that has been proposed is the so-called complete electrode model
(Polydorides, 2002a). In this thesis the complete electrode model is used as the forward model and it is
solved using the [mite element method (FEM). FEM is a feasible technique for solving partial differential
equations with complex geometries (Brenner and Scott, 1994; Miller and Henriquez, 1990; Hinton and
Owen, 1979). Lastly, the different current injection and voltage measurement strategies used in EIT are
discussed.
2.1 Physics of the Problem
In the EIT experiment, an array of electrodes, usually 16 or 32, are attached around the object, n, and a
small alternating current of the order of a few rnA with a frequency between 1 kHz-IOO kHz is applied to a
subset of electrodes and the resulting voltages are measured on the remaining electrodes. In the following,
the physical models for ElI with different boundary conditions are derived. A mathematical model of the
problem is derived from Maxwell's equations
aB
VxE=--
at
aD
VxH=J+-
at
(I)
(2)
Where V x is the curl operator, .!. partial derivative with respect to time, E is the electric field, B is the
at
magnetic induction, H is the magnetic field, J electric current density and D the dielectric displacement.
Moreover, in a linear isotropic medium the following relations are valid
D=£E
(3)
B=,uH
(4)
J=aE
(5)
In EIT it is also common practice to further simplify Maxwell's equations. The quasi-static approximation is
usually used. The injected currents are assumed to be "slowly varying," so the time-dependence is
neglected. Since the frequencies of injected alternating currents are small the quasi-static approximation is a
valid assumption (Vauhkonen, 1997). In the quasi-static approximation the time derivative in equations (1)
and (2) are zero. Moreover, the electric field is conservative and the conduction currents are dominant with
respect to the displacement currents, such that
VxE=O
(6)
4
V'xH=J
(7)
Since V x E = 0 , there is electrical potential
u such that
E=-V'u
(8)
Taking the divergence of both sides of equation (7) gives
V' . J = V' . (V' x H) = 0
(9)
Using equation (5)
V' . (IE = 0
(10)
Finally using equation (8), it gives
V'·(aVu) = 0
(11)
This is the equation that recovers the electric potential u inside the body Q. In order to account for electrode
interactions within the object, a reasonable and appropriate model should be developed. In the following
section the boundary conditions, i.e., the electrode models are briefly discussed.
2.2 Electrode Models
In this section a brief description of the four most common electrode models are discussed. This thesis will
utilize the Complete Electrode Model which is considered to be the most accurate description of the
electrode and body interaction.
Continuum Electrode Model
The Continuum Model is the simplest of the models used in EIT. The model assumes that there are no
electrodes facing the boundary of the object, but the model assumes that the current density,
j, is a
continuous function on the entire boundary of the object. In this case, the relation
J .nl = -J .nl
inside
outside
on an
is valid. Here n is the normal vector to the boundary of the object,
an.. Furthermore, using equation 8
E =
-V'u , equation (12) can be recast as a Neumann boundary condition:
(I au = -J . n ==
j
an
(12)
(13)
where
j is the negative nonnal component of injected current. Equation (13) together with equation (11) is
called the
continuum model.
5
Gap Model
The gap model assumes discrete electrodes on the surface of an object. The injected
currentj is represented
as
. au I[
}=CY-=-
an A
. au 0
}=CY-=
an
on &[ ,
1= 1,2, ... ,
L
(14)
where A is the area of the electrode,
II is injected current into the
rth electrode and L is the number of
electrodes. The shunting effect of the electrodes and contact impedance are ignored in both the continuum
and the gap models.
Shunt Model
The shunt electrode model refines the Gap model by taking into account the shunting effect of the electrode
(Holder, 2005). It assumes the potential on the electrode is constant. Also the model states that the net
current density through the surface of electrode has to be equal to the total injected current
f
au
CY-ds =11
an
one,
,1=1,2, ... ,L
(15)
on 8/ ,
I = 1,2, ... ,
L
where
VI is the measured voltage on the
rth electrode. This model ignores the contact impedances of
electrodes (Somersalo
et ai, 1992). The model is completed by defming an arbitrary choice of ground to
ensure existence and uniqueness of the result
(16)
Complete Electrode Model
The Complete Electrode Model is a refmement of the shunt electrode model in which the contact impedance
between the electrodes and the object is taken into account (Paulson
et al., 1992; Somersalo
et al., 1992).
The contact impedance layer exists between the surface of the metal electrode and the object either in
medical applications or in phantom studies with ionic solutions. The effect of the contact impedance
between the electrode-surface interface is a voltage drop when the voltages are measured on the current
carrying electrodes (Hua
et aI., 1991). This modifies the shunting effect so that the assumption of constant
voltage under the electrode is no longer valid. If both the shunting effect of the electrodes and the contact
impedance are taken into account the
complete electrode model is obtained. This model has improved the
6
accuracy of the simulated voltages to be comparable to the precision of the data acquisition system (Cheng
et al., 1989). In this model the equation (15) is replaced by following equation
au
u+zzcr-=Vz oncz ,I =
1,2, ... ,L
au an
fcr- ds =1/
on
c/ ,
1= 1,2, ... ,
L
an
(17)
Where z/ is the contact impedance between the l'th electrode and object. The complete electrode model
consists of the following equations and the boundary conditions
v .
(oVu) = 0
on n
f
au
a-ds =1{
an
on
c[ ,
1= 1,2, ... ,
L
onC{ ,1=1,2, ... ,L
L
on aQ\UC[
1=1
(18)
(19)
(20)
(21)
Equation (I9) states that the surface measured voltages on each electrode consist of the voltage on the
boundary surface underneath that electrode plus the voltage dropped across the electrode impedance.
Equation (20) states that the integral of the current density over the electrode surface is equal to the total
injected current to that electrode. Finally, equation (21) means that there is no current entering or leaving the
boundary of object on the inter-electrode gap. In addition, the following two conditions for the conservation
of charge and an arbitrary choice of ground are needed to ensure existence and uniqueness of the result:
(22)
(23)
2.3 Numerical Solution of the Forward Problem
In order to solve the inverse problem in EIT, the physical model
(complete electrode mode!), which is called
the forward part, has to be solved. An analytical solution for the forward problem is only available in simple
cases, such as circular geometries with homogenous conductivity distributions (Fuks
et al., 1991; Pidcock
et
al., 1995a). However, in most practical situations especially medical applications, the geometry of object is
sufficiently complex that numerical methods are required. The most popular and feasible numerical method
7
for solving the partial differential equations with complex geometries is the finite element method (FEM). In
comparison to other methods, such as the finite difference method (FDM), the complex geometries can be
implemented relatively easily with FEM. In contrast, the implementation of complex boundary conditions is
not easy when using [mite difference methods. Furthermore, in FEM the unknown functions can be
approximated inside each element by linear or higher order polynomials. In FEM, the geometry of the object
is subdivided into a [mite number of discrete elements. The elements are triangles in two dimensions (2D)
and tetrahedra or hexahedra in three dimensions (3D). Each element consists of nodes and faces.
Figure 1
shows a 3D mesh constructed of tetrahedra with electrodes attached on the surface of a cylinder.
The procedure in FEM for the solution of partial differential equations starts to form the so-called
variational form, or weak form, of the original equation. In this procedure, the residual of the original partial
differential is formed first. Then this residual is multiplied by a test function
vex) and integrated over the
domain. The next step is to find an approximate solution of the variational problem. The solution
u(x) of the
10
8 ....
6
-E
(.)
N 4
2
0
5
Y(cm)
X(cm)
FIG. 1: 3D meshes for a cylinder with circular electrodes attached on surface. The mesh was
generated with NETGEN mesh generator.
partial differential equation is represented as a linear combination of piecewise polynomial interpolation
functions
8
N
U(X) =
Luitpj(x)
j""l
h
{
Ion vertex
i
were
tp.=
1
0 otherwise
(24)
Where tpj (x),
i = 1,2, ... ,
N are the piecewise linear basis functions which fonn the basis of the discrete
solution space,
Ui are the value of potential at i'th vertex, and N is the number of the vertices. In the weak
fonnulation of FEM, the test functions v are also selected to be the same functions as used to approximate
the potential,
N
vex) =
L VjlPj(X)
j""l
(25)
In the following, the weak fonn for the complete electrode model is derived. Much of the derivation in the
following section was interpreted from (Holder 2005; Polydorides, 2002a; Vauhkonen, 1997). By
multiplying equation (11) by test functions v and by integrating over the domain n
fVV'. (oVu) dV = 0
n
is obtained. Using the vector derivative identity
V' .
(vaVu) =
aVu . V'v +
vV'.(oVu)
Substitution of (27) into (26) yields
fV'· (voVu) dV -
foVu. V'v
dV = 0
n
n
(26)
(27)
(28)
Invoking divergence (Gauss) theorem on (28) enables a fonn for the implementation of boundary conditions
J oV
U . V'v
dV =
f (oV
u . n)v
dS =
f () :: v
dS
n
an
an
(29)
The boundary integral on the right side of equation (29) can be separated into two sections, the domain
L
under each electrode, 5[ and domain between the electrodes an \ U 5[ • In this case, right side of equation
1=1
(29) can be written in following fonn
f au
f au
L
fau
()-vdS=
()-vdS+ L ()-vdS
an
an
L
an
1=1 e
an
an\Ue/
I
(30)
1=1
o
9
The first part of the integral vanishes because of the boundary condition in equation (21). Furthennore,
using the boundary equation (19) this tenn can be rearranged as
au 1
u-=-(V/ -u)
an
Zt
Now inserting equations (30) and (31) into (29) gives
L
1
foVU. VvdV =
L f-:;-(Vt -u)vdS
n
1=1 &/
1
(31)
(32)
Substituting the test function from equation (25) and the approximate solution in equation (24) into equation
(32) gives
(33)
The combination of the boundary equations (19) and (20) and using the approximate solution from equation
(24), the total current has the following fonn
1
I
N {f~qJj
dS}
I, =
I-(Vt -u)dS= I-vz dS -
L Zz
uj
z,
Zz
.-1
&,
&,
&,
1-
where Ie 11 is the electrode area. For the FEM model, a piecewise constant basis function is used to
approximate the conductivity on each simplex as
where
K
U =
L uJJj (x)
i=1
on element
i
otherwise
(34)
(35)
(36)
In this case,
Clj can be taken outside of the integral for each element. For a set of current
injection patterns,
I
K
, the forward model with the complete electrode model is formulated as a system of linear equations as
follows:
10
(37)
In equation (37), Al is an (NxN) symmetric matrix called the global stiffness matrix
K
Al (i,}) = JoVqJj
,VqJjdV= LO"k fVqJj
,VqJjdV
(38)
o
k~
Ok
The integral is taken over the volume of each element. For convenience the integral part is called local
stiffness matrix
s~ = Jv qJi . V
qJjdV
Ok
(39)
The other components of equation (37) are as follows:
(40)
(41)
i =
j
i,
j =
1,2, .. .L
(42)
otherwise
IK is the vector of the injected currents, and U (Nxl) is the vector of nodal voltages, and VL (LxI) is the
vector of electrode voltages. The whole first set of matrices in equation (37) is called the global conductance
matrix, A with the size of
(N +
L x
N +
L ), Assembling the global conductance matrix, A, starts first with
the construction of global stiffness matrix, At, from the local stiffness matrices and augmenting A2 to the
global stiffness matrix. Following, A3 and ~ are added to form the complete electrode model FEM for the
forward part of the problem.
The other two boundary conditions for the complete electrode model, equations (22) and (23), are needed to
impose existence and uniqueness on the forward solutions (Somersalo
et al., 1992). The matrix form in
equation (37) is solved for a set of
d driving current patterns:
(43)
The different methods used to solve the system of equations in (37) will be discussed in the next section.
11
2.4 Algorithms to Solve the Forward Model
The forward linear equation (37) can be expressed in a condensed form as:
Au=b
(44)
with u = [U VLY E
9{N+L and b = [0 IK r E
9{N+L. There are two general approaches to solve
equation (44): direct or iterative methods. The matrix A is square, sparse, symmetric and positive defmite.
Positive definiteness allows a matrix to be decomposed into two triangular factors using Cholesky
factorization. A square matrix A with non-zero pivots can be changed as the product of a lower matrix Land
an upper matrix U; this is called LU-factorization. Furthermore, since A is positive defmite, the U can be
chosen to be the transpose ofL, which is called Cholesky decomposition (Golub and Van Loan, 1996). With
Cholesky decomposition, we can solve Au = b by first computing the Cholesky decomposition A = LL T,
then solving Ly = b for y, and finally solving LTu = y for u. The Cholesky decomposition is twice efficient
as the LU decomposition. The factorization process has a computational cost 0 (n3
). The computational
benefits of iterative procedures are more profound in three dimensional cases. Iterative processes such as the
conjugate gradient (CG) method shows a cost of 0 (n2k) per iteration and requires less than n iterations to
converge (Nocedal and Wright, 1999). In iterative process such as CG, the solution is sequentially updated
at each iteration through minimization of the least square residual
min~IIAu -bll~
2
u
(45)
The CG method searches for a minimum of the above objective function by taking conjugate search
directions for every iteration and require only the computation of the gradient of the objective function. The
convergence rate of CG can be improved by a process which is called preconditioning. The application of
CG in EIT problem is discussed in detail in (Monlinari
et al., 2002; Polydories
et al., 2002b; Player
et ai.,
1999).
2.5 Current Patterns and Measurement Procedure
The common approach in EIT systems is to inject current and to measure voltages on the surface electrodes
instead of applying voltages and measuring currents. This approach reduces the effect of contact impedance
on voltage measurements due to the high input impedance of the voltage measurement system. The EIT
system can be categorized based on the number of current sources either as a single current driver system or
a multiple current driver system (Holder, 2005). In single source systems, a single current source generates a
current flow in the object using a pair of electrodes. The resulting voltages are then measured between pairs
12
of the remaining electrodes. The current source is then switched to another pair of electrodes and the voltage
measurement repeated until a complete set of measurements has been collected. In multiple current drive
systems, there are N current sources, one for each electrode. Multiple drive systems have the ability to drive
current in more than two electrodes at a time. From the voltage measurement perspective, there are two
main strategies: two-electrode and four-electrode. In the two-electrode method, the voltages are measured
from the same current injection electrodes. In the four-electrode method, the current is injected through a
pair of electrodes while the voltages are measured from another pair. The error in the voltage measurement
due to contact impedance is much lower in the four-electrode than the two-electrode method.
The most commonly used four-electrode injection approach is the
adjacent drive method also known as the
neighbouring method (Barber, 1989; Hua
et at., 1987; Webster, 1990). As it is shown in
Figure 2, current is
passed through two adjacent electrodes and the voltages differences are collected from successive pairs of
adjacent electrodes. This is repeated through the next pair of electrodes and the same procedure is repeated
for all the electrodes. From all possible voltage measurements, the numbers of independent measurements
that can be used in image reconstruction are limited to a certain number. For N electrodes, the total numbers
of adjacent measurements are N2
• In the four-electrode measurement strategy, the voltage is not measured at
a) First drive pair
"
,
.f 10 '"
b) Second drive pair
FIG. 2: Adjacent current driver. The current is passed through two neighbour electrodes and the
voltages differences are collected from successive pairs of adjacent electrodes.
a current injecting electrode thus the number of measurements is reduced to N(N -3). However, the numbers
of independent measurements are even lower due to the reciprocity theorem which states that reversing the
voltage measurement and current injection electrodes would give an identical value of the impedance
(Geselowitz, 1971; Webster, 1990). Therefore only N (N - 3)/2 of the measurements are independent.
13
In the adjacent drive method, most of the injected current travels near the boundary of the object. The
method is therefore very sensitive to conductivity variations near the surface of the object. This will reduce
the sensitivity of the voltage measurements to the conductivity changes at the centre of the object.
Another common four-electrode method is the so-called
opposite method (Avis and Barber, 1994). In this
approach current is injected through electrodes that are 180
0
apart while voltage differences are measured
with respect to one reference electrode adjacent to the current electrode. In the opposite method, the current
distributions are more uniformly distributed and hence results in good sensitivity. The opposite pattern is
less sensitive to conductivity variations close to the boundary in comparison to the adjacent injection. The
opposite pattern offers fewer current injections than can be applied for the same number of electrodes in the
adjacent strategy.
Techniques to devise an optimal current pattern in EIT have been studied by several researchers. Seagar
(1983) studied the optimal placing of a pair of point drive electrodes on a disk in order to maximize the
measured voltage differences between a homogeneous background and an offset circular anomaly.
Optimizing the measured voltages with respect to the conductivity distributions has been addressed by
means of
distinguishability. Distinguishability was defined (Isaacson, 1986) as the ability to distinguish
between two conductivity distributions 0' 1 and 0'2 up to a precision of E if there exist a current pattern I with
11111=1 for which
(46)
where R(O') is the nonlinear function associated the resistivity to the boundary voltage measurements. Gisser
et at (1987) defined the optimal current pattern for distinguishability as an optimization problem such that
IIR(al)I -
R(a2)III
I opt = arg max
11111
(47)
The optimal current pattern that maximizes equation (46) is calculated as the eigenvector of the operator
(R(al ) -
R(a2 ))2 associated with the largest eigenvalue. The next best current pattern is the one associated
with the second largest eigenvalue. All the eigenvectors for which the distinguishability reaches the
accuracy threshold of the measurement system can be considered. The concept of optimal current in terms
of distinguishability has been further studied in (Koksal
et aI., 1995; Cheney and Isaacson, 1992).
As an example of an optimal current pattern to distinguish a central circular anomaly embedded in a larger
circular object are the so-called
trigonometric current patterns (Isaacson, 1986). In this pattern, current is
injected on all electrodes and voltages are measured on all electrodes.
14
j
COS(k 27r/)
1= 1, ... ,
L k = 1, ... ,
L / 2
Ilk =
NL 27r1
sin«(k - -) -)
1= 1, ... ,
L k =
L / 2 + 1, ... ,
L-1
2 N
(48)
From a practical perspective, the application of such an optimal current requires multiple current drivers
which make the whole system more complex and expensive. The use of multiple current sources requires
higher precision current sources than when using a single source, as unmatched gains would increase the
common-mode current, i.e. the sum of all current (Cook
et al., 1994; Zhu
et al., 1994). In contrast, in a twin
drive system any gain error will affect all measurements due to the repetition through the set of patterns. As
pointed out in (Borsic, 2002) the use of multiple drive optimal patterns does not necessarily guarantee a
better accuracy over pair drive systems. The other important consideration for medical applications
(Eytiboglu and Pilkington, 1993) is that the total amount of applied current is restricted due IEC safety
regulations (lEC, 2007), so that the concept of distinguishability should be maximized with respect to this
constraint.
3.
Image Reconstruction Algorithms
One of the most difficult and challenging aspects of EIT is solving the inverse conductivity problem. As
discussed above, the forward problem is well posed and straightforward to solve, however the inverse part is
both non-linear and ill-posed. As defined by Hadamard (Hadamard, 1902) a problem is "well-posed" if it
has three properties: a) Existence- a solution exists, b) Uniqueness - the solution is unique, and c) Stability -
the solution is continuously dependent upon the observed data. According to Hadamard's criteria, a
problem is ill-posed ifat least one of the three conditions (existence, uniqueness, stability) is not satisfied.
There are broad ranges of techniques for solving ill-posed problems with the majority of them dealing with
stability and uniqueness by implementation of a
priori information through a technique called
regularization. The EIT image reconstruction methods can be divided into two major streams: stochastic
and deterministic. The stochastic method is based on probabilistic-statistical reasoning (Kaipio
et al. , 2000;
Kaipio
et al., 2003). The deterministic approach can be subdivided to linear (Mueller
et al., 2002; Barber,
1989) and non-linear methods (Molinari
et al., 2001 ; Borcea, 2001). However in this thesis the main focus is
on the deterministic approach. Before a discussion of the methods used in EIT image reconstruction, some
general background of the linearization and computation of the Jacobian are given.
3.1 Linearization and the Jacobian
With the aid of the physical modeling from the previous section, a forward (discrete) operator F which maps
the conductivity to the boundary measurement can be given as
v = F(O") = R(CT)J
(49)
15
where V E m
LxK is a vector of the boundary voltages on the L electrodes and K is the number of current
patterns and
F(a) E
mLxK
is the (discrete) forward operator (equation 19-22). The matrix R(O") is the part
of the inverse of the matrix, A, that was used in the FEM calculations (equation 37), previously. The
measured boundary data depends nonlinearly on the internal conductivity distribution, so that there are two
different approaches to solve this type of nonlinear problem. Either the forward operator is replaced by a
linear approximation at
ao and only the differences, ~a are calculated or some iterative algorithm, such as
the Gauss-Newton is applied (Holder, 2005).
In the linearized approach in which only ~a are calculated, the forward solution is needed only for the
computation of the so-called
Jacobian. The Jacobian is defined as first derivative of the forward
operatof, J =
F'(a). To reconstruct the conductivity differences, ~a, two separate measurements are
needed. For static image reconstruction the forward problem has to be solved in each iteration. However, in
both techniques the Jacobian of the discrete forward operator,
F(a) is required. The Jacobian relates the
changes in the boundary measured voltages due to small changes in the conductivity distribution of each
element, such that ~F: ~aln H
~Vl an. The Jacobian has the following form
av/ av! av!
-- --
aa1 aa2 aaN
avd av~
av~
----
aa1 aa2 aaN
J=
avi av~
(50)
----
av/ av;
-- --
In the early days of EIT and mostly for 2-D models, the calculation of the Jacobian matrix was done in a
way that each conductivity element was perturbed by small amount,
ba, then the corresponding change on
each electrode, bV was calculated through solving the forward problem. This had to repeated fOf all current
patterns k = 1, 2 ... K, and for all elements. By dividing the change in the voltages by the change in the
conductivity, an approximation for the Jacobian elements was obtained. However, when the finite element
method is three dimensional this becomes inefficient and time consuming.
16
Yorkeyand Webster (Yorkey and Webster, 1987) proposed two approaches for calculating the Jacobian.
The fIrst approach, the so-called
standard method, is based upon the matrix equation (44). In terms of the
FEM, each column element of the Jacobian can be obtained from
(51)
Invoking the chain rule the right hand side of equation (51) :
a(A-1b) =
_A-1 aA A-1b
au n
au n
(52)
The only derivative that must be calculated is the derivative of the stiffness matrix which can be calculated
as
:J
(53)
where
aAI(i, j)
1
f
..
---
= - - 2
Vrpj ' V rpj dV, I, } =1,2, ... ,N
au n
u n l',n
(54)
A more effIcient approach to calculate the Jacobian was presented as the sensitivity method (Breckon,
1990). In the sensitivity method the element of the Jacobian (derivative with respect to the conductivity) is
calculated based on two different current injections:
av~
f /
k
J1k n = --= -
Vu(1 )· Vu(J )dV
aUn
.ok
(55)
Polydorides (Polydorides, 2002a) describes an implementation of this product that is the most effIcient way
from the computational perspective to calculate the Jacobian.
In the linearized approach, in which only the difference is solved, F(u) is approximated with a fIrst order
Taylor polynomial at u 0 as
F( u) = F( u 0 ) + F I ( U 0 )( u -u 0 ) + D(llu - u 0 11 2 )
(56)
If the second order terms are dropped, equation (56) can be written in the form
17
(57)
where
5V =
F(c;") -
F(ao), oa =
a -
ao, and
J =
F'(ao) . In the noniterative approach the following is
solved.
minllJ
oa - 5V11
oa
(58)
where 8V is a difference between two successive boundary measurements. However, this minimization
tends to fail as the problem is ill-posed, so that the problem has to be regularized.
3.2 Regularization
Since the EIT problem is inherently ill-posed, the minimization process of the weighted residual is doomed
to failure (Holder, 2005). Due to the unstable nature of the problem, even small changes in the boundary
measurements (contaminatation with noise) can eventually result in wide changes in the final solution.
Under these circumstances, the final solution is no longer unique as there would be more than one solution
that can fit the data up to the measurement noise level. In this case,
regularization is an approach that
attempts to restore uniqueness and stabilize the solution.
A wide range of regularization techniques have been applied to EIT including truncated singular value
decomposition, maximum entropy, and a number of generalized least squares schemes including Tikhonov
regularization methods (Cheney and Isaacson, 1995). The most widely used regularization method is the
Tikhonov regularization. Tikhonov regularization imposes prior information regarding the solution. The
prior term is incorporated into the solution as an additional term in the least squares minimization. Based
upon the linear approximation in the previous section, instead of minimizing IIJ
Oa - 5V11 ' one minimizes an
expression of the form:
arg min
{IIJ 60" - b'V112 +
ailL 60"11
2
}
60"
(59)
where L is the regularization matrix, and the prior term is represented by IILoal12
which penalizes the
minimized function, along with a regularization parameter, a., which controls the amount of prior
information imposed over the whole minimization process. From the point of view of optimization, this is a
quadratic minimization problem that is guaranteed to have a unique solution for a. > O. The identity matrix
and the matrices corresponding to the first and second difference operators are the most common
regularization matrix that have been used in Tikhonov regularization (Hua
et ai., 1988, Woo
et ai., 1993).
The identity matrix is used in the classic Tikhonov regularization. The solution to equation (59) is calculated
from the regularized inverse
18
(60)
For large enough a the regularized solution is stable for noisy measured data. The most common approaches
to determine a are the L-curve (Lawson
et aI. , 1974), Generalized Cross Validation (Wahba, 1977), and
Discrepancy Principle methods (Morozov, 1984).
3.3 Iterative Method
Iterative methods that have been used in EIT reconstruction are mostly based on the modified Gauss-
Newton algorithm with Tikhonov regularization (Woo
et ai., 1993; Hua
et ai., 1991; Ruan
et ai., 1996; Woo
et a., 1992). The Tikhonov regularized form for the non-linear problem ofEIT can be written as
(61)
where V meas are the measured boundary voltages, F(
a) is the forward operator which simulates the voltage
measurements from an object with conductivity distribution
a , L is the regularization matrix, and
ao is the
initial estimate of background conductivity, a is a regularization factor. The non-linear solution of equation
(61) is solved iteratively using a linearized step at each iteration. The equation (61) can be solved
recursively via Newton's method as follows
0'
n+l = 0'
n +!10' n
!1O'n =-H-1VrJ)
= (JT J +
aLT L)-l [JT (V -
F(O'n)) -
aLT LO'J
(62)
where H is the approximation for the Hessian of the objective function. In this case, the regularization
parameter,a, is held constant during whole iteration process. In contrast, Levenberg (Levenberg 1944) and
Marquardt (Marquardt 1963) used a modified version of the Newton method called the damped Newton
method. In this way, a diagonal matrix, AI is added in addition to the regularization term, aLTL, but the A is
reduced to zero as a solution is achieved. When the current iteration is far from the solution, the AI term
dominates and the step direction approaches the steepest-descent direction. As
A ----j- 0 the effect of the AI
term diminishes and the direction is almost the same as the Gauss-Newton direction, which is implemented
for the [mal iterations.
It is obvious that the injected currents are not restricted only to the electrode plane, therefore the currents
spread out all over the three-dimensions. The 2D assumption for current distributions has been used since
the early days of EIT. This assumption simplifies the problem and demands less computational time. The
main difficulty with 3D models is computation. The number of [mite elements is so large for three
dimensional complex geometries compared to 2D models. 3D reconstruction algorithms used in medical
19
applications can be found in Gobel
et al (Gobel
et al.,1992), Metherall
et al (Metherall
et aI.,1996),
Polydorides (Polydorides, 2002a), Blue
et al (Blue
et aI., 2000) and Molinari
et al (Molinari
et al., 2002 ).
3D models need more efficient algorithms due to the large storage and computation time. The iterative
Gauss-Newton method is suitable for small-scale EIT problems. However in large 3D problems the large
number of [mite elements needed for sufficient accuracy makes the Gauss-Newton method unsuitable. In
this thesis, an efficient three-dimensional reconstruction algorithm is proposed.
An outline of the thesis contents is listed here:
Chapter (ii) (Paper I) This chapter describes a trust region implementation for image reconstruction of
conductivity changes in EIT. The main contribution of this chapter is the development of a fast and
reliable method for 3D EIT models which is demonstrated to be as good as or better than the existing
methods while being stable and repeatable. The dogleg method is presented for this objective. The
dogleg algorithm approximates a Levenberg-Marquardt step within the trust region of the model
function with a quadratic model. The comparison of Levenberg-Marquardt and dogleg algorithms is
presented using the reconstructed images. This comparison of two techniques suggests the
implementation of dogleg method could result in the reduction of the execution time to less then 50% of
that of Levenberg-Marquardt algorithm without any quantifiable loss of quality of reconstructed
images. This work is described in detail in the paper "Dogleg trust-region application in electrical
impedance tomography," by Mehran Goharian, Aravinthan Jegatheesan and Gerald R Moran, published
in the lOP Journal Physiological Measurement.
28 (2007) 555-572.
Chapter (iii) (Paper 11) The EIT image reconstruction is a non-linear and ill-posed inverse problem, which
requires regularization to ensure a stable solution. Traditionally the regularization matrix has been
selected to be the identity matrix, a diagonal matrix or an approximation of first or second order
differential operators. Such matrices guarantee stability but enforce smoothness in the inverse solutions
thus eliminating the possibility of recovering non-smooth solutions. In this paper, we propose a novel
approach to build a subspace for regularization using a spectral and spatial multi-frequency analysis
approach. The approach is based on the construction of a subspace for the expected conductivity
distributions using principal component analysis. It is shown via simulations that the reconstructed
images obtained with the proposed method are better than with the standard regularization approach.
The advantage of this technique is that prior information is extracted from the characteristic response of
an object at different frequencies and spatially across the finite elements. This work is described in
detail in the paper "A Novel Approach for EIT Regularization via Spatial and Spectral Principal
Component Analysis" by Mehran Goharian, Mark-John Bruwer, Aravinthan Jegatheesan,
Gerald.R.Moran, and John.F.MacGregor, published in the lOP Journal Physiological Measurement.
28
(2007) 1001-1016.
20
Chapter (iv) (Paper III) In this chapter a new approach for regularization of the ill-posed EIT problem is
presented. This approach is based on the trust region subproblem (TRS), which uses L-curve maximum
curvature criteria to find a regularization parameter. A comparison of the TRS method with conjugate
gradient least squares (CGLS) for an experimental phantom study is presented. Our results show that
both methods converge to the same point on the L-curve when the noise level is known. The TRS
algorithm has the advantage that it does not require any knowledge of the norm of the noise. This work
is described in detail in the paper "Regularization of EIT Problem Using Trust Region SubProblem
Method" by M.Goharian, M. Soleimani , A. Jegatheesan, and G.R. Moran , published in conference
proceeding of the 13th International Conference on Electrical Bioimpedance combined with the 8th
Conference on Electrical Impedance Tomography, Vienna, Austria, 2007 (Paper number : 72665).
Chapter (v) (Paper IV) EIT is an imaging modality which reconstructs images of the electrical properties
distribution of an object from boundary measurements. A series of electrodes are attached to a subject
in a transverse plane. These are linked to a data acquisition unit, which outputs data to a Pc. By
applying a series of small currents to the attached electrodes on surface of body a set of potential
difference measurements can be made from non-current carrying pairs of electrodes. Although the
physiological electrical conductivity contrast is very large, and acquired images of high definition
should be achievable, this has not been the case to date. This chapter introduces a novel EIT system
developed based on a digital signal processor (DSP) architecture which allows flexibility in design
options such as an arbitrary shape of signal generation and flexibility in the communication between
system components and data collection. The prototype system provides 48 electrodes with continuous
frequency selection from 0.1 kHz to 125 kHz. A software based phase-sensitive demodulation
technique is used to extract amplitudes and phases from the raw measurements. Signal averaging and
automatic gain control are also implemented in voltage and phase measurements. System performance
was validated using a Cardiff-Cole Phantom and a saline filled cylindrical tank. This EIT system offers
image reconstruction of both conductivity and permittivity distributions in three dimensions. This work
is described in detail in the paper "A DSP Based Multi-Frequency Electrical Impedance Tomography
System," by Mehran Goharian, Aravinthan Jegatheesan and Gerald R Moran, submitted to the lOP
Journal Measurement Science and Technology (2007).
Chapter (vi) (Paper V) In this paper, we discuss and present the results of a series of phantom experiments
to show the baseline imaging performance of our designed EIT system. The conductivity and
permittivity images for a piece of cucumber were reconstructed using four different approaches: dog-
leg, principal component analysis (PCA) based approach, Gauss-Newton, and difference imaging. In
the case of a single frequency approach, the dog-leg technique proved to be robust by converging
towards solutions with lower systematic error. However, in the case of multi-frequency analysis, a
PCA-based approach demonstrates a substantial improvement over using a Gauss-Newton technique at
21
a single frequency in terms of systematic error reduction. Our EIT system recovered the value of the
electrical conductivity and the location of the targets at the correct order of magnitude. This work is
described in detail in the paper "3D Multi-Frequency Electrical Impedance Tomography Imaging in
Phantoms" by Mehran Goharian, Aravinthan Jegatheesan, Mark-John Bruwer, Gerald.R.Moran, and
John.F.MacGregor, submitted to the lOP Journal Physics in Medicine and Biology (2007).
Chapter (vii) (Conclusions) The thesis results are summarized and, based on them, suggestions are made for
the possible directions of future work.
Appendix (Paper VI) This section presents the design and fabrication of a novel polymer based phantom
which is called polyvinyl alcohol (PV A)-based cryogel (PV A-C). PV A-C consists of water and the
polyvinyl alcohol that have been mixed to form a solution which is converted into a solid hydrogel. The
PV A-C phantoms are suitable for medical imaging such as MRI and Ultrasound with simulated normal
tissues and lesions. This section also presents the effect of radiation on the elastic stiffness, electrical
and MRI properties of PVA-C. The benefits of PVA-C based phantom are include: long-term stability
(electrical properties of tissue specimens change with time); the ability to perform long imaging
procedures without having to be concerned with subject motion or discomfort; and also to test and
develop novel medical procedures. This work is described in detail in the paper "Modifying the MRI,
elastic stiffness and electrical properties of polyvinyl alcohol cryogel using irradiation" by Mehran
Goharian, Gerald R. Moran, Kyle Wilson, Colin Seymour, Aravinthan Jegatheesan, Michael Hill,
R.Terry Thompson, and Gordon Campbell, published in the Nucl. Instr. and Meth. B (2007),
doi:l0.1016/j.nimb.2007.04.11 1.
22
References
Alessandrini, G., Rondi, L., 1998. Stable determination of a crack in a planar inhomogeneous conductor.
SIAM .J Math Ana1.30, 326-340.
Avis, N.J., and Barber, D.e., 1994. Image reconstruction using non- adjacent drive configuration. Physiol.
Meas. 16, A153- A160.
Bagshaw, A.P., Liston, A.D., Bayford, R.H., Tizzard, A., Gibson, A.p., Tidswell, A.T., Sparkes, M.K.,
Dehghani, H., Binnie, e.D., and Holder, D.S., 2003. Electrical impedance tomography of
human brain function using reconstruction algorithms based on the finite element method.
Neurolmage. 20(2), 752-764.
Barber, D.e., 1990. Quantification in impedance imaging. Clin Phys Physiol Meas, Suppl A. 11,45-56.
Barber, D.e., 1989. A review of reconstruction techniques for electrical impedance tomography. Med Phys.
16, 162-169.
Barber, D.e., Brown, B.H., 1984. Applied potential tomography. J Phys E: Sci Instrum. 17,723-733.
Blue, R.S., Isaacson, D., and Newell, lC., 2000. Real-time three-dimensional electrical impedance imaging,
Physiol Meas. 21, 1-12.
Borcea, L., 2001. A nonlinear multi grid for imaging electrical conductivity and permittivity at low
frequency. Inverse problems,
17 (2) 329-359.
Borsic, A., 2002. Regularisation Methods for Imaging from Electrical Measurements. PhD Thesis. Oxford
Brooks University. 2002.
Breckon, W. 1992. The problem of anisotropy in electrical impedance tomography. In Proc 14th Int Conf
IEEE Eng Med BioI Society. 1734-1735.
Breckon, W.R., 1990. Image reconstruction in electrical impedance tomography. PhD Thesis.Oxford
Polytechnic, UK.
Brenner, S.c., and Scott, L.R., 1994. The Mathematical Theory of Finite Element Methods. Springer.
Brown, B.H., Sinton, A.M., Barber, D.e., Leathard, A.D., and McArdle, F.l, 1992. Simultaneous display of
lung ventilation and perfusion on a real-time EIT system. In Proc 14th Int ConfIEEE Eng Med
BioI Society. 1710-1711.
Calderon, A.P., 1980. On an inverse boundary problem. In Seminar on Numerical Analysis and its
Applications to Continuum Physiscs. W.H. Meyer and M.A. Raupp, Brasilian Math. Society,
Rio de Janeiro, 65-73.
Cheng, K.S., Isaacson, D., Newell, le., and Gisser, D.G., 1989. Electrode models for electric current
computed tomography. IEEE Trans Biomed Eng. 36, 918-924.
Cheney, M., and Isaacson, D., 1995. Issues in Electrical Impedance Imaging, IEEE Computational Science
and Engineering. 2, 53-62.
Cheney, M., and Isaacson, D., 1992. Distinguishability in impedance imaging. IEEE Trans Biomed Eng, 39,
852-860.
Cherepenin, V.A., Karpov, A.Y., Korjenevsky, A.V., Komienko, V.N., Kultiasov, Y.S., 2002. Three-
dimensional EIT imaging of breast tissues: system design and clinical testing. IEEE Trans.
Med. Imaging 21(6), 662- 67.
Cook, R.D., Salunier, G.J., Gisser, G.J., Goble, J.C., Newell, lC., and Isaacson, D., 1994. ACT3: High-
speed, high-precision electrical impedance tomograph. IEEE Trans.Biomed. Imag. 41, 7l3-
722.
23
Dickin, F., Wang, M., 1996.Electrica1 resistance tomography for process tomography. Meas Sci Tech. 7,
247-260.
Dijkstra, A.M., Brown, B.H., Leathard, A.D., Harris, N.D., Barber, D.e., and Endbrooke, D.L., 1993.
Review clinical applications of electrical impedance tomography. J Med Eng & Technol, 17,
89-98.
Doerstling, B.H., 1995. A 3-D reconstruction algorithm for the linearized inverse bound- boundary value
problem for Maxwell's equations. PhD thesis. Rensselaer Polytechnic Institute, Troy, New
York.
Eyuboglu, B.M., and Pilkington, T.e., 1993. Comments on distinguishability in electrical impedance
tomography. IEEE Trans. Biomed. Eng. 40, 1328- 1330.
Eyuboglu, B.M., Brown, B.H., Barber, D.e., 1989. In vivo imaging of cardiac related impedance changes,
IEEE Eng Med BioI Mag. 8,39-45.
Fuks, L.F., Cheney, M., Isaacson, D., Gisser, D.G., and Newell, lC., 1991. Detection and imaging of
electric conductivity and permittivity at low frequency. IEEE Trans Biomed Eng, 38, 1106-
1110.
Geddes, L.A., and Baker, L.A., 1967. The specific resistance of biological material: a compendium of data
for the biomedical engineering and physiologist. Med .Biol. Eng. 5, 271-293.
Geselowitz, D.B., 1971. An application of electrocardiographic lead theory to impedance plethysmography.
IEEE Trans. Bio. Eng.18, 38-41.
Gibson, A., Bayford, R.H., Holder, D.S., 2000. Two-dimensional finite element modelling of neonatal head,
Physiol Meas.21 , 45-52.
Gisser, G., Isaacson, D., Newell, lC., 1987. Current topics in impedance imaging, Clin Physiol Meas. 8,
38-46.
Glidewell, M., and Ng, K.T, 1995. Anatomically constrained electrical impedance tomography for
anisotropic bodies via a two-step approach. IEEE Trans Med Imaging. 14, 498-503.
Gobel, J.e., Cheney, M., and Isaacson, J., 1992. Electrical Impedance Tomography in three dimensions,
Appl Comput Electromagn Soc l 7, 128-147.
Golub, G.H., and Van Loan, C.F., 1996. Matrix Computations. Johns Hopkins University Press.Baltimore.
Griffiths, H., Zhang Z., 1989.A dual-frequency electrical impedance tomography system, Phys Med BioI.
34, 1465-1476.
IEC, 2007. 60601-1 Medical Electrical Equipment: I-II. General Requirements for Basic Safety and
Essential Performance-Collateral Standard: Electromagnetic Compatibility- Requirements
and Tests 3rd edn
Isaacson, D., 1986. Distinguishability of conductivities by electric current computed tomography. IEEE
Trans Med Imag. 5, 91-95.
Hadamard, 1., 1902. Boundary value problems for partial differential equations and their physical
interpretation. Bull. Univ. Princeton, vol. 13, no. 49.
Hansen, P.e., 1998. Rank-Deficient and Discrete Ill-Posed Problems: Numerical Aspects of Linear
Inversion, SIAM, Philadelphia.
Harris, N.D., Brown, B.H., and Barber, D.C., 1992. Continuous monitoring of lung ventilation with
electrical impedance tomography. In Proc 14th Int Conf IEEE Eng Med BioI Society, 1754-
1755.
24
Hinton, E., and Owen, D.R., 1979. An Introduction to Finite Element Computations. Pineridge Press
Limited, Swansea, U.K.
Holder, D.S., 2005. Electrical Impedance Tomography: Methods, History and Applications. Institute of
Physics.
Hua, P., Woo, E.J., and Webster,
l .G ., 1991. Iterative reconstruction methods using regularisation and
optimal current patterns in electrical impedance tomography. IEEE Trans Med Imag. 10, 621-
628.
Hua, P., Webster, J.G., and Tompkins, W.J., 1988. A regularised electrical impedance tomography
reconstruction algorithm.Clin. Phys. Physiol.Meas,. 9 (Suppl. A) 137- 141.
Hua, P., Webster, J.G., and Tompkins, W.J., 1987. Effect of the measurement method on noise handling and
image quality ofEIT imaging. In Proc 9th Int ConfIEEE Eng Med BioI Society. 1429-1430.
Kaipio, J.P., Aku, S., Somersalo, E., and Haario, H., 2003. Statistical inversion approach for optimizing
current patterns in EIT. 3rd World Congress on Industrial Process Tomography. Banff,
Canada, September 2-5.
Kaipio, J.P., Vauhkonen, M., Somersalo, E., and Kolehmainen, V., 2000. Statistical inversion and Monte
Carlo sampling methods in electrical impedance tomography. Inverse problems. 16, 1487-
1522.
Koksal, A., and Eyiiboglu, B.M., 1995. Determination of optimum injected current patterns in electrical
impedance tomography. Physiol Meas, 16, A99-AI09.
Kunst, P.W.A., Vonk, Noordegraaf, A., Hoekstra, O.S., Postmus, P.E, de Vries P.M., 1998. Ventilation and
perfusion imaging by electrical impedance tomography: A comparison with radionuclide
scanning, Physiol Meas.19, 481-490.
Lawson, C.L, and Hanson, R.J., 1974. Solving least square problems Prentice-Hall, Englewood Cliffs, NJ.
Levenberg, K., 1944. A method for the solution of certain problems in least squares, Quart. Appl. Math. 2,
164-168.
Malmivuo, 1., and Plonsey, R., 1995. Bioelectromagnetism. Oxford university press, New York.
Marquardt, D., 1963. An algorithm for least squares estimation of nonlinear parameters SIAM J. Appl.
Math. 11, 431-4l.
McArdle,
F.l ., Brown, B.H., Pearse, R.J., Barber, D.C., 1988. The effect of the skull of low-birth weight
neonates on applied potential tomography imaging of centralized resistivity changes. Clin Phys
Physiol Meas. Suppl A, 9, 55-60.
Metherall, P., Barber, D.C., Smallwood, R.H., and Brown, B.H., 1996. Three dimensional electrical
impedance tomography, Nature. 380, 509-12.
Miller, C.E., and Henriquez, C.S., 1990. Finite element analysis of bioelectric phenomena. Crit. Rev. in
Biomed. Eng. 18, 207-233.
Molinari, M., Cox, S.J., Blott, B.H., and Daniell, G.J., 2002. Comparison of algorithm for non-linear inverse
3D electrical tomography reconstruction. Physiol.Meas.23, 95-104.
Molinari, M ., Cox,
S.l., Blott, B.H., and Daniell, G.J., 2001. Efficient non-linear 3D Electrical Tomography
reconstruction 2. 2nd World Congress on Industrial Process Tomography ed
Hannover,Germany.
Morozov,V.A., I 984.Methods for solving incorrectly posed problems. Springer, New York,Berlin,
Heidelberg.
25
Mueller, J.L., Siltanen, S., and Isaacson, D., 2002. A direct reconstruction algorithm for electrical
impedance tomography," IEEE Trans. Med. Imaging. 21(6) 555-559.
Murphy, D., Burton, P., Coombs, R., Tarassenko, L., and Rolfe, P., 1987. Impedance imaging in the
newborn. Clin Phys Physio} Meas, Suppl A. 8 ,131-140.
Newell, lC., Edic, P.M., Larson-Wiseman, lL., Ren, X., and Danyleiko, M.D.,1996. Assessment of acute
pulmonary edema in dogs by electrical impedance tomography. IEEE Trans Biomed Eng.
43
,133-138.
Newell, lC., Isaacson, D., Cheney, M., Saulnier, G.1., Gisser, D.G., Goble J.c., Cook, R.D., and Edic,
P.M.,1992. Impedance images of the chest. In Proc 14th Int Conf IEEE Eng Med BioI
Society.1752-1753.
Nocedal, 1., and Wright,
S.l. , 1999. Numerical optimization. Springer series in operations research.
Springer-Verlag. New York.
Paulson, K., Breckon, W., and Pidcock, M., 1992. Electrode modelling in electrical impedance tomography.
SIAM lournal on Applied Mathematics. 52(4), 1012- 1022.
Pidcock, M.K., Kuzuoglu, M., and Leblebicioglu, K., 1995a. Analytic and semi-analytic solutions in
electrical impedance tomography: 1. Two-dimensional problems. Physiol Meas. 16, 77-90.
Pidcock, M.K., Kuzuoglu, M., and Leblebicioglu, K., 1995b. Analytic and semi-analytic solutions in
electrical impedance tomography: II. Three-dimensional problems. Physiol Meas. 16, 91-110.
Player, M.A., Van Weereld, J., Allen, A.R., and Collie, D.A., 1999. Truncated-Newton algorithm for three-
dimensional electrical impedance tomography. Electronics Lett.35, 2189-219l.
Polydorides, N., 2002a. Image reconstruction algorithms for soft-field tomography. PhD Thesis.
UMIST.UK.
Polydorides, N., Lionheart, W.R.B., and McCann, H., 2002b.Krylov subspace iterative techniques: on
detection of brain activity with electrical impedance tomography. IEEE Trans.Med.Imaging.
21 , 596-603.
Ruan, W., Guardo, R., and Adler, A., 1996. Experimental evaluation of two iterative reconstruction methods
for induced current electrical impedance tomography. IEEE Trans Med Imaging. 15, 180-187.
Sadleir, R.J., Fox, R.A., VanKann F.1., and Attikiouzel, Y., 1992. Estimating volumes of intra-abdominal
blood using electrical impedance imaging. In Proc 14th Int Conf IEEE Eng Med BioI Society.
1750-175l.
Seagar, A.D., 1983. Probing with low frequency electric current. PhD. Thesis. University of Canterbury,
Christchurch, Newzeland.
Smallwood, R.H., Mangnall, Y.F., Leathard, A.D., 1994. Transport of gastric contents. Physiol Meas. 15,
175-188.
Smallwood, R.H., Nour, S., Mangnall, Y., Smythe, A.,and Brown, B.H.,1992. Impedance imaging and
gastric motility. In Proc 14th Int ConfIEEE Eng Med BioI Society. 1748-1749.
Somersalo, E., Cheney, M., and Isaacson, D., 1992. Existence and uniqueness for electrode models for
electric current computed tomography. SIAM J Appl Math. 52, 1023-1040.
Sylvester, 1. 1990. An anisotropic inverse boundary value problem. Commun Pure Appl Math. 43, 201-232.
Vauhkonen, M., 1997.Electrical Impedance Tomography and Prior Information, PhD thesis, department of
Applied Physics. Kuopio University.
Wahba, G.,1977. Ill-posed problems: Numerical and statistical methods for mildly, moderately and severely
ill-posed problems with noisy data. SIAM l .Numer.Anal. 14, 651-667.
26
Webster, J.G., 1990. Electrical Impedance Tomography. Adam Higler.
Wexler, A., Fry, B., and Neuman, M.R., 1985. Impedance-computed tomography algorithm and system.
Appl. Opt. 24, 3585-3992.
Woo, E.J., Hua, P., Webster, J.G., Tompkins,W.J., 1992.Measuring lung resistivity using electrical
impedance tomography. IEEE Trans Biomed Eng. 39, 756-60.
Woo, E.J, Hua, P., Webster, J.G., and Tompkins, W.J., 1993. A robust image reconstruction algorithm and
its parallel implementation in electrical impedance tomography. IEEE Trans. Med.
Imaging.12, 137-146.
Yorkey, T.J., and Webster, J.G., 1987. A comparison of impedance tomographic reconstruction algorithms.
Clin Phys Physiol Meas. Suppl A, 8, 55-62.
Zhu, Q.c., McLeod, C.N., Denyer, C.W., Lidgey, F.J., and Lionheart, W.R.B., 1994. Development of a real-
time adaptive current tomography. Physiol Meas. 15, A37-A43.
27
CHAPTER II
Paper I
The following paper examines the feasibility of implementing a trust-region based algorithm, the so-called
"dogleg," for the purpose of image reconstruction of conductivity distributions in electrical impedance
tomography. The dogleg algorithm approximates a Levenberg-Marquardt step within the trust region of the
model function with a quadratic model. In this paper a comparison between Levenberg-Marquardt and
dogleg algorithms is presented for different objects using simulated data.
The work presented in this chapter is the first paper in a series of papers which implements new algorithms
for image reconstruction of EIT. The simulation work presented in this paper was performed by me under
the supervision of Dr. Moran. The manuscript was written by me and edited by Dr. Moran and Mr.
J egatheesan.
(Reprinted from Physiol. Meas. 28 (2007) 555- 572 , Mehran Goharian, Aravinthan Jegatheesan and Gerald
R Moran, Dogleg trust-region application in electrical impedance tomography, Physiol.Meas.28 (2007)
555-572, Copyright (2007), with permission from lOP. http://www.iop.org/EJ/pm). Abstract is available at
http://www.iop.org/EJ/toc/0967-3334/28/5
28
lOP PUBLISHING
PHYSIOLOGICAL M EASUREMENT
Physiol. Meas. 28 (2007) 555- 572
doi: 1 0.1 088/0967-3334/28/5/009
Dogleg trust-region application in electrical
impedance tomography
Mehran Goharian1
, Aravinthan Jegatheesan2 and Gerald R Moran1,3
1 Medical Physics and Applied Radiation Sciences, McMaster University, Hamilton, Ontario,
Canada
2 McMaster School of Biomedical Engineering, McMaster University, Hamilton, Ontario, Canada
3 Diagnostic Imaging, Hamilton Health Sciences, Hamilton, Ontario, Canada
E-mail: goharim@mcmaster.ca. morang@HHSC.CAand jegatha@mcmaster.ca
Received 17 December 2006, accepted for publication 27 March 2007
Published 20 April 2007
Online at stacks.iop.org/PM/28/555
Abstract
This paper presents a trust-region implementation for image reconstruction
of conductivity changes in electrical impedance tomography. A dogleg
trust-region algorithm is applied in different cases to detect abnormalities.
The dogleg algorithm approximates a Levenberg-Marquardt step within
the trust region of the model function with a quadratic model. The
comparison of Levenberg- Marquardt and dogleg algorithms is presented using
the reconstructed images. This comparison of two techniques suggests the
implementation of the dogleg method could result in the reduction of the
execution time to less than 50% of that of the Levenberg-Marquardt algorithm
without any quantifiable loss of quality of reconstructed images.
Keywords: electrical impedance tomography (ElT), trust region, Levenberg-
Marquardt, dogleg, image reconstruction
(Some figures in this article are in colour only in the electronic version)
1. Introduction
Electrical impedance tomography (ElT) is an imaging technique that calculates the
conductivity
a (x, y, z) of an object from boundary voltage measurements. A series of
electrodes is attached to the boundary of an object under investigation, and alternating current
(or Voltage) is passed to a pair of electrodes. The induced voltages on the surface of the object
are measured through the non-current carrying electrodes. These boundary measurements
are used to calculate the internal conductivity distribution based on mathematical models that
relate the internal conductivity distribution to the boundary measurements. The calculation
of the boundary voltages when the injected currents and the conductivity maps are available
0967-3334/07/050555+18$30.00 © 2007 lOP Publishing Ltd Printed in the UK
555
29
556
M Goharian
et at
is called the
forward problem. Cheng
et al (1989) and Somersalo
et al (1992) have discussed
different models for EIT. The complete electrode model is the most accurate model (Vauhkonen
et at 1999). In this model the contact impedances between the electrodes ane the object and
the shunting effect of the electrodes are incorporated into the EIT mathematical model. The
most common approach to solve the forward problem is the finite element method (FEM).
The FEM is a feasible technique for solving the partial differential equation~; governing the
EIT problem using complex geometries and boundary conditions (Miller and Henriquez 1990,
Brenner and Scott 1994).
The EIT image reconstruction problem is considered to be an ill-posed irrverse problem
as defined by Hadamard (Groetsch 1993). Different techniques for solving inverse problems
in EIT have been proposed. One of the most common approaches uses iterati'/e least squares
to minimize the difference between the measured and the calculated boundar) voltages using
FEM. In order to obtain a stable solution for the ill-posed EIT inverse problem, the optimization
process used in the reconstruction has to use so-called regularization methods. A regularization
method stabilizes the solution by adding extra terms and constraints to the system and thereby
reduces the influence of noise. If these 'extra' constraints can be extracted from the actual
physical systems and used as an input, this would be referred to as
a prieri information.
Since the small singular values always correspond to high frequency components in the
image, regularization therefore acts on the image as a low-pass filter to remove them. The
most common types of regularization are the identity, and first- and second-order difference
matrices (Vauhkonen 2004). The regularized form of the EIT inverse problem can be written
as
min/(u) = ~IIV -
F(a)112 + ~a1lWa112
a
(1)
where
u is the true conductivity,
F(u) is the conductivity to voltages map given by FEM,
V is
the measured voltages,
W is the given regularization matrix and
a is the regularization factor.
Equation (1) can be linearized using a Gauss-Newton (GN) type optimization technique. The
recurrence relation for minimizing the function is
where
i3.an =
-H-I V /
= (JT J +aWTW)-I[JT(V -
F(un )) -aWTWun ]
(2)
H is the Hessian approximation for
f (a) in the GN method.
J is the Jacobian of the mapping
operator
F( a), and Y
n is called step length. The step length determines how far from the current
iteration one moves (in the direction of V
f) for the next iteration. The perlormance of the
GN method is highly dependent upon the condition number of the J ocobian matrix, especially
when the singular values of the Jacobian decay gradually to zero. In the case of the Hessian
approximation, the condition number of
JT J is double that of
J, so a small perturbation
(noise) from the measured value can cause instability in the solution, which may cause the line
search to fail to converge. If the Hessian is not positive definite and symmetric then this may
cause the line search to fail as well. The Gauss-Newton method has an oscillatory behaviour
during iterations and is not as robust as the steepest descent approach (Jacoby
et aI1972). The
steepest descent computes the gradient direction followed by a one-dimensional search that
approaches the minimum point of the objective function iteratively.
In this paper we use two different approaches instead of the GN to solve the nonlinear
equation in (1). The first approach is the
Levenberg-Marquardt (LM), and the s,!cond approach
is the dogleg (DL) method. In both cases the line search is replaced with a trust region strategy.
30
Dogleg trust-region application in electrical impedance tomography
557
2. Theory and algorithms
2.1. Levenberg-Marquardt (LM) method
Levenberg (1944) and Marquardt (1963) used a modified version of the Newton method called
the damped Newton method. In this case the step vector,
him, at each iteration is defined by
the following equation:
Cv2 f(x) +
AI)hlm = -\1
f(x).
(3)
This modification of the Hessian of the objective function has several effects:
(a) For A > 0 the modified Hessian is positive definite which ensures a descending direction.
(b) When the current iteration is far from the solution, the AI term dominates and the step
direction approaches the steepest-descent direction.
(c) As A ~ 0 the effect of AI term diminishes and the direction is almost the same as the
Gauss-Newton direction, which is implemented for the final iterations.
In the case of nonlinear least-squares problems, such as equation (1), the iteration has the
following form (Vauhkonen 2004):
ak+l =
ak +
J3.k
tlk =
(JT J + AI +
aWTW)-I[JT (V -
F(an» -
aWTWan ].
(4)
In this way the diagonal elements of
JT J are modified in a process called damping, where
A is the damping term. For example, if the damping term is small in the present iteration
that implies that the step will be calculated by the Gauss-Newton approach. In this way,
the LM technique has the ability to switch between slow convergence steepest descent when
it is far from the minimum and fast quadratic convergence when it is close to the solution.
An effective approach (Nielsen 1999) for updating the damping teIm is used in this work.
If we assume
f (x) :
ffi"n ~ ffi" is our objective function and just the first two terms of the
Taylor-series expansion of
f(x) are used to create the
model function at iteration
Xk (Nocedal
and Wright 1999), we have:
(5)
wherefk =
f(Xk), \1
fk = \1
f(Xk),P is the direction of the next step vector andBkis asymmetric
matrix. The updating is controlled using a
gain factor formula (5) which is used to make sure
that the actual decrease in the objective function is higher than that predicted by the quadratic
model in equation (5).
P
==
f(x) -
f(x +
p) > 8
L(O) -
L(p)
(6)
where 8 is a positive number. Using the
model function from equation (5) the denominator of
the gain factor in equation (6) is calculated as follows:
L(O) -
L(p) =
-hTmJT f -
!h!mJ
T
Jh1m
(7)
where
him is the step vector in the LM method. The higher orders in this calculation are
eliminated. On each iteration, if
p is close to one it means that the
model function is a good
representation of the objective function, and then A can be decreased for the next iteration.
The damping term,
AI, has the ability to deal with the rank deficiency of the Jacobian, and
therefore handles the singularity of the
JT J matrix. The following
criteria were used to specify
an endpoint in order to teIminate the algorithm (for the objective function in equation (1»,
Ilglloo ~ 81
g = \1
!(a) =
JT(V -
F(a» -
aWTWa
(8)
31
558
M Goharian
et at
Algorithm:
Input: objective function!; voltage measurements vector
V; starting values for conductivity 0'0;
Output: updated conductivity vector a
k:==O; v:==2; a= O'c;
Ilinitializationil
A:= J(a)T J(a)+aWTW; g:= J(a)T (V -F(a»-aWTWa;
stop := (11g 1'"-, ~ £1);
A:= r * max
i=I ... ,m (Aii) ;
while
(stop:ttrue) and
(k<kma)() do
k:=k+l;
Solve (A+Al)8o=-g;
//cholesky//
if 118all"" $ (£2110"11 + £2) then
stop:=true;
else
anew :=0'+80';
Compute p by (6-7)
if p>O then
else
A:== J(O"l J(O") +
aWTW ;
g :=
J(a/ (V -
F(a)) -
aWTWO";
stop:= (1Igll~ ~ £1);
1
A:= A * max{-,
1-(2p-l)3}; v:==2;
3
A :=
A * v ;
v:=v*2;
end if
end if
end while
Figure 1. The complete Levenberg-Marquardt algorithm.
where £1 is a small positive number defined by the user. The algorithm stops when the
following
criterion is met:
(9)
where £2 is a small positive number defined by the user. The LM algorithm is presented in
pseudo code in figure 1, and more detail can be found in Madsen
et al (2004). The starting
value for A is chosen based on the maximum element of
(JT J +
aWTW) as per Nielsen
(1999).
2.2. Dogleg method
The dogleg (Powell 1970) method combines the Gauss-Newton and the method of steepest
descent (similar to LM technique) in solving unconstrained optimization problems. In this
situation, this happens explicitly through the control of the radius of a
trust region. With the
trust region approach, an approximation of the objective function,
f, called a
model function,
Lk, same as equation (5) is constructed to represent the behaviour of the objective function
at each iteration (Nocedal and Wright 1999). The iteration steps are determined to minimize
the
model function inside the trust region. The size of the trust region is detennined based on
32
Dogleg trust-region application in electrical impedance tomography
559
the performance of algorithm in the previous iteration. The following subproblem needs to be
solved to obtain the search direction,
P (Nocedal and Wright 1999)
min
{Lk(Xk +
p) =
fk +
pT\1 fk +
!pT BkP}
P E
Rn
(10)
subject tollpli2 ( ~k
where ~k is called the trust region radius. The 11-112 is the L2-norm. The vector
p* is a
globally optimized solution to problem (10) for a scalar J.. ~ 0 if the following conditions are
satisfied (Nocedal and Wright 1999):
(i)
(B +
J..I)p* = -\1
fk
(ii) J..(~ - Ilplb) = 0
(11)
(iii)
(B +
AI) is positive semidefinite.
The candidate step is the solution to the constrained subproblem in equation (10). The
necessary and sufficient conditions for the solution are given in equation (11). The radius
of the trust region has a critical effect on the convergence of the optimization. It is selected
based on the success of the
model function (equation (5» in predicting the behaviour of the
objective function during each iteration. If the model correctly predicts the objective function
behaviour then the radius is increased to allow for bigger steps. In the case of failure to predict,
the radius is decreased and the subproblem (equation (10» is solved again using the reduced
radius. More and Sorensen (1983) suggested solving a nonlinear equation by the Newton
method using the sparse Choleski decomposition of IIplA)1I = �� with(B +
J..l)p(J..) = -\1
f.
This method requires two to three Choleski decompositions per iteration. There are other
methods which are based on the minimization of the subproblem in equation (10) in the
T
two-dimensional subspace containing the Cauchy step,
Psd =
-(g; 19)9 and the Newton step
pgn = -
B-1 g. The most common technique is the dogleg (DL) method which was suggested
by Powell (1970), where
p =
Pgn if
IIPgn II :::;; ~ and
p = CIP~dll
)Psd if
IIPsd II ) ~. In the
remaining case,
P is a convex combination of
Pgn and
Psd such that
IIpll = ~. The advantage
of the dogleg method is that it requires only one Choleski decomposition per iteration. It
should also be mentioned that in the case of the EIT least -squares problem, the matrices g and
B are given as follows:
g =
JT (V -
F(a» -
aWTWa
B =
J(a)T J(a) +
aWTW.
(12)
(13)
The dogleg step can be selected based on Powell's suggestion through the above-mentioned
algorithm listed as equation (14). The dogleg algorithm is presented in pseudocode in figure 2
and more details can be found in Madsen
et at (2004).
if
IIPsdll ~ fi then
Pdl := CP~dll)
P,d:
else if
I\Pgn 1\ :::;; ~ then
Pdt :=
Pgn;
else
Pdt:= Psd +
f3(pgn -
Psd);
Where
f3 is chosen to meet
I/Pdlli = ~
end
33
(14)
560
M Goharian
et al
Algorithm:
Input: objective functionf, voltage measurements vector V, starting values for conductivity 00
Output: updated conductivity vector, 0
K:=O; 0:= 0'0; ~:= ~o' GN:=true;
//initialization//
B:= l(erl l(er) +
aWTw; g :=
l(erl (V -
F(er) -
aWTWer;
stop = (1lglloo ~
Cl) ;
while
(stop*true) and
(k<kmax) do
k:=k+l;
P
'- IIgl12
g'
sd .- IIlgl12
'
ifGN=true then
P
.
B-1 •
gn'= -
g,
GN:=false;
end if
Compute
Pdl by (14)
/J.O":=
Pdl;
if II~all ~ (£211011 +
C2) then
stop:=true;
else
Onew :=0+ ~o;
Compute p by (6-7)
if p>O then
a:=anew ;
B:= l(al 1(0") +aWTW; g:= l(erl (V -
F(a)) -aWTWer;
stop:= «llgt ::; £1) or ellfex)t ::; £3»);
GN:=true;
if p>O.75 then
~:= max{~,3*llpdlll}
else if p<O.25 then
~:= ~/2;
end if
stop := II~erll ~ (c211erll +
c2) ;
end if
end if
end while
Figure 2. The complete Powell's dogleg algorithm.
2.3. Forward model
In order to solve the
inverse problem in EIT, the forward problem must be solved. When
the conductivity distribution and injected currents are known, the
forward solution gives the
estimated measured voltages on the electrodes. The most common model for the forward
34
Dogleg trust-region application in electrical impedance tomography
561
problem is the complete electrode model (Somersalo
et al 1992, Cheng
et al 1989). The
equations for the complete electrode model are:
V' .
(aVa) = 0,
au
u +z[a- =
U[,
an
f Cf
au dS =
Iz
,
an
a
au = 0,
an
XEn
(15)
x E
e[, I = 1,2, ...
,L
(16)
X E
ez, I = 1, 2, ... ,
L
(17)
(18)
where Q is the object,
u is the scalar potential distribution,
n is the outward unit normal of
the boundary
an, z[ is the electrode contact impedance,
U[ is the potential on each electrode,
II is the injected current,
L is the number of electrodes and
e[ denotes the lth electrode. For
existence and uniqueness of the solution the law of conservation of charge and a reference
voltage has to be added to the model
L
Lit =0
(19)
[=1
L
LU[=O'
(20)
[=1
In the finite element method a linear system is formed
Au =
f
(21)
where
A E
rn-(N+L)x(N+L) is the master matrix which has integrations over the elements and
over the boundary elements,
b E
rn-(N+Lxl) is the vector in which the
N first elements are
the voltages in each node and the last
L are the voltages on the electrodes, and
j(N+Lxl) is a
vector in which the first
N elements are zero and the last
L elements consist of the injected
currents. More details regarding assembling and solving the above-mentioned matrices can
be found in Polydorides (2002). The forward and inverse parts of the image reconstruction
code used in this study were implemented in MATLABTM (The Mathworks, Natick, MA). The
forward solution is an adaptation of the code provided on the EIDORS website (Polydorides
and Lionheart 2002).
3. Methods
The LM and dogleg algorithms were tested on simulated case studies. The simulated geometry
is a cylinder with both height and diameter equal to 3 (arbitrary length units). This cylinder
is surrounded by 48 electrodes in three vertically-stacked rings with each ring composed of
16 electrodes. The length and width of each electrode are 0.29 cm and 0.42 cm, respectively.
There is 0.3 cm gap between each electrode. The geometry was created and discretized
into tetragonal finite elements using FEMLABTM (COM SOL 2004). The total number of
tetrahedrons was 7149 with 1575 nodes. The geometry is displayed in figure 3 with the
electrodes indicated. Different cases were studied. In the first study, a single abnormality is
35
562
3
2.5
2
N 1.5
0.5
o
-1
Figure 3. Simulated geometry with 48 electrodes.
3
2.5
2
N
1.5
0.5
o
x
y
-1 .5
Figure 4. Tumour geometry for EIT simulation in the case of two abnormalities.
M Goharian
et at
assigned at different locations within the volume of the cylinder, while in the second case, two
abnormalities (of approximately equal volume) are located in different positions inside the
volume. The two abnormalities are depicted in figure 4. The sizes of the objects were selected
to be 1 % of the total number of finite elements with 71 finite elements belonging to each
object. The adjacent current injection pattern and adjacent voltage measurement ('adjacent'-
'adjacent') were implemented and in total 2160 pairs of measurements were obtained with
the 48 electrodes on three rings. Different Gaussian noise levels were added to the voltage
36
Dogleg trust-region application in electrical impedance tomography
563
z=1.20
z=1.0
z=0.80
z=0.60
z=0.40
z=0.20
Figure S. Original conductivity distribution for the case of one abnormality at (-0.5, -0.5, 0.5)
at different heights,
z, of cylinder for 30% changes.
measurement. The standard deviation of the noise was chosen to be 0.5 and 1 % of the mean
absol ute measured voltage level. The highest value of the added noise (1 %) is a relatively
large noise level for an EIT measurement system (Eyiiboglu
et al 1994). The regularization
parameter,
a, determines the degree of smoothing against the noise in the reconstructed
image. There are more robust approaches for choosing the value of
a, such as L-curve,
generalized cross validation (GCV) and the discrepancy principle (Hansen 1992, Johnston and
Gulrajani 2000). For instance, the L-curve method balances the semi-norm of the regularized
estimated vector against the norm of the residual vector. Graham and Adler (2006) introduced
and evaluated some objective hyperparameter selection techniques for EIT. The value of the
regularization parameter for this study was empirically chosen to be 5 x 10-5
• All simulations
in this study were performed on an AMD Opteron 175 running Linux with 4 GB of RAM
(Random Access Memory). In any iterative method, the initial guess for the conductivity
distribution has to be determined. In this study, a homogeneous estimate for the conductivity
distribution was used. The best homogeneous conductivity distribution was computed based
on Kao
et al (2006).
4. Results
4.1. Performance for single abnormality
In this study, an increased conductivity (by 30% above background) was assigned to a small
region containing 71 elements where the position of an abnormality was varied from the side
to the centre and from the bottom to the top of the cylinder. Figure 5 shows the original
conductivity distribution when the abnormality is located at (-0.5, -0.5, 0.5). The X,
Y
coordinates of the cylinder varied from [-1.5, 1.5] and height from [0, 3]. Figures 6 and 7
show the reconstructed images for the LM and dogleg methods with 1 % noise respectively.
Even in the case of 1 % noise the image quality is sufficient to detect the abnormality in both
methods. It can be noted that both methods reconstruct almost the same image since both
the LM and dogleg methods work with combinations of Gauss-Newton and steepest-descent
37
564
M Goharian
et al
z=1.20
z=1.0
z=0.80
z=O 60
z=0.40
z=0.20
Figure 6. Reconstructed images with dogleg method with I % noise for the case of one abnormality
at (-0.5, -0.5, 0.5) at different heights,
z, of cylinder.
II
z=1.20
z=1.0
z=0.80
,. IS! 1
11
,o~~ L'
Oc,
,]1
z=0.60
z=0.40
z=0.20
Figure 7. Reconstructed images with LM method with 1 % noise for the case of one abnormality
at (-0.5, -0.5, 0.5) at different heights,
z, of cylinder.
method. In other simulations the abnonnality moves from the bottom to the top of cylinder
while the X and
Y positions of the abnonnality stay in the centre of the cylinder (0, 0, 0.5)
to (0, 0, 1.4) and (0, 0, 2.5). In interest of space only reconstructed image with the dogleg
method for 1 % noise level for object at (0, 0, 0.5) is shown in figures 8 and 9. Both algorithms
perfonn worse for an object at the centre particularly for the 1 % noise level. In the adjacent
method the current is passing mostly in the outer part of imaged object. Therefore the adjacent
method is very sensitive to conductivity variation close to the boundary and insensitive to the
central conductivity variation (Dickin and Wang 1996).
38
Dogleg trust-region application in electrical impedance tomography
565
z=1.20
z=1.0
z=0.80
z=0.60
z=0.40
z=0.20
Figure 8. Original conductivity distribution for the case of one abnonnality at (0, 0, 0.5) at different
heights,
z, of cylinder for 30% change.
z=1 20
z=1.0
z=0.80
z=0.60
z=0.40
z=0.20
Figure 9. Reconstructed images with dogleg method with 0.5% noise for the case of one
abnonnality at (0, 0, 0.5) at different heights,
z, of cylinder.
The last study for a single abnormality investigated a small change in conductivity (10%
higher than background) with an abnormality at the bottom of the cylinder. The dogleg
reconstructed images with a 0.5% noise level are shown in figure 10. This figure is an
indication of the ability of the method to detect an abnormality, even with this small variation.
4.2. Performance with two abnormalities
Two cases were used in this part to check the distinguishability of the two algorithms with two
abnormalities. In the first case, there are two small regions containing 71 tetrahedron elements
39
566
M Goharian
et al
z=1.20
z=1.0
z=0.80
:r
,
'r
':
-- t
I·
..
.~
I
•
:~
• c
z=0.60
z=0.40
z=0.20
Figure 10. Reconstructed images with dogleg method with 0.5% noise for the case of one
abnormality at (-0.5, -0.5, 0.5) at different heights,
z, of cylinder with 10% conductivity change.
z=1.0
z=0.80
z=0.60
z=0.40
z=0.20
z=0.10
Figure 11. Original conductivity distribution for the case of two abnormalities both with 30%
conductivity change at (-0.5, -0.5, 0.5) and (1, 1, 0.5) at different heights,
z, of cylinder.
which are positioned diagonally at the same plane at the bottom of the cylinder (-0.5, -0.5,
0.5) and (1, 1, 0.5). In the first study the same conductivity changes of 30% higher than the
background are assigned. The original distribution for this case is shown in figure 11. The
reconstructed image with the dogleg method in the case of 1 % noise is shown in figures 12. In
the second case, the conductivity changes
(+30% and -50%) with respect to the background
are defined in the two regions at the same plane. Figures 13 and 14 show the original and
reconstructed conductivity. In both cases even with 1 % noise the abnormalities can be clearly
40
Dogleg trust-region application in electrical impedance tomography
567
z=1.0
z=0.80
z=0.60
z=0.40
z=0.20
z=0.10
Figure 12. Reconstructed conductivity images with dogleg for 1 % noise for the case of two
abnormalities both with 30% conductivity change at (-0.5, -0.5,0.5) and (1,1, 0.5) at different
heights,
z, of cylinder.
z=1.0
z=0.80
z=0.60
I'
) .
. ·:It
.1) 4
z=0.40
z=0.20
z=0.10
Figure 13. Original conductivity distribution for the case of two abnormalities at (-0.5, -0.5,
0.5) and (1, 1, 0.5) at different heights,
z, of cylinder one with 30% higher and the other one with
50% lower than background.
identified. In the last study the abnormality with 30% higher conductivity than the background
is defined in a region towards the top of the cylinder at the position (0. 5, 0.5, 1.5) and the
abnonnality with -50% is located towards the bottom of the cylinder at (-0.5, -0.5, 0.5).
The original and reconstructed images for this case are shown in figures 15 and 16.
41
568
5. Discussions
M Goharian
et at
z=1.0
z=O.80
z=O.60
.0
z=0.40
z=O.20
z=O.10
Figure 14. Reconstructed conductivity images with dogleg for the case of two abnormalities at
(-0.5, -0.5, 0.5) and (1, 1, 0.5) at different heights,
z, of cylinder at I % noise one with 30%
higher and the other one with 50% lower than background.
z=2.0
.1 .
z=1.20
z=1.80
z=1.0
~
~::,_L-
z=1.60
z=O.40
Figure 15. Original conductivity distribution for the case of two abnormalities at (-0.5, -0.5,
0.5) and (0.5, 0.5, 1.5) at different heights,
z, of cylinder one with 30% higher and the other one
with 50% lower than background.
Both algorithms proved to be robust under the tested conditions by converging towards
solutions. Even at the highest tested level of noise (1 %), the algorithms produced reconstructed
conductivity maps in which the abnonnalities were visually separable from background and
42
Dogleg trust-region application in electrical impedance tomography
569
z=2.0
z=1.80
z=1.60
.r'
z=1.20
z=1.0
z=0.40
Figure 16. Reconstructed conductivity images with dogleg for the case of two abnormalities at
(-0.5, -0.5,0.5) and (0.5, 0.5, 1.5) at different heights,
z, of cylinder with 1% noise, one with
30% higher and the other one with 50% lower than background.
Table 1. False positive and negative percentages of finite elements that were determined based on
an image intensity threshold that ensures 95% of tumour finite elements are correctly classified.
False positive False negative
DU
1.35
0.58
LMa
1.12
0.85
DLb
18.91
1.2
LMb 13.76
0.96
a Two objects at (-0.5, -0.5, 0.5) and (0.5, 0.5, 1.5) with conductivity changes of +30% and
-50% at 0.5% noise leveL
b Two objects at (-0.5, -0.5, 0.5) and (0.5, 0.5, 1.5) with conductivity changes of +30% and
-50% at 1 % noise level.
centred at the expected location. However, the presence of noise in measurements proved to
have similar negative affects on both algorithms. A visual examination of the reconstructed
images reveals that as noise increases, the number and sizes of non-existent protuberances also
increase (see figures 7 and 9). A statistical analysis of the reconstructed conductivities also
demonstrates the affects of noise on the reconstruction process. A quantitative comparison
of image reconstruction using two methods is presented in table 1. Table 1, a comparison of
reconstructions from the two algorithms with 0.5% and 1 % levels of noise, exemplifies the
results, which was typical of all tested cases. The comparison was performed by setting a
threshold on image intensity defined to include 95% of all pixels belonging to tumour tissue.
All background finite elements with image intensity higher than this level will indicate the
percentage that can be considered as false positives. The results show that as the level of
noise increases, the number of false positives increases drastically while the number of false
negatives was not significantly affected. The dogleg algorithm is slightly more sensitive to
noisy environments, where at 1 % level of noise, it has 37% more false positives than the
Levenberg-Marquardt algorithm. The level of misclassification with 1 % noise would not be
43
570
M Goharian
et at
Table 2. Perfonnance comparison of LM and DL methods: initial value of objective function,
!(a), final value of objective function
f(a), objective function and Jacobian calculations, total
number of iteration/linear system calculation, execution time in minutes.
Initial value of
f(a)
Final value of
f(a)
DL
LM
DL
LM
Case 1 3.08 x 10-4 3.08 X 10-4 4.8 x 10-6 4.8 x 10-6
Case 2 1.6 x 10-5
1.6 x 10-5
4.9 x 10-7 4.9 x 10-7
Fun/lac
DL LM
10/5 10/10
10/5 10/10
Iter/sys. solved
DL LM
10/5 10/10
10/5 lO/tO
Exec. time
DL LM
26 57
25 55
Case 1: two objects at (-0.5, -0.5,0.5) and 0, 1,0.5) with conductivity changes of +30% and -50% at 1 % noise
level.
Case 2: one object at (-0.5, -0.5, 0.5) with conductivity changes of +30% at 1 % noise level.
acceptable in most applications, however modem EIT systems report a level of lower than 1 %,
and therefore this should not prove to be a problem (McEwan
et al 2006).
The placement of the abnormalities within geometry plays a very significant role in
the ability of an algorithm to accurately reconstruct an image. Abnormalities which were
located near the surface of the geometry and closer to the electrodes were reconstructed with
much higher accuracy, greater contrast against the background and had more clearly defined
boundaries. Test cases in which the abnormality was closer to the surface were also less
sensitive to the non-existent protuberances due to noise as seen in figure 6. This was equally
true of both the dogleg and the Levenberg-Marquardt reconstructions. The abnormalities
located in the centre were also reconstructed to be much larger in size than actuality with both
algorithms. However these issues are common to all EIT reconstruction algorithms and have
to do with the very nature of EIT. Changes from abnormalities located in the centre cause a
much smaller boundary voltage change than a voltage change located near the surface (Dickin
and Wang 1996). This affect is amplified by the fact that the comparison tests used only
adjacent firing patterns. The majority of the current would not have travelled to the centre of
the cylinder and the abnormalities located in the centre would not have been probed to the
extent as an abnormality located close to the surface. The lack of clear information in voxels
located in the centre means that the regularization smoothes over conductivity changes located
in the centre of the object. The use of optimal current patterns (Lionheart
et al 2001, Goble
and Isaacson 1989), or more measurements would compensate for the lack of information
from the centre, but this could be said of any EIT reconstruction algorithm and does not go
towards demonstrating the viability of the dogleg and LM algorithm in EIT reconstruction.
A statistical comparison of the two methods is given in table 2. This comparison shows
the total number of evaluations of the objective function, the Jacobian, and solving the linear
system of equations. As can be seen, the DL algorithm requires fewer objective and linear
system of equations evaluations. It has already been mentioned that LM needs to solve
equation (3) on each iteration to determine the step. This means that when the LM step fails,
the algorithm needs to solve equation (3) with an increased damping factor, A. In other words,
every failed step results in an unproductive effort. The advantage of the dogleg method is that
it does not require finding the Gauss-Newton step in the case of failure. Once the Gauss-
Newton step has been found, the algorithms can solve the subproblem (10) for a different
radius,
L\. Reducing the number of times that the algorithm needs to solve the Gauss-Newton
step has a crucial effect on the performance of the algorithm. For many optimization problems,
computing the gradient of the objective function requires more computing time than computing
the function value, and computing the Hessian sometimes requires far more computing time
and memory than computing the gradient. As shown in table 2, the total execution time for
the dogleg method is less than 50% of that for the LM method. If the size of the mesh were to
increase, the benefits of using the dogleg algorithm would be even greater. Figure 17 shows
44
Dogleg trust-region application in electrical impedance tomography
571
.10-3 ,..--_ _--,-_ _ _- , _ _ _-,-_ _ _,.---_ _--, _ _ _--,-_ _ _- , _ _ _-,--_ _----,
\
~ LM I
-t-DL
10~~--~---~---~---~--~---~---~---~---~
1
3
6
9
10
Iteration number
Figure 17. Convergence rate change for objective function versus iteration for the case of two
objects at (-0.5, -0.5,0.5) and (1, 1,0.5) with conductivity changes of +30% and -50% at 1 %
noise level.
the convergence rate change for objective function versus iterations for the case of two objects
and is indicative of convergence rates for all test cases. As can be seen, the algorithms are able
to converge after only five iterations, with very minimal change at further iterations. If further
reduction in computation time is required, the number of iterations can be reduced with only
a small sacrifice in accuracy.
6. Conclusion
This paper has introduced the Powell's dogleg algorithm as an alternative method to
Levenberg-Marquardt for solving the EIT inverse problem. As demonstrated in all the
simulations, the dogleg method produces similar solutions to those of LM with less
computational effort. The introduction of a trust region method to the EIT problem has
additional benefits. It can be used to solve EIT problems in a constrained optimization
formulation.
Acknowledgments
The authors would like to thank NSERC for funding this project (Moran-NSERC Discovery).
References
Brenner S C and Scott L R 1994
The Mathematical Theory of Finite Element Methods (Berlin: Springer)
Cheng K S, Isaacson D, Newell J C and Gisser D G 1989 Electrode models for electric current computed tomography
IEEE Trans_ Biomed. Eng_ 36918-24
45
572
M Goharian
et al
COMSOL 2004 The FEMLAB Reference Manual COMSOL Burlington, MA, USA
Dennis J and Schnabel R 1996
Numerical methods for unconstrained optimization and nonlinear equations Classics
in Applied Mathematics (Philadelphia, PA: SIAM)
Dickin F and Wang M 1996 Electrical resistance tomography for process tomography
Meas. Sci. Technol. 7247-60
Eyiiboglu B M, Pilkington T C and Wolf P D 1994 Estimation oftissue resistivities from multiple-electrode impedance
measurements
Phys. Med. BioI. 39 1-17
Goble J and Isaacson D 1989 Optimal current patterns for three-dimensional electric current computed tomography
Proc. 11
th Annual Int. Conf. IEEE Engineering in
Medicine and Biology Society vol 2 pp 463-4 (Seattle, WA)
Graham B M and Adler A 2006 Objective selection ofhyperparameter for EIT
Physiol. Meas. 27 565-79
Groetsch C W 1993
Inverse Problems in the Mathematical Sciences (Braunschweig: Vieweg)
Hansen P C 1992 Analysis of discrete ill-posed problems by means of L-curve
SIAM Rev. 34 561-80
Jacoby S L S, Kowalik J 5 and Pizzo J T 1972
Iterative Methods for Nonlinear Optimization Problems (Englewood
Cliffs, NJ: Prentice-Hall)
Johnston P Rand Gulrajani R M 2000 Selecting the comer in the curve approach to Tikhonov regularization
IEEE
Trans. Biomed. Eng 47 1293-6
Kao T J, Kim B S, Isaacson D, Newell J C and Saulnier G J 2006 Reducing boundary effects in static EIT imaging
Physiol. Meas. 27 S13-23
Levenberg K 1944 A method for the solution of certain problems in least squares
Quart. App/. Math. 2 164-8
Lionheart W R B, Kaipio J and McLeod C N 2001 Generalized optimal current patterns and electrical safety in
EIT.
Physiol. Meas. 2285-90
Madsen K, Nielsen Hand Tingleff 0 2004
Methods for Non-Linear Least Squares Problems 2nd edn IMM,
Technical University of Denmark) available at http://www2.imm.dtu.dk/pubdb/views/edoc_download.php/
3215/pdf/imm3215.pdf
Marquardt D 1963 An algorithm for least squares estimation of nonlinear parameters
SIAM J. Appl. Math. 11431-41
McEwan A, Romsauerova A, Yerworth R, Horesh L, Bayford R and Holder D 2006 Design and calibration of a
compact multi-frequency EIT system for acute stroke imaging
Physiol. Meas. 27 S199-21O
Miller C E and Henriquez C S 1990 Finite element analysis of bioelectric phenomena
Crit. Rev. Biomed. Eng. 18
207-33
More J J and Sorensen D C 1983 Computing a trust region step
SIAM J. Sci. Stat. Comput. 4 553-72
Nielsen H B 1999 Damping parameter in Marquardt's method. IMM, DTU.
Report IMM-REP-1999-05
Nocedal J and Wright S J 1999 Numerical optimization
(Springer Series in Operations Research) (New York:
Springer)
Paulson K, Breckon Wand Pidcock M 1992 Electrode modelling in electrical impedance tomography
SIAM J. Appl.
Math. 52 1012-22
Polydorides N 2002 Image reconstruction algorithms for soft-field tomography
PhD Thesis UMIST, UK
Polydorides N and Lionheart W R B 2002 Matlab toolkit for three-dimensional electrical impedance tomography:
a contribution to the electrical impedance and diffuse optical reconstruction software project
Meas. Sci.
Technol.13 1871-83
Powell M J D 1970
A Hybrid Method for Non-Linear Equations Numerical Methods for Non-Linear Algebraic
Equations ed P Rabinowitz (London: Gordon and Breach) pp 7-114
Somersalo E, Cheney M and Isaacson D 1992 Existence and uniqueness for electrode models for electric current
computed tomography
SIAM 1.
Appl. Math. 52 1023--40
Vauhkonen P J 2004 Image reconstruction in electrical impedance tomography
PhD Thesis University of Kuopio,
Kuopio, Finland
Vauhkonen P J, Vauhkonen M, Savolainen T and Kaipio J P 1999 Three-dimensional electrical impedance tomography
based on the complete electrode model
IEEE Trans. Biomed. Eng. 461150-60
46
CHAPTER III
Paper II
The following paper is a continuation of a series of papers for new image reconstruction algorithms for
EIT. In this paper, we introduce a novel approach for the regularization of the EIT problem using a spectral
and spatial multi-frequency analysis approach. The approach is based on the construction of a subspace for
the expected conductivity distributions using principal component analysis. The advantage of this technique
is that prior information is extracted from the characteristic response of an object at different frequencies
and spatially across the finite elements.
The SS-MIA algorithm which is presented in this paper is a modification of a suggested algorithm in PhD.
thesis of the second author of the paper (Bruwer). The simulation work in this paper was performed by me
under the supervision of Drs. Moran and Macgregor. The manuscript was written by me and edited by Drs.
Moran and MacGregor.
(Reprinted from Physiol. Meas. 28 (2007) 1001-1016, Mehran Goharian, Mark-John Bruwer, Aravinthan
Jegatheesan, Gerald R Moran, and John F MacGregor, A novel approach for EIT regularization via spatial
and spectral principal component analysis , Physiol. Meas.
28 (2007) 1001-1016, Copyright (2007), with
permission from rop.
http://www.iop.orglEJ/pm).
Abstract is available at:
http://www.iop.orglEJ/abstractl0967-3334/28/9/003
47
lOP PuBLISHING
PHYSIOLOGICAL MEASUREMENT
Physiol. Meas. 28 (2007) 1001-1016
doi:10. 1 088/0967-3334/28/9/003
A novel approach for EIT regularization via spatial
and spectral principal component analysis
Mehran Goharianl
, Mark-John Bruwer2, Aravinthan Jegatheesan3,
Gerald R Moranl and John F MacGregor2
1 Medical Physics and Applied Radiation Sciences, McMaster University, Hamilton, Ontario,
Canada
2 Department of Chemical Engineering, McMaster University, Hamilton, Ontario, Canada
3 School of Biomedical Engineering, McMaster University, Hamilton, Ontario, Canada
E-mail: goharim@mcmaster.ca
Received 15 February 2007, accepted for publication 7 June 2007
Published 21 August 2007
Online at stacks.iop.org/PM/28/l001
Abstract
Electrical impedance tomography, EIT, is an imaging modality in which the
internal conductivity distribution of an object is reconstructed based on voltage
measurements on the boundary. This reconstruction problem is a nonlinear
and ill-posed inverse problem, which requires regularization to ensure a stable
solution. Most popular regularization approaches enforce smoothness in the
inverse solution. In this paper, we propose a novel approach to build a
subspace for regularization using a spectral and spatial multi-frequency analysis
approach. The approach is based on the construction of a subspace for the
expected conductivity distributions using principal component analysis. It is
shown via simulations that the reconstructed images obtained with the proposed
method are better than with the standard regularization approach. Using this
approach, the percentage of misclassified finite elements was reduced up to
twelve fold from the initial percentages after five iterations. The advantage
of this technique is that prior information is extracted from the characteristic
response of an object at different frequencies and spatially across the finite
elements.
Keywords: electrical impedance tomography, peA, multi-spectral analysis,
regularization
(Some figures in this article are in colour only in the electronic version)
0967-3334/07/091001+16$30.00 © 2007 lOP Publishing Ltd Printed in the UK
1001
48
1002
M Goharian
et of
1. Introduction
In electrical impedance tomography the internal electrical conductivity distribution is
reconstructed from electrical measurement on the surface of the object under study. Arrays
of electrodes are attached to the object surface and an alternating current is passed through at
least a single pair of electrodes. The induced voltages are measured on the remaining surface
electrodes. There are many potential applications of EIT in medical imaging (Cheney
et al
1999) and industrial process monitoring (Dickin and Wang 1996). A more comprehensive
review of the theory and application of EIT is provided in a recent text by Holder (2005).
This type of problem, in which the property under investigation is recovered from
boundary measurements, is a nonlinear inverse problem and is severely ill posed. The ill
conditioning stems from conductivity distributions for which the voltage measurements are
extremely small and influenced by noise. Such unreliable measurements make the solutions
unstable. To obtain stable and meaningful solutions regularization is used. There are two
primary regularization classes: the penalty method and the projection method. In the penalty
method a penalty term is added to the original problem to 'penalize' noise contributions to the
solution. A popular example of this class is Tikhonov regularization. The penalty methods
regularize ill-posed problems by inserting terms to filter out unwanted components in the
solutions. These unwanted properties in the solutions originate from the noisy measurements.
Tikhonov regularization achieves this by means of adding a penalty term to the objective
function.
In the projection method the solution is projected to a subspace to remove the unwanted
components (Jacobsen 2004). The truncated singular value decomposition (TSVD) is an
example. Another class of regularization is the so-called hybrid method which is combination
of projection and penalty methods (Hanke and Hansen 1993). The implementation of
regularization is equivalent to introducing
a priori information in the reconstruction process,
which may affect the ability of detecting sharp variations in the conductivity distribution (in
EIT for example). The Tikhonov regularized form for the nonlinear problem of EIT can be
written as
min
<I>(a) =
{IIV -
F(a)1l2 +a \lLaI12
}
(1)
(J
where
a is the conductivity distribution,
V are the measured voltages,
F (a) is a functional map
between the conductivity distribution and the voltage measurements,
L is the regularization
matrix, and
a is a regularization factor. The regularization factor tries to balance the amount of
prior information that is used in the final solution. The most common approaches to determine
a are the L-curve (Lawson and Hanson 1974), generalized cross validation (Wahba 1977) and
discrepancy principle methods (Morozov 1984). Equation (1) can be solved recursively via
Newton's method as follows:
(2)
/).an =
-H-1\1<1>
= (JT J +
aLT L)-l[JT (V -
F(an) -
aLT Lan]
(3)
where
H is the approximation for the Hessian of the objective function, and
J is the Jacobian
for the mapping function
F(a). Differential operators (typically, first and second order) are
the most common penalty functions that have been used in Tikhonov regularization (Hua
et al
1988, Woo
et al 1993).
In the NOSER algorithm (Cheney
et al 1990) a diagonal matrix was suggested for
regularization which chooses
L =
J JT J. Smooth solutions are expected when using these
49
A novel approach for EIT regularization via spatial and spectral principal component analysis
1003
regularization matrices. Vauhkonen
et al (1997) introduced the basis constraints method
(BCM) in which the conductivity distribution was considered to be a linear combination of
preselected basis functions; Wi, that is,
a = L~l Ci
Wi, Ci E m are the weight functions
for the basis functions. A representative ensemble of possible conductivity distributions was
created based on prior information on the structures and conductivities of the object. The
basis functions are created from the principal eigenvectors of the ensemble. The drawback of
the BCM method is that it is biased towards a particular set of solutions. In another study,
Vauhkonen
et at 1998 introduced the subspace regularization method (SSRM) to overcome
the prior bias problem of BCM. The SSRM method weights the solution towards the subspace
which is spanned by a representative ensemble of orthonormal vectors. Vauhkonen suggests
calculating a regularization matrix whose null space is Sw in the following way:
L =
I -
WWT
(4)
where
I is the identity matrix and
W is the matrix holding all orthonormal column vectors,
Wi, i = 1, ... ,
M. If
a lies in the subspace
Sw (spanned by the orthonormal column vectors,
wd then
WWT a =
a (since
WWT gives the orthogonal projection of
a onto the subspace
Sw). Thus, the implementation of
L gives,
La =
(l -
WWT)a =
a -
WWT a = 0 (if
a
lies in
Sw) and such
a are favored. In Vauhkonen
et al (1998) it was shown that the SSRM
technique is able to reconstruct conductivities which do not match the prior information. This
is due to the fact that SSRM only draws the solution towards the subspace instead of forcing
it, unlike BCM. Vauhkonen (1997) built a subspace from an ensemble of 81 possible thoracic
conductivity distributions, using the anatomical structure of the tissue being imaged from MRI
images, and published conductivity values.
In this paper we proposed subspace regularization using a similar approach to SSRM.
However, in our work it is assumed that the anatomical structure and conductivities are not
available
a priori to build the subspace. We propose a novel approach to build the subspace
for regularization using spectral and spatial multi-variate analysis techniques developed by
Bruwer (2006). The paper is organized as follows. Section 2 develops the principles and
implementation of the proposed regularization. Section 3 introduces the simulation-based
studies used to demonstrate the feasibility of this technique. Section 4 presents and discusses
the results.
2. Theory and algorithms
2.1. Multi-frequency peA modeling of the EIT problem
The electrical conductivity of mammalian tissue depends on the tissue structure at the cellular
level and it is known to change with frequency (Schwan 1957). The conductivity of excised
tissue samples can be modeled via the Cole-Cole model (Cole and Cole 1941). The Cole-Cole
model has been directly implemented in the EIT image reconstruction algorithms by using
data across a range of frequencies (Mayer
et al 2006). Constructing the model in this way
is often referred to as parametric imaging. The premise of parametric imaging is that the
collection of data simultaneously across a range of frequencies will decrease the ill posedness
of the problem.
A similar multiple-frequency approach has been applied in microwave image
reconstruction (Fang
et al 2004). Fang
et al implemented simultaneous multiple-frequency
measurement data to reconstruct a single image of tissue property dispersion. The high-
contrast objects were successfully reconstructed through the multiple frequency-dispersion
reconstruction algorithm (MFDR) for which a single-frequency algorithm failed to reconstruct
50
1004
M Goharian
et al
an image. Image quality was improved by adding more frequencies in the reconstruction
model. In another approach a multispectral constrained approach has been implemented
for near-infrared frequency domain tomography for the purpose of creating images of total
hemoglobin, water, oxygen saturation and scatter parameters (Srinivasan
et a/200S). In this
method the spectral constraint of the chromophore and scattering spectra was implemented
directly in the algorithm, thus reducing the parameter space from 12 images to 5 parametric
images. In addition, the spectral constraint approach reduced the noise in the reconstructed
chromophore concentration, especially those in water and scatter power images. The
usefulness of the multi-frequency approach has also been investigated in the bioluminescence
optical tomography (BLT) technique. For bioluminescence imaging, it is important to be able
to correctly localize and quantify the three-dimensional (3D) distribution of bioluminescent
sources inside a small animal to study various molecular and cellular activities. Dehghani
et al (2006) used multi wavelength emission data in spectrally resolved bioluminescence optical
tomography to improve the accuracy of the reconstructed image. The bioluminescent source
location in the reconstructed image was detected to be within 1 mm of the original source
location by using six wavelength bands.
In a related work (Bruwer 2006) it was shown how principal component analysis (peA)
can be used as an alternative, non-parametric approach to account for the frequency dependence
of tissue conductivity in EIT image reconstruction. peA perfonns a mathematical procedure
that transfonns a number of correlated variables into a smaller number of un correlated variables
called
principal components. The first principal component contains as much of the variability
in the data as possible, and each succeeding component accounts for as much of the remaining
variability as possible. Low-order components often contain the most important aspects of the
data.
peA implementation on multi-frequency BIT data (in this study conductivity distribution)
can effectively compress highly correlated data and project it onto a reduced subspace via a
linear combination of the original conductivity data. In this way only a few reduced dimensions
are needed to describe all of the significant variation in the multi-frequency EIT data. In order
to implement multi-frequency analysis via peA on the EIT data, the original image data must
be fonnatted into a matrix. The following procedure is perfonned. Each finite element in the
image has conductivity values over the range of discrete frequencies. Each such characteristic
spectrum is considered as a row vector and these row vectors concatenated vertically to form
a matrix, X, of dimension
(M x
N), where
M is equal to the number of finite elements and
N
is the number of discrete frequencies.
To find the first principal component from the original conductivity data, X, the following
optimization problem should be solved:
max '11 = pi
XT XPI
PI
S.t. pi PI = 1
(5)
The solution to (5) gives PI, the eigenvector with the largest eigenvalue (AI) of the covariance
matrix,
XTX, scaled to unit length. To find all principal component vectors the additional
constraint should be added to the optimization problem in (5) which states that they are
orthogonal to each other:
max'l1i == pi
XT XPi
i = 1, ... ,
k
PI
S.t. pi Pi = 1
(6)
pi Pj = 0
j = 1, ... , (i - 1),
y
i ~ 2
51
A novel approach for EIT regularization via spatial and spectral principal component analysis
1005
The principal component vectors,
PI, P2,' .. ,
Pk, are concatenated into the loading matrix,
P. For a given observation, Xj, the projection onto the peA model gives a set of 'scores', ti,
through
tT =
xi P. The structure in the original conductivity data matrix, X, is modeled as
A
X =
Ltipr +E =
TpT +E
i=l
(7)
where
A is equal to or less than the rank of X,
T is the score matrix
(M x
A), pT is the
loading matrix
(A x
N), and
E is the residual matrix
(M x
N). See Wold
et al (1987) for more
details on the peA method. In this way, the frequency-dependent structure in the original
conductivity data, which excludes most of the unstructured noise, can be calculated using only
the
A dominant principal components via
-
T
X=TP =X-E.
(8)
Each finite element has its own unique conductivity spectrum with respect to all the frequencies,
referred to here as its
spectral signature. This method of image analysis is called 'multivariate
image analysis' (MIA) (Ge1adi and Grahn 1996).
2.2. Simultaneously extracting spatial and spectral structure
Tumor tissues have high vasculature leading to higher water content when compared to
n ann al tissue. The high rates of necrosis in tumor tissue also lead to breakdown of cell
membranes, lowering the resistance to low frequency current transmitted through tissue.
These physiological changes manifest as changes in electrical properties. It has been shown
that the dielectric constant and conductivity for cancerous tissue (for instance breast tissue)
is three or more times greater than that of nonnal tissue (Surowiec
et al 1988). Therefore
it would be possible to identify and cluster finite elements in a 3D EIT reconstruction based
on the detected electrical properties. If the EIT mesh was spatially finer than the underlying
physiological traits, it would be expected that the electrical properties of neighboring finite
elements would exhibit similar electrical properties. In Bruwer (2006) a new method which
combines the spatial correlation structure with the frequency correlation has been developed.
This method is called 'spatial-spectral MIA', or 'SS-MIA'.
The proposed SS-MIA approach is as follows. EIT data from multiple frequencies are
gathered and an image is reconstructed independently for each frequency. The reconstructed
conductivities are then rearranged into a matrix, X, where the number of rows is given by the
number of finite elements in the 3D mesh and the number of columns is given by the number
of frequencies for which data were reconstructed. The
ith row in this matrix would represent
conductivities across the range of frequencies for the
ith finite element of the mesh. peA is
applied to the matrix of multi-frequency EIT data, X, to get dominant score vectors (matrix
T
in equation (7», where the
ith row vector
ofT, tT is the compressed spectral signature ofthe
ith
finite element. For the
ith finite element, the finite elements,j = 1, ...
,ki' are selected which
have their centers of mass within a specific radius,
r, from the center of mass of the
ith finite
element. The bounding volume is a sphere in three dimensions and a circle in two dimensions.
The size of radius acts as a balance between image noise reduction and image smoothing.
Thereafter, the vector distance between each spectral signature,
tj, of the
jth selected finite
element and that of the
ith finite element
(ti) is calculated as follows:
dij =
L (tjl -
til)2,
l
j = 1, ...
,kj,
52
(9)
1006
M Goharian
et at
where
I = 1, ... ,
A, A being the number of principal components in the peA model of
equation (7). These distances are concatenated into the vector, dj =
[dn ...
djk]T. This
analysis is repeated for all finite elements in the model, giving a set of these spectral distance
vectors,
D = {di ,
i = 1, ... ,
M}, where M is the total number of finite elements. The set
D has
two important features. First, the lengths
(kj ) of the vectors, dj , vary due to the variable number
of finite elements within a radius,
r, for each given finite element. The second feature of the
set,
D, is that the order of the elements in each vector, dj , varies based on mesh geometry. The
following approach is implemented to deal with the variable vector lengths. First, the greatest
vector length,
K = max{kj,
i = 1, ... , M}, in the set,
D, is determined. Then, for each vector,
dj ,
aj =
floor(K/kj), and hj =
K -
kjaj are calculated. If hj = 0
(ki divides exactly into K) then
ddj is formed simply by concatenating dj with
aj - 1 copies of itself. If
hi > 0 then the first
(kjai) elements of ddi are filled as described above, while a different approach is needed to fill
the last
hi elements. A matrix,
Si. is formed in which the first column is the vector, Sj, and the
second column the vector di. The elements of Sj are the spatial distances from the surrounding
finite elements to the centre finite element (between centers of mass). Next, the rows of Si are
sorted so that the elements of the first column are in ascending order. Finally,
hi equally spaced
elements from the second column of Si are selected and used to fill the last
hi elements of ddi.
This approach has the effect of finding a subset of finite elements in a manner which does not
introduce bias, but also samples the range of radial distances (in the spatial image-domain)
evenly. The above operation ensures that the new set,
DD = {ddi,
i = 1, ... ,
M}, has vectors
of equal length,
K. A key feature of this approach is that variation of finite element density
throughout the whole volume is handled automatically. The elements of the augmented vector
ddi are sorted in ascending order to remove the random arrangement. This procedure creates a
'spectral-distance spectrum' for each finite element. Finite elements whose spectral signature
is consistent with those of its immediate spatial neighborhood will have a spectral-distance
spectrum comprised of mostly small values. By contrast, if the spectral signature of a given
finite element is very different from its spatial neighborhood, its spectral-distance spectrum
will show many large values.
The last step is to perform peA analysis on the combination of the compressed spectral
signature,
ti, for each finite element and its associated spectral-distance signature, ddj • The
peA is performed on an augmented matrix, Xaug =
[T Q], where
T is the matrix of scores
from the original peA model (cf equation (7)), and
Q is a matrix formed by vertically
concatenating the row vectors, dd? A block scaling on
T and
Q is applied so that they have
equal weight in the peA model (Eriksson
et al 20a 1). Figure 1 shows a flowchart for the
SS-MIA algorithm.
2.3. Image reconstruction procedure with SS-MIA
The Newton-type method (2) was used for the solution of the nonlinear problem (1). In
this study a three-step approach was implemented for image reconstruction. The proposed
procedure is as follows:
(1) Assuming a uniform conductivity, take the first step in the iterative Newton method for
the inverse problem using equations (2) and (3) simultaneously to reconstruct the first
estimate of the conductivity distribution for each discrete frequency.
(2) Generate the matrix, X, by concatenating the individual conductivity vectors, and perform
peA analysis to estimate X =
T pT.
(3) Each column of X gives the updated conductivity estimated at the corresponding
frequency. These estimated conductivity distributions are used for the next iteration.
53
A novel approach for EIT regularization via spatial and spectral principal component analysis
1007
Reconstruct EIT conductivity
at multiple frequencies and
arrange results into column
vectors
..
Concatenate conductivity
column vectors into matrix X
.!.
Perlorm PCA on matrix X to
find matrix T
L
Repeat for each finite element i
.....
Find all finite elements whose
centre-of-mass is within
radius rfrom the centre-of-
mass of element
i
+
For each element within
radius, calculate spectral
distance using (7) and
arrange into vector matrix d,
+
Calculate k, the vector
lengths for all d,
I ..
Find the greatest vector of
greatest length K
L
Repeat for each finite element i
y
Find the number of times the
vector can be completely
filled with copies of Itself
using aj=ftoor(K/kl}
~
Find the difference in number
of elements from the
maximum number of
elements using b,=K-k,·a,
ifb,=O
I
If b, > 0
y
•
~ Concatenate d, with a,-1
I
copies of itself
Create a copy of d, and
arranged according to spatial
distance from element I.
Select b, number of elements
equally spaced apart and add
spectral values to d,
+
Concatenate d, with a,-1
copies of Itself
I
•
Sort d, according to spatial
distance from finite element
i
I
".
Create matrix 0 by vertical
concatenatIOn all d"
i=1 .. M as
row vectors
~
Create Xaug=[T OJ
Figure 1. Flowchart for SS-MIA algorithm.
54
1008
2.5
N 1.5
1
.
0.5
o
1.5
Figure 2. Simulated EIT geometry with electrodes.
M Goharian
et af
1.5
(4) Generate Xaug as described in section 2.2 and perform PCA analysis to generate X aug (the
SS-MIA approach). The first A columns of this matrix correspond to the updated score
matrix,
T. Equation (8) is then implemented to obtain the updated X =
T pT.
(5) The regularization matrix,
L, can be updated using equation (4). In this case the
a priori
subspace,
W, for a given frequency would simply be the corresponding column of the
updated X from SS-MIA, normalized to unit length.
(6) Petform the next step of Newton's method (equations (2) and (3)) with the updated value
of
L. This is done for each frequency.
(7) Return to step 2 and repeat until convergence or the maximum iteration number is reached.
The EIT image reconstruction code used in this study was implemented in MatLab™.
The forward solver is an adaptation of the code provided on the EIDORS website (Polydorides
and Lionheart 2002).
3.
Methods
The proposed SS-MIA regularization technique was tested on simulated cases. The simulated
geometry is a cylinder with both height and diameter equal to 3 cm. There are 48 electrodes
in three rings with each ring composed of 16 electrodes connected peripherally around the
cylinder. The geometry was created and discretized into tetragonal finite elements using
FEMLABTM. The total number of tetrahedral for the inverse part was 7149 with 1575 nodes.
The number of tetrahedrons used for forward computations was 8125, in order to avoid what
is referred as an
inverse crime (Lionheart 2004). The geometry is displayed in figure 2 with
the electrodes indicated. Two tumors are positioned at (0.5, 0.5, 1.5) and (-0.5, -0.5, 0.5)
respectively within the tissue volume. These are depicted in figure 3. The tumor volumes were
calculated to be 0.24 and 0.25 cm3 for upper and lower tumors respectively, and the number
of finite elements belonging to each tumor was approximately 1 % of the total number of finite
elements.
55
A novel approach for EIT regularization via spatial and spectral principal component analysis
Figure 3. Tumor geometry for EIT simulation.
3 1----------------------------~
+ Normal tissue
~ 2.51 -t- Tumor tis~ue
E ()
en
oS
f 1.51
15
-g 11
o
<.) 0.5~~--&e--~e-~e .... e~-~e~-~e~e~et-<t
10.1
Frequency, [MHz)
1.5
Figure 4. Breast tissue conductivity variation verses frequency (Surowiec
el al 1988).
1009
The opposite current injection pattern and adjacent voltage measurement ('opposite'-
'adjacent') were implemented. Different Gaussian noise levels were added to the measured
voltages. The standard deviation of the noise was chosen to be either 1 % or 3% of the mean
absolute measured voltage.
The conductivities of the tumor and the background as a function of frequency were chosen
to be the same as the experimental data presented by Surowiec
et al (1988) for cancerous and
normal breast tissue, respectively. Ten discrete frequencies were chosen from 10 KHz to
1 MHz in this study. The conductivity frequency dependence is shown in figure 4. To
reduce the size of the simulation and thus reduce the computational task only the real part of
conductivity was simulated.
The regularization parameter,
ex, determines the degree of smoothing against the noise
in the reconstructed image. More robust approaches, such as L-curve, generalized cross
validation (GCV), and the discrepancy principle (Hansen 1992, Johnston and Gulrajani 2000),
56
1010
04
0.35
0.3
$
.g,0 25
.~
Ol 0.2
c
:c
~ 0.15
..J
0.1
0.05
2
4 6 8
Discrete frequency
~o
M Goharian
et al
1
Figure 5. Loading weights for first principal component of the MIA peA model.
can be used to detennine the optimal value for the regularization factor. However, in this
study, the hyperparameter,
a, was chosen heuristically to be 5 x 10-5
.
There are other parameters for the MIA and SS-MIA approaches that need to be set
for the simulation study. Since the conductivity variables all have the same units, the data
X for the PCA model was not scaled. In implementing the SS-MIA method to update the
regularization matrix on each iteration the data in Xaug were mean centered and block scaled.
This is performed by multiplying each variable in a given block by the additional scaling
factor
1/,Jk, where
k is the number of variables in that block (Eriksson
et al 2001). The other
parameter to be chosen was the radius,
r, of the bounding volume. The value was chosen
empirically to balance image noise reduction against image smoothing. A value of
r = 0.8
was found empirically to be reasonable.
A quantitative comparison between the three methods was performed. In this case a
threshold on image intensity was defined to include 95% of all pixels belonging to tumor
tissue. All background finite elements with image intensity higher than this level will indicate
the percentage that can be considered to be false positives.
4. Results
4.1. Image reconstruction with MIA and SS-MIA peA model
The proposed image reconstruction procedure in section 2.3 was applied to the simulated
geometry. Cross validation (Eriksson
et al 2001) was used to find the most significant
principal components for MIA and SS-MIA. The loading vector corresponding to the most
significant component for MIA is plotted in figure 5. It is obvious from the loading
vector graph that this component captures the variation in conductivity over all frequencies
(cf figure 4). Figure 6 plots the loading weights for SS-MIA. The loading weights are related
to the spectral-distance spectrum (ddj). The spectral signature consists of just one variable (tii
from the MIA model). Its associated loading weight is 0.78 (the first element in figure 6-off
the scale). Comparing the first element with the loading weights for the other elements in
figure 6 indicates a strong correlation between the spectral signature of a given finite
element and the distances in the spectral domain between it and the finite elements in its
spatial neighborhood. The physical interpretation of this is that finite elements with a high
conductivity have a very different spectral signature from their neighbouring finite elements.
57
A novel approach for EIT regularization via spatial and spectral principal component analysis
0.04 -
0.035
0.03
!II
~O. 025
.~
C) 0.02
t:
'6
~ 0.015
...J
0.01
0.005
o
--'-
100
200
300
400
5cio 600
70c)
Spectral-distance spectrum
1011
Figure 6. Loading weights for first principal component of the SS-MIA peA model. The model
shows the correlation structure between spectral signature and its relationship with its spatial
neighbors. The spectral signature comprises just one variable. Its associated loading weight is
0.78 (the first element-offthe scale).
z=2.5
z=1 .70
z=1 .50
z=1 .30
z=1 .00
z=O.50
Figure 7. Original conductivity distributions for the geometry illustrated in figure 2 shown as a
series of 2D slices (at lowest frequency 10 kHz).
This is sensible since tumors have high conductivity and also tend to have a smaller volume
relative to the total tissue volume.
Figure 7 shows the original conductivity distribution as series of 20 slices at different
levels for the geometry with two abnormalities in figure 3. Figures 8 and 9 show the
reconstructed images for 1 % and 3% noise level, respectively, using both MIA and SS-
MIA to update the conductivity and regularization matrices, respectively.
4.2. Image reconstruction with only SS-MIA peA model
In this case the inverse problem was solved using the Gauss-Newton (equations (2)-(3))
method with the regularization matrix was updated using the SS-MIA regularization
procedures described in the previous sections. MIA was not used to update the conductivity
matrix. Figure 10 shows the reconstructed images for 3% noise.
58
1012
M Goharian
et al
z=2.5
z=1 .70
z=1 .50
z=1 .30
z=1 .00
z=0.50
Figure 8. Reconstructed conductivity distributions for the case of 1 % noise level with both MIA
and SS-MIA shown as a series of 2D slices. Images reconstructed after five iterations (at lowest
frequency 10 kHz).
z=2.5
z=1.70
z=1.50
z=1.30
z=1 .00
z=0.50
Figure 9. Reconstructed conductivity distributions for the case of 3% noise level with both MIA
and SS-MIA shown as a series of 2D slices. Images reconstructed after five iterations (at lowest
frequency 10kHz).
4.3. Image reconstruction with total variation method
For the comparison with the MIA and SS-MIA methods the total vanatlOn (TV) as a
regularization approach was used. Total variation assumes that the data set is blocky and
discontinuous (Borsic 2002). Figures 11 and 12 show the reconstructed images for 1 % and
3% noise level using TV regularization.
5. Discussion
Figure 13 shows a quantitative comparison for the case of 3% added noise level. This shows
that the percentage of misclassified finite elements is reduced up to twelve fold from the
initial percentage after five iterations. SS-MIA has been presented as a new approach for
59
A novel approach for EIT regularization via spatial and spectral principal component analysis
1013
z=2.5
z=1 .70
z=1 .50
z=1 .30
z=1 .00
z=O.50
o
•. :
Figure 10. Reconstructed conductivity distributions for the case of 3% noise level witj1 only
SS-MIA shown as a series of 2D slices. In this case the regularization matrix is updated in each
iteration based on the SS-MIA PCA model. Images reconstructed after five iterations (at lowest
frequency 10 kHz).
z=2.5
z=1 .70
z=1 .50
z=1 .30
z=1 .00
z=O.50
Figure 11. Reconstructed conductivity distributions for the case of 1 % noise level with TV
regularization shown as a series of 2D slices. Images reconstructed after five iterations (at lowest
frequency 10kHz).
regularization through multi-frequency EIT imaging. It is clear, both from comparing the
images in figures 8 and 9 with figures 11 and 12, and from the metrics in figure 13, that the
combined SS-MIA and MIA approach improves the reconstructed images over the traditional
approach of a regularized image reconstruction at a single frequency. The benefit of this
technique is that prior information is extracted from both the frequency and spatial behavior
of the object under study. As can be seen in figure 13, the false positive percentage for TV
regularization at the fifth iteration is higher than the third iteration which is an indication that
the algorithm did not converge.
Figure 10 indicates that the implementing solely the SS-MIA approach produced a more
uniform image (compared to the combination of MIA and SS-MIA) in the non-tumor regions,
60
1014
M Goharian
et al
z=2.5
z=1 .70
z=1 .50
z=1 .30
z=1 .00
z=O.50
Figure 12. Reconstructed conductivity distributions for the case of 3% noise level with TV
regularization shown as a series of 2D slices. Images reconstructed after five iterations (at lowest
frequency 10 kHz).
80
c 70
o
~60
~
~SO
U III
'E40
'0
~30
~
~ 20
Q)
0..'0
o
72.4
12.B
rn MIA & SS-MIA
nnSS-MIA
eHV
47.15
iteration number
Figure 13. False positive percentage of finite elements that was determined based on an image
intensity threshold that ensures 95% of tumor finite elements are correctly classified.
but also smoothes the image over the tumors as well. The main characteristic of the SS-MIA
approach is to extract the key spatial correlation from the given geometry. This correlation
structure is shown by equation (9). Finite elements which show a large conductivity value
(i.e., belonging to tumors) seem to be surrounded by many more finite elements that have a
different conductivity. This is due to the fact that the tumor size is small, so that moving far
enough in any direction will intersect normal tissue with a lower conductivity. In contrast, the
normal tissue with low conductivity tends to be surrounded by many finite elements with the
same conductivity.
6. Conclusions
The traditional EIT regularization approach at a single frequency needs to include
a priori
information which could be anatomical structure from another modality or the noise level
61
A novel approach for EIT regularization via spatial and spectral principal component analysis
1015
of the measurement system. Combined MIA and SS-MIA have been presented as an
alternative approach. The advantage of this approach is that the principal components and then
regularization matrix are determined from the statistical nature of the multi-frequency data, a
property that is called the' spectral-distance spectrum', rather than any
a priori assumptions
about the model structure. Since PCA can transform input data from the high-dimensional
space into a low-dimensional subspace while retaining most of the intrinsic information of the
input data, it is an effective tool to extract information from the multi-frequency data.
Acknowledgment
The authors would like to thank the reviewers of this paper for their valuable critical comments.
References
Borsic A 2002 Regularisation methods for imaging from electrical measurements
PhD Thesis, Oxford Brookes
University, England,UK
Bruwer M J 2006 Process systems approaches to diagnostic imaging and identification
PhD Thesis, McMaster
University, Canada
Cheney M, Isaacson D and Newell J C 1999 Electrical impedance tomography
SIAM Rev. 41 85-101
Cheney M, Isaacson D, Newell J C, Simske S and Goble J 1990 NOSER: an algorithm for solving the inverse
conductivity problem
Int. 1.
Imag. Syst. Techno!. 2 66-75
Cole K S and Cole R H 1941 Dispersion and absorption in dielectrics: I. Alternating current characteristics
1. Chem.
Phys. 9341-51
Dehghani H, Davis S C, Jiang S, Pogue B W, Paulsen K D and Patterson M S 2006 Spectrally resolved bioluminescence
optical tomography
Opt. Lett. 31 365-7
Dickin F and Wang M 1996 Electrical resistance tomography for process tomography
Meas. Sci. Technol. 7 247--60
Eriksson L, Johansson E, Kettaneh-Wold N and Wold S 2001
Multi-and Megavariate Data Analysis: Principles and
Applications (Tvistevagen UmeA, Sweden: Umetrics Academy)
Fang Q, Meaney P M and Paulsen D P 2004 Microwave image reconstruction of tissue property dispersion
characteristics utilizing multiple frequency information
IEEE Trans. Microw. Theory Tech. 52 1866-75
Geladi P and Grahn H 1996
Multivariate Image Analysis (Chichester, UK: Wiley)
Hanke M and Hansen P C 1993 Regularization methods for large-scale problems
Surv. Math. Ind. 3 253-315
Hansen P C 1992 Analysis of discrete ill-posed problems by means of L-curve
SIAM Rev. 34 561-80
Holder D S 2005
Electrical Impedance Tomography: Methods, History and Application (Bristol, UK: Institute of
Physics Publishing)
Hua P, Webster J G and Tompkins W J 1988 A regularised electrical impedance tomography reconstruction algorithm
Clin. Phys. Physiol. Meas. (Suppl. A) 9137-41
Jacobsen M 2004 Modular regularization algorithms
PhD Thesis Technical University of Denmark, DTU
Johnston P Rand Gulrajani R M 2000 Selecting the comer in the L-curve approach to Tikhonov regularization
IEEE
Trans. Biomed. Eng. 47 1293-6
Lawson C L and Hanson R J 1974
Solving Least Square Problems (Englewood Cliffs, NJ: Prentice-Hall)
Lionheart W R B 2004 reconstruction algorithms: pitfalls, challenges and recent development
Physiol. Meas.
25125-42
Mayer M, Brunner P, Merwa Smolle-Juttner R, Maier A and Scharfetter H 2006 Direct reconstruction of tissue
parameters from differential multi frequency EIT
in vivo Physiol. Meas. 27 S93-S101
Morozov V A 1984
Methods for Solving Incorrectly Posed Problems (New York: Springer)
Polydorides N and Lionheart W R B 2002 A matlab toolkit for three-dimensional electrical impedance tomography:
a contribution to the electrical impedance and diffuse optical reconstruction software project
Meas. Sci.
Technol
13 1871-83
Schwan H P 1957 Electrical properties of tissue and cell suspensions
Adv. Bioi. Med. Phys. 5147-208
Srinivasan S, Pogue B W, Jiang S, Dehghani H and Paulsen K D 2005 Spectrally constrained chromophore and
scattering NIR tomography provides quantitative and robust reconstruction
Appl. Opt. 44 1858-69
Surowiec A J, Stuchly S S, Barr J R and Swamp A 1988 Dielectric properties of breast carcinoma and the surrounding
tissues
IEEE Trans. Biomed. Eng. 35 257-63
Vauhkonen M 1997 Electrical impedance tomography and prior information
PhD Thesis Kuopio University, Finland
62
1016
M Goharian
et al
Vauhkonen M, Kaipio J P, Somersalo E and Karjalainen P A 1997 Electrical impedance tomography with basis
constraints
Inverse Problems 13 523-30
Vauhkonen M, Vadasz D, Kaipio J P, Somersalo E and Karjalainen P A 1998 Tikhonov regularization and prior
information in electrical impedance tomography
IEEE Trans. Med. Imaging 17285-93
Wahba G 1977 Ill-posed problems: numerical and statistical methods for mildly, moderately and severely ill-posed
problems with noisy data
SIAM J. Numer. Anal. 14 651-67
Wold S, Eshensen K and Geladi P 1987 Principal component analysis
Chemometr. Intel!. Lab. Syst. 2 37-52
Woo E J, Hua P, Webster J G and Tompkins W J 1993 A robust image reconstruction algorithm and its parallel
implementation in electrical impedance tomography
IEEE Trans. Med. Imaging 12 137-46
63
CHAPTER IV
Paper III
The following paper presents a new approach for the regularization of the ill-posed EIT problem through a
quadratic constrained least square approach. The proposed approach is based on the trust region
subproblem (TRS), which uses L-curve maximum curvature criteria to find a regularization parameter. In
this technique the trust region radius changes during each iteration which forced the algorithm toward the
correct radius that corresponds to the elbow, the point of maximum curvature on the L-curve.
The simulation and experimental work presented in this paper was performed by me under the supervision
of Dr. Moran. The manuscript was written by me and edited by Drs. Moran and Soleimani.
This manuscript has been accepted for publication in the proceedings of the 13th International Conference
on Electrical Bioimpedance combined with the 8th Conference on Electrical Impedance Tomography,
Vienna, Austria.
64
Regularization of EIT Problem Using Trust Region SubProblem
Method
M.Goharianl
, M. Soleimane, A. Jegatheesan3, and G.R. Moran 1.
4
IMedical Physics and Applied Radiation Sciences, McMaster University, Hamilton, Canada
2William Lee Innovation Centre, School of Material, Manchester University, England
3School of Biomedical Engineering, McMaster University, Hamilton, Canada
4Diagnostic Imaging, Hamilton Health Sciences, Hamilton, Ontario, Canada
ABSTRCT
The electrical impedance tomography CElT) problem is to image electrical proper-
ties, such as conductivity and permittivity, in the interior of a body given mea-
surements of voltages at the boundary of the object. In this paper we present a new
approach for regularization of the ill-posed EIT problem. This approach is based
on the trust region subproblem (TRS), which uses L-curve maximum curvature cri-
teria to find a regularization parameter. A comparison of the TRS method with
conjugate gradient least squares (CGLS) for an experimental phantom study is pre-
sented. CGLS is an efficient technique when the norm of measured noise is exactly
known. Our results show that both methods converge to the same point on the L-
curve when the noise level is known. The TRS algorithm has the advantage that it
does not require any knowledge of the norm of the noise.
Keywords: Regularization, EIT, Image reconstruction, Trust region subproblem, CGLS.
1.
INTRODUCTION
In electrical impedance tomography (EIT) the internal electrical conductivity can be reconstructed from
voltage measurements on the surface of an object under study. In the case of complete and noiseless
boundary measurements, the EIT problem is known to have a unique solution. In practice however, the
measured data is noisy and incomplete. Hence, in this situation it is difficult to obtain a satisfactory
solution from the nonlinear and ill-posed EIT problem. So it is necessary to use regularization tech-
niques because of the ill-posed nature of the problem.
Each regularization method employs a special parameter known as the regularization parameter, to
control the effect of the noise on the solution. The nature of the regularization parameter is different for
each method. For instance, in Tikhonov regularization the penalty parameter acts as a regularization
parameters and in iterative approaches the number of iterations will serve as the regularization parame-
ter. In the EIT problem, the relation between perturbations in the internal conductivity distribution and
65
the perturbations in the boundary measurements is nonlinear. However, this nonlinear relation can be
formulated using a linearized form. If in the continuous region,
1 is the Frechet derivative of the poten-
tial, U, with respect to conductivity, cr, for the nonlinear mapping
CF H
U(CF) , then the linearized form
has following form
Jx=8U
(1)
where x E
iRn is the perturbation in conductivity distribution that maps to the differential measure-
ments
OU E
9lm
, and
J E 91
mxn is the lacobian matrix. Due to the ill-posed nature of problem, x cannot
be recovered from (1). In practice, the most popular technique is to apply the conjugate gradient me-
thod to the normal equations associated with problem (1). In this way a least square conjugate gradient
(CGLS) approach to solve the inverse problem can be implemented in a similar fashion to a least
square problem
minllJx - OUl12
(2)
x
CGLS creates components in the direction of singular vectors related to large singular values at the
early stages of iterations. While components associated with small singular values will be effective at
later iterations this means that the number of iterations acts as a regularization parameter. Thus, the
regularized solution would be obtained by stopping the iteration before the unwanted components add-
ed to the current solution.
The success of CGLS strongly depends on knowing when to stop the iteration, which is a difficult task
in itself. On the other hand, we could use CGLS with the Tikhonov regularization approach, which is
called a damped least squares problem. The convergence of this technique also depends on both a good
choice of the damping parameter and a preconditioner.
In this paper, we introduce the regularization problem as a quadratic constrained least squares problem.
It has been shown that this approach is equivalent to Tikhonov regularization (Elden, 1977). In the area
of optimization this problem is a special case of trust-region methods, which is known as the trust-
region subproblem (Bjorck,1996). A comprehensive study regarding different approaches for solving
the regularization problem as a trust-region subproblem is given in (Rojas, 1998).
We adopted the recently developed method for solving the trust-region subproblem in the regulariza-
tion case (Grodzevich, 2004) for the EIT problem. In this method the regularized solution was obtained
using a parameterized trust region approach to estimate the region of maximum curvature of the L-
curve.
66
2. TRUST REGION SUBPROBLEM
2.1 Structure of Problem
We formulate the EIT regularization problem as the following quadratically constrained least squares
problem
R c,. := min
IIJX -
bUll 2
s·t.lxII2$~
(3)
with ~>O.Using the method of Lagrange multipliers it was shown (Elden, 1977) that this formulation is
equivalent to Tikhonov regularization. This means that any solution to (3) with the value of ~ is equiv-
alent to solving the following Tikhonov regularization with parameter a
(4)
It gives
x a =
xc,.. The solutions for (3) are the same as the following problem, which can be formulated
by squaring both objective and constraint in (3),
1]
c,.:= min
Q(x):= x
T
Hx -
2G
T
x
s.t. "x"~ S ~2
(5)
where
H =
JT J and G =
JT 8U . The optimization problem in (5) is called the trust region subproblem
(TRS). The solution for the regularization problem (3) is found by solving the TRS (5) sequentially.
The TRS can be used to form the L-curve,
L(J,OU):::: {(log(A),logIIJxl'! -oul\): A > o}
(6)
The regularization parameter, ~ can be found through the point of maximum curvature, or the elbow,
on the L-curve (Hansen, 1992). Therefore, the trust region radius, 6. needs to change iteratively to steer
the algorithm to the elbow of the L-curve. A feasible vector
x· =
Xc,. is a solution to (5) if and only if
(Sorensen, 1997):
(H -
A*J)Xc,. = G
(H-A*I)~O
A: (lIxc,.ll~ - d
2
) = 0
),* sO
(7)
for a Lagrange multiplier
,.1* =
,.16.. These conditions relate the Tikhonov regularization (4) to the TRS
(5). A solution
x* to the TRS (5) is a solution to (4) corresponding to
a
2 = -26.
(8)
67
Conversely any solution Xa to (4) for a, solves (5) for Ll2
= jjxa II~ . So the objective function value in
(3) and ~ represent a unique point on the L-curve.
2.2
Creation of L-Curve With TRS
From (Stem and Wolkowicz, 1995) we know the strong Lagrangian duality is satisfied for TRS so that
TRS can be formulated as an unconstrained concave maximization problem,i.e.
17fl = max
minL(x, A)
A:::;O x
where
L (x, A) =
X T Hx - 2 G
T X +
A (~2 -llxll ~ ) represents Lagrangian of TRS. Define
(9)
D(t) :::
(t - G
T), k(t):= (Ll2 +
1),.1, J
(D(t)) -
t
(10)
-G
t
where the AJ(D(t» is the smallest eigenvalue for D(t). Then the unconstrained dual problem for TRS is
as follows
17 fl = max
kef)
f
(11)
The L-curve is formed using ~ in TRS as a parameter and finding the residual for the corresponding
optimal X.c\ • The L-curve can be formed using any of t ,~, and Afl . They are interchangeable and used
to parameterize the regularization problem to give points on the L-curve. These parameters are related
to each other as follows:
t =
At. +BUT J(JT J -At.1)-1 JT 8U
~2 =
c5[!T J(JT J -
A
fl
1)-2 JT &J
At. = Al
(D(t»
More details regarding derivation of the above parameters are given in (Grodzevich, 2004).
3. EXPERIMENTAL MEASURMENT SETUP
3.1 Measurements on phantom
(12)
To validate the proposed approach, experiments were conducted on a cylindrical phantom of 10cm
height and 5cm radius containing 48 circular stainless steel electrodes in three rings with each ring
composed of 16 electrodes connected peripherally around the cylinder. The radius of each electrode
was 0.3 cm and the gap between electrodes was 1.35 cm. The phantom was filled with 0.2M saline and
connected to a newly designed EIT system. Our EIT system currently consists of 48 channels operating
at multiple frequencies from 100Hz-125 KHz that are continuously selectable. A digital signal proces-
sor (DSP) is used to control the operation of every module. In order to convert measured raw data into
68
amplitude and phase information, we have implemented a PC-based phase-sensitive detection, which
uses lock-in amplifiers.
The test objects were a metal rod 6 cm tall with diameter 1.5 cm and a plastic rod with 5.5 cm height
and 1.3 cm diameter. These are high contrast objects when used in the saline solution background.
Figure 1 depicts the objects and the overhead view of the electrode geometry.
Fig.l A cylindrical phantom with two test objects.
The metal rod was 1.5 cm from electrode 1 and the plastic rod at the same depth near electrode 9. The
adjacent current injection pattern and adjacent voltage measurement ('adjacent'-'adjacent') were im-
plemented and in total 2160 pairs of measurements were obtained. The frequency of the sine wave
current was chosen for this test to be 125 KHz with a peak current of 4 rnA. Two sets of measurements
without and with test objects were performed. Each set of measurement was repeated fifty times for the
purpose of standard deviation calculation. The calculated standard deviation was used as an estimation
of measurement noise, which was used to determine the stopping point for the CGLS approach.
3.2 Summary of Results
Figures 2-3 show the background-subtracted images for two TRS and CGLS methods. The metal rod
appeared as red while the plastic rod appears blue.
Figure 4 shows the different points on L-curve for
both methods. The best possible solution that was obtained using the Tikhonov method is also shown in
the graph.
69
z=8.0
z=5.0
z=4.0
z=3 .0
z=2.0
z=1 .0
Fig.2 Reconstructed images with TRS at different heights of phantom.
z=8.0
z=5.0
z=4.0
z=3.0
z=2.0
z=1.0
Fig.3 Reconstructed images with CGLS at different heights of phantom.
70
-3
7.5X 10
7
6.5
~ 6
Ci
o
=-5.5
C1J
:J
" .~
5
Q::
4.5
4
*
0.5
o
*
* 0
* 0
*
* **. • ~
CGLS solution
TRS solution
~ + __
Tlkhonov solution
*
1
1.5
2
2.5
3
solution, I og(ll xJl)
Fig.4 L-curve with TRS, CGLS, and best solution with Tikhonov.
Eventhough the curve is not strongly L-shaped, both horizontal and vertical parts are distinguishable.
The TRS algorithm was able to locate the elbow and this was close to the best solution provided by the
Tikhonov approach.
Figures 5 shows the reconstructed images with CGLS where the level of mea-
surement noise was set to be 30% of the calculated standard deviation.
Figure 6 shows the associated
L-curve for CGLS.
4. DISSCUSION
The results show that the TRS was able to reconstruct images with the same performance as CGLS.
The two test objects were clearly distinguishable at correct position at all different heights as indicated
in
Figure 2. The TRS algorithm was able to follow points on the curvature and locate the elbow as
shown in
Figure 4. The TRS selected solution (marked with +) was close to the best possible Tikhonov
solution.
CGLS is known to be a fast and robust approach that has been used for regularization of ill-posed prob-
lems. It requires very low memory. The disadvantage of CGLS is the behaviour is called semiconver-
gence , which requires a knowledge of exactly when to terminate the iteration. The stopping criterion
for CGLS is based on the discrepancy principle. It terminates when the residual is smaller than a prede-
fined level. This level is set based on the norm of measurements noise. In this paper the norm of the
measurement noise was estimated from the standard deviation of repeated measurements. When the
norm of the noise was set equal to the standard deviation of the measurements, the CGLS converged to
71
the same result as TRS (Figure.3). The CGLS results at each iteration were shown as circles above the
TRS points in Figure 4.
z=8.0
z=5.0
z=3.0
z=2.0
Fig.S Reconstructed images with CGLS at30% standard deviation.
-3
7.5 x 10
7
6.5
~
6
~
'5) 5.5
.Q
~ 5
'0
'iii
~ 4.5
*
o
*
* 0
* 0
* *,0 0
4 TRS solution-- -"') +
3.5
o
o
* 0 0
000 0
z=4.0
z=1 .0
o 000.
CGlS soluli n----... )
0.5
1.5
2
~5
3
15
4
4.5
solulion,log<llxlI)
Fig.6 L-curve with TRS and CGLS at 30% standard deviation.
5
The advantage of the TRS method compared to CGLS is that the former does not require any know-
ledge of the noise for its process. As it shown in Figures 5-6 when the norm of noise is estimated to be
72
30% of standard deviation of measurements the CGLS required many more iterations to terminate with
a much worse final result (Figure.5).
An additional advantage of the constrained least square approach in comparison to the Tikhonov me-
thod is that the physical properties of the problem could be used to estimate the norm of constraint 1:1.
5. CONCLUSIONS
We have successfully implemented a quadratically constrained least square approach to the EIT regula-
rization problem. This approach was solved using a parameterized trust region to estimate the region of
maximum curvature of the L-curve. The experimental results have proven the feasibility and benefit of
this technique for the EIT problem.
6. ACKNOWLEDGMENT
Authors from McMaster University would like to thank for NSERC for support for this project via a
discovery grant awarded to Dr.G.R Moran.
73
REFERENCES
Bj6rck, A., 1996. Numerical Methods for Least Squares Problems. SIAM, Philadelphia
Elden, L., 1977. Algorithms for the regularization of ill-conditioned least squares problems. BIT
17,134-145
Hansen, P.C., 1992. Analysis of discrete ill-posed problems by means of the L-curve. SIAM Rev.,
34(4),561-580
Grodzevich, 0., 2004. Regularization using a parameterized Trust-Region subproblem. M.Sc. Thesis,
Department of Combinatorics and Optimization, University of Waterloo, Canada
Rojas, M., 1998. A Large-Scale Trust-Region Approach to the Regularization of Discrete Ill-Posed
Problems. Ph.D. Thesis, Technical Report TR98-19, Department of Computational and
Applied Mathematics, Rice University, Houston
Sorensen, D.C., 1997. Minimization ofa large-scale quadratic function subject to a spherical constraint.
SIAM Journal on Optimization, 7(1), 141-161
Stem, RJ., and Wolkowicz, H., 1995. Indefinite trust region subproblems and nonsymmetric eigenva-
lue perturbations. SIAM J. Optimiz., 5(2), 286-313
74
CHAPTER V
Paper IV
The following paper describes the design, implementation, and testing of a 48-channel multi-frequency EIT
system. The operating frequency for the designed EIT system ranges from 0.1 to 125 kHz. The system
specifications are comparable with the existing EIT systems with the capability of 3D measurement over
selectable frequencies. A phase-sensitive detection is implemented to recover the absolute phase and
amplitude of the measured signal over all channels for 3 dimensional (3D) measurements. The designed
EIT system offers image reconstruction of both conductivity and permittivity distributions in 3D.
The work presented in this paper was performed by me under the supervision of Dr. Moran. Kenrick Chin
has helped me in the design and fabrication of the hardware part of the EIT system. The majority of the
user control software was developed by Aravinthan Jegatheesan. The manuscript was written by me and
edited by Dr.Moran.
(Manuscript submitted electronically to lOP journal of Measurement Science and Technology, Manuscript
# MST/2549911PAP)
75
A nsp Based Multi-Frequency Electrical Impedance Tomography
System
Mehran Gohariani
, Aravinthan Jegatheesan 2, Kenrick Chini
, and Gerald.R.Moranli ,3
IMedical Physics and Applied Radiation Sciences, McMaster University, Hamilton, Ontario, Canada.
2McMaster School of Biomedical Engineering, McMaster University, Hamilton, Ontario, Canada
3Diagnostic Imaging Hamilton Health Sciences, Hamilton, Ontario, Canada
ABSTRACT
This paper describes the design of a novel multi-frequency EIT system which provides
a flexible mechanism for addressing up to 48 electrodes for imaging conductivity and
permittivity distributions. A waveform generator based on a digital signal processor is
used to produce sinusoidal waveforms with the ability to select frequencies in the range
of 0.1 to 125 kHz. A software based phase-sensitive demodulation technique is used to
extract amplitudes and phases from the raw measurements. Signal averaging and
automatic gain control are also implemented in voltage and phase measurements.
System performance was validated using a Cardiff-Cole Phantom and a saline filled
cylindrical tank. The signal-to-noise ratio (SNR) using saline tank was greater than 60
dB and the maximum reciprocity error less than 4% for most frequencies. The
common-mode rejection ratio (CMRR) was nearly 60 dB at 50 kHz. Image
reconstruction performance was assessed using data acquired through a range of
frequencies. This 1;:IT system offers image reconstruction of both conductivity and
permittivity distributions in three dimensions.
Keywords: Electrical Impedance Tomography (EIT), Multi-frequency EIT, Hardware design.
1. INTRODUCTION
The electrical impedance of mammalian tissue varies with frequency and is known to depend on the tissue
structure at the cellular level (Schwan, 1957). The electrical impedance characteristics of different tissues
in vitro have been measured by different researchers (Gabriel 1996; Foster and Schwan, 1989; Pethig and
Kell, 1987). Mainly because of the wide variation in the electrical characteristics of tissue, and hence the
potential image contrast, electrical impedance measures have the potential to be used as clinical diagnostic
tools for tissue pathology (Blad 1996; Crile
et a11992, Griffiths
et al ,1994) monitoring ischemia (Ristic
et
al ,1997), and distinguishing between normal and cancerous tissue (Blad
et al ,1996).
Electrical impedance tomography (EIT) is an imaging technique that reconstructs images based on the
electrical conductivity of tissue in response to an injected alternating electrical current. The first successful
tomographic impedance image were reported by Barber and Brown (Barber and Brown, 1986) using the
Sheffield Mark 1 system (Brown, 1987) at a single frequency. Single frequency EIT gives limited
76
information regarding the electrical characteristics of tissue. Injecting alternating current over a range of
frequencies on the other hand would allow a more comprehensive probe of the tissue properties. The
spectroscopic electrical response of tissue can be characterized by creating a Cole-Cole plot (Cole 1940;
Cole and Cole 1941).
There are several EIT systems designed by different research groups which are operating in different
frequency ranges. The EIT system in (Yerworth
et al ,2002) covers frequencies between 225Hz and 77 kHz
which was designed for imaging brain function. It was designed based on previous 16 electrodes UCH
Mark 1 a but has the capacity to address for 64 electrodes. Halter
et al 2004 built an EIT system operating
between 10kHz and 10 MHz. It has 64 electrodes and was designed for breast imaging. McEwan
et al
2006 has reported multi-frequency EIT system called UCLH Mk 2.5 for acute stroke imaging operating
between 20Hz-256 kHz. Wilson et al 2001 has reported the Sheffield Mk3.5 EIT system which covers 30
frequencies between 2 kHz and 1.6 MHz. Their system has eight electrodes with an adjacent drive/receive
electrode data acquisition.
The work reported here introduces a novel EIT system developed based on a digital signal processor (DSP)
architecture which allows flexibility in design options such as an arbitrary shape of signal generation and
flexibility in the communication between system components and data collection. The prototype system
provides 48 electrodes with continuous frequency selection from 0.1 kHz to 125 kHz. A phase-sensitive
detection is implemented to recover the absolute phase and amplitude of the measured signal over all
channels for 3 dimensional (3D) measurements. The designed EIT system offers image reconstruction of
both conductivity and permittivity distributions in 3D. In the following section the details of the hardware
design are presented. The results outline the implementation and the performance specifications (See
methods) of the EIT system.
2. HARDWARE SYSTEM OVERVIEW
In this section the major block components of the EIT system will be described. The system was built using
a modular design concept. It consists mainly of a DSP based current source, data acquisition boards and a
multiplexing unit. All of these modules are assembled on a main board, which can communicate with a
computer through a serial port. A graphical user control interface was written using Lab View
(www.ni.com). The EIT system is shown in
Figure 1. In the following sections, the design details, the
implementation of the EIT system, and the performance benchmarks of the system signal-to-noise,
reciprocity error, common-mode-rejection, and measurement accuracy, are presented.
77
FIG. 1: Multi-frequency EIT system.
2.1 Com puter
The controlling computer consists of a Windows based 1.7-GHz Pentium-IV with 256 Mbytes of RAM
which contains several expansion cards to house the data acquisition boards.
2.2 Waveform Generator
The waveform generator uses a DSP controller to achieve flexible waveform geometry and high speed
performance. A block diagram of the DSP controller is shown in
Figure 2. The function of the DSP
controller and waveform generator is to generate any number of predetermined current or voltage pulses of
any given shape (sine or pulse in present design), duration, frequency, and repetition rates. Parameters
which determine the full specifications of the desired injection pulse are transferred from the host personal
computer (PC) via a standard RS-232 serial communications port. This provides the maximum flexibility
without the need to modify the hardware. In a future design, the serial interface to the host PC will most
likely be replaced with a universal serial bus (USB) interface.
There are different ways to implement the sine-wave generator on a DSP such as look-up table,
interpolation, polynomials, etc. to the present design generates a sine wave using a seventh-order Taylor
78
polynomial approximation. The DSP generates a sine wave of a desired frequency for the first quadrant.
Trigonometric identities are then applied to fill in the sine wave over the entire cycle, i.e., from 0 to 21[. The
firmware for the DSP is developed using CodeWarrior from Metrowerks (www.metrowerks.com). Code is
written using both C and assembly language. A block diagram and schematic are shown in Figures 2-3.
The DSP used is the DSP56F807PY80 from Freescale Semiconductor (http://www.freescale.com)
(formerly Motorola Inc). The main features of the DSP56F807 are:
• 16-bit Digital Signal Processor with RlSC (reduced instruction set computer) architecture
• clock frequency of 80MHz
• 60K words of flash program memory
• 8K words of random access memory (RAM)
2.2.1 Digital to Analog Converter (DAC)
The DAC is a 16-bit converter, AD669, from Analog Devices (http://www.analog.com). with a settling
time of 2.5J.!s for one least significant bit (LSB) step. With a bipolar output range of -IOV to + 1 OV, one
LSB is equivalent to a voltage step of300J.!V.
2.2.2
Amplifier
The purpose of the output amplifier is to convert the voltage output of the DAC to a low impedance output
voltage or a high impedance current output as desired by the particular experiment. The amplifier chosen is
a monolithic power operational amplifier LM675 from National Semiconductor (www.national.com). The
device operates at a maximum supply voltage of 60V, capable of outputting 3A. The maximum slew rate is
8V/ J.!S. We chose this amplifier because it has wide bandwidth and low input offset voltage.
2.2.3
Current Driver
The DSP is used as a waveform generator to produce the sinusoidal (or pulse) voltage signal with the
frequency range from 100 Hz to 125 kHz. The voltage signal is then converted to current by a precision
balanced voltage controlled current source (VCCS), based on the enhanced Howland topology (Bertemes-
Filho
et aI, 2000). The enhanced Howland current source uses a single operational amplifier with both
negative and positive feedback. This design has only a single active device and the ability to adjust the
output resistance. The schematic diagram of the enhanced Howland circuit is shown in Figure 4. The
output current of the Howland circuit, IL is given in Eq. (1):
(1)
79
Trigger
--~
RS-232
Comm Port
DSP56F807
DSP
16-bit DAC
FIG. 2: Block Diagram ofDSP based waveform generator.
AMP
The output current is independent of the load resistance. The output resistance is given in Eq. (2)
R =
R1R4b (R3 +R4)
out
R2R3 -
RI (R
4a
+R4b
)
-~ Output
(2)
In the actual circuit, the total impedance is limited due to the presence of output stray capacitance. The
presence of stray capacitance degrades the output impedance and hence the circuit performance. The
presence of both positive and negative feedback paths may cause the enhanced Howland circuit to oscillate
under large capacitive loads at high frequencies. To avoid this effect, a 10 pF capacitor was added in
parallel with the negative feedback element (Cook
et ai, 1994).
2.3 Current Injection Multiplexer
Forty-eight channels are implemented as sixteen channels on each of three separate custom printed circuit
boards. Each board can select sixteen electrodes. The output of the balanced current generator can be
injected to any pair of electrodes using multiplexers.
For maximum flexibility, it is desirable to be able to select any pair of electrodes which will be the sources
of the injected EIT current or voltages. This selection must be controlled from the host PC. In a typical
experiment, the injection electrodes are first selected and voltages are recorded from a set of the remaining
electrodes. The injection then proceeds to a different pair of electrodes and the measurements are repeated.
The electronic analog multiplexer chosen is ADG406 from Analog Devices (www.analog.com).This is an
integrated circuit providing a 16-channel analog multiplexer which is capable of selecting anyone of
sixteen channels at switching speeds up to 2MHz. Typical resistance in the ON state is 50n. Leakage
current in the OFF state is typically lnA.The selection of the pair of injection electrodes is controlled by the
host PC. Since a serial communication link has already been established with the waveform generator DSP,
an auxiliary microcontroller unit
(MeV) is used to control the ADG406 analog multiplexers. A block
diagram of this arrangement is shown in
Figure 5.
80
FB-92
93
XTAL EXTAL-
IRQA lea
I ROO tN--
DSP56F807PY80
DSP
A15 17
MAX232
A14 16
R50232
--'1 R2
RD2 9
160 RXOO
A13 15
R 50232 to HOST
~ 12
TD2 10
159 1XDO
A12 14
A11
13
A10 12
p.s
11
PJ3
10
GPIOB7 47
GPIOB6 46
+5V
LTC1487
GPI085 45
~1k
R50485
GPIOB4 44
+RS485 6
RD 1
56
GPIOB3 43
RXD1
330 -RS485 7
T04
55 lXD1
GPlOB2 42
- W +5V 8_
RE 2
54 GROOS
GPIOB1 41
1k
GND 5_
OEkJ
GPIOBO 40
R50485 LAN
+
nsp Schematic Diagram
FIG. 3: Schematic diagram of DSP waveform generator.
81
+5V~
500kHz
GN~
XTALOSC
-Cl
6
[23
CS LllA.C
-
GNO~
AOO69
L 1
DAC
7 OB15
8 OB14
9 0813
10 CE12
11 CE11
12 CE10
13 CE9
14 CE8
15 CE7
16 CE6
17 CB5
18 CE4
19 CE3
20 CE2
21 CE1
22 CEO
28
:£L- 51 [
26
51
25
1k
~121
10k
-
10k
'Ii
2
+ TLD7
f20P
Rl
Input ~""""--
R3
FIG. 4: A modified Howland current driv~r.
R2
R5
LOAD
A local area network (LAN) is created to communicate between the master DSP and any number of slave
MCUs. The physical layer of the LAN is based on a 3-wire RS-485 interface standard operating at 9600
baud. The auxiliary
MCV selected is an ATtiny2313 from Atmel (www.atmel.com).This design is chosen
because it has a number of advantages: It can be expanded to accommodate any number of multiplexers;
The multiplexer components can be physically located a substantial distance from the master DSP; and, the
ATtiny2313, besides being low power, can be placed into STOP mode in order avoid emitting any digital
RF radiation which may degrade sensitive analog measurements.
RS-232
Cornm Port .0(-
DSP56F807
DSP
FIG. 5: Block diagram of multiplexers
RS-485 IAN AThny2313
MCU r---I~
82
AIXJ406
Analog Multiplexer
Injection
CurrentN oItage
16Cbannels
+5V 20
+12V 1
19 51
GND 10
-12V 27
20 52
GND 12
ADG406BN
21 53
ATMEL
MUX
22 S4
ATtiny2312
23 55
MCU
24 56
25 57
A 28 0
26 sa
11 59
LTC1487 '
10 510
RS-485
9 511
+RS485 6_
RD 1
2 RO
PB3 15
14 A3
8 512
-RS485 7_
TO 4
3 TO
PB2 14
15
A2
7 513
+5V 8
RE 2
8 PD4
PB1 13
16 A1
6 514
GND 5
DE 3
P03
PBO 12
17 AO
5 515
WAKEUP 6 PD2
P06 11
18 EN
516
+12V 1
19 51
-12V 27
20 52
GND 12
ADG406BN
21 S3
MUX
22 S4
23 S5
-RS485 10
9 B
24 S6
+RS485 8
7 A
WAKEUF 6
5 TRIGGER
25 57
+5V 4
3 -12V
B 28 0
26 sa
GND 2
1 +12V
11 S9
1 O-pin ribbon cable MUX connector
10 510
9 511
PB7 19
14 B3
8 512
PB6 18
15 B2
7 513
PBS 17
16 81
6 514
PB4 16
17 SO
5 515
PD5 9
18 EN
516
MUX Schematic Diagram
FIG. 6: Schematic diagram of
MeV based multiplexers
83
The command protocol used to communicate between the host PC and the DSP was developed to be
efficient using error-checking procedure. The same protocol is used to communicate between the master
DSP and each slave
MCV. Any command send from the PC to the DSP, intended for any MCU, is simply
relayed to the LAN and hence to all MCUs. Both communications links operate at 9600 baud. In a future
design, the primary link with the PC will be replaced with a USB port which is capable of operating at
much higher rates. A schematic diagram of the multiplexers is shown in
Figure 6.
2.4 Data Acquisition
Analog-to-digital
(AID) signal conversion is performed using a PCI-6251 data acquisition board (DAQ)
from National Instruments (NI) (www.ni.com).This card can sample up to 16 inputs at a rate of up to 1.25
MHz with 16-bit resolution. It has seven programmable input ranges (��lOO mV to
��10 V) per channel. The
analog input of the PCI-6251 is equipped with a differential programmable gain input amplifier (PGIA).
The PGIA precisely interfaces to and scales the signal passed to the analog-to-digital converter (ADC)
regardless of source impedance, source amplitude, DC biasing, or common-mode noise voltages. The PGIA
helps to minimize the digitization noise on the receive channel so that the digitized voltage represents a
sizeable proportion of the full range. A large first-in-first-out (FIFO) buffer holds data during acquisitions
to ensure that no data is lost. The DAQ board can handle multiple
AID conversion operations with direct
memory access (DMA), interrupts, or programmed input-output
(UO). PCI62S1 DAQ card is selected
because of its high resolution and fast sampling rate. The schematic diagram of DAQ board is shown in
Figure 7.
2.5 Software
Software was written to control all components of the EIT system including the waveform generator,
multiplexer,
AID DAQ boards and signal post processing. The software is divided into four components,
i.e., data capture module, data processing module, control and communication module, and graphical user
interface. F or the data capture module, an event-driven module is developed using the Lab VIEW
(www.ni.com) event structure. The module can respond synchronously to the external trigger coming from
DSP. This is important as it synchronizes the data collection between DAQ and DSP. The data processing
module processes the raw data collected by the data capture module. The raw data is converted into
amplitude and phase information using phase sensitive detection (see next subsection). This module
contains different components such as: Fast Fourier Transform (FFT), filtering, and mathematical
calculation. The timing and communication between DSP and DAQ boards are handled in the control and
communication modules. A graphical user interface was built which allows the user to select the desired
current pattern, injection frequency, measurement pattern, current amplitude, and number of waves to
average. The converted data is then passed to a program written in Matlab v.6.S (Mathworks Inc.) which
reconstructs the interior conductivity and permittivity distributions.
84
2.5.1
Phase-Sensitive Detection
In order to convert measured raw data into amplitude and phase information, a virtual implementation of
phase sensitive lock-in amplifier was programmed. A lock-in amplifier takes a reference signal and a noisy
input signal and extracts only that part of the noisy input signal whose frequency matches the reference
signal. The virtual lock-in amplifier consists of three distinct portions; reference signal detection, signal
mixing and signal filtration. A block diagram of the lock-in amplifier is presented in Figure 8. A reference
signal generated by the hardware is fed into an internal synthesizer based on a Fast-Fourier Transfonnation
(FFT) that generates sine and cosine waves consistent with the reference signal. The two synthesized
waves are then mixed with the measured voltage signal coming from the DAQ board. A low-pass filter
following the mixer rejects the remaining high frequency components and retains the DC component which
is proportional to the amplitude of the signal component and the cosine of its phase relative to the phase of
the internal reference. The same process of mixing and filtering is then applied using a cosine wave. Using
the combined results from the two processes, the amplitude and phase of the measured voltage signal can
be determined.
3. METHODS
3.1. System Specifications
In this section a number of tests of system performance characteristics are summarized including common
mode rejection ratio (CMRR), reciprocity, signal-to-noise ratio (SNR) and precision. The testing was
performed using either a Cardiff-Cole phantom, or a saline filled cylindrical tank - as outlined below.
3. 1. 1.
The Cardiff Cole Phantom
The Cardiff-Cole phantom (Griffiths, 1995) is a resistor-capacitor wheel phantom designed to test and
compare multi-frequency EIT systems. The phantom is representative of a cylinder of homogeneous
conductor with complex impedance. The design generates identical complex impedance distributions to
each drive configuration. The design used here has 32 drive and 32 receiver electrodes interleaved. There
are three rings of jumpers inside the phantom. The Cole equation parameters can be changed using the first
set of two jumpers (11&J2). These parameters are specified in terms of characteristic frequency (fc), and
the ratio of the limiting values of resistance at high and low frequencies,
R 0/
Roc; . The third set of jumpers
(13) shorts out simulated electrode impedances.
3.1.2. Saline Cylindrical Tank
A cylindrical phantom was built from acrylic plastic (Plexiglass®). The phantom is 10cm tall with a radius
of Scm. It contains 48 circular stainless steel electrodes in three rings with each ring consisting of 16
85
electrodes located radially around the cylinder. The radius of each electrode was 0.3 cm and the gap
between electrodes was 1.35 cm.
16
Analogue
MUX
Input range
selection
L.....-L_O;_I:_e~_ss---'k ADe
I I~_F_IF_O_ ..... ~::'
FIG.7: Block diagram ofNI-6251 DAQ board (Adapted from M series user manual-371022G-Ol ,
(www.ni.com ).
ASin(OJst + ¢)
Reference
Signal
1= Acos¢
Low-Pass-
I-------t~
Filter ...... - - - -
Q= Asin¢
Low-Pass-
~-----t~
Filter
1------
FIG.8: Block diagram ofFFT-based phase sensitive detection.
86
3.2. Reciprocity Error Measurement
Four-point measurements are used in EIT instrumentation to remove the effects of electrode/skin interface
impedance. Under ideal circumstances a current injected through electrodes 1 and 2 on an arbitrary region
will induce a voltage V12, at electrodes 3 and 4. If the current and voltage sensing electrodes are reversed,
a voltage V34 will be recorded at electrodes 1 and 2. These two voltage measurements should be identical
based on reciprocity theory (Geselowitz, 1971). In practice, these two measurements will be different due
to difference in the resistance and capacitance of the different electrode and electrode/skin interfaces.
Reciprocity error is defined as
RE =
(V12 -
V34 ) x 100%
V12
(3)
The reciprocity error experiment was conducted on the Cardiff phantom (13s removed). It was averaged
over 50 frames of image data. The reciprocity error for each electrode combination was calculated and the
largest error was reported for each frequency.
3.3. Common-Mode Rejection Ratio Measurement
The common-mode rejection ratio (CMRR) shows the ability of the voltage measurement system
(differential amplifier in this case) to reject any signal common to both input electrodes. Under ideal
conditions the output voltage of a differential amplifier is related to two input voltages as VO=Ad(V+-V-),
where Ad is the differential gain. However, in reality the output will be as Vo=A.ieV+-V-)+O.5*As(V++V-),
where As is the common-mode gain. The CMRR was measured using the same approach as has been done
previously using two circuit configurations (Yerworth
et at, 2003). The following equation was used to
calculate the EIT system CMMR (Brown
et ai, 1999)
(4)
where Vern is the common mode voltage and V in is the differential voltage.
3.4. Signal-to-Noise Ratio
The signal-to-noise ratio (SNR) is defined as the ratio of the mean and standard deviations (Yerworth
et al
,2002).It was calculated for 50 sets of data from the cylindrical tank containing 780 ml of 0.2M saline. A
maximum 4-mA peak-to-peak sine wave current was injected. A full scale range of �� 1 V for the DAQ
board was used for all frequencies (0.1-125 kHz).
87
in
~
a:::
a:::
:E
0
100
90 -
80
70
60
50
40
30
20
10
0
0
I-+-system CMRRi
, ~ DAQ-CMRR !
10 20 30 40 50 60 70 80 90 100 110 120 130
Frequency (KHz)
FIG.9: CMRR variation verses frequency for EIT system and DAQ board.
20 ~----------------------------------------------------------------~
18
~
16
g 14
Q)
~
12
(J
o
is. 10
��u Q)
0:::
8
E
::s 6
E
Ox
C'G
4
:E
2
o +---~----~----r---~----~--~----~----r----r----~--~----'---~
o
10
20
30
40
50
60
70
80
90 100 110 120 130
Frequency (KHz)
FIG. 1 0: Maximum reciprocity variation verses frequency.
89
3.5.
Measurement Accuracy
The system's accuracy was estimated by measurement on a network of known resistors (480 Q) and
capacitors (68nF) in parallel. Resistors and capacitors with �� 1 % tolerance were used. The data was
collected at selected frequencies and the amplitude and phase of measured voltages were determined
through phase-sensitive detection. The resulting voltages were compared with expected voltages. The
accuracy at each frequency was reported as the deviation of measured voltages from the theoretical values
as a percentage of the full-scale voltage measurement. The full-scale measurement is based on the 16-bit
��lOVDAQ.
3.6. 3-D Images of Test Objects
To evaluate the ability of the system to reconstruct images, experiments were carried out on the cylindrical
phantom with a few test objects. The first set of test objects was a metal rod 6 cm tall with diameter
1.5 cm
and a plastic rod with 9 cm height and 2 cm diameter. The metal rod was placed 1.2 cm from electrode 1
and the plastic rod at the same depth near electrode 9. The objects are depicted in Figure 13a. The second
set of test objects were the metal rod and a piece of cucumber - depicted in Figure 14a. In this case, the
cucumber was placed close to electrode 3 and the metal rod at the same depth near electrode 11 The
adjacent current injection pattern and adjacent voltage measurement ('adjacent'-'adjacent') were
implemented and 2160 pairs of measurements were obtained in total. Eleven different frequencies, 1 kHz,
2.5kHz, 5kHz, 10kHz, Is.62skHz, 25 kHz, 31.35kHz, 50kHz, 62. 5kHz, 83.3kHz, and 125 kHz were
chosen for this test with a peak-to-peak current amplitude of 4 rnA. Two sets of measurements without and
with test objects were performed. The conductivity and permittivity images were reconstructed using the
difference approach.
4. RESULTS
4.1. CMRR
Figure 9 shows the CMRR variation with frequency. At low frequencies «15 kHz), the CMRR was 80 dB
and decreased below 60 dB at high frequencies (>50 kHz). It decreased further to 50 dB at 125 kHz.
CMRR values measured were close to the values of the instrumentation amplifier of DAQ board alone.
4.2. Reciprocity
The variation of maximum reciprocity error for each electrode combination is shown in Figure 10. It shows
a maximum error of <3% for low to mid frequencies (lkHz-30kHz) and less than 5% up to 100 KHz
frequencies. At the highest frequency the reciprocity error is increased.
88
4.3.
Signal-to-Noise Ratio
The SNR results are presented in Figure 11. It shows a mean SNR of approximately 62 dB for "opposite-
adjacent" configuration over all three sets of electrodes rings. The SNR reaches 60 dB if one set of ring is
used in "adjacent-adjacent" configuration. The SNR perfonnance of the system in "adjacent-adjacent"
configuration using three sets of rings is more than 45 dB. This SNR is independent of the frequency.
4.4.
Measurement Accuracy
Figure 12 shows the accuracy of the system. Data was collected at selected frequencies and the measured
voltages were compared to theoretical values. The values in figure 12 do not represent actual errors in the
measured data, but rather the error relative to the applied current. To make these conversions, the
percentages in the graph must be multiplied by the ratio of the full-scale current (20 rnA) to the applied
peak current (4 mA)-a factor of 5. The system's accuracy reported values for our system are comparable
to those reported by the Dartmouth group (Hartov
et ai, 2000).
4.5.
Saline Tank Imaging Results
From Figure 13 the positions of different targets are clearly visible, that is, they all had more clearly
reconstructed conductivity distribution images. The EIT system performance has also been validated by the
ability to reconstruct the conductivity and pennittivity distribution images of the cucumber in the
cylindrical tank in the 3-D difference image approach (Figure 14).
5. DISCUSSION
5.1.
Summary of Results
We have designed, implemented, and tested a 48-channel multi-frequency EIT system. The system
specifications are comparable with the existing EIT systems with capability of 3-D measurement and over
selectable frequencies.
The most important parameter of the EIT data collection system is the variation of SNR over the range of
operating frequencies. The SNR measurements were perfonned on a 10-cm diameter tank filled with O.2M
saline. The mean SNR profile of measured values is shown in Figure 11. Three different sets ofSNR were
reported. Measurement from only one set of rings was performed for the purpose of comparison with SNR
values reported by other experimenters. The SNR from only one ring approached 60dB independent of
frequency. This SNR value is higher than 40 dB reported for Sheffield MK3.5 (Wilson
et al ,2001) and also
MK3a which was 55dB (Lu, 1995b). In comparison with the UCLH Mark lb our SNR profile shows less
frequency dependence and higher value than the 50dB reported (Yerworth
et ai, 2002). The specification of
60 dB is equivalent to a root-mean-square (nns) noise of less than 0.1 % of the signal level. In the case of
the "adjacent-adjacent" arrangement for three-ring setup the results show variation between 41-45 dB
90
65
60
55
50
45
_ 40
en
~ 35
~ 30
CI)
25
20
15
10
5
0
0
10 20 30
1---Three-ring adjacent injection
...... Single ring- adjacent injection
~ Three-ring opposite injection
40 50 60 70 80 90 100 110 120 130
Frequency (kHz)
F IG.ll: Mean SNR variation with frequency for three different cases; three-ring with adjacent-adjacent
Three-ring with opposite-adjacent and single ring with opposite-adjacent configuration.
1
0.9
0.8
-~
0
C=68 nt, R=480 ohms
; 0.7
ns
0
0.6
I/)
:::l
't; 0.5
...
0 ... ...
(I)
(I)
0.4
>
;;
0.3
ns
Cii
0::
0.2
0.1
0
0
20
40
60
80
100
120
140
Frequency (KHz)
F IG.12: Measurements accuracy variation with frequency relative to full-scale ofDQA measurement.
91
which is in a good range for 3-D imaging. In this case, measurements close to the drive electrode pair
produced higher amplitude than those further away on a different ring.
Common mode errors are shown in Figure 9. The common mode errors are at reasonable levels for a
mUltiplexing EIT system (Boone
et at, 1996, table 6). The common mode errors are considered the major
limitations of EIT which are more critical issues in the case of a multi-frequency EIT system (Boone and
Holder, 1996). There are various sources of common mode signals, such as unequal drive or receive
electrode impedances and the relative position of the drive and receive electrodes (Riu
et at, 1995). In
addition, the current source driver is not perfectly balanced and hence any imbalanced current may pass
through the differential amplifier of the DAQ and generate a common mode signal. To reduce this effect
the voltage measurement needs to use separate signal grounds by means of separate power supplies
(Jennings
et at, 2001). In practice, stray capacitances will cause coupling between independent grounding
signals and thus limiting the cancellation of the common mode signal using this approach.
The overall CMRR of the system was close to that of the differential amplifier in the DAQ. The amplifier
has a CMMR of 90 dB at 1 kHz and 60 dB at 50 kHz. The overall CMRR is greater than 70 dB at 20 kHz
and greater than 50 dB at 125 kHz. There are several factors that degrade the CMRR - the finite common
mode rejection (CMR) of the instrumentation amplifier inside the DAQ and unequal drive or receive
electrode impedances. To reduce the CMRR error, common mode voltage compensation can be
implemented. (Rosell
et at, 1992). Common mode compensation reduces the common mode signal coming
from the measurement system and hence improves the overall CMRR of system (Smith
et at, 1995). In
order to see an improvement, a DAQ with an improved CMRR would need to be purchased or designed.
The measured reciprocity errors for our EIT system are comparable with those reported for the UCLH
Mark Ib (Yerworth
et at, 2002). However, our reciprocity errors increased at the higher frequencies rather
than at both low and high frequencies as reported for the UCLH Mark 1 b.
The present system needs safety consideration based on IEC60 1 for clinical application. IEC60 1 regulation
requires that the injected current must be less than 1 00 ~A rms for frequencies less than 1 kHz and less than
0.1 rnA / kHz up to a maximum of lOrnA for frequencies above 1 kHz.
The system's ability to reconstruct conductivity and permittivity distributions has been validated with the
cylindrical tank. The conductivity and permittivity reconstructed images in Figures 13-14 corresponded to
the true position of the different objects. The system performance enabled clear images from different
objects to be reconstructed in the tank studies. In future, there are some additional planned improvements in
the system design.
92
(a)
z=8.0
z=5.0
z=4.0
z=3.0
z=2.0
z=1 .0
(b)
93
G
G
0
CD
0
0
Plastic
@ CD
Rod
Rod
@
@
@
@
@
(C)
FIO.13: a): Metal and plastic rods inside of saline filled cylindrical tank. b): reconstructed conductivity
images of metal and plastic rods as a series of 2-D slices using subtraction from empty tank at
50 KHz, metal rod is in red and plastic rod is in blue color c): The original position of two
objects relative to electrodes in lower part of phantom. The metal rod was placed 1.2 cm from
electrode 1 and the plastic rod at the same depth near electrode 9.
5.2.
Further Work
The system needs further hardware improvements and attention to safety factors in order to be suitable for
clinical trials. The current prototype is for laboratory testing purposes only and hence safety features must
be developed in anticipation of actual clinical testing on humans. We are planning development of a next
generation multi-frequency system with parallel data acquisition operating at frequencies up to IMHz. In
addition, we are planning to use high resolution 24-bit delta-sigma analog-to-digital converters with no
amplification to digitize the electrode voltages directly and digital subtraction to remove the common-mode
signal. In order to achieve at least 20-bit accuracy this requires the output impedance of the current driver
to be greater than Zo ~ (220
-1)Z Lmax' where ZLmaz is the maximum load impedance (Holder, 2005). In
this situation, the output impedance must be greater than Ion which presents quite a challenge. In our
present design we have chosen to use a multiplexed system as it provides flexibility in selecting electrodes
and performing different measurement with different electrode patterns. In the future design we will
eliminate the multiplexing approach in order to reduce or eliminate the common mode error, as with the
Sheffield Mark 3.5 system (Wilson
et aI, 2001).We will also employ test objects which more closely
simulate human electrical properties (Ooharian
et al ,2007a).
94
(a)
z=8.0
z=s.O
z=4.0
z=3.0
z=2.0
z=1.0
(b)
95
z=8.0
z=5.0
z=4.0
z=3.0
z=2.0
z=1.0
(c)
G
0)
CD
8
G
8
~
Ro d
@
@
@
(d)
Fig.l4: a): cucumber and rod inside of saline filled cylindrical tank, b): reconstructed conductivity images
of cucumber and rod using subtraction from empty tank as a series of 2-D slices at 125 KHz,
metal rod is in red and cucumber is in blue color c): reconstructed permittivity images of
cucumber and rod using subtraction from empty tank as a series of 2-D slices at 125 KHz, metal
rod is in blue and cucumber is in red color d): The original position of two objects relative to
electrodes in lower part of phantom. The cucumber was placed close to electrode 3 and the
metal rod at the same depth near electrode 11.
96
6. CONCLUSION
We have designed an EIT system capable of collecting 3-D measurements and capable of reconstructing
complex impedance distributions. System performance has been checked with a saline tank and a Cardiff
Cole Phantom. Multi-frequency image reconstruction algorithms will be implemented (Goharian
et at
2007b) to determine the viability of this approach as a reliable imaging technique and suitability for clinical
application. There are several EIT systems currently being used for clinical or industrial purpose but most
of them are for 2-D image reconstruction. We believe that this design offers 3-D measurement and also
capable of reconstructing 3-D complex impedance distributions.
7.
ACKNOWLEDGMENT
Authors would like to thank for NSERC for support for this project via a discovery grant awarded to
Dr.G.R Moran.
97
REFERENCES
Barber, D.C., Brown, 8.H., 1984. Applied potential tomography. J Phys E: Sci Instrum. 17 ,723-733
Bertemes-Filho, P., Brown, B.H., Wilson, A.J., 2000. A comparison of modified Howland circuits as
current generators with current mirror type circuits. Physiol. Meas. 21 ,1-6
Blad, 8., 1996. Clinical application of characteristic frequency measurements: Preliminary in vivo study
Med. BioI. Eng. Comput. 34 ,362-365
Blad, B., Baldetorp, 8., 1996. Impedance spectra of tumor tissue in comparison with normal tissue: A
possible clinical application for electrical impedance tomography. Physiol. Meas 17, A I 05-
AI15
Boone, K.G., Holder, D.S., 1996. Current approaches to analogue instrumentation design in electrical
impedance tomography. Physiol. Meas. 17,229--47
Brown, B.H., Seagar, A.D., 1987. The Sheffield data collection system. Clin Phys Physiol Meas. 8 (suppl
A),91-7
Brown, B.H., Smallwood, RH., Barber, D.C., Lawford, P.V., Hose, D.R, 1999. Medical Physics and
Biomedical Engineering (Bristol: Institute of Physics Publishing)
Cook, R.D, Saunier, G.1., Gisser, D.G., Goble, J., Newel, le., Isaacson, D., 1994. ACT3: a high-speed,
high-precision electrical impedance tomography. IEEE Trans. Biomed. Eng. 41, 713-22
Cole, K.S., 1940. Permeability and impermeability of cell membranes for ions Symp. Quant. BioI. 8, 110-
22
Cole, K.S., Cole, RH., 1941. Dispersion and absorption in dielectrics. J. Chern. Phys. 9,341-51
Crile, G.W., Hosmer, H.R, Rowland, A.F., 1992. The electrical conductivity of animal tissues under
normal and pathological conditions Amer.J. Physiol. 60, 59-106
Foster, K.R., Schwan, H.P., 1989. Dielectric properties of tissues and biological materials: a critical review
Crit.Rev. Biomed. Eng. 17,25-104
Gabriel, c., Gabriel, S., Corthout, E., 1996. The dielectric properties of biological tissues-I: Literature
survey Phys. Med. BioI. 41,2231-2249
Gabriel, S., Lau, RW., Gabriel, c., 1996. The dielectric properties of biological tissues-II: Measurements
in the frequency range 10 Hz to 20 GHz .Phys. Med. BioI. 41, 2251-2269
Geselowitz, D.B., 1971. An application of electrocardiographic lead theory to impedance plethysmography
IEEE Trans. Biomed. Eng. 18,38--41
Goharian, M., Moran, G.R,Wilson, K., Seymour, C., Jegatheesan, A., Hill, M., Thompson, R.T., Campbell,
G., 2007a. Modifying the MR!, Elastic Stiffness and Electrical Properties of Polyvinyl
Alcohol CryogeJ Using Irradiation, Nucl.Instr. and Meth. in Phys. Res. B, in press.
Goharian M, Bruwer MJ, 1egatbeesan A, Moran GR, MacGregor 1F., 2007b. A Novel Approach for EIT
Regularization Via Spatial and Spectral Principal Component Analysis. Physiol. Meas, in
press
Griffiths, H., Jossinet, J., 1994. Bioelectric tissue spectroscopy from multifrequency EIT.Physiol. Meas.
Suppl. 2A, 15, 29-35
98
Griffiths, H., 1995. A Cole phantom for EIT Physiol. Meas. 16 (suppl) A29-A38
Halter, R., Hartov, A., Paulsen, K.D., 2004. Design and Implementation of a high frequency electrical
impedance tomography system Physiol. Meas. 25, 379-390
Hartov, A., Mazzarese, R., Reiss, F., Kerner, T., Osterman, S., Williams, D., Paulsen, K.A., 2000. A multi-
channel continuously-selectable multi-frequency electrical impedance spectroscopy
measurement system IEEE Trans.Biomed. Eng. 47, 49-58
Holder, D.S., 2005. Electrical Impedance Tomography: Methods, History and Application. lOP Bristol,
UK.
Jennings, D., Schneider, I.D., 2001. Front-end architecture for a multi-frequency electrical impedance
tomography system Med &Bio!.Eng.&Comput 39,368-374
Lu, L., Brown, B.H., Barber, D.C., Leathard, A.D., 1995a. A fast parametric modeling algorithm with the
Powell method. Physio!. Meas. 16, A39-47
Lu, L., 1995b. Aspects of an electrical impedance tomography spectroscopy (EITS) system PhD Thesis
University of Sheffield
McEwan, A., Romsauerova A, Yerworth R, Horesh L, Bayford R, Holder D., 2006 Design and calibration
of a compact multi-frequency EIT system for acute stroke imaging Physiol. Meas. 27, S199-
S210
Pethig, R., Kell, D.B., 1987. The passive electrical properties of biological systems: their significance in
physiology,biophysics and biotechnology Phys. Med. BioI. 32, 933-70
Ristic, B., Kun, S., Peura, R., 1997. Muscle tissue ischemia monitoring using impedance spectroscopy:
Quantitative results of animal studies Proc. IEEE Engineering in Medicine and Biology
Society. 2108-2111
Riu,
P.I., Rosell, 1., Lozano, A., Palbis-Areny, R., 1995. Multi-frequency static imaging in electrical
impedance tomography: partl instrumentatiion requirement Med&Biol.Eng.&Comput 33,
784-792
Rosell, 1., Riu, P., 1992. Common-mode feedback in electrical impedance tomography
Clin.Phs.Phsio.Meas 13, AII-14
Schwan, H.P., 1957. Electrical properties of tissue and cell suspensions Advances in Biological and
Medical Physics. 5, 147-208
Smith, R. W.M., Fresston, LL., Brown, B.H., 1995. A real- time electrical impedance tomography system
for clinical use-design and preliminary results IEEE Trans.Biomed. Eng. 42, (2) 133-140
Wilson, A.1., Milnes, P., Waterworth, A.R., Smallwood, R.H., Brown, B.H., 2001. Mk3.5: a modular,
multi-frequency successor to the Mk3a EIS/EIT system Physio!. Meas. 22, 49-54
Yerworth, R.I., Bayford, R.H., Cusick, G., Conway, M., Holder, D.S., 2002. Design and performance of
the UCLH mark 1 b 64 channel electrical impedance tomography (EIT) system, optimized for
imaging brain function Physiol. Meas. 23, 149-158
Yerworth R.I, Bayford RH, Brown B, Milnes P, Conway M, Holder DS., 2003. Electrical impedance
tomography spectroscopy (EITS) for human head imaging Physiol. Meas. 24,477-489
99
CHAPTER VI
Paper V
The following paper discusses the performance tests for the proposed algorithms under experimental
conditions. The performance of three proposed algorithms was compared to standard BIT image
reconstruction algorithms for different objects. In the case of the multi-frequency analysis the peA-based
approach represents a substantial improvement over the Gauss-Newton technique at a single frequency
analysis in terms of systematic error reduction.
The work presented in this paper was performed by me under the supervision of Dr. Moran. The
manuscript was written by me and edited by Drs. Moran, MacGregor and Bruwer.
(Manuscript submitted electronically to Physics in Medicine and Biology, Manuscript #
PMB/2587671PAP1161553)
100
3D Multi-Frequency Electrical Impedance Tomography Imaging in
Phantoms
Mehran Goharianl
, Aravinthan Jegatheesan2
, Mark-John Bruwer3
, John.F.MacGregor3
, Gerald.R.Moranl
,4,
!Medical Physics and Applied Radiation Sciences, McMaster University, Hamilton, Ontario, Canada.
2McMaster School of Biomedical Engineering, McMaster University, Hamilton, Ontario, Canada
3Department of Chemical Engineering, McMaster University, Hamilton, Ontario Canada
4Diagnostic Imaging Hamilton Health Sciences, Hamilton, Ontario, Canada
ABSTRACT
We have recently built and reported a DSP-based multi-frequency (0.1 kHz to 125
kHz), 48 channel electrical impedance tomography (EIT) system. In this paper, we
present the results of a series of phantom experiments to show the baseline imaging
performance of our EIT system. Conducting and nonconducting objects with various
widths were positioned at different distances from a Plexiglas tank edge (10 cm
diameter) filled with 0.2 M saline solution. The results suggest that the lcm diameter
conductor can be detected while placed near to the centre of the tank, using the
difference imaging approach. The conductivity and permittivity images for a piece of
cucumber were reconstructed using four different approaches: dog-leg, principal
component analysis (PCA), Gauss-Newton, and difference imaging. In the case of the
multi-frequency analysis, the PCA-based approach provided a substantial improvement
over the Gauss-Newton technique in terms of systematic error reduction. Our EIT
system recovered a conductivity value of 0.08 Sm'! for the 0.07 Sm·l piece of
cucumber (14% error).
Keywords: Electrical Impedance Tomography (EIT), Image reconstruction, multi-frequency EIT
1. INTRODUCTION
Electrical impedance tomography (EIT) is an imaging technique that reconstructs images of an object's
electrical properties such as conductivity and permittivity. An array of electrodes is connected on the
surface of the object to be imaged. A small current is directed into a subject through a pair of these
electrodes; the current pathway inside the object is determined by the internal conductivity, pennittivity and
the geometry of the subject. A number of different schemes are possible for current excitation, for instance
trigonometric, adjacent and opposite drive pairs. The most commonly used injection current approach is the
adjacent drive method also known as the
neighbouring method (Barber, 1989; Hua
et aI., 1987). In this
approach current is passed through two adj acent electrodes and the voltages differences are collected from
successive pairs of adjacent electrodes. Another common method is the so-called
opposite method (Avis
and Barber, 1994). In this approach current is i~ected through electrodes that are 180
0
apart while voltage
differences are measured with respect to one reference electrode adjacent to the current electrode.
101
Any local perturbation in the conductivity distribution affects the current pathway which eventually
transfers this effect to the boundary voltages. The relationship between these boundary voltages and
intemallocal changes is highly nonlinear, that is large changes in the interior conductivity or permittivity of
the body give rise to very small changes in the surface voltage measurements (Breckon and Pidcock, 1987).
The aim of EIT is to use these limited numbers of surface voltage measurements to reconstruct the
conductivity and permittivity distribution inside the object. Therefore, the EIT image reconstruction
problem is known to be severely ill-posed. There are two main routes to reconstruct an image in EIT either
through linearization of the problem (Cheney
et ai, 1990) or using an iterative process (Paulsen and Jiang,
1997).
EIT has a wide range of applications in both industrial process control (Brown, 2001) and in clinical
environments. EIT has been used for different medical applications such as head imaging (McEwan
et al
2006 ,Holder 200 I), lung imaging (Cherepenin
et al ,2002; Mueller
et ai, 2001), and breast imaging (Soni
et al ,2004; Kerner
et ai, 2002; Cherepenin
et ai, 2001; Osterman
et al ,2000). EIT has been used as a
diagnostic tool to map the functional activity inside the human body (Holder, 2005).
We have recently designed and constructed a multi-frequency EIT system that is capable of continuously
selecting the operating frequency from 0.1 to125 KHz (Goharian
et ai, 2007a). This instrument is capable
of reconstructing both the conductivity and the permittivity distributions at multiple frequencies. This paper
presents the evaluation of the imaging capability of this new EIT hardware system. In the first section of
this paper, a brief overview of the hardware and software is presented. In the following sections an outline
of the experimental imaging protocol is explained. The results are presented as the maximum detectable
depth as a function of object width, which shows the overall EIT system performance. Three different
image reconstruction methods were used to recover conductivity distributions for a piece of cucumber.
2.
OVERVIEW OF EIT HARDWARE SYSTEM
A block diagram of our EIT system is shown in
Figure 1. This EIT system currently supports 48 electrodes
but it has the capability to expand up to 256. In this system for the purpose of multi-frequency EIT and to
improve measurement accuracy, a digital signal processor (DSP) based system was designed using a
modular structure. It includes mainly the following modules: a multi-frequency DSP based current source,
a micro-controller based multiplexing system, a data acquisition module (DAQ), a software based phase
sensitive detection system, and a user friendly real-time controlling software. The high speed DSP controls
all aspect of the communication and processes between each module. The main specifications of each
module are as follows:
a) DSP waveform generator and current driver: The purpose of the DSP controller and waveform generator
is to create current pulses of any given shape (sine curve or pulses in present design), with the flexibility to
102
change the duration, frequency, and repetition rates. The generated pulses are passed through the digital-to-
analogue converter (DAC) and an amplifier. A precision balanced voltage controlled current source
(VCCS) is used to convert the voltage signal to current.
b) Micro-controller multiplexer: The output of the balanced current generator can be directed to any pair of
electrodes using the multiplexer system. Three separate custom printed circuit boards each consisting of
sixteen channels were used to switch the forty-eight electrodes. This can be expanded to accommodate any
number of multiplexers in future designs. An auxiliary microcontroller unit (MCU) is used to control the
selection of a pair of driving electrodes. A local area network (LAN) is established to communicate
between the master DSP and any number of slave MCUs.
Software
Control
Program
~________ ~~~2~mrr,uni~!I~________ ~
PCI COMlmunlcaoo, ,, _ _ _ __
_ __
_
1.rigg"',---1
Data
Acquisition
Device
16-bit
ADC
FIG. 1: Hardware block diagram of Multi-frequency EIT
Digital Signal
Processor
~
~
Analog Injection Signal
Amplifier
c) DAQ module: The digital signal conversion is handled by the PCI-6251 data acquisition board (DAQ)
from National Instrument (NI). There are seven programmable input ranges (��100 mV to ��10 V) per
channel in the PCI-6251. The DAQ board has a differential programmable gain input amplifier (PGIA)
which helps to eliminate the digitization noise on the receive channel.
103
d) Phase sensitive detection: In order to obtain amplitude and phase infonnation from the raw data, a PC
based phase sensitive detection mechanism is used. One advantage of using a lock-in amplifier is its ability
to reduce noise, i.e. to improve the signal-to-noise ratio of the signal to be measured. The main component
of the lock-in amplifier is a phase sensitive detector. A lock-in detector takes a reference signal and a noisy
input signal and uses a phase-sensitive detector to extract only that part of the output signal whose
frequency and phase match the reference. The software implementation of lock-in-amplifier has advantages
such as: it reduces the physical size of hardware system, it is more accurate than its hardware
implementation, and it much easier to run over a very wide frequency range in software demodulation.
The hardware perfonnance has been reported in more detail elsewhere (Goharian
et aI, 2007a), where it is
concluded that our EIT system performance is competitive with other recently reported multi or single
frequency systems. The signal-to-noise ratios of greater than 60 dB are achieved on saline phantom
experiments. This SNR value is higher than 40 dB reported for Sheffield MK3.5 (Wilson
et aI, 2001) and
also MK3a which was 55dB (Lu, 1 995a).The maximum reciprocity error achieved is less than 4% for most
frequencies. The measured reciprocity errors for our EIT system are comparable with those reported for the
UCLH Mark Ib (Yerworth
et ai, 2002).The measured resistance and capacitance values of a parallel RC
load are detennined to better than 0.2% in most cases that indicates good accuracy of our system. In the
following section, we provide an outline of the experimental imaging protocol we have adopted.
3.
METHODS
3.1 Experiment Protocol
The perfonnance of the EIT hardware and software were evaluated on a designed cylindrical phantom of
10cm height and 5cm radius made from Plexiglas. The phantom has three rings with each ring composed of
16 embedded electrodes. Each stainless steel electrode has a radius of 0.3 cm and the gap between
electrodes was l.35 cm. The phantom was filled with 0.2M saline and connected to the multi-frequency
EIT system. Eleven different frequencies, 1 kHz, 2.5kHz, 5kHz, 10kHz, 15.625kHz, 25 kHz, 31.3 5kHz,
50kHz, 62.5kHz, 83.3kHz, and 125 kHz were chosen for this test with a peak-to-peak current amplitude of
4 rnA. These frequencies were chosen arbitrary and could have been selected anywhere between 0.1 kHz
and 125 kHz which is the range of operating frequency of our system. The adjacent current injection
pattern and adjacent voltage measurement ('adjacent'-'adjacent') were implemented and 2160 pairs of
measurements were obtained in total. The first set of high-contrast objects were 6 cm tall metal rods of
various diameters (l.2cm, 1.5cm). The nonconducting test objects were 6cm tall plastic (insulating)
cylinders of the same diameters. These test objects were positioned at different depths from the phantom
edge close to electrode 1 (1 cm, 2 cm, 3 cm and 4 cm). The figures will be oriented such that the actual
104
position of electrode 1 appears at the 3 O'clock position. Figure 2b shows the relative position of objects
with respect to electrode 1 in the cylindrical phantom.
In another experiment a piece of cucumber was inserted inside the phantom holder to evaluate the
capability of system to reconstruct both conductivity and permittivity distributions of this test sample at
multiple frequencies. Cucumber was chosen as it is cellular structure should display impedances closer to
what would be expected in vivo. An Agilent (Agilent Technologies, Mississauga, ON) 4294A Impedance
analyzer was used together with an Agilent 16452A liquid test fixture to measure the conductivity and
permittivity of the piece of cucumber. The impedance analyzer operates in the frequency range of 40Hz to
110MHz; however the liquid test fixture only operates accurately up to 30MHz. The measurements were
performed only between 1 kHz and 1 MHz to match the frequency range of interest in EIT experiments.
Using the impedance analyzer, the complex impedance (resistance R, and reactance X), the complex
admittance (G and B), and the capacitance were measured between 1 kHz and 1 MHz. The measurement
procedure is the same as what was implemented in (Goharian
et ai, 2007d).
3.2
Image Reconstruction Algorithm
Two sets of measurements without and with test objects were performed. The conductivity and permittivity
images were reconstructed using the difference approach. The absolute images of cucumber were
reconstructed using the Dogleg technique (Goharian
et ai, 2007b) at a 125 kHz frequency. The difference
method was applied to reconstruct the conductivity and permittivity changes over the frequency range (Lu
et ai, 1995b). In this case the data over multiple frequencies was collected and the images were
reconstructed based on differences between each frequency verses a reference saline tank at the same
frequency. The instrumentation errors, mostly due to stray capacitance, have a large impact on the absolute
imaging accuracy. In the frequency difference approach, the instrumentation errors can be eliminated to a
considerable degree through the use of a difference measurement. A multi-frequency approach using
Principal Components Analysis (PCA) was used to model the frequency dependence of tissue electrical
properties (Goharian
et ai, 2007c). The method is based on Multivariate Image Analysis (MIA) and our
approach is called "spatial-spectral MIA", or "SS-MIA" (Goharian
et ai, 2007c). The PCA approach
models the multi-frequency variation that is directly observable in the data. PCA is a technique used to
project multidimensional data sets to lower dimensional subspaces in order to build a regularized model
that can easily handle the highly correlated frequency-dependent conductivity variations.
4.
RESULTS
In this study the three-dimensional reconstructed images are displayed as a series of two-dimensional
image slices along the z-axis. Figures 2a and 3a indicate the difference images of the metal rod 1 cm and 4
cm from electrode 1 at 125 kHz respectively. The electrode one position is at the 3 o'clock position, which
is shown pictorially in Figures 2b and 3b. When the object is 1 cm from the edge, one can localize the
105
z=8.0
z=5.0
z=3.0
z=2.0
(a)
(b)
a
~
z=4.0
z=1.0
FIG.2:
a) Difference conductivity images of a metal rod with 1.2 cm diameter at 1 cm far from electrode 1
b) the actual position of metal rod relative to electrodes.
106
z=8.0
z=3.0
z=5.0
z=2.0
(a)
t.::\
~)
(b)
z=4.0
z=1.0
FIG.3: a) Difference conductivity images ofa metal rod with 1.2 cm diameter at 4 cm far from electrode 1
b) the actual position of metal rod relative to electrodes.
107
z=8.0
z=5.0
z=4.0
z=3.0
z=2.0
z=1 .0
(a)
(b)
108
@
PlastIc
Rod
@
(C)
~
~
FIG.4: a): Difference conductivity images of a metal rod with 1.2 cm diameter at 1 cm far from electrode
1 and 1.5 cm diameter plastic rod at same depth from electrode 9, both rods were 6 cm tall. b)
Metal and plastic rods inside of cylindrical phantom, c): The original position of two objects
relative to electrodes in lower part of phantom. The metal rod was placed 1 cm from electrode 1
and the plastic rod at the same depth near electrode 9.
target in difference images. The electrode artifacts were clearly more pronounced in the situation where the
target is almost at the centre of the tank (object 4cm from electrode 1). One can still however detect the
presence of the target within the difference image while the object is 4 cm from the edge of tank. We
defined this type of blurry difference image as the maximum discernible depth. Figure 4a shows the
difference images for two targets separated radially by 1 cm (one from electrode 1 and one from electrode
9) at 125 kHz. The metal rod (conductive) was 1.2 cm in diameter and the plastic rod (non-conductive) was
1.5 cm. The positions of the conductive and nonconductive targets are indicated in Figures 4b and 4c.
Figure 5a illustrates the position of the cucumber inside the phantom holder. Whereas figure 5b is a plot of
the measured conductivity of the cucumber, showing that the conductivity of the cucumber has frequency
dependency. The four different techniques which were used to reconstruct the conductivity variation of
cucumber are the dogleg, PCA (SS-MIA), Gauss-Newton, and difference imaging. Figure 5c shows the
reconstructed images using the dogleg technique. Figure 5d shows pictorially the placement of the
cucumber near electrode 3. A comparison between images reconstructed using the SS-MIA approach and
the Gauss-Newton approach is shown in Figure 6. As can be seen, the SS-MIA was able to remove most of
the uncorrelated frequency-dependent systematic error and build images based on the frequency-dependent
conductivity variations of the cucumber. Although at the height ofz=2cm and z=3cm, it might be said that
the cucumber is detected, clearly the Gauss-Newton reconstruction is not of the same quality. Figure 7
shows the conductivity and permittivity images of the cucumber using the difference imaging approach.
109
5.
DISCUSSION
It is shown in (Kerner
et aI, 2000) that the electrode artifacts are much less at the higher frequency; hence
EIT system can achieve to maximum resolution at the highest operating frequencies. A large number of
experiments were performed and images have been reconstructed. The findings are reported in terms of
maximum detectable depth as a function of object width. The best results (highest sensitivity) with the
designed EIT system were obtained when positioning the conductive and non-conductive objects near the
tank edge. The object at 4 cm from the edge of the tank can still be detected within the difference image.
Since the injected current pattern was adjacent therefore it is expected that the detection will be less
sensitive at the centre of the tanle In the adjacent injection pattern the sensitivity would be best at the edges,
as most of the current will pass though the edges of the tank. A better choice to detect objects in the centre
of the tank is to implement opposite current injection patterns since more of this current will pass through
the centre of the tank. Unfortunately in medical applications we do not always know the position of the
objects we are looking for. The better option would be using the optimal current pattern which has been
studied in (Koksal
et at., 1995; Cheney and Isaacson, 1992). Although we are exploring optimal firing
patterns, the purpose ofthis work is to investigate the limits of this prototype.
Figure 4a shows the background-subtracted images for the metal and plastic rods. The metal rod appeared
as red while the plastic rod appeared blue. The two test objects were clearly distinguishable at the correct
position at all different heights as indicated in
Figure 4a. It seems that the size of the plastic rod in the
conductivity image is almost the same as the actual size up to height, z=5cm. The aim of the multi-
frequency EIT approach is to detect the frequency related variation in electrical conductivity that appears in
biological tissues. The data collection of EIT hardware system contributes to the systematic measurement
errors (SchIappa
et aI, 2000).The hardware systematic error might not be identical for all drive/receive
electrodes combinations, hence the errors will produce spurious features on reconstructed images.
110
(a)
0.25
-0.2
E
-
en
>0.15
->
-~ 0.1
"'C
c:
0
00.05
0
1
10
100
1000
Frequency (kHz)
(b)
111
z=8.0
z=3.0
z=5.0
z=2.0
(c)
CD
(d)
z=4.0
z=1.0
8
FIG.5: a) A piece of cucumber inside of saline tank, b) conductivity variation verses frequency for
cucumber from 1kHz up to IMHz c) Reconstructed conductivity images of cucumber using
Dogleg approach at 125 KHz. The cucumber was close to electrode 3, d) position of cucumber
relative to electrodes.
112
z=8.0
z=5.0
z=4.0
0.08
0.07
0.06
0.05
0.04
0.03
0.02
0.01
z=3.0
z=2.0
z=1.0
(a)
z=8.0
z=5.0
z=4.0
z=3.0
z=2.0
z=1 .0
(b)
FIG.6: a): Reconstructed conductivity images of cucumber using PCA approach at 125 KHz. The
cucumber was close to electrode 3. Eleven different frequencies were used for multi-frequency
analysis b): Reconstructed conductivity images of cucumber using Gauss-Newton approach at
125 KHz.
113
z=8.0
z=s.O
z=4.0
z=3.0
z=2.0
z=1.0
(a)
z=8.0
z=s.O
z=4.0
z=3.0
z=2.0
z=1.0
(b)
FIG.7: a): Reconstructed conductivity images of cucumber using difference method at 125 KHz. The
cucumber was close to electrode 3. b): Reconstructed pennittivity images. The cucumber was
close to electrode number 3.
114
These spurious objects can be confused with actual changes in the target objects. A common approach to
eliminate the systematic errors is a subtraction of data images from calibration saline tank data at different
frequencies. It is shown that (SchIappa
et aI, 2000, Griffiths, 1994) the systematic errors are load
dependent, so that the subtraction from the calibration tank data will not remove all the errors when the
system is used for different patients or objects. Load dependence is also referred to as a 'transimpedance
error'. A comprehensive review for errors in multi-frequency EIT system is given in McEwan
et a12007.
In the static imaging approach, a single set of data is used to reconstruct images, consequently this
approach will not eliminate or reduce the systematic error which could occur in reconstructed images. In
this study a PCA based approach is used to reduce multidimensional data sets (capturing multi-frequency
data over series of frequencies) to lower dimensions. The advantage of the PCA-based approach is that it
can remove the frequency-dependent systematic error from the true frequency-dependent conductivity
variation of the obj ect under study.
It is clear, from comparing the images in Figure 5 and 6, that our PCA-based approach represents a
substantial improvement over using a Gauss-Newton image reconstructed at a single frequency to
distinguish between background and the object (here a cucumber).
The dog-leg algorithm has proven to be robust by converging towards solutions as it is shown in Figure 5c.
The dog-leg algorithm is introduced as an alternative method to Levenberg-Marquardt for solving the EIT
inverse problem in single frequency imaging (Goharian
et aI, 2007a). For the case of a cucumber, the
conductivity values reconstructed in the images using both the dog-leg and PCA based approaches can be
shown to make sense physically. For the cucumbers at 125 kHz, our EIT system recovered a conductivity
value of 0.08 Sm-
I
whereas the measurement of the cucumber directly using the HP impedance analyzer
recorded a value of 0.07��O.004 Sm-1
• A graph of measured conductivity variation verses frequency for a
piece of cucumber is shown in Figure 5b. Clearly, our system can recover the conductivity value and
location of the object fairly accurately. Although tissue conductivity changes have been identified on the
images, which were reconstructed using the Gauss-Newton approach, the systematic errors were significant
enough to severely distort these images.
The complex impedance distributions of a cucumber were clearly reconstructed within the tank in the three-
dimensional using the difference method (Figure 7). The conductivity and permittivity reconstructed
images in Figure 7 corresponded to the true position of the cucumber. These images demonstrate the
ability of our EIT system to correctly reconstruct both the conductivity and the permittivity distributions in
three dimensions.
115
6.
CONCLUSION
The great advantage of our EIT system is its versatility. There are several EIT systems currently being used
for clinical or industrial purposes but most of them are for 2-D image reconstruction. Our EIT system
proved that is capable to reconstruct 3-D complex impedance distributions. The results of the saline
phantom experiments clearly show the feasibility of a PCA based method of detecting an object using a
multi-frequency approach. PCA can transform high-dimensional data into a low-dimensional subspace
while retaining most of the intrinsic information of the input data. It is also an effective way to eliminate or
reduce the frequency-dependent systematic error from the true frequency-dependent conductivity variation
of the object. For multi-frequency EIT systems that are operating at higher frequency up to 1 MHz, the
PCA based approach is necessary technique to reduce systematic error as the systematic errors increase at
the high frequency.
7.
ACKNOWLEDGMENT
Authors would like to thank for NSERC for support for this project via a discovery grant awarded to
Dr.G.R Moran.
116
REFERENCES
Avis, N.J., and Barber, D.C., 1994. Image reconstruction using non-adjacent drive configuration. Physio!.
Meas.
16, AI53-AI60.
Barber, D.C., 1989. A review of reconstruction techniques for electrical impedance tomography. Med Phys.
16, 162-169.
Breckon, W.R., and Pidcock, M.K., 1987. Mathematical aspects of impedance imaging CHn. Phys. Physiol.
Meas. 8 (Suppl. A) 77-84
Brown, B.H., 200 l. Medical impedance tomography and process impedance tomography: a brief review
Meas. Sci.Technol.
12, 991-6
Cherepenin, Y., Karpov, A., Korjenevsky, A., Komienko, Y., Kultiasov, Y.U., Mazaletskaya, A.,
Mazourov, D., and Meister, D.,2002. Preliminary static EIT images of the thorax in health and
disease Physiol. Meas. 23, 33--41
Cherepenin, Y., Karpov, A., Korjenevsky, A., Kornienko, Y., Mazaletskaya, A., Mazourov, D., and
Meister, D., 2001. A 3D electrical impedance tomography (EIT) system for breast cancer
detection Physiol. Meas. 22, 9-18
Cheney, M., and Isaacson, D., 1995. Issues in Electrical Impedance Imaging, IEEE Computational Science
and Engineering. 2,53-62
Cheney, M., Isaacson, D., Newell, J., Simske, S., and Goble, J., 1990. NaSER: an algorithm for solving the
inverse conductivity problem Int. J. Imaging Syst. Technol. 2, 66-75
Griffiths, H., 1994. Modelling of systematic errors in dual-frequency EIT Innov. Tech. BioI. Med.
15, 116-
22
Goharian, M., Jegatheesan, A., and Moran, G.R., 2007a. A DSP based Multi-Frequency electrical
impedance tomography, submitted to Measurement Science and Technology.
Goharian, M., Jegatheesan, A., and Moran, G.R., 2007b. Dogleg trust-region application in electrical
impedance tomography Physiol. Meas. 28,555-572
Goharian, M., Bruwer, MJ., Jegatheesan, A., Moran, G.R., and MacGregor, J.F., 2007c. A novel approach
for EIT regularization via spatial and spectral principal component analysis PhysioI. Meas in
press
Goharian, M.,Moran, G.R.,Wilson,K., Seymour, C., Jegatheesan, A., Hill, M, Thompson, R.T., Campbell,
G., 20007d. Modifying the MRI, elastic stiffuess and electrical properties of polyvinyl alcohol
cryogel using irradiation. Nucl. Instr. and Meth. B. in press.
Holder, D., 2005. Electrical Impedance Tomography: Methods, History and Applications, lOP Publishing,
London, UK
Holder, D.S., 1992. Electrical impedance tomography (EIT) of brain function Brain Topogr. 5,87-93
Hua, P., Webster, J.G., and Tompkins, WJ., 1987. Effect of the measurement method on noise handling
and image quality of EIT imaging. In Proc 9th Int Conf IEEE Eng Med BioI Society.1429-
1430.
Kerner, T.E., Hartov. A., Soho, S.K., Poplack, S.P., and Paulsen, K.D., 2002. Imaging the breast with EIS:
an initial study of exam consistency Physiol. Meas. 23, 221-36
Kerner, T.E., Williams, D.B., Sunshine, Osterman, K.S., Reiss, F.R., Hartov, A, and Paulsen, K.D.2000.
Electrical impedance imaging at multiple frequencies in phantom. Physiol. Meas. 21, 67-77
117
Koksal, A., and Eytiboglu, B.M., 1995. Determination of optimum injected current patterns in electrical
impedance tomography. Physiol Meas, 16, A99-AI09.
Lu, L., 1995a. Aspects of an electrical impedance tomography spectroscopy (EITS) system PhD Thesis
University of Sheffield
Lu L, Brown BH, Barber DC, Leathard AD., 1995b. A fast parametric modeling algorithm with the Powell
method. Physiol. Meas. 16,
A39--47
McEwan, A., Romsauerova,
A., Yerworth, R., Horesh, L., Bayford, R., and Holder, D., 2006. Design and
calibration of a compact multi-frequency EIT system for acute stroke imaging, Physio!.
Meas, 27, S199-S210
McEwan, A., Cusick, G., and Holder, D.S., 2007. A review of errors in multi-frequency EIT
instrumentation Physio!. Meas. 28, S197-S215
Mueller, J.L., Isaacson, D., and Newell, J.C., 2001. Reconstruction of conductivity changes due to
ventilation and perfusion from EIT data collected on a rectangular electrode array Physiol.
Meas. 22,97-106
Osterman, K.S., Kerner, T.E., Williams, D.B., Hartov,
A., Poplack, S.P., and Paulsen, K.D., 2000.
Multifrequency electrical impedance imaging: preliminary in vivo experience in breast
Physio!. Meas. 21, 99-109
Paulsen, K.D., and Jiang, H., 1997. An enhanced electrical impedance imaging algorithm for hyperthermia
applications J. Hyperthermia. 13, 459-80
Schiappa, 1., Annese, E., and Griffiths, H., 2000. Systematic errors in multi-frequency EIT. Physiol. Meas.
21,111-8
Soni, N.K., Hartov, A., Kogel, c., Poplack, S.P., and Paulsen, K.D., 2004. Multi-frequency electrical
impedance tomography of the breast: new clinical results Physio!. Meas. 25, 301-314
Wilson,
A.1., Milnes, P., Waterworth, A.R., Smallwood, R.H., Brown, B.H., 2001. Mk3.5: a modular,
multi-frequency successor to the Mk3a EISIEIT system PhysioL Meas. 22, 49-54
Yerworth, RJ., Bayford, R.H., Cusick, G., Conway, M., Holder, D.S., 2002. Design and performance of
the UCLH mark 1 b 64 channel electrical impedance tomography CElT) system, optimized for
imaging brain function Physiol. Meas. 23, 149-158
118
Chapter VII
Summary and Future Work
The key results of this work are summarized in this chapter. This work consists of two main parts. The first
part is related to the development of 3D image reconstruction algorithms for single and multi-frequency
EIT. The performances of the proposed algorithms are confirmed with simulated objects. The second part
relates to the design of novel multi-frequency hardware and the performance testing of the hardware using
a designed phantom. The proposed algorithms are finally tested under experimental conditions using the
designed hardware.
8.1 Development of 3D Algorithms for Image Reconstruction of EIT
This portion of the work is discussed in chapters II, III, and IV of this thesis. In this work, 3D
reconstruction methods for EIT have been proposed. Chapter II has introduced the dogleg algorithm as an
alternative method to the Levenberg-Marquardt for solving the EIT inverse problem. The main aim was to
compare the efficiency and the computational time of these two different methods. It was found that the
dogleg algorithm and the Levenberg-Marquardt converged to the same solution because the functional to
be minimized is same in both methods. The dogleg method combines the Gauss-Newton and the steepest
descent explicitly through the control of the radius of a trust region (Nocedal and Wright, 1999). The
dogleg algorithm method is preferable to Levenberg-Marquardt method because of its significantly shorter
computational time (Goharian
et at, 2007a). The three dimensional EIT problem especially in medical
imaging requires a finite element model with a large number of elements for complex shapes (such as a
breast or a body). In this case the number of finite elements will be so large that it may increase
computational time and storage requirements such that it might cause the computer to run out of memory
during the calculation. Either parallelizing the problem onto several processors, or introducing more
efficient algorithms will be the solution to 3D image reconstruction. Based upon the work presented in this
thesis, it would seem that the dogleg method would be one of the most efficient techniques for the 3D
models.
Chapter III has introduced a novel multi-frequency approach for 3D image reconstruction. The use of
Tikhonov based regularization techniques is equivalent to introducing a priori knowledge to the
reconstruction process. Such approaches guarantee the stability but enforce smoothness in the inverse
119
solution thus eliminating the possibility of non-smooth solutions (Borsic, 2002). The standard EIT
regularization approach at a single frequency needs to include a priori information (Vauhkonen, 1997).
This could be anatomical structure information from another modality or the noise level of the
measurement system for example. In this chapter, we propose a novel approach to build a subspace for
regularization using a spectral and spatial multi-frequency analysis approach. The approach is based on the
construction of a subspace for the expected conductivity distributions using principal component analysis
(PCA). The advantage of this approach is that the principal components and then the regularization matrix
are determined from the statistical nature of the multi-frequency data, a property that is called the
"spectral-distance spectrum," rather than any a priori assumptions about the model structure (Goharian
et
ai, 2007b). The simulations results proved that the reconstructed images obtained with the proposed
method are reliable even at high level of noise in compare to the standard regularization approach. A
quantitative comparison between the proposed multi-frequency approach and standard single frequency
regularization approach proved the ability of algorithm in reduction of the percentage of misc1assified finite
elements up to twelve fold from the initial percentages after five iterations. The biggest advantage of this
technique is that prior information is extracted from the characteristic response of an object at different
frequencies and spatially across the finite elements.
In chapter IV we presented a new approach for the regularization of the ill-posed EIT problem. A
comparison of the proposed technique with the conjugate gradient least square (CGLS) is presented. We
have successfully implemented a quadratically constrained least square approach to the EIT regularization
problem. The proposed approach is based on the trust region subproblem (TRS), which uses L-curve
maximum curvature criteria to find a regularization parameter. This approach was solved using a
parameterized trust region to estimate the region of maximum curvature of the L-curve. We implicitly
changed the trust region radius during each iteration. This forced the algorithm toward the correct radius
that corresponds to the elbow, the point of maximum curvature on the L-curve. In this chapter a comparison
of the TRS method with CGLS for an experimental phantom study was also presented. CGLS is a fast and
robust technique for regularization of ill-posed problems. The disadvantage of CGLS is that requires a
knowledge of exactly when to terminate the iteration process. It terminates when the residual is smaller
than a predefined level. This level is set based on the norm of measurements noise. Our results show that
both methods converge to the same point on the L-curve when the noise level is known. The big advantage
ofTRS algorithm is that it does not require any knowledge of the norm of the noise.
One possible attractive application ofTRS technique could a combination of PC A based approach (chapter
III) with TRS. In this way, TRS can be used as solver for inverse part of EIT, which use L-curve
maximum curvature criteria to find a regularization parameter. At the same time the multi-frequency PCA
approach can be used to regularization matrix from the characteristic response of an object at different
120
frequencies. This could leads to a conclusion that combining both two approaches may result in a better
algorithm that can provide a reliable and a fast way to locate a regularized solution in the absence of prior
knowledge about the final solution.
8.2 Design and Construction of Multi-Frequency Hardware System
The second part of thesis is composed of two chapters namely V and VI. The Chapter V discuses the
design, implementation; and testing of a 48-channel multi-frequency EIT system. There are several EIT
systems designed by different research groups which operate over different frequencies. The system
specifications are comparable with existing EIT systems with the capability of 3-D measurements over
selectable frequencies. We have decided to use multiplexing in our designed system to obtain flexibility in
addressing electrodes and performing experiments with differing approaches. Even though this may
increase the common mode error. A multiplexer design gives flexibility to arrange a four-terminal
measurement that is the best technique to minimize the effect of contact impedance. In four-terminal
measurements two electrodes are used to drive current and two to measure voltage. This method in
principle eliminates series impedances and contact-resistance errors. The electrode impedance on the
voltage measurement becomes negligible as the input impedance of the differential amplifiers is so high
(Richards, 1996). The receive channels used to measure the electrode voltages are also multiplexed in the
existing design. This does not mean however that future design of this prototype may not have an
increasingly parallel acquisition to improve speed. From an image reconstruction perspective, the multiple
sources seem to be a better choice. The multiple drive approach however needs to be accurately matched
for a large number of current drivers in order to reduce common-mode currents. The measured average
signal-to-noise ratio (SNR) using a saline phantom was greater than 60 dB across most frequencies and this
figure was independent of the frequency of the injected current. Our system SNR value is higher than the
Sheffield MK3.5 (Wilson
et ai, 2001) and also MK3a which was 55dB (Lu, 1995). In comparision to the
UCLH Mark 1 b our system shows less frequency dependence the on SNR and higher value than the 50dB
reported by others (Yerworth
et ai, 2002). The common mode errors are at reasonable levels for a
multiplexing EIT system (Boone
et ai, 1996). The overall CMRR of the system was close to that of the
differential amplifier in the DAQ.
A software based phase-sensitive demodulation technique is used to extract amplitudes and phases from the
raw measurements. This type of design gives much more flexibility in comparison to hardware based phase
sensitive demodulation; which has a limited range of operating frequency.
121
There are several EIT systems currently being used for clinical or industrial purposes but most of them are
for 2-D image reconstruction.
Our hardware offers 3-D measurement and also capable of reconstructing 3-
D conductivity and permittivity distributions that has been validated with the cylindrical tank experiments.
The chapter VI discusses the testing of the proposed algorithms with phantom experiments. A piece of
cucumber positioned inside of saline phantom. The performance of four different approaches: dog-leg,
principal component analysis (PCA), Gauss-Newton, and difference imaging were checked in respect of
reconstructed conductivity images. The results show that the PCA-based approach presents a substantial
improvement over the Gauss-Newton technique in terms of systematic error reduction. Our EIT system
recovered a conductivity value of 0.08 Sm-1 for the 0.07 Sm-l piece of cucumber (14% error).
8.3 Future Research for Multi-Frequency EIT
Although this initial prototype has exceeded our initial expectations, our research and testing has
demonstrated some key areas where modifications could be made to improve the system:
• Our designed system needs further hardware improvements and attention to safety factors in order
to be applicable for clinical setup. With respect to patient saftey, for frequencies above 1 kHz, the
injected current is limited by approximately 100 J..lA x f, where f is the frequency of injected
current in kHz. (lEC, 2007). From the safety point of view any dc voltage should not pass to the
injection electrode to avoid the skin irritation or burns. This can be done through dc blocking
capacitors. The disadvantage with this approach is that any DC error in current driver output will
accumulate charge in these capacitors which eventually appears as transient voltage error during
switching electrode from 'drive' to the 'sense' mode in multiplexing system (McEvan
et at ,
2007). A better design is to implement a sensing circuit to detect dc currents (Cherepenin
et ai,
2001).
• The next generation of our designed multi-frequency system should include parallel injections and
data acquisition. In this approach multiple current sources can be used along with no amplification
to digitize the electrode voltages directly and digital subtraction to remove common-mode signal
in data acquisition system. Even though this type of design removes multiplexing system and
associate parasitic impedance however it requires designing higher precision current sources
(Holder, 2005). The high precision current sources will keep the common-mode currents, i.e.
summation of all injected currents, as small as possible.
122
• A novel cryogel-polymer was introduced in the Appendix which offers the opportunity to meet the
requirements for phantoms in order to validate 3D reconstruction algorithms and measurement
systems. In future a simple tissue-mimicking phantom that mimics the electrical conductivity of
the normal tissue as well as benign and malignant lesions should be developed. This will give a
pre-clinical trial for testing the hardware and image reconstruction performance that will be a
crucial step in the development of a breast cancer detection tool. The operating frequencies of our
hardware system should extend up to IMHz to enable better conductivity and permittivity image
reconstruction. In a real experimental data set, the combined conductivity and permittivity data
will show a more complex dependence on frequency. In this case the proposed multi-frequency
approach can be used to assess and distinguish between tissues in different functional states.
• The novel PCA-based method has been developed in chapter three is able to extract the
frequency-dependent structure in multi-frequency data sets. It was proved through simulation and
experimental case studies that the PCA based approach (termed "SS-MIA") achieves a substantial
improvement in clarity of EIT images over standard EIT approaches. One possible attractive
application could be a combination of peA based approach (chapter III) with TRS approach
(chapter IV). The TRS can be inserted into SS-MIA approach as a solver for inverse part ofEIT.
In this way, both the regularization parameter and regularization matrix can be simultaneously
extracted and built from the characteristic response of an object at different frequencies. The TRS
use L-curve maximum curvature criteria to find a regularization parameter and the multi-
frequency PCA approach can build the regularization matrix from the multi-frequency data sets.
This could lead to a reliable and fast approach without imposing any prior assumption about the
final solution.
• One promising current injection pattern from image reconstruction perspective is the optimal
current injection strategy that produces the most uniform sensitivity over the whole volume of
study (Cheng
et al., 1990). This current pattern has some drawbacks which need to be considered
in practical application. Optimal current patterns methods need as many current drivers as there
are electrodes. In this case voltages are measured from the same current-carrying electrodes. To
avoid the contact impedance effect in this approach the electrodes that are so-called compound
electrodes can be used (Hua
et al., 1993). The compound electrode has different parts for the
current injection and voltage measurements which are not in connection with each other. An
example of an optimal current pattern is so-called
trigonometric current patterns (Isaacson, 1986).
This type of current can be easily implement using the DSP based design.
123
• Increasing the number of the electrodes in order to collect much more information has the
potential to improve the image reconstruction precision. For example, if 256 electrodes are used
there would be 256 excitations for adjacent current strategy. Each excitation pair yields 254
measurements this therefore results in total 65,024 dependent and independent voltage
measurements. This extension would enable to approach the imaging of objects of 5 mm in size.
Obviously this approach needs a very fast data collection system, high speed DSP, and USB link
to transfer high data rate.
124
References
Boone, K.G., Holder, D.S., 1996. Current approaches to analogue instrumentation design in electrical
impedance tomography. Physiol. Meas. 17, 229-47
Borsic, A., 2002. Regularisation Methods for Imaging from Electrical Measurements. PhD Thesis. Oxford
Brooks University.
Bruwer, MJ., 2006. Process systems approaches to diagnostic imaging and identification PhD Thesis.
McMaster University, Canada.
Cheng, K.S., Simske, S.l., Isaacson, D., Newell, J.e., and Gisser, D.G., 1990. Errors due to measuring
voltage on current-carrying electrodes in electric current computed tomography. IEEE Trans
Biomed Eng, 37,60-65.
Cherepenin, V., Karpov, A., Korjenevsky, A., Kornienko, V., Mazaletskaya, A., Mazourov, D., and
Meister, D., 2001. A 3D electrical impedance tomography (EIT) system for breast cancer
detection Physiol. Meas.
22 , 9-18
IEC, 2007. 60601-1 Medical Electrical Equipment: I-II. General Requirements for Basic Safety and
Essential Performance-Collateral Standard: Electromagnetic Compatibility-Requirements
and Tests 3rd edn
Isaacson, D., 1986. Distinguishability of conductivities by electric current computed tomography. IEEE
Trans Med Imag. 5, 91-95.
Goharian, M., legatheesan, A., and Moran, G.R., 2007a. Dogleg trust-region application in electrical
impedance tomography Physiol. Meas. 28, 555-572
Goharian, M., Bruwer, MJ., Jegatheesan, A., Moran, G.R., and MacGregor, J.F., 2007b. A novel approach
for EIT regularization via spatial and spectral principal component analysis Physiol. Meas.28,
1001-1016.
Holder, D., 2005. Electrical Impedance Tomography: Methods, History and Applications, lOP Publishing,
London, UK
Hua, P., Woo, E.l., Webster, J.G., and Tompkins, W.l., 1993. Using compound electrodes in electrical
impedance tomography. IEEE Trans Biomed Eng. 40,29-34.
Lu L 1995., Aspects of an electrical impedance tomography spectroscopy (EITS) system PhD Thesis
University of Sheffield
McEwan, A., Cusick, G., and Holder, D.S., 2007. A review of errors in multi-frequency EIT
instrumentation Physiol. Meas. 28, S197-S215
Nocedal, 1., and Wright, S.J., 1999. Numerical optimization (Springer Series in Operations Research) (New
York: Springer)
Richards,1., 1996.1nnovations in impedance testing. Medical Device & Diagnostic Industry
Vauhkonen, M., 1997.Electrical Impedance Tomography and Prior Information, PhD thesis, department of
Applied Physics. Kuopio University.
Yerworth, RJ, Bayford, R.H., Cusick, G., Conway, M., Holder, D.S., 2002. Design and performance of the
UCLH mark 1 b 64 channel electrical impedance tomography (EIT) system, optimized for
imaging brain function Physiol. Meas. 23, 149-158
Wilson, A.l., Milnes, P., Waterworth, A.R., Smallwood, R.H., and Brown, B.H., 2001. Mk3.5: a modular,
multi-frequency successor to the Mk3a EIS/EIT system Physiol. Meas. 22,49-54.
125
Appendix A
Paper VI
The following paper presents a novel cryogel-polymer design that offers the opportunity to meet the
requirements for phantoms in order to validate 3D reconstruction algorithms and measurement systems.
The purpose of this work was to study the effect of radiation on the elastic stiffness, electrical and MRI
properties of polyvinyl alcohol (PV A)-based cryogel (PV A-C).
The work presented in this paper was performed by me under the supervision of Drs. Moran, Campbell and
Thompson. The experimental data for 1.89T and 3T were collected by Kyle Wilson. The data were
analyzed by me. The manuscript was written by me and edited by Drs. Moran, Thompson and Campbell.
(Reprinted from Nucl. Instr. And Meth. B , this article is in press as: Mehran Goharian, Gerald.R. Moran,
Kyle Wilson, Colin Seymour, Aravinthan Jegatheesan, Mike Hill, Terry.R.Thompson, Gord Campbell,
Modifying the MR!, elastic stiffness and electrical properties of polyvinyl alcohol cryogel using irradiation.
Nucl. Instr. and Meth. B. (2007) in press. doi:l0.l016/j.nimb.2007.04.111, Copyright (2007), with
permission from Elsevier)
126
Available online at www.sciencedirect.com
ScienceDirect
Beam Interactions
with Materials & Atoms
ELSEVIER
Nuclear Instruments and Methods in Physics Research B xxx (2007) xxx-xxx
www.elsevier.com/locate/nimb
Modifying the MRI, elastic stiffness and electrical properties
of polyvinyl alcohol cryogel using irradiation
Mehran Goharian a, Gerald R. Moran a,b,c,*, Kyle Wilson c, Colin Seymour a,
Aravinthan Jegatheesan a, Michael Hill a, R.Terry Thompson h,c, Gordon Campbell c,d
a
Medical Physics and Applied Radiation Sciences, McMaster University, 1280 Main St. W, Hamilton, Ontario, Canada L8S 4KI
b
Imaging Division, Lawson Health Research Institute, London, Ontario, Canada
C
Department of Medical Biophysics, University of Western Ontario, London, Ontario, Canada
d
Integrated Manufacturing Technologies Institute, Medical Devices, NRC, London, Ontario, Canada
Available online
Abstract
The aim of this work was to study the effect of radiation on the elastic stiffness, electrical and MRI properties of polyvinyl alcohol
(PVA)-based cryogel (PVA-C). The PVA-C samples were irradiated with a 60CO y-source, at 2.18 x 106 Rads. The indentation measure-
ments (an indication of elastic stiffness) reduced by about 14.6% for PVA-3C and 5.7% PVA-6C after irradiation, indicating that the
material became harder/stiffer. It was found that MRI relaxation times provide an alternative and non-destructive method to evaluate
the radiation effect on PV A-C. The
TI of PV A-C that had undergone three freeze thaw cycles decreased with irradiation by 10%, 25% and
35% at 1 T, 1.89 T and 3 T respectively. The
TI ofPVA-C that had undergone six freeze thaw cycles decreased with irradiation by 18%,
15% and 11% at I T, 1.89 T and 3 T respectively. The
T2 ofPVA-C decreased with irradiation only at I
T, however this change is hypoth-
esized to be due to the interaction of two spin pools in the gel. The electrical conductivity
(a) and permittivity constant
(e) of the unir-
radiated and y-irradiated PV A-C samples were measured at different frequencies in the range 40 Hz to 1 MHz. The results demonstrated
that the conductivity increased with irradiation by 50% for PVA-3C (three freeze thaw cycles) and 75% for PVA-6C (six freeze thaw
cycles) at frequencies greater than 1 KHz.The permittivity decreased with irradiation up to 25% for 3C and 35% for 6C at frequencies
less than 1 KHz.
© 2007 Elsevier B.V. All rights reserved.
Keywords: PVA-C; Irradiation; Conductivity; MRI relaxation times
1. Introduction
A PV A-cryogel (or PV A-C) is a hydrogel generated by
freezing ( - 20��C) and thawing (+20 ��C) an aqueous PV A
solution [1]. The cryogel molecular organizational structure
changes during the freeze thaw cycles as hydrogen bonds
form between water and the hydroxyl groups on the PYA
molecules. As the number of freeze-thaw cycles (FTCs)
increases the degree of hydrogen bonding increases. A
* Corresponding author. Address: Medical Physics and Applied Radi-
ation Sciences, McMaster University, 1280 Main St. W., Hamilton,
Ontario, Canada L8S 4Kl. Tel.: +1 905 525 9140x26887.
E-mail address:morang@mcmaster.ca (G.R. Moran).
0168-583X1$ - see front matter © 2007 Elsevier B.V. All rights reserved.
doi: 1O.1016/j.nimb.2007.04.l11
PV A-C gel phantom has been successfully used to mimic
organs and tissues in MR and ultrasound imaging studies
[2]. Polymer crosslinking is a possible mechanism to
explain the hardening of PYA cryogels. Three models have
been proposed in an effort to explain the mechanisms
inducing the solidification and the crosslinking. They
involve hydrogen bonding, polymer crystallite formation,
or a liquid-liquid phase separation process [3]. The factors
affecting the PV A -C electrical and mechanical characteris-
tics are the number of freeze thaw cycles, PV A concentra-
tion and additive materials.
The methods for determining elastic stiffness include
rheological (storage and loss moduli), materials tensile test-
ing, indentation and ultrasound propagation. Wan et al. [4]
Please cite this article in press as: M. Goharian et aI., Modifying the MR1, elastic stiffness and electrical properties ... , Nucl. 1n8tr. and
Meth. B (2007). doi:l0.l016/j.nimb.2007.04.l11
M. Goharian et al. J Nuc!. Instr. and Meth. in Phys. Res. B xxx (2007) xxx-xxx
Lowed a statistically significant difference among 1-6
eeze thaw cycles for materials tensile testing. The indenta-
::m method was chosen in this study as the measure of
astic stiffness since it can provide a single measurement
ith sufficient sensitivity to compare the relative changes
~ elastic stiffness for each sample.
Campbell et al. have demonstrated that the electrical
~operties of PV A-C can be changed by adding water sol-
ble salts [5]. The measurement of electrical properties is
Ie of the most convenient methods for studying the poly-
ler structure [6]. These properties are affected by the struc-
lre of polymer, type of doping agent, concentration and
.anufacturing process [7]. In an earlier work, we demon-
rated that there are significant changes in the electrical
roperties of PVA-C that undergoes 3FTCs as compared
, 6FTCs [27].
In this same manuscript we had presented MRI results
f plain 3C and 6C (15%) PYA-C. We measured PVA-C
Lat had undergone three or six freeze thaw cycles at three
~fferent fields, 1 T, 1.89 T and 3 T. It was determined that
.} increased linearly with field (419 ms, 755 ms and
)69 ms for 3C and 391 ms, 602 ms and 929 ms for
6C).
was also found that
T2 remained relatively constant
~ound 50-60 ms.
Other groups have performed similar measurements.
or example Lukas et al. [8] have determined the MRI
:laxation times at 1.5 T of a 10% PVA-C over 1-2 freeze
law cycles. They estimate for example a
T} of approxi-
lately 1050 ms and a
T2 of 150 ms at 35��C for 10%
VA-C that has undergone two freeze thaw cycles. Surry
al. [2] also measured MRI relaxation times of 100/0
V A-C, however they measured PV A-C that had under-
)ne from 1 to four freeze thaw cycles. For three freeze
rcles for example, they measured a
T} of approximately
W ms and a
T2 of 115 ms at 1.5 T.
It has been demonstrated that the PVA-C properties
lake it suitable for the construction of medical imaging
rrantoms [2,9,10]. PVA-C polymer offers the opportunity
, meet the phantom requirements suitable for impedance
Iaging techniques, such as electrical impedance tomo-
'aphy (EIT) [11]. Tissue mimicking phantoms made from
VA-C have been used: for MRI to simulate normal
ld diseased tissue [12], for temperature dosimetry in
[RI [8] and for simulations of ultrasound guided breast
.opsy [13].
Past studies have examined the effect of radiation
duced polymerization (cobalt-60 [14,15J , electron beam
6,17]) of PVA-C hydrogels. Degradation rather than
'osslinking occurs with PYA alone. However, PVA-C
ydrogels can be crosslinked by irradiation. It has been
lown in rheological experiments that both the elastic
~operties of the storage (G') and the loss (G") moduli
.creased with absorbed radiation dose from 0 to 30 kGy
8]. There has not been any study to elucidate the com-
ned affects of freeze thawing and radiation processes on
le properties (electrical impedance, magnetic susceptibility
Id elastic) of PYA-C.
Ionizing radiation is an efficient tool for sterilization. A
great part of single-use medical products are sterilized by
this technique [19]. In order to use irradiation to sterilize
PVA-C phantoms, it is crucial that the changes in the mate-
rial properties, in particular the elastic stiffness, the conduc-
tivity and the MRI properties be characterized.
The characterization of these changes also allows the
possibility of using radiation to either "customize" the
properties of a phantom, or to use PV A as a radiation
dosimeter. Two new Fricke gel systems using PV A hydro-
gel and PVA-C have been described [9]. Both systems
showed excellent linear coefficients
(R2
= 0.99), signifi-
cantly lower diffusion coefficients and allowed storage for
longer periods before radiation exposure. The authors con-
cluded that the PV A based gels were a significant improve-
ment over previous Fricke gel systems.
The long term goal of this research is to fabricate PVA-
C gels to be utilized in medical devices and to be used as
tissue mimicking phantoms for imaging research. A
requirement of the imaging studies is that the phantom
constructions be visible with sufficient contrast, when
viewed with both MRI and impedance imaging. The pur-
pose of the present work is to determine the changes
induced in the elastic stiffness, the MRI properties and
the electrical properties of the cryogel, resulting from
y-
irradiation.
2. Materials and methods
2.1. PVA-C manufacture
In this study 15% PVA-C specimens were prepared that
had undergone either three or six freeze thaw cycles. The
details of the manufacture have been detailed elsewhere
[27J. Samples will be designated "PVA-3C" for 15%
PVA-C with three freeze thaw cycles or "PVA-6C" for
15% PVA-C with six freeze thaw cycles.
2.2. Elastic stiffness (indentation)
An indentor was devised that applies a minor load of
0.1 N and a major load of 0.71 N to the flat surface of a
PVA-C specimen through a 2.5 mm diameter hemispherical
foot. The relative depth of penetration between the major
and minor loads was recorded for 10 measurements on
each specimen. This method is similar to the method
described in ASTM D 1415 - Standard Test Method for
Rubber Property - International Hardness [20]. However
the major load applied according to ASTM 01415 was
judged to be too severe for the cryogels used in this study
thus the major load of 0.71 N was used. A smaller depth
of penetration indicates a harder/stiffer material.
2.3. MRI
MRI was performed on all of the samples at field
strengths of 1 T (OrthoOne, ONI corp, Wilmington MA,
Please cite this article in press as: M. Goharian et al., Modifying the MRI, elastic stiffness and electrical properties ... , Nucl. Instr. and
Meth. B (2007), doi:l0.l016/i.nimb.2007.04.111
M. Goharian et al. /
Nuc!. Instr. and Meth. in Phys. Res. B xxx (2007) xxx-xxx
3
www.onicorp.com). 1.89 T (referred to as 2 T here after)
(Magnex magnet, SMIS console) and 3 T (lMRIS 3 T,
Innovative Magnetic Resonance Imaging Systems, Winni-
peg, Canada). The spin-lattice relaxation time,
TJ and
the spin-spin relaxation time,
T2 were measured at each
field strength.
For
TI measurements a standard 2D inversion recovery
sequence was repeated with different inversion times (TI).
The imaging parameters were as follows: 6 em FOV,
64 x 64 matrix, TR = 4000 ms, TE = 13.8 ms, TI = 20,
60, 100, 200, 400, 600, 800 and 1000 ms at 1 T. At 2 T
and 3 T, the same parameters were used except TR =
5500 ms and TI = 20, 60, 100, 200, 300, 400, 500, 600,
700, 800, 1000, 2000 and 3000 ms. Signal intensity was
determined in a cross-sectional region of interest in each
image for each of the PVA-C samples. The
Tl recovery
can be described by the equation:
S(TI) = So(1 - 2exp( - TI/Tt)),
(1)
where S(TI) is the signal intensity of the selected ROI of the
image and So is the signal at time TI = 0. This equation can
be rearranged to form:
TI/T\ = -In[I/2 - S(TI)/(2So)]
(2)
from which
T\ can be obtained by a linear regression where
both So and
TJ are fitted parameters.
T2 was measured at 1 T using a 2D spin echo sequence
with: a 6 cm FOV, 64 x 64 matrix, TR = 4000 ms and
TE = 25, 50, 100, 200, 400 and 500 ms at 1 T. The same
parameters were used at 2 T and 3 T except TR = 5500 ms.
The peaks of the echoes were fit to the equation:
S(TE) = So exp( -
TE/T2) ,
(3)
where S(TE) is the signal intensity of the selected ROI and
So is the signal intensity at TE = 0. In a previous manu-
script [27J we had postulated the presence of two different
T2 components. At 1 T, there was not sufficient resolution
in the time domain to fit to two exponentials. However at
2 T and 3 T the data was fit assuming a bi-exponential
behaviour.
2.4. Impedance analysis
An Agilent (Agilent Technologies, Mississauga, ON)
4294A Impedance analyzer was used together with an Agi-
lent 16452A liquid test fixture to measure the conductivity
and permittivity of the PYA-C. The details of the proce-
dure have been outlined elsewhere [27].
The conductivity,
(J , was calculated using
(4)
where
t is the spacing between the electrodes,
A is the elec-
trode area (radius = 38.0 mm) and
Rp is the (parallel) resis-
tance measured.
The permittivity,
e, was calculated using:
(5)
where Cp is the parallel capacitance measured in the liquid
fixture and
eo = 8.85e- 12 C2/(N m2
) is the permittivity of
free space.
2.5. Irradiation
The PVA-C samples were irradiated with the calibrated
60Co y-rays source in the McMaster Nuclear Reactor. The
source was a cylindrical Cobalt-60 of approximately 5600
Curies (single 'pencil'). The samples were arranged equidis-
tant upon a circle of 15.2 cm from the centerline of the
source. The high dose rate was determined to be 51.4 x
104 Rad/h. The total delivered dose was 2.18 x 106 Rads.
2.6. Statistics
Multiple electrical measurements were made to deter-
mine the average and standard deviation. Standard
t-tests
were performed to determine if measurements before and
after irradiation were different at the
p < 0.05 level.
3. Results
The results of the elastic stiffness (indentation) proce-
dure are summarized in Fig. 1.
The MRI relaxation time measurements are summarized
in Figs. 2 and 3. For unirradiated samples as is shown in
Fig. 2 there is an increase in
Tl with field strength while
T2 is less sensitive to field strength. This behaviour has been
documented before [27]. Fig. 3 shows the
TI and
T2 of irra-
diated PVA-3C and PVA-6C samples measured at the dif-
ferent field strengths. y-irradiation appears to have
decreased the
T] of the PV A samples at all field strengths.
At 1 T, the
T2 data fit a single exponential well whereas at
2 T the data was better fit by a bi-exponential behaviour.
The higher
T2 component (300 ms) is not reported. At
3 T the
T2 data fit a single exponential relaxation.
0.9 ~----,--------,-------,--==____ --,
0 .8 +-------::-::~---I------+-----_+_---___l
0.7
II)
II)
41
:E 0.6
.. II)
.!::! 0 .5
'1i)
~ 0.4
W
0.3
0 .2
0.1
o
No.of freez-thaw cycles
Fig. 1. Variation of the elastic stiffness (indentation) of unirradiated and
irradiated PVA-3C and PVA- 6C.
Please cite this article in press as: M. Goharian et aI., Modifying the MRI, elastic stiffness and electrical properties .", Nucl. Instr. and
Meth. B (200'1), doi:.lO_1016/i.nimh.200J_04.111
M. Goharian et al. /
Nuc!. Instr. and Meth. in Phys. Res. B xxx (2007) xxx-xxx
Relaxation Times Pre-irradiation
1100
900
Ii)
..��.. 700
Q)
E
.. c: 500
.2
10 )(
co 300
Q)
D::
100
-100
MRI Field Strength (Tesla)
:;'ig.2. Variation of
TJ and
T2 values for unirradiated PVA-3C and PVA-
iC for various field strengths.
Relaxation Times Post-irradiation
1100
900
Ii)
..��..
714
(I)
700
E
j::
c:
500
0
.. co
)(
co
Q) 300
D::
100
-100
MRI Field Strength (Tesla)
;ig.3. Variation of
Tl and
T2 values for irradiated PVA-3C and PVA-6C
or various field strengths.
The electrical properties are summarized in Figs. 4-6.
<ig. 4 shows the conductivity spectra for the PVA-C sam-
)les before and after y-irradiation in the frequency range of
·0 Hz to 1 MHz. Fig. 4 indicates that there is a significant
lecrease in the conductivity with freeze thaw cycles for
mirradiated samples and the conductivity has significantly
ncreased with radiation. The resistance plot (Fig. 5) shows
.n equivalent change in behaviour as conductivity. The
rradiated PVA-C has lower resistance than the unirradi-
.ted case.
In Fig. 6, the permittivity is plotted versus frequency.
~ote that in this plot, the error bars fall within the data
ymbols. There is a significant decrease in the permittivity
or both PV A-3C and PVA-6C at the lower frequencies.
~or frequencies greater than approximately 1 KHz, this
lifference in the permittivity of the samples becomes less.
'. Discussion
The results indicate a change in the elastic stiffness, the
~RI and the electrical properties of PV A -C between three
~
~
>-
'5
~
;:,
"0
c
0
U
0.1
0.09
0.08
0.07
0.06
0.05
0.04
0.03
0.02
0.0~01
t::,
No Rad Plain-3C
No Rad Plain-6C
~
Rad Plain-3C
Rad Plain-6C
1~3
1' ,
Frequency (Hz)
Fig. 4. Variation of electrical conductivity for unirradiated and irradiated
PVA-3C and PVA-6C versus frequency.
350
l!! Q)
-a> 300
"E
E
.c
Q, 250
Q)
o
c
co
1i) 200
'(ii
Q)
0::
~ 150
·c
U
Q)
[ij 100
b.
No Rad Plaln-3C
v No Rad Plain-6C
... Rad Plain-3C
Rad Plain-6C
Frequency (Hz)
Fig. 5. Variation of electrical resistance for unirradiated and irradiated
PVA-3C and PVA-6C versus frequency.
8 X 10·
6
2
A
No Rad Plain-3C
No Rad Plain-6C
• Rad Plain-3C
• Rad Plain-6C
Frequency (Hz)
Fig. 6. Variation of permittivity for unirradiated and irradiated PVA-3C
and PV A-6C versus frequency.
Please cite this article in press as: M. Goharian et aI., Modifying the MRI, elastic stiffness and electrical properties ... , Nucl. Instr. and
Meth.
(200J)' doi:lO.1016/i.nimh.2007.04. 11 1
M. Goharian et al. I Nuc!. Instr. and Meth. in Phys. Res. B xxx (2007) xxx-xxx
5
and six freeze thaw cycles as well as before and after ')i-radi-
ation. For the unirradiated specimens the indentation in
the PVA-6C is significantly less than PVA-3C, indicating
that the PVA-6C specimen is harder/stiffer. The difference
is statistically significant difference
(P« 0.001). This sup-
ports the tension measurements reported by Rosiak et al.
[21J showing that PVA-6C was stiffer than PVA-3C. The
irradiated specimens for both the PVA-3C and PVA-6C
were harder/stiffer by 14.6% and 5.70% respectively, than
the unirradiated samples. The difference (pre and post radi-
ation) for the PVA-3C was confirmed statistically
(P« 0.001), but the difference for PVA-6C was not statis-
tically significant
(P> 0.05).
The increased hardness/stiffness after irradiation may be
explained by an increased crosslinking induced by the irra-
diation. Since PVA-6C has undergone more entanglements
than PVA-3C resulting from the freeze thaw process, the
smaller change in stiffness after irradiation (which was
not significant) may be attributable to fewer available sites
for crosslinking.
Both the spin-lattice
(Td and the spin-spin
(T2) relaxa-
tion times can provide information regarding the nature
and frequency of molecular motion occurring within mate-
rials. The spin-lattice relaxation time in a polymer is inde-
pendent of molecular weight and is mainly determined by
the main chain motion. In contrast, the spin-spin relaxa-
tion time decreases with increasing molecular weight and
hence viscosity [22].
T2 is normally determined by the slow-
est motions occurring within the sample [23,24]. This will
depend also upon the radiation dose, since a polymer
may be crosslinked using ionizing radiation. The number
of crosslinks is proportional to the radiation dose, rand
the efficiency of the process which is stated as a G-value,
G(X) [25]. The G-value is defined as the number of reac-
tions or events (in this case number of crosslinks) per
100 e V of absorbed energy. The initial effect of radiation
is to increase the molecular weight,
Mw,
Mw = M~/(I - c5),
(6)
where M~ is the initial molecular weight and
b is the num-
ber of crosslinks units per weight average molecule [22].
The gelation dose,
r g, is reached when c5 = 1. One crosslink
involves two crosslinked units, one on each molecule. Any
energy absorption from radiation of (0.624 x 1020
rg) eVg-1
produces
No/2Mw crosslinks per gram, where
No is Avoga-
dro's number and therefore the gelation dose is as follows
[22]:
(7)
G-values for many reactions are in the range from 0.1 to 10,
with 1 being a typical value for crosslinking pure polymers.
The spin-lattice relaxation in a polymer has been
described by a single exponential relaxation time. In the
unirradiated and low molecular weight case, because of
the great flexibility of the polymer chain and because pro-
tons are in rapid motion, the spin-spin relaxation decay
has a single exponential behaviour (fast exchange). For
the irradiated polymer, particularly above the gelation
dose, the spin-spin relaxation decay can be described by
a double exponential function of the form [21,22J
A(t)/A(O) =
!e-t
/ T2L + (1 -
!)e-t
/ T2S
,
(8)
where
A(t)/ A(O) represent the normalized ( at
t = 0) signal
amplitude and
T2L and
T 2S refer to the long and short
relaxation times respectively and
f is the fraction of protons
relaxing with
T2L. As the dose increases
f decreases and
more of the signal decay comes from the short component.
The dose at that point were the short component starts to
decrease depends on the initial molecular weight of the
polymer [25]. The crosslinked network can decrease the
T2 which reflects the more restricted nature of translational
chain motion of the network structure. So the spin-spin
relaxation is sensitive to the crosslinked network structure.
At the highest dose,
Tl [25], which is most sensitive to
molecular motions occurring at frequencies near to the res-
onance frequency can be reduced. At higher fields a larger
decrease in
Tl relaxation time is expected which was indeed
observed in MRI measurements at 2 T and 3 T.
Note that at I T, the
T2 data fit a single exponential well,
indicating that if two water pools (one hydration layer and
one free mobile water) are present, that there is not suffi-
cient resolution in the time points to separate these compo-
nents. At 2 T however the data was better fit by a bi-
exponential behaviour. A higher
T2 component (300 fiS)
is postulated to be due to more mobile water, whereas
the second component (approx 50 IDS) is due to water asso-
ciated with the PYA-C. At 3 T the
T2 data fit a single expo-
nential relaxation yielding a
T2 of approximately 50 ms.
This fits a single exponential probably because the mobile
water component has a larger
T2 at 3 T and was not ade-
quately resolved at the lower TE values.
T2 changed little following radiation except at 1 T. It is
postulated that this is resulting from the I T measurement
arising from an average of a shorter (50 ms) water compo-
nent associated with the PVA-C with a longer (150 ms)
more mobile water component. As discussed above, if the
signal corresponds to a recovery such as that indicated in
Eq. (8), then radiation will not modify the
T2's of these
components, but rather the weightings of these compo-
nents. Thus at I T following radiation, a
T2 decrease from
96 ms to 70 ms for the 3C sample for example indicates an
increase in the weighting of the short component as
expected.
The data for conductivity shows that the conductivity
increases with ')i-radiation. It is speculated that this is due
to radiation induced polymeric crosslinking. This cross-
linking reduces hydrogen bonding between the PV A and
water, creating energetic free electrons, ions and radicals,
which are able to migrate through the polymer, eventually
changing the electrical conductivity. A clear description for
conduction mechanism of disordered polymers has not
been adapted to the data due to complexity of their struc-
ture [26].
Please cite this article in press as: M. Goharian et aI., Modifying the MRI, elastic sti.ffness and electrical properties ... , Nucl. Instr. and
Meth. B (2007), doi:IO.1016/j.nimb.2007.04.111
M. Goharian et at. I Nuc!. Instr. and Meth. in Phys. Res. B xxx (2007) xxx-xxx
Conclusion
The elastic stiffness, MRI relaxation times and electrical
operties have been measured for unirradiated and irradi-
~d PV A-C samples (three or six freeze thaw cycles). All of
:!se techniques illustrated clear differences due to radia-
In exposure on the PYA cryogel polymer. Irradiation of
T A-C samples increased their elastic stiffness by 14.6%
r PVA-3C and by 5.7% for PVA-6C. The
T} relaxation
lles were decreased by irradiating the samples. The
T}
PVA-C that had undergone three freeze thaw cycles
creased with irradiation by 10%, 25% and 35% at 1 T,
~9 T and 3 T respectively. The
Tl of PVA-C that had
!dergone six freeze thaw cycles decreased with irradiation
18%, 15% and 11% at 1 T, 1.89 T and 3 T respectively.
le
T2 of PV A-C decreased with irradiation only at 1 T,
Iwever this change is hypothesized to be due to the inter-
tion of two spin pools in the gel.
Irradiation increased the conductivity of PVA-C by up
50% for 3C and 75% for 6C at frequencies greater than
K..Hz. This can be explained from the basis that radiation
juces crosslinking and therefore ions and free radical are
rmed which can be trapped in the bulk of material. The
rmittivity was also decreased with irradiation up to
% for 3C and 35% for 6C at frequencies less than 1 KHz.
:knowledgements
We acknowledge funding support from an NSERC
ant awarded to G.R. Moran and we thank Multi-Mag-
tics Incorporated (MMI) for the use of MRI time on
e 2 T and 3 T MRI facilities in London, Ontario. We
ank Mr. W. Wells at the Integrated Manufacturing Tech-
,logies Institute who prepared the PYA samples. We are
ry grateful to the staff of the McMaster University Nu-
clear Reactor for the use of the gamma cell. We particu-
larly thank Rob Pasuta for irradiating the samples, for
his time and for his insight.
References
[1] Y. Mori et ai., J. Mater. Sci. 32 (1997) 491.
[2] K.J.M. Surry et a!., Med. Phys. BioI. 49 (2004) 5529.
[3] S.R. Stauffer, N.A. Peppas, Polymer 33 (1992) 3932.
[4] W. Wan et a!., J. Biorned. Mater. Res. (Appl. Biomater.) 63 (2002)
854.
[5] G. Campbell et aI., Technical Report IMTI-CTR-077, NRC, 2003.
[6] P.D. Garrett, D.T. Grubb, J. Polym. Sci. Part B: Polym. Phys. 26
(1988) 2509.
[7] A. Shehap et aI., J. Appl. Polym. Sci. 68 (1998) 687.
[8] L.A. Lukas et aI., Magn. Reson. Med. 46 (2001) 1006.
[9] K.C. Chu, B.K. Rutt, Magn. Reson. Med. 37 (1997) 314.
[10] I. Mano et aI., Magn. Reson. Med. 3 (1986) 921.
[11] J. Wtorek et aI., Ann. N.Y. Acad. Sci. 873 (1999) 520.
[12] K. Chu et al., Phys. Med. BioI. 45 (2000) 955.
[l3] K.J.M. Surry et aI., Med. Image Anal. 6 (3) (2002) 301.
[14] M.A. Khaled et al., Egypt. J. Sol. 26 (1) (2003).
[15] B.J. Lambert, J.M. Hansen, Radiat. Phys. Chem. 52 (1998) 11.
[16] J. Bray, E. Merrill, J. Appl. Polym. Sci. 17 (1973) 3779.
[17] N. Peppas, E. Merrill, J. Polym. Sci.: Polym. Chem. Ed. 14 (1976)
441.
[18] C. Zainuddin et al., J. Biomater. Sci. Polym. Ed. 13 (9) (2002) 1007.
[19] A. Danno, J. Phys. Soc. Jpn. 13 (7) (1958) 722.
[20] ASTM D 1415, Standard Test Method for Rubber Property,
International Hardness, ASTM International, West Conshohocken,
PA,USA.
[21] J. Rosiak et aI., Rad. Phys. Chern. 32 (1988) 793.
[22] A. Charlesby, R. Folland, Rad. Phys. Chern. 15 (1980) 393.
[23] R.N. Ibbett, NMR Spectroscopy of Polymers, Chapman& Hall, 1993.
[24] P .A. Mirau, A Practical Guide to Understanding of the NMR of
Polymers, John Wiley & Sons, 2005.
[25] R. Folland, A. Charlesby, Rad. Phys. Chern. 10 (1977) 6l.
[26] J. Berkowitch et al., J. Polym. Sci. 25 (111) (1957) 490.
[27] Moran et al. Magnetic Resonance in Medicine, submitted for
publication.
-lease cite this article in press as: M. Goharian et at., Modifying the MRI, elastic stiffness and electrical properties ... , Nucl. Instr. and
tleth. B (2007), doi:lO.1016/j.nimb.2007.04.111
Stephanie AIviS/IOPP
22106/2007 08: 16
To PermissionsJlOPP@IOPP
cc
Stephanie Alvis
Publishing Administrator
bee
Subject Fw: Permission letter
-- Forwarded by Stephanie
Alvis/IOPP on 221061200708:15--
"rnehran-
<goharim@mcmaster.ca>
To <pmea@iop.org>
cc
22106/2007 05:24
Subject Permission lelter
Dear Stephanie Alvisj
I am completing a Ph.D. thesis at McMaster University, Hamilton, Canada
entitled r 'New hardware and software design for three dimensional
Electrical
Impedance Tomography' '. I would like your permission to reprint the
following journal articles in my thesis.
l-Mehran Goharian, Aravinthan Jegatheesan. and Gerald R. Moran, "Dogleg
trust-region application in electrical impedance tomography", Physiol.
Meas. 28 (2007) 555-572.
2-Goharian M, Bruwer M J, Jegatheesan A, Moran G R, and MacGregor J F I 'A
Novel Approach for BIT Regularization Via Spatial and Spectral Principal
Component Analysis I I
in press. (PMEA/244261/PAP/161553)
Please note that I am the principle author of these works which represents
research performed during my graduate studies. I am also requesting that
you
grant me an irrevocable, non-exclusive license to McMaster University and
the National Library of Canada to reproduce this material as part of my
thesis. Proper acknowledgement of your copyright of the reprinted material
will be given in the thesis.
If these arrangements meet your approval, would you please provide me with
an appropriate letter for this request?
Sincerely,
Mehran Goharian
PERMISSION TO REPRODUCE AS REQUESTED
JS GIVEN PROVIDED THAT
(a)
tbo ssn6ant eftR8 alfl~effB) is obtained
(b) t~e source of the material including author/editor,
tItle, date and publisher is acknowledged.@
rop Publishing Limited
Dirac House
Temple Back
~:/;,
(I
~;!~~~L
~ .ff.~"'/
... ~.v.~ ~.:::! •• _::::C
~
On_
..
Da
Rlgh~&;P;-;;j~~' '------
~r~ ~~ ~
c<J ~ "~~~ ~6
_ _ al> .ak-
'WWW.~"I'.~~~ ~
~'>L. .;~\:.t..~c.!.. M, k...s...
~~
oZ:;;
ov.... ~ ~~l...e ~
..l.Lo -D~~rh,."....&.i~
,rn __ '_ ..rl, ~fto..,...J.LoA~
From: Rights and Permissions (ElS) [mailto:Permissions@elsevier.com]
Sent: Tuesday, July 17, 2007 7: 12 AM
To: mehran
Su bject: RE: Req uest
Dear Dr Goharian
We hereby grant you permission to reproduce the material detailed below at no charge in your
thesis subject to the following conditions:
1. If any part of the material to be used (for example, figures) has appeared in our publication
with credit or acknowledgement to another source, permission must also be sought from that
source. If such permission is not obtained then that material may not be included in your
publication/copies.
2. Suitable acknowledgement to the source must be made, either as a footnote or in a reference
list at the end of your publication, as follows:
"Reprinted from Publication title, Vol number, Author(s), Tit]e of article, Pages No., Copyright
Owner, Copyright Elsevier (or appropriate Society narne)(Year)"
3. Your thesis may be submitted to your institution in either print or electronic form.
4. Reproduction of this material is confined to the purpose for which permission is hereby given.
S. This permission is granted for non-exclusive English rights only. For other languages please
reapply separately for each one required.
6. This inc1udes permission for the Library and Archives of Canada to supply single copies, on
demand, of the complete thesis. Should your thesis be published commercially, please reapply
for permission.
Yours sincerely
Marion Moss
Senior Rights Assistant
Elsevier Ltd
The Boulevard
Langford Lane
Kidlington
Oxford
OXSIGB
Tcl:+441865843280
Fax: +44 1865 853333
E-mail: m.moss@elsevier.com
-----Original Message-----
From: mehran [mailto:goharim@mcmaster.ca]
Sent: 16 July 2007 18:18
To: Rights and Permissions (ELS)
Subject: Request
Dear Clare;
I sent an email dated June 22 regarding a permission to reprint the following journal article in my
thesis.
Goharian M, Moran G R,Wilson K, Seymour C, Jegatheesan A, Hill M, Thompson R T, and
Campbell G, "Modifying the MRI, Elastic Stiffness and Electrical Properties of Polyvinyl Alcohol
Cryogel Using Irradiation", Nucl. Instr.and Meth. in Phys. Res. B, ( in press).
As I am completing a Ph.D. thesis at McMaster University, Hamilton, Canada entitled "New
hardware and software design for three dimensional Electrical Impedance Tomography". I need
to get the permission letter to be able to submit my thesis. Would it be possible to process my
request ASAP?
I am also requesting that you grant me an irrevocable,
non-exclusive license to McMaster
University and the National Library of Canada to reproduce this material as part of my thesis.
Proper acknowledgement of your copyright of the reprinted material will be given in the thesis.
Sincerely,
Mehran Goharian
I FM BE Proceedings
ct~~ IFMBE
IFMBE CONFERENCE
13th International Conference on Electrical Bio-Impedance & 8th Conference on
Electrical Impedance Tomography
August 29 - September 2,2007, Graz, Austria
COPYRIGHT FORM
Please return this signed form to the Conference Secretariat by fax at: +43-316-873-7890
(one signed copyright form is f~quired for one paper)
Please type or print in block letters.
Co-authors: M
Sde.iio1an,· #
A. Je(jd~eesan
J aad at:. ml2211
ID and title # 92-66
5" - i~{,( 0.Yi2a-tjC,~1 () f
E<T Pr;;biel1-'1 .
of the paper:
i< Sil"ft -(list ret}.'Ct fI1
S;vdo Pro lier?q O1e-iJzoc:i
COPYRIGHT TRANSFER
The undersigned hereby assigns the rights of publishing the above work to the International
Federation for Medical and Biological Engineering (IFMBE). The work, if accepted for
publication, will be published in the IFMBE Proceedings Series. The undersigned (on behalf
of all authors) hereby represents and warrants that the work is original and that he/she/they
is/are the author(s) of the work.
Authors retain all proprietary rights in any process, procedure, or article of manufacture
described in the work. Authors may reproduce or authorize others to reproduce the above
work, material extracted verbatim from the above work, or derivative works for the author's
personal use or for company use without a request for permission from the IFMBE.
IMPORTANT! Submission of papers implies the author's willingness to register at the
conference and present the paper. After a paper is accepted for presentation, one of the
authors must complete a registration form and pay the appropriate fee before the paper can be
published in the Proceedings. One author's registration will guarantee publication of one
accepted paper.
Date
Hart,! 2-5,
2..c;P 7:. Sjgnature_+!i ____ ·--\;;c1or,v:!,~~.r-<...a~L,'l__ _______ _