Home > New Hardware and Software Design for Electrical Impedance Tomography

New Hardware and Software Design for Electrical Impedance Tomography

Page 1
NEW HARDWARE AND SOFTWARE DESIGN FOR ELECTRICAL IMPEDANCE TOMOGRAPHY

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

Page 28
jCOS(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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

Page 46
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)

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

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

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

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

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

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

Page 53
566 M Goharian et al
z=1.20 z=1.0 z=0.80
:r
,
'r
': -- t

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

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

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

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

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

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

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

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

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

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

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

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

Page 65
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)

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

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

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

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

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

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

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

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

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

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

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

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

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

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

Page 80
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 = 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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

Page 96
+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

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

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

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

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

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

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

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

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

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

Page 106
(a) z=8.0 z=5.0 z=4.0 z=3.0 z=2.0 z=1 .0 (b)
93

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

Page 108
(a)
z=8.0 z=s.O z=4.0 z=3.0 z=2.0 z=1.0
(b) 95

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

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

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

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

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

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

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

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

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

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

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

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

Page 121
z=8.0 z=5.0 z=4.0 z=3.0 z=2.0 z=1 .0
(a) (b)
108

Page 122
@
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

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

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

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

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

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

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

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

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

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

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

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

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

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

Page 136
• 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

Page 137
• 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

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

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

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

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

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

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

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

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

Page 146
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~

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

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

Page 149
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__ _______ _

Set Home | Add to Favorites

All Rights Reserved Powered by Free Document Search and Download

Copyright © 2011
This site does not host pdf,doc,ppt,xls,rtf,txt files all document are the property of their respective owners. complaint#nuokui.com
TOP