UFDC Home  myUFDC Home  Help 



Full Text  
PAGE 1 1 COMPENSATION OF SHAPE CHANGE AR TIFACTS AND SPATIALLYVARIANT IMAGE RECONSTRUCTION PROBLEMS IN ELECTRICAL IMPEDANCE TOMOGRAPHY By SUNGHO OH A DISSERTATION PRESENTED TO THE GRADUATE SCHOOL OF THE UNIVERSITY OF FLOR IDA IN PARTIAL FULFILLMENT OF THE REQUIREMENTS FOR THE DEGREE OF DOCTOR OF PHILOSOPHY UNIVERSITY OF FLORIDA 2009 PAGE 2 2 2009 Sungho Oh PAGE 3 3 To my family in Korea PAGE 4 4 ACKNOWLEDGMENTS To begin with, Dr. Rosalind Sadleir de serves all my gratitude for giving me the researchopportunity in the field of El ectrical Impedance Tomography. If it were not for our casual meeting in the late summer of year 2004, all this work would not have existed. I feel that I have never expressed my gratitude towards her in a proper way. Now I would like to take this opportunity to say that I am ge nuinely thankful for the enduri ng support she provided with me and my research. My transition from electrical engineering to biomedical engineering would have been impossible without Dr. Hans van Oostrom. Hi s presence has been tremendously helpful throughout my graduate studies. I also want to mention that hi s current role as a graduate coordinator is highly appreciated. With Dr. Huabei Jiang in my dissertation committee, I always tried to think outside the box. Fundamental mathematics in image reconstructi on problems are essentially the same. I actively searched for relevant literatures in areas other than Electrical Impedance Tomography, hoping that the methods developed in my resear ch could be useful in other areas. Exchanging ideas with people in other discip lines is one of the rewarding aspects of biomedical research. My extern al committee member Dr. Charlene Krueger gave me wonderful opportunities to work in her fetal/ prenatal research projects. I tr uly enjoyed working with her in the General Clinical Research Center and Neonatal Intensive Ca re Units at Shands hospital. From Dr. Paul Carney I learned that there is nothing wrong with holding on to my dreams. Dr. Mathew Saxonhouse provided me with an inspir ational opportunity in the Neonatal Intensive Care Unit. I owe what I am to my family in Seoul, Kor ea. Life alone and away from home made me realize how much I love my fa mily. I love my parents who put up with their most stubborn and PAGE 5 5 quirky child. I love my brother and sister who now have become big pillars of our family. I also want to express gratitude to Alrajhi family in Florida for their friendship and for being there whenever I felt like there was a hole in the chest. As far as I know, words are the lowest medi um in teaching. My best guiding light always has been the lives of my grandparents and ancestors I hereby swear that I will always do my best to meet your expectations and not to become a disgrace to you as your son, your student, and your friend. PAGE 6 6 TABLE OF CONTENTS page ACKNOWLEDGMENTS...............................................................................................................4 LIST OF TABLES................................................................................................................. ..........9 LIST OF FIGURES.......................................................................................................................10 ABSTRACT...................................................................................................................................13 CHAPTER 1 INTRODUCTION..................................................................................................................16 Electrical Impedance Tomography......................................................................................... 16 Clinical Applications.......................................................................................................... ....17 Rationale..........................................................................................................................17 Advantages......................................................................................................................18 Continuous Patient Monitoring....................................................................................... 19 Quantitative EIT..............................................................................................................19 Limitations.................................................................................................................... ...20 About This Dissertation..........................................................................................................21 Main Goal........................................................................................................................21 Organization of the Dissertation...................................................................................... 21 2 FORWARD PROBLEM IN EIT............................................................................................23 Electrode Layout.....................................................................................................................23 The Full Array.................................................................................................................23 Hemiarray........................................................................................................................23 Electrode Arrays in 3D EIT............................................................................................. 24 Measurement Method.............................................................................................................25 Adjacent Electrode Configuration................................................................................... 26 Opposite Electrode Configuration................................................................................... 26 The Forward Problem............................................................................................................ .27 3 IMAGE RECONSTRUCTION IN EIT.................................................................................. 30 Inverse Problem................................................................................................................ ......30 Time Difference Imaging................................................................................................30 The Sensitivity Matrix.....................................................................................................31 Inverse Model.................................................................................................................. 32 Regularized Inverse............................................................................................................ ....33 Illposed Nature of the Inverse Problem.........................................................................33 Singular Value Decomposition........................................................................................ 33 Condition Number........................................................................................................... 35 PAGE 7 7 Singular Images............................................................................................................... 36 Truncated Singular Value Decomposition......................................................................37 Weighted Minimum Norm Method................................................................................. 38 Tikhonov Regularization.................................................................................................40 Lcurve............................................................................................................................41 Quantity Index................................................................................................................. 42 4 ESTIMATING ELLIPTICAL SHAPE CHANGE OF AN OBJECT BOUNDARY USING JOUKOWSKI TRAN SFORMATION IN EIT......................................................... 43 Background.............................................................................................................................43 Method....................................................................................................................................45 The Joukowski Transformation and Elliptical Boundary Shape Change........................ 46 Boundary Shape Estimation............................................................................................ 51 Computer Si mulation....................................................................................................... 53 Phantom Experiment....................................................................................................... 54 Results.....................................................................................................................................58 Simulated Data................................................................................................................ 58 Phantom Data.................................................................................................................. 60 Discussion and Conclusion.....................................................................................................61 5 NORMALIZATION OF SPATIALLYVAR IANT IMAGE RECONSTR UCTION PROBLEM IN EIT USING SYST EM BLURRING PROPERTIES..................................... 64 Background.............................................................................................................................64 Method....................................................................................................................................65 Point Spread Function and the Blur Matrix..................................................................... 65 Normalizing Terms..........................................................................................................67 Pixelwise Scaling........................................................................................................... 67 Weighted Pseudoinverse................................................................................................68 WMNM Normalization................................................................................................... 68 Normalized Sensitivity Matrices..................................................................................... 69 Reconstruction Matrix after Normalization ..................................................................... 69 Computer Si mulation....................................................................................................... 70 Phantom Experiment....................................................................................................... 71 Measuremen t System....................................................................................................... 72 Results.....................................................................................................................................72 Ideal Reconstructions...................................................................................................... 72 Simulated Data................................................................................................................ 73 Phantom Data.................................................................................................................. 74 Discussion...............................................................................................................................79 Conclusion..............................................................................................................................82 PAGE 8 8 6 IMAGE RECONSTRUCTION BY WEIGHTED PSEUDOINVERSION METHOD IN HEMIARRAY EIT................................................................................................................. 83 Introduction................................................................................................................... ..........83 Method....................................................................................................................................83 Weighted PseudoInversion.............................................................................................83 Singular Value Decomposition........................................................................................ 84 Computer Si mulation....................................................................................................... 84 Results.....................................................................................................................................86 Noiseless Simulated Measurement.................................................................................. 86 Noisy Simulated Measurement........................................................................................89 Localization of Anomaly.................................................................................................89 Discussion and Conclusion.....................................................................................................89 7 FUTURE WORK.................................................................................................................... 92 Estimation of Shape Change in 3D......................................................................................... 92 Normalization Method in EIT Image Reconstruction............................................................ 92 Test for Anomalies of Different Conductivities.............................................................. 92 Test for Multiple Anomalies........................................................................................... 92 Normalization in 3D........................................................................................................93 APPENDIX A MATHEMATICAL FORMULATION.................................................................................. 94 Least Squares Solution......................................................................................................... ..94 TSVD Solution (Full Rank)....................................................................................................94 Weighted Minimum Norm Method........................................................................................ 95 Tikhonov Regularization........................................................................................................95 Joukowski Transformation of a Unit Circle........................................................................... 96 B MATLAB/COMSOL SCRIPTS............................................................................................. 98 Solve.m...................................................................................................................................98 Solve_batch.m.................................................................................................................. ....102 Measure.m...................................................................................................................... ......103 Measure_all_8.m..................................................................................................................103 Check_sol.m.................................................................................................................... .....104 Sensitivity.m.................................................................................................................. .......104 Reconstruct.m.......................................................................................................................107 LIST OF REFERENCES.............................................................................................................108 BIOGRAPHICAL SKETCH.......................................................................................................116 PAGE 9 9 LIST OF TABLES Table page 11 Conductivity values (Sm1) of mammalian tissues............................................................ 18 21 The complete sequence of measurements using 8electrode adjacent electrode configuration.................................................................................................................. ....26 22 The complete sequence of measuremen ts using 8electrode opposite electrode configuration.................................................................................................................. ....27 31 Properties of the sensitivity matrices (S ) of 8electrode full array and hemi array topologies...........................................................................................................................36 41 Estimated c (percentile) from the m easurements made on cylindrical and elliptical phantoms I and II.............................................................................................................. .62 51 Condition numbers of the nor ma lized system matrices..................................................... 69 52 Index assignment for va rious norma lization methods....................................................... 70 53 RSD of QI (511) calculated from the id ea l measurement reconstructions (at full rank, k = 20) of eight electrode cases................................................................................. 72 54 Maximum deviation error () of the relative QI, calculated from the simulated measurement reconstructions ( k = 16) of eight electrode cases......................................... 73 55 Maximum deviation error of the relative QI (), calculated from phantom reconstructions of eight electrode full array and hemiarray topologies............................. 74 61 Condition numbers of sensitivity matrices with and without weighting........................... 90 PAGE 10 10 LIST OF FIGURES Figure page 11 Schematic of a typical EIT system.....................................................................................16 21 Illustration of 8electrode full array applied on the boundary of a disk object.................. 23 22 Illustration of 8electrode hemiarray applied on the boundary of a disk object................ 24 23 Mesh generation for a 2D disk model with 8electrode full array and a disk anomaly. .... 28 24 Forward solution for the model in Figure 23.................................................................... 29 31 An inverse model composed of 344 elements................................................................... 32 32 Singular value spread of sensitivity matrices for 8electrode full array and hemi array.... 34 33 Singular images of the full array se nsitivity matrix (8electrodes).................................... 36 35 TSVD reconstructions for FEM mode ls containing a single anomaly.. ............................ 38 36 The l curve example in Tikhonov regularization.. ............................................................. 41 41 Extent of shape change..................................................................................................... .44 42 A 2D disk model after three different degrees of elliptical shape change (c = 0.1, 0.2 and 0.3) under the Joukowski transforma tion.................................................................... 46 43 Inverse Joukowski transformation of an elliptical doma in with 16 equidistant boundary electrodes...........................................................................................................47 44 Measurements made on the disk model............................................................................. 48 45 Reconstruction image from the measurem ent on a disk object containing an anomaly at the center. ................................................................................................................. ......49 46 Absolute values of maximum and mi nimum of reconstruction amplitude ( ) for various degrees of shape change c (%)..............................................................................50 47 QI vs. the degrees of shape change ( Joukowski parameter : c (%))................................... 50 48 Summary of our method to estimate boundary shape........................................................ 51 49 The boundary elements that were used for the template (fil led with color)......................52 410 Template vs Joukowski parameter ( c = 0 ~ 50%)............................................................. 52 PAGE 11 11 411 Elliptical phantoms simulating two different degrees of shape change ( c = 0.1 and 0.2)........................................................................................................................... ..........55 412 Prototype KHU Mark1 with 16 channels applied to a cyli ndrical phantom...................... 56 413 A malleable cylindrical phantom with 16 bar electrordes................................................. 57 414 Reconstructed images for a single anomaly on five different locations. ........................... 58 415 Cost function curve for estima ting the amount of shape change. ...................................... 59 416 Estimation errorbar plot for various amounts of shape change ( c ). Measurements were noisy with additiv e W GN of 20dB SNR................................................................... 59 417 Reference measurements from homogene ous me dium for four different amounts of shape change......................................................................................................................60 418 Elliptical boundary shape estimation by the Joukowski parameter c (%)......................... 61 419 Area of the domains for the elliptical shape change of various degrees ( c )...................... 63 51 Overview of the ePack system........................................................................................... 71 52 Eight electrode full array images ( k = 17) reconstructed from phantom measurem ents made for a bloodlike anomaly (100mL) in four different locations on the negative axis (rA/rP = 0, 0.25, 0.5, and 0.75)................................................................................. 75 53 Eight electrode hemiarray images ( k = 16) and relative QI curves reconstructed from a bloodlike anom aly (100mL) using methods A1, B1, B4 and C5.................................. 78 54 Radial curves of relati ve QI calculated using TSVD (A1) and norma lization method (C5)....................................................................................................................................80 61 Singular images of se nsitivities ma trices (S and SP )......................................................... 85 62 Comparison of singular value spectra S SP and SW in eightelectrode hemi array case.....................................................................................................................................86 63 Images of a single anomaly at 15 different locations reconstructed from noiseless m easurements using A) TS VD, B) WMNM and C) the WPI in eightelectrode hemiarray case ( k =20).......................................................................................................87 64 Images of a single anomaly at 15 di fferent locations reconstructed from noisy m easurements (AWGN of 30dB) using A) TSVD, B) WMNM and C) the WPI in eightelectrode hemiarray case ( k =20)..............................................................................88 65 Localization of a single anomaly at various positions from simulated noisy m easurements (AWGN of 40dB SNR).............................................................................. 90 PAGE 12 12 66 Localization of a single anomaly at various positions from simulated noisy m easurements (AWGN of 30dB SNR).............................................................................. 91 PAGE 13 13 Abstract of Dissertation Pres ented to the Graduate School of the University of Florida in Partial Fulfillment of the Requirements for the Degree of Doctor of Philosophy COMPENSATION OF SHAPE CHANGE AR TIFACTS AND SPATIALLYVARIANT IMAGE RECONSTRUCTION PROBLEMS IN ELECTRICAL IMPEDANCE TOMOGRAPHY By Sungho Oh May 2009 Chair: Rosalind J. Sadleir Major: Biomedical Engineering Electrical Impedance Tomography (EIT) is an imaging modality in which electrical conductivity within the object is estimated from surface voltage measurements. Conductivity images in EIT can provide clinically signifi cant information, since conductivity changes are closely related with physiological changes inside the body. Clinical applications of EIT are noninvasive, costeffective and simple to apply to patients. Despite the above advantages, performance of EIT systems is affected by th e uncertainties pertai ning to patients body shape change and spatial variability of the image reconstruction problem. In this research, we identify these uncertainties to be the major sources of reconstruction errors. In EIT image reconstruction, the surface shape of an image object is often assumed to be known. In clinical environments, shape informa tion is not always available. Discrepancies between the assumed and actual shapes can result in errors that may have clinical significance. We suggest an algorithm that estimates domain shap es for the use in 2D EIT. We investigated elliptical boundary distortions of a unit disk object as changes from circular to elliptical geometries, defined using the J oukowski transformation. Boundary sh apes of a real domain were then estimated as ellipses after investigating th e spatial characteristics of image artifacts caused PAGE 14 14 by shape changes. Our method was tested with boundary voltage measurements obtained using a full array electrode layout from elliptical simulation and phantom models containing a small disk anomaly at various positions. We found that the proposed method c ould estimate elliptical shape changes with relatively small error. The EIT image reconstruction problem is spatia llyvariant, meaning that the same anomaly placed at different locations within an imag e plane may produce different reconstruction signatures. Correcting errors du e to this spatial variability should improve reconstruction accuracy. We present methods to normalize the spatiallyvariant image reconstruction problem by equalizing the system Point Spread Function (PSF). In order to equalize PSF, we used blurring properties of the system derived from the sensitivity matrix. We compared three mathematical schemes: PixelWise Scaling (PWS), Weighted PseduoInversion (WPI) and Weighted Minimum Norm Method (WMNM) to normalize reconstructions. The Quantity Index (QI), defined as the integral of pixel values of an EIT conductivity image, was considered in investigating spatial variability. The QI values along with reconstructed images are presented for cases of 2D full array and hemiarray electrode topologies. We found that a less spatiallyvariant QI could be obtained by appl ying normalization methods to conventional regularized reconstruction methods such as Truncated Si ngular Value Decomposition (TSVD) and WMNM. The normalization methods were tested with boundary voltage measurements obtained from simulation disk models containing a smaller disk anomaly, and cylindrical phantom models with anomalies of various volumes placed at various locations within the electrode plane. For anomalies of the same volume, QI error caused by spatial variability was reduced the most among the tested methods when WMNM normaliz ation was applied to WMNM regularized reconstructions for both hemiarray and full array cases. PAGE 15 15 The use of the blurring properties was further investigated in hemiarray EIT, where the electrodes cover only one half of the object boundary. Boundary measurements are relatively not sensitive to the conductivity anomaly that lies far away from electrodes, and the anomaly may be invisible or undetected in the images reconstructed using conventional methods. We propose a WPI method to enhance sensitivity in the region distant from the electrodes. The method was tested with data obtained from a 2D circular object. A smaller disk anomaly was varied in location within the object. The WPI method detect ed anomalies with relatively small errors for the hemiarray case. PAGE 16 16 CHAPTER 1 INTRODUCTION Electrical Impedance Tomography Electrical Impedance Tomography (EIT) is a noninvasive im ag ing modality that displays conductivity distribution with in an object. In EIT, a set of el ectrodes is positioned on an objects surface, input currents are a pplied, and the resulting surf ace voltages are measured. The measurements are dependent on the objects in ternal conductivity distribution. Given the knowledge of input current patterns and the objec ts shape, crosssectio nal conductivity images can be reconstructed from th e measurements (Figure 11). Figure 11. Schematic of a typical EIT system. PAGE 17 17 EIT conductivity images can be useful in many areas, mainly clinical, industrial and geophysical. An emerging area is the molecular medicine where realtime imaging of cell electroporation may become f easible using EIT (Davalos et al. 2004). Some examples of industrial and geophysical a pplications include: Imaging of fluid flows in pi pes (Dickin and Wang 1996), Imaging material distribution within the process vessels (Heikk inen et al. 2006), Crack detection (Alessan drini and Rondi 1998), Detecting air bubbles in pipe lines (Ljaz et al. 2008), Detecting the free surface be tween the liquid and air (T ossavainen et al. 2004), Imaging heterogeneities of sandy sa mples (Borsic A et al. 2005), Groundwater studies (Nobes 1996), Detecting corrosion faults in metallic plates (Vilhunen et al. 2002). In this dissertation, we focused on the clinical aspect of EIT. The main clinical subjects that EIT have demonstrated applicability in are: Gastrointestinal tract (gastric emptying and acid secretion), Lungs (regional ventilation), Brain, Pelvis (pelvic blood volume) (H older 2005, Morucci and Rigaud 1996). Now let us discuss clinical applic ations of EIT in more details. Clinical Applications Rationale The human body is m ade of tissues well contrast ed in electrical c onductivity (Barber and Brown 1984). Therefore, an imaging modality such as EIT producing distributional maps of conductivity within th e body may find useful clini cal applications (Dijkstra et al. 1993). Table 11 shows that the conductivity values of mammalian tissues vary widely with tissue types, and tissues such as blood and cerebrospinal fluid have smaller conductivity values compared to fat and bone tissues. PAGE 18 18 Table 11. Conductivity values (Sm1) of mammalian tissues (Barber and Brown 1984) Tissue Conductivity Cerebrospinal fluid 1.54 Plasma 1.52 Blood 0.17 White matter 0.15 Grey matter 0.35 Fat 0.040.05 Lung (expirationinspiration) 0.040.14 Bone 0.01 Human arm (longitudinal) 0.42 Human arm (transverse) 0.15 Skeletal muscle (longitudinal) 0.670.8 Skeletal muscle (transverse) 0.040.06 Cardiac muscle (longitudinal) 0.170.63 Cardiac muscle (transverse) 0.020.24 The conductivity characteristics of pathological tissues differ from those of normal tissues (Walker et al. 2000, Sha et al. 2002, Keshtkar and Keshtkar 2007, Shini and Rubinsky 2007). Therefore, EIT conductivity images are anticipat ed to provide useful screening methods for breast cancer (Cherepenin et al. 2001, Soni et al. 2004), prostate cancer (Halter et al. 2007) or dental caries (Pretty I A 2006). Advantages There are several advantages of EIT that ma kes it especially a ttr active to clinical application. EIT is completely noninvasive innocuous (Ghahary and Webster 1989), and simple to apply. It can easily be made por table using either a laptop computer or a PDA (Tang et al. 2006). It can be made as a handheld probe and be applied wh erever it is needed (Kao et al. 2006). Additionally, integration with wi reless technology is plausibl e (McEwan and Holder 2007). EIT systems can be developed to promptly acquire data and produce images (Barber and Seagar 1987) with relatively low cost Their imaging systems can produce realtime images and are cheaper than other imaging modalities such as MRI or CT. PAGE 19 19 Continuous Patient Monitoring The advantages me ntioned above make EIT sy stem attractive as a continuous patient monitor. The systems can be used to monitor the physiological events that cause detectable conductivity changes continuously, since they can be made to produce realtime conductivity images (Edic et al. 1995, Lionheart et al. 1997). We categorized the potential imaging areas where EIT has been considered as a con tinuous patient monitor in the following. Brain imaging: detecting intraventricular he morrhage (Murphy et al. 1987), brain function monitoring for cerebral ischemia and epile psy (Boone and Holder 1995) and stroke detection (Romsauerova et al. 2006). Lung imaging: neonatal lung imaging (Tak tak et al. 1995), monitoring pulmonary ventilation and perfusion (Campbell et al. 1994, Adler et al 1997, Frerichs 2000), pulmonary edema (Adler et al. 1995) and emphysema (Eyboglu et al. 1995). Digestive system: monitoring gastric emptying (Mangnall et al. 1991, Nour et al. 1995). Internal hemorrhage: bleeding detection in th e pelvis (Meeson et al. 1995), bleeding rates monitoring in the peritoneum (Sadleir and Fox 2001, Xu et al. 2007), Blood flow imaging (Brown et al. 1991, Vonk Noordegraaf et al. 1997) Temperature distribution: thermal monitoring during hyperthermia (Mos kowitz et al. 1995) and hyperthermia treatment (Conway 1987). For neonatal patients, applica tions in monitoring gastric emptying, pulmonary functions and intraventricular hemorrhage we re considered promising (Murphy et al 1987). Quantitative EIT The possibility of quantitative EIT was raised at an early stage of EITs development. Exem plary applications suggested by Barber ( 1990) were gastric emptying, cardiac function, pulmonary ventilation and perfus ion. Studies that followed demons trated that estimating volume of conductive anomalies was plausible for some physiological changes. These include the changes of pulmonary liquid volume (Adler and Guardo 1996, Adler et al. 1997), bleeding PAGE 20 20 volume within the pelvic bowl (Thomas et al. 1994), bleeding volume inside the ventricles (Blott et al. 2000), and bleeding rate inside the peritoneum (Sadleir and Fox 2001). An earlier work by Pomerantz et al. (1970) suggested the use of fourelectrode transthoracic impedance measurement in monitoring in trathoracic fluid volume. For patients with significant intrathoracic fluid accumulation, consistent changes were observed in baseline measurements. In this type of impedance measurement, however, localization of the internal anomalies is not feasible. Add itionally, the measurement is highl y dependent on electrode errors and on patients anthropometric parameters that are changeable during measurements (Zlochiver et al. 2007), possibly causing false positives. These limitations can be overcome in EIT by solving for the spatial conductivity distribution. The specificity of observations can be verified by inspection of reconstructed spa tial distributions, and effects of electrode contact problems and anthropometric variations on EIT recons tructions can be reduced (Soleimani et al. 2006, GmezLaberge and Adler 2006). Limitations The number of electrod es that can be placed on the body surface is limited. As a consequence, the images produced have low spatial resolution. By in creasing the number of electrodes, image quality will generally improve. Additionally, using nonstationary electrodes, one can obtain more independent measurements using the same number of electrodes (Murphy and York 2006). However, Tang et al. (2002) found that improvements were found near the object boundary rather than in the image center wh en more electrodes are considered. Therefore, having more measurements may not directly relate to the improvements of measurement sensitivity with respect to the changes in the body center. The EIT image reconstruction problem is illposed, meaning that small measurement noise can cause large image ar tifacts (Holder 2005). The problem usually has to be regularized in order PAGE 21 21 to produce meaningful images. Regularized im ages may still produce inaccurate quantity estimation due to spatialvariability of the reconstruction problem. The effect of the spatial variability is that the reconstr uction of an anomaly varies with its position. This can lead to significant errors in quantitative EIT (Adler et al. 1997, Sadleir and Fox 1998). Correcting these errors will result in more accurate reconstructio n, and after the correction reasonable accuracy may be achievable with sma ller number of electrodes. Another important factor limiting the performa nce of EIT systems is the patients body movements (i.e. breathing) and uncertainties related with the use of elect rodes (i.e. misplacement or positional changes). For convenience, we use the term shape change whose general meaning includes all the uncertainties mentioned above (Oh et al. 2008). Shape change problems are almost unavoidable in clinical practices, a nd the image artifacts are hard to remove by conventional regularization methods, producing significant errors (Adler et al. 1996). About This Dissertation Main Goal We aim to i mprove reconstruction accuracy by reducing errors caused from the above identified problems in EIT image reconstruction. Our main goal is to compensate for the effects of shape change and spatial variability in EIT image reconstruction Organization of the Dissertation This dissertation is organized as follows. In Chapter 1, we give a brief introduction about EIT and its application s. Limitations of clinical EIT are identified, which lead us to the goal statement of this work. In Chapter 2 and 3, forward and inverse problems in EIT image reconstruction are described. In Chapter 4, we describe methods to compensate for artifacts caused by shape change, and demonstrate a method to estimate domain boundaries as ellipses. The estimated shapes are found to be effective in reducing shape change artifacts. In Chapter 5, PAGE 22 22 errors due to spatial variability in the image reconstruction problem are identified. We aim to reduce these errors by normalization using the syst em blurring properties. We finish the chapter by proposing one normalization method that show s the best performance among the ones in comparison. In Chapter 6, a blurring property is found to be useful in enhancing localization of anomalies when applied as the WPI for the hemia rray case. In Chapter 7, we conclude this work by suggesting future works. PAGE 23 23 CHAPTER 2 FORWARD PROBLEM IN EIT Electrode Layout The Full Array In the full array layout, electrodes are position ed equally spaced around a 2D object boundary. This ensures current projections to be axially symmetric within the assumed object shape. Barber and Seagar (1987) showed that clin ically useful images can be reconstructed in difference EIT using this type of 2D array. In Figure 21, use of the 8electrode full array on a 2D disk model is illustrated. The angl e between two adjacent electrodes is 45o for this model. Figure 21. Illustration of 8electrode full array applied on the boundary of a disk object. Electrodes are indexed from E1 to E8. Hemiarray The novel hemi array electrode layout was motivat ed by EIT apparatus fo r the detection of intraperitoneal hemorrhage (Sadleir et al. 2008). Patients in blunt trauma accidents, who are suspected to have developed intr aperitoneal hemorrhage, may also have injured their spines. In PAGE 24 24 such a case, lifting the body in order to place elec trodes on the back can further damage the spine. Therefore, the preferred electrode locations will be the anterior surface of the body for patients in supine position, and this is where the hemiarray can be practical. In Figure 22, use of the 8electrode hemiarray is illustrated on a 2D disk model. The angle of each pair of adjacent electrodes is 22.5o, with the exception of the angle (225o) of the bottom pair (E8, E1). Figure 22. Illustration of 8electrode hemiarray applied on the boundary of a disk object. Electrodes are indexed from E1 to E8. Variable was used to represent clockwise angular displacement away from the northpole. Axes whose values were 0o, 90o, and 180o were denoted as and respectively. Although using the hemiarray electrodes can be very practical in some cases, the current projections formed by using the hemiarray are highly asymmetric within the image plane. Typically, measurement sensitivit y is comparatively low in the region, which is further away from electrodes. Electrode Arrays in 3D EIT 3D expansion of the full array can be useful when imaging objects have near cylindrical shapes. Frerichs et al. (19 99) used three layers of 16electrode full arrays placed on the chest to PAGE 25 25 image ventilation. Dehghani et al. (2005) used four layers of 16electrode full arrays to image human breast. The linear electrode array suggested by Powell et al. (1987) has an advantage in which it can be placed on the desired part of the body su rface. Kotre (1996) suggest ed the use of linear electrode arrays placed ort hogonally for subsurface imaging. When applying electrodes in EIT, the design of electrode arrays can be clinically or practically motivated. Rectangular electrode arrays were considered for breast cancer detection (Mueller et al. 1999), in which a handheld EIT probe coul d be pushed to flatten the breast. For head imaging (as in stroke cases), a modified 1020 EEG electrode layou t is a practical choice (Tidswell et al. 2001). In the neonatal skull, fontanel st ays open for a relatively long period of time. It has higher electric conductivity than the surrounding skull bone tissues. Based on this feature, an electrode placed on t op of the fontanel could be used for current injection consistently. A previous study by Sadleir and Tang (2009) appl ied this configuration to the detection of neonatal intraventricular hemorrhage. Measurement Method Once positioned, the electrodes are used to inject curren ts and measure the resulting voltages. The sequence of obtaining measurements us ing the electrodes has to be configured in a certain way. In this chapter, we discuss two current injection methods that uses adjacent and opposite electrode pairs. Based on the superposition principle, measurem ents made using all available adjacent pairs are generally seen to be sufficient. For exampl e, the measurement made using the (E1, E3) pair will, theoretically, be identical to the sum of the measurements made using the (E1, E2) and (E2, E3) pairs (Figure 21). In practice, electrodes used for current injection are not shared for voltage PAGE 26 26 measurements in order to avoid uncertainties that can arise from the changes of contact impedances in electrodes. Adjacent Electrode Configuration In the adjacent electrod e configuration, an ad jacent pair of electrode s is chosen to drive input current. A frame of input boundary voltage measurements is made using all adjacent pairs excluding the input pair. The ne xt neighboring input pair produces the next frame of measurements and so on. This procedure is repe ated for all the possible adjacent input pairs. Table 21 shows the complete sequence of measurements using the adjacent electrode configuration in 8electrode case. The measurement is conventionally configured so that the first measurement always begins from the lo west number of the electrode indexes. Table 21. The complete sequence of measur ements using 8electrode adjacent electrode configuration. Input pair Output pairs (E1, E2) (E3, E4), (E4, E5), (E5, E6), (E6, E7), (E7, E8) (E2, E3) (E4, E5), (E5, E6), (E6, E7), (E7, E8), (E8, E1) (E3, E4) (E1, E2), (E5, E6), (E6, E7), (E7, E8), (E8, E1) (E4, E5) (E1, E2), (E2, E3), (E6, E7), (E7, E8), (E8, E1) (E5, E6) (E1, E2), (E2, E3), (E3, E4), (E7, E8), (E8, E1) (E6, E7) (E1, E2), (E2, E3), (E3, E4), (E4, E5), (E8, E1) (E7, E8) (E1, E2), (E2, E3), (E3, E4), (E4, E5), (E5, E6) (E8, E1) (E2, E3), (E3, E4), (E4, E5), (E5, E6), (E6, E7) This type of measurement method yields a total of n(n3) measurements for n electrodes. Due to reciprocity, only half the measuremen ts are independent (Malmivuo and Plonsey 1995). For example, using 8electrodes we can obtain 40 measurements in total, but only 20 of them are independent. Opposite Electrode Configuration In the opposite electrode configuration, a pair of electrodes that are diametrically opposite is selected for input current in jection. Boundary voltage m easurements are made using adjacent pairs of all combinations selected from the rest of the electrodes. The opposite electrode PAGE 27 27 configuration showed improved measurement sens itivity and SNR in domain center (Avis and Barber 1994). However, it produces fewer meas urements compared to adjacent electrode configuration. When n electrode are used, it produces n(n4) measurements in total, and only half of them are independent. This indicates that localization accuracy may improve, but image resolution may not. In cerebral imaging, the pres ence of the skull prevents the current from entering the brain, significantly decreasing measurement sensitivity and SNR. The opposite current injection was suggested by Bayford et al. (1996) for cerebral imaging. Where opposite or approximately opposite current injection was us ed for 3D cerebral imaging, more accurate localization was expected (Tidswell et al. 2001, Dong et al. 2005, Romsauerova et al. 2006). Table 22 shows the complete sequence of measurements using the opposite electrode configuration in 8electrode array case. Measur ements were configured so that the first measurement always began from the lowest number of the electrode indexes. Table 22. The complete sequence of measur ements using 8electr ode opposite electrode configuration. Input pair Output pairs (E1, E5) (E2, E3), (E3, E4), (E6, E7), (E7, E8) (E2, E6) (E3, E4), (E4, E5), (E7, E8), (E8, E1) (E3, E7) (E1, E2), (E4, E5), (E5, E6), (E8, E1) (E4, E8) (E1, E2), (E2, E3), (E5, E6), (E6, E7) (E5, E1) (E2, E3), (E3, E4), (E6, E7), (E7, E8) (E6, E2) (E3, E4), (E4, E5), (E7, E8), (E8, E1) (E7, E3) (E1, E2), (E4, E5), (E5, E6), (E8, E1) (E8, E4) (E1, E2), (E2, E3), (E5, E6), (E6, E7) The Forward Problem The forward problem in EIT is to calc ulate the in ternal electric potential from the applied current density J for a known conductivity distribution within a given object The object boundary d is assumed to be smooth. When the prescribed input J is transmitted through PAGE 28 28 d, the corresponding is generated following the Lap ace equation shown in Equation 21a with the boundary condition as in Equation 21b. 0 in subject to (21a) Jnzyx ,, on d (21b) Here, n is a unit outward normal to the surface and is a function of positional variables ( x,y,z ). Ideally, forward problems can be solved analyti cally. However, analytic al solutions are not always possible for domains that have complex shape or nonuniform conductivity distribution. In such cases, the Finite Element Method (FEM) can be used to provide numerical approximations of the analytical solution. In FEM, the object domain is divided into small elements. A mesh structure composed of the tria ngular elements is displayed in Figure 23. We designed a 2D disk model containing a smaller di sk anomaly with 8elect rode full array (Comsol Multiphysics: Burlington MA). The model was disc retized to 1472 secondorder triangular finite elements. Figure 23. Mesh generation for a 2D disk model with 8electrode full array and a disk anomaly. PAGE 29 29 Then, the solution is approximated as a pol ynomial function in each element where pixel values are assumed to be constant. If the elements are sufficiently small, the finite element solution is supposed to approach the exact solution closel y. We assigned electrodes the conductivity of copper (5.997 [Sm1]). Conductivity values for background and anomaly were 1 [Sm1] and 2 [Sm1] respectively. We applied inward current density on the boundary of electrode E1 (1 [Am2]) and E2 (1 [Am2]). The rest of the boundary segments were electrically insulated. The solution (electric potential: ) of our forward problem was solved using the direct linear system solver UMFPACK (Unsymmetric Mu ltiFrontal PACKage) and displayed in Figure 24. Figure 24. Forward solution for the model in Figur e 23. Electric potential field is generated for the current density applied on electrode E1 (1 Am1) and E2 (1 Am1) in a 2D disk model (1 Sm1) containing a smaller disk anomaly (2 Sm1). PAGE 30 30 CHAPTER 3 IMAGE RECONSTRUCTION IN EIT Inverse Problem For a fixed input current pattern and an object shape, a boundary voltage measure ment V can be written as a function of conductivity distribution as in Equation 31. fV (31) The inverse problem in EIT image reconstruction is to estimate for the solution using the nonlinear relation f between and V Time Difference Imaging Reconstructing images of s mall conductivity changes from measurement changes over time is termed time difference imaging. In EIT, a sequence of difference images can be obtained in realtime, which is desirable for continuous monitoring or functional imaging applications. For a small conductivity change away from the reference conductivity distribution 0, the measurement f() can be expressed as a Taylor series as in Equation 32. !20 2 0 0 0n f f ffVn (32) By ignoring higher order terms, we obtain a linear approximation as in Equation 33. 0 0ffV (33) where f() is the reference measurement V0, and f() the Jacobian. If we move the reference measurement V0 to the left hand side of th e equation, we get Equation 34. 0 0fVVV (34) where V is the measurement change away from f() One can express f() in matrix form, which gives us a system of linear equations as in Equation 35. PAGE 31 31 S V, (35) where S is the Jacobian matrix, also termed the sensitivity matrix. Now, the nonlinear inverse problem or Equation 32 is linearized to Equa tion 35. In Equation 36, all the terms of the components in Equation 35 are laid out. N j NM M ji i N j M iV V V V 2 1 1, 1, 2,21,2 ,1 ,1 2,11,1 2 1S S S S SS SSSS, (36) where M is the total number of measurements, and N is the total number of elements (pixels) within the inverse model. In Equation 35, will be solved from measurement V for an assumed S While inversion of the linearized equation as in Equation 35 can be iterated to obtain more accurate solution, a single step inversion was shown to produce reasonable reconstructions (Cheney et al. 1990, Avis and Barber 1992). The method used in th is dissertation is also based on a single step inversion method, which will be expl ained later in this chapter. The Sensitivity Matrix For every unit elementwise anom aly, the forw ard problem can be solved to produce a set of measurements. S can be built by using these sets of measurements as the columns of S This approach necessitates solving as many forwar d problems as the total number of elements. Using sensitivity theorem or lead field theo rem by Geselowitz (1971), Murai and Kagawa (1985) could calculate S with much less computational burden. Their method requires only one forward problem solved. Each entry Sij can be calculated as in Equation 37. PAGE 32 32 j j ijdv II 0 0S (37) The i and j are measurement and element indexes respectively. Sensitivity Sij is a negative integral of the inner product of lead fields due to measurement i over a finite vo lume of element j (vj). The input and output lead fields are /I and /I respectively, for the reference conductivity distribution 0. The rows and columns of the matrix S are arranged according to the measurement index i and the element index j. Measurement index was derived from the complete measurement sequence as the one shown in Table 21.The jth column of S is socalled the ideal measurements for a unit conductivity anomaly appearing in the element j. Inverse Model We designed an inverse model which was us ed throughout this work (Figure 31). Our model is a simple structure composed of 344 elem ents. All elements are square, and their areas are the same except the ones lying on the object boundary. We calculated S for the inverse model having assumed uniform conductivity distribu tion, since accurate knowledge of internal conductivity is not easi ly available (Meeson et al. 1995). Figure 31. An inverse model composed of 344 el ements. All elements are square and have the same area excluding the ones on the object boundary. PAGE 33 33 Regularized Inverse Illposed Nature of the Inverse Problem In EIT, the total number of measurements is typically smaller than the total number of elements. The number of measurements is limite d by the number of electrodes applied. Thus, Equation 35 is typically an underdetermined problem, and S is not square. For a nonsquare matrix, the direct inverse does not exist in genera l, in which case the least squares solution can be obtained. In the least squares method, solution is estimated so as to give the best fit that minimizes the sum of squared errors of Equation 35 as in Equation 38. 2 2minS V (38) This gives us the least square s solution (39) (Appendix A). VTT1SSS (39) The inverse problem in EIT image reconstr uction is a severely illposed problem (Vauhkonen et al. 1998). As a result, small amount of noi se in the measurement can cause large image artifacts. Simple least squa res solutions as in Equation 39 will fail to produce meaningful images in EIT image reconstruction. In order to obtain reasonable images, the problem in Equation 35 has to be further regularized. Singular Value Decomposition Singular Value Decomposition (SVD) is a factor ization of a rectangular matrix that has many applications (Golub and Reinsch 1970). Here, we apply SVD to compute the pseudoinverse of matrix S. Using SVD any matrix can be represente d as a product of a unitary matrix, a diagonal matrix and another unitary matrix as in Equation 310. TUDVS (310) PAGE 34 34 where D is a diagonal matrix whose entries are termed the singular values. The rank of D is identical to the rank of S. The singular values are nonne gative numbers, arranged in a descending order. We calculated the matrices S for the adjacent electrode confi guration in the 8electrode full array and hemiarray cases. In Figure 32, the singul ar values of both matri ces are displayed as a log scale plot. Both curves decay gradually to zero indicating that our inverse problem is illposed. Figure 32. Singular value spread of sensitivity matrices for 8el ectrode full array and hemiarray. In Equation 310 U and V are unitary matrices whose inverses are identical to their transpose. The columns of U and V are termed the left and right singular vectors of S. Written in column form, U and V can be expressed as in Equation 311. Muuu21U (311a) Nvvv21 V (311a) PAGE 35 35 Here, ui (M ) and vi (M ) are ith column of U and V respectively. This column form lets us rewrite Equation 310 as in Equation 312. T N T T Mv v v uuu 2 1 21D S (312) Consequently, S is sum of the rank one matrices (uivi T) weighted by diagonal entries of D (di) as in Equation 313. N i T iiivud1S (313)Condition Number The system matrix of an illposed inverse pr oblem is badly conditioned. The condition of a matrix can be represented by the condition number. The condition number of a matrix is defined to be the ratio between its maximum and mini mum singular values as in Equation 314. min maxd d S (314) where dmax and dmin are the maximum and mini mum singular values of S respectively. Here, the condition number (S) provides a way to measure the sensitivity of S to measurement noise. If (S) is close to 1, the matrix columns are very independent. When (S) is large, the columns are nearly dependent, making the inverse problem illposed. The condition numbers in Table 21 are extr emely large, indicating that our inverse problem is severely illposed. The hemiarray matrix condition number is larger than the full array, indicating that its inverse problem is more illposed. The ranks of both matrices are 20, half the number of total meas urements, owing to the recipr ocity principle (Malmivuo and Plonsey 1995). PAGE 36 36 Table 31. Properties of th e sensitivity matrices (S) of 8electrode full array and hemiarray topologies. Full array Hemiarray Dimensions 40 40 Rank 20 20 Condition Number 2.2117 2.0218 Singular Images Columns of the right singular vector V are termed the singular images, and they form the basis of the reconstruction images (Zadehkoochak et al. 1991). These singular images can provide useful insights about how regularization determines spatial resolution in reconstruction images. Figure 33. Singular images of the full array sensitivity matrix (8electrodes). The mode number increases from left to right and top to bottom. In Figure 33, singular images of the sensitivit y matrix in 8electrode full array case are displayed. Two main observations may be made from the figure. The boundary region, which is PAGE 37 37 near the electrodes, will be described ahead of the central region, which is far away from the electrodes. Similarly, lower spatial frequency components in the image will be depicted before higher spatial frequency components. Singular images in the hemiarray case is seen to be similarly characte rized (Figure 34). In the hemiarray case, the region is furthest away from the elec trodes. It is seen that the singular images are unable to describe details of the region, indicating potential difficulties in reconstructing anomalies. Figure 34. Singular images of the hemiarray sensitivity matrix (8e lectrodes). The mode number increases from left to right and top to bottom. Truncated Singular Value Decomposition In TSVD, we attempt to get a meaningful a pproximate of the inverse solution by replacing S by a wellconditioned matrix derived from S. Based on the SVD of S in Equation 310, the pseudoinverse of S (S+) can be calculated as in Equation 315. PAGE 38 38 TUVDS (315) where D+ is a truncated inverse of D. Therefore, an estimated conductivity solution can be estimated using the pseudoinverse S+ as the reconstruction matrix. V S (316) Figure 35. TSVD reconstructi ons for FEM models containing a single anomaly. A) FEM models containing a single anomaly in 5 di fferent locations (8electrode full array). B) Corresponding images generated from TSVD reconstruction ( k = 16). C) FEM models containing a single anomaly in 5 di fferent locations (8electrode hemiarray). D) Corresponding images generated from TSVD reconstruction ( k = 16). The higher order singular values of S correspond to high spatial frequency components in the singular image. As a result of inversion, the higher order values that are small will be reciprocated and show as artifacts in reconstruction image. By way of truncation of the small values, image artifacts can be successfully s uppressed. Truncated SVD (TSVD) is such a A B C D PAGE 39 39 procedure that was proposed for discrete illpos ed problems to obtain smooth solutions (Hansen 1987). The estimated solution by TSVD can be expressed as the weighted sum of rankone matrices as in Equation 317. k i i i T iv d Vu1, (317) where k is the truncation number. The k can be chosen based on the amount of noise in measurement V (Hansen 1987). When k equals the rank of S TSVD solution is the same as the least squares solution, since S+ is equal to ST( SST)1 (Appendix A). Once the reconstruction vector is calculated, images are made by assigning estimated values to their corresponding elements (pixels) in the image plane. In Figure 35, images were produced from TSVD solutions and displayed after linear smoothing. Weighted Minimum Norm Method Sensitiv ity v alues are typically low in the regions that are further from the boundary electrodes. The FOCal Underdetermined Syst em Solver (FOCUSS) algorithm employs the strategy that power normalization of the ideal measurement should improve spatial resolution of images. It was originally proposed for solving inverse problems in Electroencephalogram (EEG) and Magnetoencephalogram (MEG) by Gorodnits ky and Rao (1997). The algorithm produces a least squares solution based on minimization of a reweighted norm of the solution (Gorodnitsky and Rao 1997). This method and its modifications were successfully adapted to EIT image reconstruction problems of various electrode to pologies: full array (Cla y and Ferree 2002, Dong et al. 2004), hemiarray (Sadleir et al. 2008), and modified 1020 syst em of EEG electrode layout (Dong et al. 2005). Clay and Ferree (2002) suggested a single step inversion of the FOCUSS algorithm in the context of EIT image reconstruction. This wa s termed the Weighted Minimum Norm Method PAGE 40 40 (WMNM) (Sadleir et al. 2008). Using WMNM, columns of S are equalized in terms of its power before pseudoinversion. The weighting terms ( wj) are defined as in Equation 318. 2/1 1 2 M i ij jSw (318) The changes made due to this equalization can be recovered by resca ling the pixels. Thus, using a diagonal weighting matrix W comprised of these terms, the final form of WMNM regularized solution can be obtaine d as in Equation 319 (Appendix A). V SWW (319)Tikhonov Regularization When a meaningful solution cannot be computed by TSVD, imposing additional constraints on the problem may resolve this problem. In the general version of Tikhonov regularization, one obtains the solution by minimizing the quantity in Equation 320. 2 2 0 2 2minarg R S V (320) where R is regularization matrix and is regularization parameter. The Newtons OneStep Error Reconstructor (NOSER) al gorithm is a well known example of this type of single step Tikhonov regularization solution with R being diag( STS ) (Cheney et al. 1990). In Equation 321, we present zeroorder Tikhonov regula rization for an identity matrix R and a zero vector 0. 2 2 2 2minarg S V (321) Larger means larger amount of regulari zation constraining the solution to have smaller norm. Algebraic form of the soluti on from Equation 321 is shown in Equation 322 (Appendix A). VIT TSSS1 (322) PAGE 41 41 Lcurve The choice of regularization param eter (or the truncation number k for TSVD) will determine the level of regularization, a nd is affected by the noise quantity in V which is generally unknown. The use of the l curve was suggested in selecting adequate regularization parameters (Hansen 1992). In the l curve, the solution norm is plotted versus the residual norm in loglog scale. The resulting curve is usua lly concave (Lshaped), hence its name. The l curve displays a tradeoff in regularization, where maxi mum curvature point is often seen as an optimal choice of a regularization parame ter (Hansen and OLeary 1993). Figure 36. The l curve example in Tikhonov regularizat ion. The corner was located at = 7.56311. PAGE 42 42 In Figure 36, l curve for Tikhonov regularization was generated from the boundary measurement made from the 8electrode FEM m odel with a central anomaly. Its corner was located at = 5.274611 using the regularization toolbox by Hansen (1994). Quantity Index Recently, es tim ation of absolute conductivity va lues has been studied using parametrical forward models, where the absolu te conductivity values producing the best fit to conductivity measurements are suggested (Brown et al. 2002, Zhang and Patterson 2005). However, design of a parametrical model can be difficult for certain body parts due to tempor al shape changes during measurement or large shape variations across patients. Here, we propose an integral of EIT images as a useful measure related with anomaly volume. The quantity defined in Equation 323 is termed the Quan tity Index (QI), which can be thought as the average of conductivity change multiplied by the number of elements N N j jjAQI1, (323) For an element (or pixel) j conductivity change and element area are denoted as j and Aj respectively. For later investigation, we define the relative QI (QI ) for the reference QI as in Equation 324. 0QIQIQI (324) Unless otherwise noted, QI for an anomaly lo cated at the domain center was used for QI0. PAGE 43 43 CHAPTER 4 ESTIMATING ELLIPTICAL SHAPE CHANGE OF AN OBJECT BOUNDARY USING JOUKOWSK I TRANSFORMATION IN EIT Background Measurements in EIT are sensitive to uncertaintie s related w ith electrod es and variations of body shapes (i.e. respiratory movements and postu ral changes). Consequen tly, artifacts due to the corresponding errors in the measurement ofte n dominate reconstruction images, which lead to significant errors in estimating anomaly volum es. These errors have to be corrected for successful EIT imaging and quantitation (Adler et al. 1996, Sadleir and Fox 2001, Zhang and Patterson 2005). Previously, it has been shown that even a slight shape mismatch in the boundary shape models can produce significant image artif acts (Breckon and Pidcoc k 1988, Kolehmainen et al. 2005). Additionally, the body is subject to temporal shape changes such as thoracic expansion (Adler et al. 1996). Shape variations of object boundaries are closely correlate d with variations in electrode locations, therefore, problems such as electrode misplacement (Barber and Brown 1988), and electrode movement (Soleimani et al. 2006) are analogous to the shaperelated problems. Uncertainties related with body shape and electrodes cause si milarly characterized image artifacts (i.e. spurious peaks near the object bounda ry) (Soleimani et al. 2006, Boyle et al. 2008). Here, we use the term shape change to describe these uncertain ties in a collective sense (Figure 41). In difference EIT, recovery of reconstructi on quality may be possible using the knowledge about the correct boundary shape. Barber and Br own (1988) showed that image artifacts due to electrode errors were effectively removed us ing the reference measurements made on the homogeneous model of the correct boundary shape. PAGE 44 44 Figure 41. Extent of shape change. Shape estimation using CT, MRI or Ultrasound (Molebny et al. 1996) is not always available. On the other hand, we can readily approximate object shapes all the time from EIT measurements, since they depend principall y on shape (Barber and Brown 1988). Breckon and Pidcock (1988) estimated domain boundary by desc ribing its shape approximately in terms of Fourier series. Kolehmainen et al. (2005) introduced a numerical me thod that finds a deformed image of the original isotr opic conductivity based on the th eory of Teichmller spaces. Information about electrode locations can roughly describe boundary shapes, and can be used to enhance reconstruction images. Kiber et al. (1990) parameterized measurements in terms of separation of electrodes, and estimated electr ode positions iteratively to approximate boundary shapes. Soleimani et al. (2006) proposed the use of an augm ented Jacobian matrix. Using the matrix, both the spatial information of electro de movement and the c onductivity change was simultaneously reconstructed. Measurement sensitivity to unit perturbation of each electrode position composed the augmented co lumns of the matrix. Heikkinen et al. (2002) conducted a PAGE 45 45 similar approach, which simultaneously recons tructed the conductivity distribution and the contact impedance of the electrodes. In this chapter, we investigate image artifacts caused by our shape change models. We model shape changes in 2D simply as elliptical shape distortions of a unit disk boundary. In 2D EIT image reconstruction, an object is often a ssumed to have circular boundary. Our motivation is that parts of human body genera lly have elliptical profiles rath er than circular. We describe elliptical shape change by the Joukowski transf ormation, which can define the relation between circular and elliptical geometri es. Based on the relation, correct reference measurements can be obtained from the solution of the circular dom ain (Lionheart 1998). We de monstrate a method to estimate shapes as ellipses defined by a si ngle geometrical parame ter of the Joukowski transformation (the Joukowski parameter ). In our method, images of artifacts due to elliptical shape changes, which we use as a template are obtained from measurements made on homogeneous elliptical domains of various eccentricities. Images reconstructed from actual measurements may then be compar ed with the template to find the ellipses that produce the best match. In order to validate our method, analysis was done with simulated measurements with and without added noise; and phantom measurements (e lliptical phantoms filled with saline solution as a background material). Data were gathered using the adjacent electrode configuration with 16electrode full array, and images were ge nerated using TSVD regul arized reconstruction method. Method In 2D EIT reconstru ctio n, the domain is often assume to be circular when the accurate shape of the model is unavailable. However, crosssections of body parts such as abdomen or PAGE 46 46 chest are noticeably more elliptical in profile th an circular. Here, we investigate shape changes simply as elliptical boundary distortions on a disk model (Oh et al. 2007). The Joukowski Transformation and Elliptical Boundary Shape Change Confor m al mapping is an anglepreserving tran sformation, meaning th at the angle between any arbitrary two vectors is preser ved after the transformation. It is useful when the difficulty of a problem can be expressed in terms of geomet ry. According to the Riemanns theorem, any object with closed boundary can be mapped co nformally into a unit circle. The Joukowski transformation can describe a conf ormal relation between elliptical and circular geometries as in Equation 41. We use the Joukowski transforma tion to model elliptical shape changes. Z cZZgW 1 )( (41) where c is termed the Joukowski parameter In Equation 41, the relation g( Z ) is conformal, and Z and W are complex numbers that represents spatial coordinates of the electrodes. Z denotes the boundary points of a unit circle and W those of an ellipse. The main axes of the ellipse are 1 + c in length horizontally and 1 c vertically (Appendix A). Therefore, c is a geometrical parameter of an ellipse determining degrees of elliptical shape changes. Figure 42. A 2D disk model after three different degrees of elliptical shape change ( c = 0.1, 0.2 and 0.3) under the Joukowski transformation. PAGE 47 47 In Figure 42, we illustrated three different degrees of changes ( c = 0.1, 0.2 and 0.3), and laid the initial disk model of no shape change over the ellipses. Reference boundary voltage measurements are in itially obtained from a disk object using equidistant boundary electrodes. After the shape change it is apparent that the object is no longer a disk, thus the initial reference measurements are incorrect. The correct reference measurement can be obtained conveniently from the forward solu tion of the disk domain. In Figure 43, due to the Joukowski transformation, locations of the boundary electrodes have changed, with electrodes near the x axis having experienced larger shifts along the circular arc than the ones near the y axis. Assuming that the electrodes had been placed equidistantly on the elliptical boundary, the correct reference measurements we re obtained from the circular model using electrodes with shifted locations by the inverse transformation ( g1( W )) of the Equation 41. Figure 43. Inverse Joukowski transformation of an elliptical domain with 16 equidistant boundary electrodes. A) Elliptical domain ( c = 0.2) with equidistant boundary electrodes. B) Disk domain with nonequi distant boundary electrodes whose locations are shifted by the inverse Joukowski transformation. Boundary voltage change (V ) after the elliptical shape change is obtained as follows. First, the reference boundary voltage ( V0) was measured on a 2D disk object (radius: rd = 1) of A B PAGE 48 48 homogeneous conductivity usi ng the 16electrode full array attached on its boundary. Then, various amounts of elliptical shape changes were applied on the disk object containing a single disk anomaly (radius: ra = 0.1) at its center. The boundary voltages ( Ve) were measured using electrodes with shifted locations (by g1( W )) on circular boundary (41). Finally, the measurement set (Ve) due to elliptical shape change was obtained by subtracting the reference V0 from Ve as in Equation 42. 0VVVee (42) In Figure 44, measurements are displayed for three degrees of elliptical shape change ( c = 0, 0.01 and 0.1). The main attribute of measuremen t errors due to shape change is sharp spurious peaks near the domain boundary. In Figure 44, peak amplitudes are observed to grow as c values increase. It is also apparent that the m easurements are completely dominated by this shape change over c = 0.1. Figure 44. Measurements made on the disk model. A) of no shape change ( c = 0), B) of elliptical shape change c = 0.01 and C) c = 0.1. A B C PAGE 49 49 Conductivity images were produced from Ve using TSVD regularized reconstruction as in Equation 43. e eV S, (43) where e is the conductivity change reconstructed from Ve. Artifacts in the images occur, because effects of the shape change on rec onstruction images are through compensatory conductivity changes within the domain (Tossavainen et al. 2004). Figure 45 shows images greatly corrupted by the shape change artif acts, whose peak amplitudes increased with c. Figure 45. Reconstruction image from the meas urement on a disk object containing an anomaly at the center. A) No shape change (c = 0), B) elliptical shape change (c = 0.1), and C) elliptical shape change (c = 0.2). Images were reconstructed using TSVD with truncation at 63. In Figure 46, the absolute maximum and minimum of e values, which represent positive and negative artifact peaks respectively, are plotted for various c values. Results of linear fits to these curves allowed several fi ndings: The increase of positive peaks was more linear, while that of negative pe aks was more rapid for increasing amount of the shape change. In order to investigate shape change in quantitative EIT, average of e (an analog of QI) was calculated and plot ted in terms of c (Figure 47). The average of e decreased as c increased, and it went below zero at about c = 0.15. Additionally, QI errors were less than 5% for shape A B C PAGE 50 50 changes whose degrees (c) were less than 0.5 (Oh et al. 2007). Therefore, shape change is seen to be a more important issue in imaging than in quantitation. Figure 46. Absolute values of maximum and minimum of reconstruction amplitude () for various degrees of shape change c (%). Figure 47. QI vs. the degrees of shape change (Joukowski parameter: c (%)). PAGE 51 51 Boundary Shape Estimation Artifacts due to shape changes are undesirabl e in reconstructed images, but they may contain concentrated information about shape changes. The main attributes of the artifacts were found to be high amplitude peaks with alternatin g signs near the object boundary (Figure 45). The TSVD reconstruction seems to have a tenden cy to push the artifacts close to the boundary. We also found that the peak magnitude of th e boundary artifacts increased consistently with c (Figure 46), while their patterns were relativel y preserved (Figure 45) Based on these findings, we suggest a method to estimate elliptical boundary shape from EIT reconstruction (Figure 48). Our method estimates the changed shape as an ellipse, defined by a single geometrical parameter c (the Joukowski parameter) as in Equation 41. Figure 48. Summary of our me thod to estimate boundary shape. In Equation 43, images of reconstructed conductivity e were obtained by TSVD. Because the artifacts due to shape changes collect mainly in the region near the boundary, we selected a group of elements (pixels) near the boundary (Figure 49), and denoted them collectively as J. Finally, a template T was defined as the reco nstruction pixel values (e) in J that depends on c as in Equation 44 (Figure 410). cJcTe, (44) PAGE 52 52 Figure 49. The boundary elements that were used for the template (filled with color). Figure 410. Template vs Joukowski parameter (c = 0 ~ 50%). PAGE 53 53 Typically in EIT, accurate know ledge of boundary shape is unavailable, but conductivity change a reconstructed by these actual measurem ents can readily be compared with T(c). Now, the shape change can be estimated by finding the c value that produces the least sum of squared errors as in Equation 45. J a ccTJ2)()( minarg (45) Method in Equation 45 often fails when the mean difference is large between T and a. In order to resolve this problem, we considered a scale factor (s) that was defined as the ratio of the maximum pixel values of a and T as in Equation 46. J a ccTJs2 1)()( minarg (46)Computer Simulation We designed a 2D disk model (rd = 1) with 16electrode full array using Comsol Multiphysics and Matlab. All electrodes had the conductivity of copper (5.997[Sm1]). The length of each electrode was 0.1, and subtended an angle of 5.7o on its circular boundary. The models were discretized to 1916 secondorder triangular finite elements, and then solved using the direct linear system solver UMFPACK (U nsymmetric MultiFrontal PACKage). The model contained a small disk anomaly (ra = 0.1rd) placed at one of five di fferent locations (relative radial distances of: 0, 0.2, 0.4, 0.6 and 0.8 from the origin) on the positive x axis. The model background conductivity was set to one, and anomalies represented a unit conductivity increase from the background. Actual measurements Ve as in Equation 42 were obtai ned using adjacent electrode configuration for the model underg oing five differe nt degrees (c = 0.1, 0.2, 0.3, 0.4 and 0.5) of the elliptical shape change describe d by the Joukowski transformation. PAGE 54 54 Phantom Experiment We constructed an acrylic cylindrical tank (rP = 14.0 cm) as shown in Figure 412. Sixteen stainless steel bar electrodes (13.0 mm 10 cm) were attached equidistantly around the tank inner surface. Image sensitivity decreases significantly when anomaly volume is outside the electrode plane (Metherall P et al. 1996), so we used vertically extended electrodes previously used in the study by Sadleir and Fox (1998). A hole was drilled on each electrode center, and the same sized holes were drilled on the tank wall at a certain height (6.5 mm) from the bottom. Bolts were inserted so that the threaded part s stick outward. They were fastened while being checked for any water leak. Finally, the bolts we re connected to 16 channels of the external measurement system. The tanks were filled with 5 L of saline solution (0.2 Sm1) as a background material. The surface level of the saline solution inside the tank was within the height of the bar electrodes. Two insulating plastic rods of radii (rA) 1.3 cm and 1.8 cm were used as anomalies I and II. Anomalies I and II were placed at 9 and 5 different locations respectively, on the positive x axis. The displacement from the phantom center was from 0 to 11.2 cm with 1.4 cm gaps for anomaly I, and with 2.8 cm gaps for anomaly II. The elliptical phantoms were design ed (Figure 411) to approximate c = 0.1 and 0.2 shape change applied to the cylindrical phantom shown in Figure 412. The lengths of semimajor and semiminor axes of the elliptical plane were 15 cm and 12.5 cm (elliptical phantom I), and 16 cm and 11.5 cm (elliptical phantom II). The c values of these ellipses were calculated to be 0.09 and 0.16 respectively. PAGE 55 55 Figure 411. Elliptical phantoms simulating two different degrees of shape change (c = 0.1 and 0.2). The 16 stainless steele bar electrode s are attached on the elliptical boundary equidistantly. Another phantom (rP = 10.0 cm) was constructed from a flexible plastic bucket (Figure 413). This phantom was filled with 4 L of saline so lution. In this phantom, shape changes could be applied by inserting flat objects beneath the base of the phantom to distort the boundary shape, thus resulting in the upper crosssections of the bucket to become more distorted than the lower. Four increasing degrees of elliptical shape change were applied in this manner. This type of shape change was deemed to approximate real body shape changes more closely than the elliptical phantoms (i.e. abdominal breathing). PAGE 56 56 Figure 412. Prototype KHU Ma rk1 with 16 channels applie d to a cylindrical phantom. Data were obtained using a prototype of KH U Mark1 system (Figure 413). KHU Mark1 is an EIT system developed in the Impedance Imaging Research Center (IIRC), Kyung Hee University (KHU), Suwon, Korea. The system allo ws flexibility for addre ssing electrodes, with a wide range of operating fre quency (10 kHz ~ 500 kHz) (Oh et al. 2007). We collected measurements at 10 kHz for i nput current applied at 1 mA. PAGE 57 57 Figure 413. A malleable cylindrical phantom with 16 bar electrordes. PAGE 58 58 Results Simulated Data Figure 414 shows reconstructed images with an anomaly at each of five different positions for an elliptical shape change (c = 0.2). The shape change affected images of a central anomaly more, and the patterns of the artifacts near the domain boundary are similar. Figure 414. Reconstructed images for a single anomaly on five different locations. A) before and B) after shape change (c = 0.2) For noiseless measurements, the results of sh ape change estimation are shown in Figure 415. In Figure 415 (A), the minimum of the cost function is seen to correspond exactly to the degree of the shape change applied (c = 0.2). We used the estimated shape information to obtain the reference measurement, which was applied to reconstruction and improved image quality greatly (Figure 415 (B)). The White Gaussian Noise (WGN) of 20dB Si gnaltoNoise Ratio (SNR) was added to corrupt these simulated measurements. Results in Figure 416 shows that the overall errors were only about 2% for shape changes up to c = 0.5 (50%). A B PAGE 59 59 Figure 415. Cost function curve for estimating the amount of shape change. Reconstructed image recovered using the estimated shape information is shown on the right. Figure 416. Estimation errorbar plot for various amounts of shape change (c). Measurements were noisy with additiv e WGN of 20dB SNR. PAGE 60 60 Phantom Data We estimated c values from the measurements made on the cylindrical and elliptical phantoms. In Table 41, the estimated values of c were averaged across anomaly locations and expressed in percentile. The results show that we succeeded in estimating c with discrepancies less than 0.04. The scale factor (s = 1.50175) was calculated from th e reconstruction of the cylindrical phantom. In the malleable phantom case, shapechange corrupted measurements (Figure 417) did not match the characteristics of the simula ted as previously shown in Figure 44. Figure 417. Reference measurem ents from homogeneous medium for four different amounts of shape change However, we confirmed the spatial charac teristics of boundary image artifacts for both were similar. Using our method, c values were estimated as in Figure 418. We do not know what the c value was for each shape change, but Figure 418 at least sh ows monotonous increase for increasing amounts of shape change. PAGE 61 61 Figure 418. Elliptical boundary shape estimation by the Joukowski parameter c (%). Four different degrees of shape change were ap plied on a malleable plastic phantom in an increasing order. Discussion and Conclusion We developed our method motivated by our obs ervation that the reconstruction process naturally divides influences of conductivity di stribution and boundary shape change. Artifacts due to boundary shape changes were observed to cause spurious boundary peaks, which were used as features to estimate elliptical shap e changes. Utilizing the reference measurement obtained from a homogeneous model of the estimat ed shape, we recovered the original image quality. We were roughly able to estimate shape cha nges from the simulated measurements with and without AWGN, and phantom measurements. Phantoms were made either to be elliptical or to be malleable so that approximately elliptical shape distortions can be applied. For simulated model of c = 0.5 shape change, estimation errors were under 5% in the pres ence of 20 dB SNR PAGE 62 62 AWGN. In phantom experiments, we could always estimate approximate values of c (Table 41 and Figure 418). Table 41. Estimated c (percentile) from the measurements made on cylindrical and elliptical phantoms I and II. Anomaly I (rA = 1.3) was positioned in nine different locations, and Anomaly II (rA = 1.8 cm) in five locations. The c values were averaged across anomaly position. Cylindrical Elliptical I Elliptical II Expected c values 0 9 16 Estimated c values (Anomaly I) 1 12 19.9 Estimated c values (Anomaly II) 1 11.8 19.6 The method by Soleimani et al. (2006) successfully reconstructed images of reduced artifacts for both conformal and nonconformal sh ape changes. However, a recent investigation by Boyle et al. (2008) found that it could not estimate conformal shape changes. Our proposed method is a templatematching algorithm and does not modify the reconstruction matrix. As long as changes of domain shapes produce similar boundary artifacts, our method is expected to always find an approximate estim ation of the changed shape. In Figure 419, areas of the ellipses transf ormed from a unit circ le via the Joukowski transformation are plotted in terms of c. The area of ellipse can be calculated as (1 + c)(1 c). For c = 0.1, the difference between the crosssecti onal areas of the original and changed models is about 1% of the original ar ea. Therefore, we regarded the Joukowski transformation to be an appropriate model for most body sh ape changes that are typically c = 0.1 and below. However, it would still be interesting to investigate realistic crosssectional area changes of the subject domain. The results from this investigation could be utilized for further improvements in image quality. PAGE 63 63 Figure 419. Area of the domains for the el liptical shape change of various degrees (c). PAGE 64 64 CHAPTER 5 NORMALIZATION OF SPATIALLYVARIANT IMAGE RECONST R UCTION PROBLEM IN EIT USING SYSTEM BLURRING PROPERTIES Background Previous works suggested quantitative EIT as a potential method to monitor volume changes inside the body. Exemplary areas include air and liquid in the lu ngs and bleeding in the peritoneum. Relevant invivo experiments were done in monitoring dynamic hyperinflation in dogs (Adler et al. 1998) and fluid volume changes in Contin uous Ambulatory Peritoneal Dialysis (CAPD) (Sadleir and Fox 2001). 1 A main challenge in quantitative EIT is the erro rs related with the spatial variability of EIT imaging systems. The image amplitude of a conduc tive anomaly varies with its position (Sadleir and Fox 1998). In monitoring changes of lung air and liquid volumes, Adler et al. (1997) reported errors that depended on ra dial positions of a same anomaly, despite exclusion of a noisy region near the object boundary. Spatial variability in quantification may not be a significant problem when the region of interest (ROI) is relatively small. In Blott et al.s study (2000) of simulated twodimensional intraventricular hemorrhage, quantitation errors due to spatial variance were likely negligible for that re ason. However, determining an ROI is not straightforward for some applications. In intraperitoneal hemorrhage a pplication, bleeding may occur over a large fraction of the imaged domain (Sadleir and Fox 1998). Therefore, reducing effects of the spatial variabil ity on estimating quantities is de sired to improve accuracy of quantitative EIT reconstruction. 1 This chapter was published in Physiological Measurement, volume 30, pages 275289 in 2009. The work was supported in part by NIH grant R01EB002389, and by the US Army Medical Research and Materiel Command under Award No. w81xwh0710591, both to RJS. PAGE 65 65 In order to reduce variations in quantity estimates, Thomas et al. (1994) used a method that scaled logresistivity images pixelwise. A st udy by Sadleir and Fox (1998) demonstrated a postreconstructive filtering method th at combined pixelwise scaling and conformal transformations in the context of a linearized backprojection method. The problem of a spatially variant imaging system was investigated in detail by Wheeler et al. (2002) in the context of resolu tion measures. In order to iden tify the resolution of EIT images, spatially invariant measures were sought. In th eir comparative study, an area integral measure turned out to be the most r obust against spatial variance. We used reconstruction blur an alysis in order to reduce erro rs caused by spatial variance. We demonstrate that normalizing the Point Spre ad Function (PSF) by its integral can result in improvements. CohenBacrie et al. (1997) used variance of this PSF as penalty term in Tikhonov regularization. In our case, we used integral s of PSF functions as normalizing terms. Three mathematical frameworks of normalization we re compared: PixelWise Scaling (PWS), Weighted PseudoInversion (W PI) and Weighted Minimum Norm Method (WMNM) for twodimensional eight electrode full array and hemiarray reconstructions. Data were obtained from numerically simulated twodimensional disk mode ls containing simulated disk anomalies placed at various locations in the image plane and a cylindrical saline phantom containing bloodlike anomalies in a less conducting background. Method Point Spread Function and the Blur Matrix In this chapter, we introdu ce the point spread function ( PSF) in order to investigate blurring in EIT image reconstruction. Consider a unit conductivity change in an element (pixel) j within the image. This unit change can be expre ssed in a vector form as an entry of 1 in the jth position and zeros elsewhere ([0 ... 0 1 0 0]T) (1N), whose weighted sum is a vector form of PAGE 66 66 the ideal image. The change in boundary voltage measurements subject to this unit conductivity change is the jth column of the sensitivity matrix S (MN). Ideally any measurement can be represented as a weighted sum of sensitivity matrix columns. Therefore, the column space of S can be termed the ideal measurement space. By reconstructing from the ideal measurements, a blurred version of the ideal element (pixel ) image is obtained, which we term the PSF in Equation 51. T j00100 B (51) where B (NN) is the blurring matrix, and j is the PSF vector for an anomaly in the jth element. The column space of B is a blurred version of the ideal images. The blur matrix B is defined as a product of the sensi tivity matrix and the reconstruction matrix. Therefore, it can be easily obtained for any reconstructio n matrix that is expressed in algebraic form. In Equation 52a and Equation 52b, definitions of B are displayed in the c ontext of TSVD and WMNM regularizations respectively. SSB (52a) SSWWB (52b) It should be noted that one can easily obtain B for other regularization methods such as Tikhonov regularization in Equati on 322. The definition of B for Tikhonov regularization is shown in Equation 53. SSSSB T TI1 (53) Using SVD, we can estimate basis images that determine the system blurring B. SSB (54a) T TUDV UDV (54b) PAGE 67 67 TDVVD (54c) T kV 00 0I V (54d) where Ik (kk) is the identity matrix. It is clear from Equation 54d that columns of the right singular matrix V can be seen as eigenimage vectors. It is also clear that the blurring depends on the truncation number k. Normalizing Terms In order to normalize QI, the column sum of B, defined as qj in Equation 55, was considered. N i ij jq1B (55) The normalization matrix used throughout this paper is defined as Q (NN) as in Equation 56, which is a diagonal matrix whose entries are qjs. Nq q q00 00 00 002 1 Q (56)Pixelwise Scaling Thomas et al. (1994) used a method that scaled logre sistivity images pixelwise in order to reduce variation in quantity estimates. We can also apply pixelwise scaling (PWS) to our case in which conductivity change is to be reconstructed. In Equation 57, PWS is presented in matrix form. V SQ1 (57) A more advanced method was proposed by Sa dleir and Fox (1998), who demonstrated a postreconstructive filtering method that co mbined pixelwise scaling and conformal PAGE 68 68 transformations. In their work, the maximum spatial error in quant itation from phantom experimental results decreased from 30% to 6%. However, a large region near the object boundary (normalized radius > 0.75) was excluded from quantification. Weighted Pseudoinverse The normalization that we use here is wei ghted pseudoinversion (WPI). The normalized reconstruction can be obtained by weighting columns of the sensitivity matrix by terms in Q prior to reconstruction as in Equation 58. V SQ (58) If SQ is a fullrank matrix, Equation 58 will become identical to Equation 57. In this vein, we speculate that WPI should normalize the reconstruction in a sim ilar way to PWS, with an additional advantage of usi ng the truncated pseudoinverse. WMNM Normalization In the work of Oh and Sadleir (2007), a W P I method was shown to decrease spatial variance of the QI. However, the normalized reconstruction in Equation 58 did not produce reasonable images. To obtain meaningful images, one can borrow the mathematical framework of WMNM, which can be done by multiplying Q postreconstructively. The final form is shown in Equation 59. V SQQ (59) It was also noticed that the WMNM framew ork can be used to introduce normalization effects in some cases. Sadleir et al. (2008) used a combined met hod of WMNM regularization in Equation 319 and WMNM normalization in Equation 59 in the eight electrode hemiarray case, which produced normalized QIs and images. Th e matrix equation of this normalized reconstruction method is shown in Equation 510. PAGE 69 69 V SWQQW (510) Normalized Sensitivity Matrices The pseudoinverted part of the reconstr uction m a trix can be understood as pseudoinversion of a matrix product of the normalization matrix ( Qt or Qw) and the sensitivity matrix ( S ). Matrices Qt and Qw were calculated using diff ering definitions of the B matrix as shown in Equation 52a and Equation 52b respectively. Cha nges of the system matrix condition number owing to this matrix multiplicati on are presented in Table 51. Table 51. Condition numbers of the normalized system matrices System matrix Condition Number (full array) Condition Number (hemiarray) S 2.2117 2.0218 SW 2.0017 3.21 18 SQt 6.3717 3.9317 SWQw 2.3117 2.5118 SWQt 2.0517 4.1017 Table 51 shows that matrix conditions can improve depending on the choice of the normalization matrix. For a full array, W tended to decrease the condition number, while Qt and Qw did not. For the hemiarray, Qt tended to decrease condition number, while W and Qw tended to increase them. Therefore, when noisy meas urements are to be reconstructed, using a normalization matrix that decreases c ondition number could be beneficial. Reconstruction Matrix after Normalization For convenience, each n o rmalization method was given an index as in Table 52. Normalization method C5 defined in the last row of the Table 52 is based on WMNM reconstruction and WMNM normalization as for B4. However, normalization matrix Qt was used instead of Qw. PAGE 70 70 Computer Simulation Twodim e nsional forward models of a disk (rd) containing a single internal anomaly (ra) at various locations were designed and solved usin g Comsol Multiphysics and Matlab. We tested Table 52. Index assignment for various normalization methods. The methods are represented by reconstruction matrices. Th e normalization matrices Qt and Qw were calculated using blur matrix definition in (5 2a) and (52b) respectively. Methods A1 ~ A4 are based on TSVD reconstruction (316), methods B1 ~ B4 and C5 are based on WMNM reconstruction (319). Index Normalization Method Reconstruction Matrix A1 None S+ A2 PWS Qt 1S+ A3 WPI (SQt)+ A4 WMNM Qt( SQt)+ B1 None W ( SW)+ B2 PWS Qw 1W ( SW)+ B3 WPI W ( SWQw)+ B4 WMNM QwW ( SWQw)+ C5 WMNM QtW ( SWQt)+ two types of models with anomaly locations as shown in Figure 21 and Figure 22. The full array model had eight boundary electrodes placed equidistantly (Figure 21). The hemiarray model had eight electrodes placed on the anteri or boundary only (Figure 22). All electrodes had the conductivity of copper (67 Sm1). The length of each electrode was 0.1 relative to the disk radius, subtending an angle of 5.7o on the disk perimeter. The models were discretized to 1374 secondorder triangular finite elements, and then solved for boundary voltage values subject to adjacent input current patterns using the direct linear system solver UMFPACK. The anomaly had a radius 0.1rd, and was centered at locations with relative radii (ra/rd) of 0, 0.2, 0.4, 0.6, 0.8 from the origin. For hemiarray cases the anomaly position was varied in angles ( ) from 0o to 180o with 5o increment. In full array cases only = 90o was used. The model background PAGE 71 71 conductivity was set to one, and anomalies repr esented a unit conductivity increase from the background. Phantom Experiment A cylindrical phantom (rP = 14.0 cm) was filled with saline solution (5 L). Vertical bar electrodes (width 13.0 mm and length 101.5 mm) were attached equidistantly around the phantom boundary to create an approximately two dimensional field pattern. An insulating plastic rod of radius 1.4 cm (0.1rP) was used as an anomaly, and it was moved along the and axes (at = 0o, 90o, and 180o), with the anomaly center pl aced successively at relative radial displacements (rA/rP) of 0 (center), 0.25, 0.5, and 0.75, ten locations in total. This procedure was repeated for different volumes (50, 100, 150, 200, and 250 mL) of anomaly having conductivity similar to that of blood (0.67 Sm1) placed in the background of saline with a conductivity of 0.2 Sm1. Figure 51. Overview of the ePack system (picture taken from Tucker et al. (2008)) PAGE 72 72 Measurement System W e used the ePack, 8 channel EIT system to collect measurements (Tucker et al. 2008). The ePack is an integrated EIT system targeted at identifying intraperitoneal hemorrhage as a result of blunt trauma. Its composed of expandable belt of bar electrodes, current source, DSPbased measurement system and imaging software (Figure 51). It provides flexibility in measurement by allowing various el ectrode configurations that can be easily set up from imaging the software. The ePack measurement system take s only 350 ms to collect one complete set of measurements, and its error range is within 0.2 % for a nominal 50 resistance (Tucker et al 2008). The device was operated at a 2 mA constant current and 62.5 kHz. Results Ideal Reconstructions Ideal m easurem ents were obtained by indi vidually perturbing all elements by unit conductivity increase. They are effectively columns of the sensitivity matrix. Their corresponding QI values were calculated for al l the measurements. We compared normalization methods in terms of Relative STD ( RSD : absolute value of variation expressed as standard deviation divided by the average) calculated using all QI values as in Equation 511 (Table 53). )(/)()( QIaverage QISTDQIRSD (511) Method A3 resulted in an RSD decrease in all cases, while C5 resulted in an overall decrease of RSD except along the hemiarray axis. Table 53. RSD of QI (511) calculated from the idea l measurement reconstructions (at full rank, k = 20) of eight electrode cases. Normalization Method Full array ( = 0o) Hemiarray ( = 0o) Hemiarray ( = 90o) Hemiarray ( = 180o) A1 0.1478 0.1749 0.1734 0.9579 A2 0.3294 0.4000 0.2586 0.4829 A3 0.0628 0.1023 0.1019 0.9095 A4 0.0878 0.1053 0.1790 0.9846 PAGE 73 73 B1 0.1124 0.5440 0.8739 0.1946 B2 1.1087 5.8956 2.9424 16.7774 B3 2.4142 0.8479 0.2373 0.5036 B4 0.2141 0.2729 0.0936 0.8172 C5 0.0672 0.3413 0.0184 0.6702 Simulated Data Sim u lated measurements were created from the computer models as described in the method section. We investigated relative QI (QI) (%) with the central anomaly as the reference. Maximum deviation errors () away from the central QI value (QI0) were calculated as in Equation 512 and displayed in Table 54. 1 QIMax (512) Method A3, which was shown to perform most st ably in the case of ideal measurements, also reduced for simulated measurements. However, it failed to decrease for hemiarray axis. Method B4 showed overall decrease of for hemiarray cases, but its increased for the full array case. Method C5 produced results sim ilar to the ideal measurement case. Its decreased in all cases except for the hemiarray axis. Table 54. Maximum deviation error () of the relative QI, calculated from the simulated measurement reconstructions (k = 16) of eight electrode cases. Truncation number k was chosen based on the l curve investigation Normalization Method Full array ( = 0o) Hemiarray ( = 0o) Hemiarray ( = 90o) Hemiarray ( = 180o) A1 0.7239 0.3134 0.7409 0.9539 A2 2.7983 2.2928 0.5365 0.5759 A3 0.2437 0.0671 0.8284 0.9075 A4 0.4398 0.3078 0.8478 0.9574 B1 0.2740 3.9886 2.8486 0.7287 B2 1.0667 4.9794 5.2909 2.2658 B3 0.3379 6.2294 2.2075 0.9253 B4 0.9985 0.1382 0.5224 0.8945 C5 0.2239 0.5944 0.6117 0.7948 PAGE 74 74 Phantom Data Phantom m easurements were made as previously described in the met hod section. In Table 55, values are displayed fo r various volumes of the bloodlike anomaly. The values were calculated from QI values calculated from ten different anomaly positions. Table 55. Maximum deviation error of the relative QI (), calculated from phantom reconstructions of eight elec trode full array and hemiarray topologies. In the full array, the anomaly was placed on 4 different locations on the negative axis with relative displacement from the domain centre (rA/rP = 0, 0.25, 0.5, and 0.75). In the hemiarray, the anomaly was placed on 10 different locations on the and axes (rA/rP = 0, 0.25, 0.5, and 0.75). Truncation number k was chosen by investigation of the l curve. Full array Anomaly Volume Normalization Method 50mL ( k = 16) 100mL ( k = 17) 150mL ( k = 17) 200mL ( k = 17) 250mL ( k = 17) A1 0.3344 0.3799 0.3048 0.2631 0.2377 A2 0.8131 0.8741 0.6714 0.6055 0.5685 A3 0.1570 0.1634 0.1517 0.1173 0.1017 A4 0.2266 0.2588 0.2129 0.1760 0.1562 B1 0.1896 0.1918 0.1831 0.1492 0.1318 B2 1.6966 9.6980 17.6351 5.3715 3.5571 B3 0.1961 0.3818 0.1346 0.1869 0.1777 B4 0.4055 0.4456 0.3568 0.3126 0.2852 C5 0.1513 0.1582 0.1483 0.1138 0.0980 Hemiarray Anomaly Volume Normalization Method 50mL ( k = 15) 100mL ( k = 16) 150mL ( k = 16) 200mL ( k = 16) 250mL ( k = 17) A1 0.5084 0.6845 0.4759 0.4213 0.4219 A2 0.5319 0.2658 0.1438 0.2028 0.2872 A3 0.6003 0.6216 0.5373 0.5155 0.5536 A4 0.5923 0.7090 0.5181 0.4819 0.4940 B1 0.6883 0.4502 0.2519 0.4558 0.3855 B2 3.0771 4.5643 8.2285 6.4299 1.3628 B3 0.7296 0.8014 0.7127 0.8695 0.9601 B4 0.4587 0.4695 0.3345 0.3314 0.4013 C5 0.3072 0.2578 0.2660 0.2253 0.2609 Our results in Table 55 show that method A2 worked better for anomalies of larger volumes in hemiarray case. B4 showed a decrease for hemiarray case, and values tended to PAGE 75 75 increase as the anomaly volume increased. This may have been due to our assumption of reconstructions being the weighted sum of systems PSFs. Method C5 decreased values for anomalies of all volumes in both full array and hemiarray cases. Cross examination of Tables 53, 54 and 55 allow us make some findings. For the full array, the spatial variability (in terms of relative STD of QI or ) decreased for A3, A4, B1 and C5. A3 and C5 seem to be best candidates since they resulted in the most decrease in most cases. Figure 52. Eight electr ode full array images ( k = 17) reconstructed from phantom measurements made for a bloodlike anomaly (100mL) in four different locations on the negative axis (rA/rP = 0, 0.25, 0.5, and 0.75). A) Images reconstructed using normalization methods A1, A3, A4, B1 and C5 are compared. B) Images reconstructed using A1, B1 and C5 were inve stigated in terms of the anomalys blur radius (resolution) and positi on error. Resolution was calcula ted as a square root of the ratio between area of the half amplitude set and the actual anomaly size. Position was estimated as a centre of mass of the half amplitude set. A B PAGE 76 76 In Figure 52 (A), methods A3 and A4, which successfully reduced spatial variability of QI, resulted in blurred reconstr uction of a central anomaly compar ed to reconstructions using B1 and C5, which produced more compact central anom aly. It should be noted that B1 and C5 are methods that are based on WNMN regularizati on. Method B1 blurred im ages of the anomaly near the boundary. C5 did not produc e a blurred reconstruction of th is location. Thus, C5 appears to be the most reasonable method for full array. In Figure 52 (B), selected methods (A1, B1 and C5) were compared in terms of resolution and pos ition error. In order to estimate resolution, we followed definition of the blur radius suggested by Adler and Guardo (1996). Elements that have values greater than half the maximum element va lue were chosen as the Half Amplitude (HA) set. Then, the resolution was calculated as a square root of the ratio between area of the HA set ( AHA) and the domain area ( Ao) as in Equation 513. o HAA A resolution (513) The position of an anomaly was located as the centre of mass of HA set in reconstruction images as in Equation 514. m m m mmp position (514) Here, pm is the position vector (xm,ym) within the domain. Small and uniform resolution values are regard ed to be desirable indicating less blurring. In Figure 52, there are overall enhancement of resolution when using B1 and C5. In terms of position error, C5 shows great improvement except when the anomaly was close to the domain boundary. PAGE 77 77 For the hemiarray case, normalization could be made particularly effective for anomalies on the axis. However, this was obtained at the expense of increasing variability along the axis. Even though B4 and C5 were most effective for and axes anomalies, they were not effective for axis anomalies (Table 54 and 55). PAGE 78 78 Figure 53. Eight electrode hemiarray images (k = 16) and relative QI curves reconstruc ted from a bloodlike anomaly (100mL) using methods A1, B1, B4 and C5. Relative QI values were calculated with the central anomaly (rA/rP = 0) as a reference. All relative QI curves were scaled the sa me (0.5~2) and overlaid on top of each corresponding reconstruction image. The anomalies were positioned at ten different locations on the and axes (rA/rP = 0.25, 0.5, and 0.75). PAGE 79 79 Investigation for various anomaly volumes (50, 100, 150, 200, and 250mL) used in our phantom experiments (Table 55) revealed that B4 and C5 consistently reduced spatial variability. In Figure 53, methods that reduced variability are comp ared for 100mL anomaly in ten locations. Anomaly images reconstructed using B1 were comp aratively less blurry than A1, B4, and C5. In reconstruction of anomalies on the axis, B1 reduced the spa tial variability. Method B4 performed well in terms of reducing spatial va riability of QI, but did not produce a good image quality, as shown in Figure 53 (B4). Method C5 resulted in reduced variability producing images of improved qualities. However, images of anomalies on the axis were still unclear. In Figure 54, performance of our proposed method C5 was demonstrated for normalizing QI in full array and hemiarray topology. The relative QI values were plotted versus radial axis. The uniformity of the relative QI curves increased, when C5 was used. The improvement was consistent for all axes considered (the and axes. Discussion Im age reconstruction in EIT is a non linear a nd illposed inverse pr oblem. Uncertainties caused by these properties prevent EIT images from having the highresolution typically found in MRI and CT imaging. However, one may extr act a unique property (e lectric conductivity) of an imaged object in EIT, which is potentially valuable in clinical monitoring. EIT image reconstruction is a spatiallyvariant estimation pr oblem. In this paper, we attempted to reduce spatial variance of QI. Spatial variance is hard to deal with because its properties also depend on the regularization method chosen. Therefore, we sought a general approach that can reduce spatial variance in this paper. We proposed met hods that normalize full array and hemiarray EIT reconstructions using blurring properties calculated directly from the sensitivity matrix as a normalization term. This term is dimensionless, and unit of reconstructed quantities were PAGE 80 80 preserved after normalization. We found that the spatial variability observed was more dependent on the regularization scheme than the inverse model. Additionally, we found that patterns of the spatial variability were similar for different forward models. We chose the model that had elements of similar areas in order to simplify QI calculation. Figure 54. Radial curves of relative QI calculated using TSVD (A1) and normalization method (C5). The QI values are obtained from reconstructions for the full array (k = 17) and hemiarray (k = 16). The anomaly (100mL) was positioned on various locations. The curves were plotted for the A) full array , B) hemiarray C) hemiarray , and D) hemiarray axes. PAGE 81 81 We implemented normalization of spatiallyvar iant EIT image reconstruction using three different methods: pixelwise s caling (PWS), weighted pseudoi nversion (WPI), and weighted minimum norm method (WMNM). We found in our comparison that PWS methods generally worked poorly in reducing spatial variance. Employing WPI methods (A3 and A4) worked well for the full array, but not for the hemiarra y. The WMNMregularized reconstruction (B1) produced images of improved contrast in bot h full array and hemiarray cases. While B1 succeeded in reducing spatial variability for th e full array, it did not reduce variability for hemiarray cases. By using the WMNM framework, but with a differently defined normalization matrix as in C5 ( WQt) rather than B1 (W ), Sadleir et al. (2008) succeeded in decreasing spatial variability of QI. C5 also produced images that were generally superior to those produced using A1. In our investigation, we confirmed the results of Sadleir et al. (2008) in terms of QI. Further, we demonstrated that C5 performed fairly well for the full array as well. The images in Figure 53 (C5) were reconstructed from phantom meas urements, while the images shown in Sadleir et al. (2008) were reconstructed from noiseless simula ted measurements. Our images did not match the images in Sadleir (2008) closely for anomalies in the axis. Even though the condition number of the system matrix of B1 ( SW) was larger than that of C5 (SWQt) (Table 3), anomaly images reconstructed by B1 appeared less diffuse than those by C5 over all regions including the axis. We tested our methods for single anomalies loca ted at various positions within the domain. We have yet to determine if our method works well for multiple anomaly cases. Our normalization scheme was based on 2D and we obtai ned the blurring matrix from the sensitivity matrix of a 2D inverse model. We applie d the 2D normalization method successfully to cylindrical phantom data acquired using bar electr odes. The blur matrix can easily be obtained PAGE 82 82 from a 3D sensitivity matrix, and there is restriction on impl ementing our normalization method in 3D. Overall, we found that the best normalized QI values were obtained by applying a WMNM normalization method to WMNM regularization (C5) for both full array and hemiarray topologies. Our comparison showed that C5 was generally superior to other methods, even though there were test cases where it did not perform be tter than other candidates. Conclusion In this chapter, we com p ared various normalization methods in terms of QI. We demonstrated that some of these methods successfully reduced spatial vari ability of QI for the measurements obtained in numerically simulations and phantom experiments. Our investigation included many different approaches to resolv e spatial variability problem in EIT image reconstruction. Although method C5 generally pe rformed well, there was no single method that performed significantly better over both electrod e topologies. However, we have showed that discrete use of this t ype of normalization method should result in increased quantitation stability in EIT reconstruction. Further testing of these methods is necessary for more realistic geometries in 2D and 3D. Formulated in a straightforward manner, the methods should be simple to implement without great computational burden. Thus, they should be easily applicable to various image reconstruction problems that are sp atiallyvariant. We anticipate that this type of normalization method could also be used for highly as ymmetric image reconstruction problems. PAGE 83 83 CHAPTER 6 IMAGE RECONSTRUCTION BY WEIGHTED PSEUDOINVERSION METHOD IN HEMIARR AY EIT Introductio n In hem i array EIT, the input electrodes are placed only on one half of the object to be imaged. Measurement sensitivity varies de pending on the anomaly location, and some reconstructed anomalies were either invisibl e or erroneous when utilizing conventional regularization methods. For example, in the posterior region ( region) of the image domain, anomalies were almost invisible when TSVD rec onstruction was used (Figure 53). Therefore, methods to enhance image quality needs to be investigated further for 8electrode hemiarray EIT. We investigated Weighted PseudoInverse (WPI) method to improve image quality and anomaly localization. The weighting coefficients were composed of system blurring property, which were directly calculated from the sensitivity matrix. The method was tested for a 2D simulated model with circular boundary, containi ng a smaller disk anomaly located in various positions. The measurements obtained from the model were contaminated with AWGN of various amounts. The effect of the WPI method was studied, and compared with plain TSVD and WMNM reconstruction methods. We found that image quality was enhanced by using the WPI method, which resulted in improve ment in anomaly localization. Method Weighted PseudoInversion The m o tivation in our method is that impr oved images may be obtained by weighting reconstruction prior to pseudoinversion. We define entries of a diagonal weighting matrix P as in Equation 61. PAGE 84 84 1 1 2 N i ij jpB, j = 1, 2, N, (61) where pj values are reciprocals of square sum of B s columns. The advantages of considering properties of the blur matrix B as weighting coefficients are that they are dimensionless and can be precalculated from the sensitivity matrix S Then, we use WPI framework to modify the problem to reconstruct the image as in Equation 62. V SP (62)Singular Value Decomposition Using SVD, we can estimate basis (singular) images that de termine system blurring. In Figure 61, we compared the singular images of S (TSVD reconstruction) and SP (WPI reconstruction). Singular images of S fail to describe region, which may explain why anomalies in the region are invisible for TSVD reconstructi on. On the contrary, singular images of SP can roughly describe the entire region. In order to investigat e the effect of the WP I method, it may be worthwhile to examine the singular value spectra of the sensitivity matrices with and without wei ghting. In Figure 62, singular value spread of S (TSVD reconstruction), SW (WMNM reconstruction) and SP (WPI reconstruction) are compared. We found that SPs singular values decrea se the most rapidly as the mode number increases. Computer Simulation For our investigation, we used the hemiarra y forward models described in Chapter 5. Anomaly position was varied with respect to relative radius (r = 0, 0.2, 0.4, 0.6 and 0.8); and angle displacement ( 0o to 180o with 5o increment) (Figure 22). PAGE 85 85 Figure 61. Singular images of sensitivities matrices (S and SP). A) S and B) SP in eightelectrode hemia rray case. Singular modes from 1 to 20 are displayed from left to right a nd top to bottom. Rank of both matrices was 20. A B PAGE 86 86 Figure 62. Comparison of singular value spectra S, SP and SW in eightelectrode hemiarray case. Results Noiseless Simulated Measurement In Figure 63, images reconstructed from noi seless simulated measurements are displayed. Using the conventional TSVD reconstruction, we could estimate anomaly locations reasonably well for and regions (Figure 63 (A)). However, anomalies in the region diffused away as they moved further away from the electrodes. Using WMNM these anom alies could be roughly recovered (Figure 63 (B)). However, a relativ ely large artifact p eak appeared in the (r = 0.8) region, and errors in localizati on are seen for anomalies in the (r = 0.6 and 0.8), and (r = 0.6 and 0.8) regions. The WPI method shown in Equation 62 produced images of improved qualities (Figure 63 (C)) compared to th e TSVD and WMNM rec onstruction methods. PAGE 87 87 Figure 63 Images of a single a nomaly at 15 different locations reconstructed from noiseless m easurements using A) TSVD, B) WMNM and C) the WPI in eightelectrode hemiarray case (k =20). The first, second and third columns of each figure display anomalies on the and axes respectively. A B C PAGE 88 88 Figure 64 Images of a single a nomaly at 15 different locations r econstructed from noisy measurem ents (AWGN of 30dB) using A) TSVD, B) WMNM and C) the WPI in eightelectrode hemiarray case (k =20). The first, second and third columns of each figure display anomalies on the and axes respectively. A B C PAGE 89 89 Noisy Simulated Measurement For noisy 3D phantom measurements, k = 12 was chosen afte r the inspection of lcurves. Here, we used k = 12 for the simulated measurements corrupted with AWGN of 30dB. Figure 64 shows that TSVD and WMNM reconstructions produced inaccurate images, while the WPI method as in Equation 62 produced images of acceptable qualities, whic h shows a potential for reasonable localization. Localization of Anomaly The images were created from the measurem ents corrupted with AWGN of various SNR. Based on the images, anomalies were localized as follows. First, image pixel values that were greater than a half of the maximum pixel value were chosen as the half amplitude set. Then, center of mass of the half amplitude set was take n as the anomalys location. In Figure 65 and 66, we present the localization re sults of a single anomaly from th e measurements corrupted with AWGN of 40dB and 30dB SNR respectively. Th e discrepancy between the original and reconstructed locations is illustrated by arrows. The mean and variance of localization errors were about 0.07 and 1.436, and 0.11 and 6.174 respectively. Discussion and Conclusion The EIT image reconstruction is a spatiallyvar iant problem. In 8elect rode hemiarray EIT, some reconstructed anomalies were either invi sible or erroneous when plain TSVD and WMNM reconstruction methods were used. In this chapter, we considered the WPI method for image reconstruction. In the WPI method, sensitivity matrix S was multiplied with a weighting matrix P before pseudoinversion of the inverse problem as in Equation 62. By using the WPI framework, we expect to benefit from regularizing effect of the truncation, which reduces uncerta inties of high spatial frequency PAGE 90 90 components. By using the WPI met hod, our results showed that th e errors in the images were reduced. Consequently, anomaly localization beca me reasonable with relatively small errors. Figure 65. Localization of a single anomaly at various positions from simulated noisy measurements (AWGN of 40dB SNR). Imag es were reconstructed by the WPI (k = 13). Anomaly positions varied as a function of radius (r = 0, 0.2, 0.4, 0.6 and 0.8) and angle (: from 0o to 180o with 5o increment). Arrows were used to illustrate localization errors. In Table 61, condition numbers of sensitivity matrices with and without weighting are presented. In 8electrode hemiarray case, matrix SP had larger condition number than S and SW. Therefore, reconstruction by SP may be more susceptible to measurement noise. Table 61. Condition numbers of sensitivity matrices with and without weighting. Reconstruction matrix Condition number S+ 9.1117 W(SW)+ 4.8817 (SP)+ 1.1520 PAGE 91 91 Figure 66. Localization of a single anomaly at various positions from simulated noisy measurements (AWGN of 30dB SNR). Imag es were reconstructed by the WPI (k = 12). Anomaly positions varied as a function of radius (r = 0, 0.2, 0.4, 0.6 and 0.8) and angle (: from 0o to 180o with 5o increment). Arrows were used to illustrate localization errors. This investigation is limited in that only a single anomaly case was considered. Besides, we tested the method for AWGN. However, noise in the realistic measurements may be differently characterized. The method should be tested further, and its performance may be further improved by combining it with other regularization techniques. The method is straightforward to implement and should be widely applicable to other image reconstruction problems. We anticipate th at the method may be useful for any highly asymmetric image reconstruction problems. PAGE 92 92 CHAPTER 7 FUTURE WORK Estimation of Shape Change in 3D Object shapes are often assumed to be either cy lindrical or spherical in 3D EIT. If we can define a conformal relation between the real and the assumed shapes in 3D, a similar approache to the one we introduced in Chapter 4 can be implemented to estimate the shape change. The stereographic projection is a know n conformal mapping that projects a sphere onto a plane. We anticipate that this conformal ma pping would play an important role in estimating shape changes in 3D EIT. Normalization Method in EIT Image Reconstruction Test for Anomalies of Di fferent Conductivities Various normalization methods in EIT image reconstruction were introduced in Chapter 5, and applied to reducing QI erro rs caused by spatial variabilit y of the image reconstruction problem. The methods were tested for various volumes of an anomaly that has the same conductivity, which was regarded to simulate bleeding, perfusion or respiration. Testing for various conductivity values is s uggested for future work, simulating temporal changes (aged or malignant tissues) or situational changes (i.e. e ffect of delivered drugs) of tissues electrical properties. Test for Multiple Anomalies The normalization methods introduced in Chapter 5 is suggested to be tested for multiple anomaly cases. In lung imaging, it is beneficial to simultaneously monitor the liquid and air volumes, which can be seen to be anomalies of two different conduc tive properties. In the case of internal bleeding, this model can simulate bleeding from multiple sources. PAGE 93 93 Normalization in 3D The normalization methods proposed in Chapter 5 are generic and simple to implement. We tested the methods to 2D full array and hemi array electrode configurations. However, it can be readily tested for the 2D and 3D measurements obtained fr om highly asymmetrical objects. Potential areas include: subsurface imaging (K otre 1996), imaging usi ng a rectangular array (Mueller et al. 1999), and 3D cerebral imaging (Tidswell et al. 2001, Sadleir and Tang 2009). PAGE 94 94 APPENDIX A MATHEMATICAL FORMULATION Least Squares Solution We will demonstrate how to obtain least squares solution of the inverse problem of the Equation 35. The least squares problem is an optimization problem, where one estimates the solution by an approximate that minimizes the re sidual norm as in Equation A1. 2 2minSV (A1) First, we take the derivative of the quantity with respect to as in Equation A2. 2 2SV (A2) Then, we get Equation A3 and Equation A4. VTSS (A3) VT TSSS (A4) By putting this quantity equal to zero, we get Equation A5. VT T SSS (A5) Finally, by inverting the square matrix STS, we can obtain the least squares solution expressed in algebraic form in Equation A6. VT TSSS1 (A6)TSVD Solution (Full Rank) Full rank TSVD solution is equivalent to the l east squares solution. For full rank truncation, D+ equals D1. Thus, the TSVD solution becomes as in Equation A7. VTUVD1 (A7) The least squares solution from Equation A6 is calculated using pseudoinverse of S as in Equation A8. PAGE 95 95 T TSSS1 (A8) If we apply SVD to S, we get Equation A9. T T T T TUDV UDV UDV 1 (A10) T T T TUDV UDV VDU1 (A11) T T TUDVVVD1 2 (A12) T TVDUVVD2 (A13) TUVD1 (A14) Finally, this is identical to the fu ll rank TSVD regularization inverse of S. Weighted Minimum Norm Method The method is motivated by that S can also be written as S(WW1). V S (A15) V 1WWS (A16) V 1WSW (A17) Then, the weighted sensitivity matrix SW is inverted instead of S as in Equation A18. V SW W1 (A18) Finally, we get the final form of th e WMNM solution as in Equation A19. V SWW (A19) Therefore, the overall change is modi fication of the reconstruction matrix S+ to W(SW)+. Tikhonov Regularization We will demonstrate how to obtain an al gebraic form of the solution by Tikhonov regularization. In the standard zeroorder Tikhonov re gularization, one estimates the solution by minimizing the weighted sum of residual norm and solution norm as in Equation A20. PAGE 96 96 2 2 2 2minarg S V (A20) Taking the derivative of the quantity with respect to gives us Equation A21. 2 2 2 2S V (A21) VTSS (A22) VT TSSS (A23) By putting the derivative equal to zero, we get Equation A24. VT T S SS (A24) VT T SISS (A25) Finally, inverting the square matrix ( STS +I ), whose condition can be improved by choice of we get the solution as in Equation A26. VT TSISS1 (A26) Joukowski Transformation of a U n it Circle Joukowski transformation from Equation 41 is rewritten here. Z cZW 1 (A27) Complex number Z in Cartesian coordinates can be expressed as in Equation A28. jyxZ (A28) where x and y are real numbers. In polar form, sin cos jrZ, (A29) where r is radial distance from Z to the origin, and is the radian. Similarly, complex number W can be written as in Equation A30. PAGE 97 97 jVUW (A30) For unit circle r = 1, Z = cos( ) + j sin( ). After the Joukowski tr ansformation, points on the unit circle will be W as in Equation A31. sin cos sin cosj c j W (A31) sin cos sin cos j cj (A32) sin cos sin cos jcj (A33) sin1cos1 cjc (A34) sin1cos1 cjc (A35) From Equation A35, we get U and V as a function of (Equations A36 and A37). cos1 c U (A36) sin1 cV (A37) U and V satisfy Equation A38. 1 112 2 c V c U (A38) Therefore, W represents an ellipse whose semimajor and semiminor axes are 1+ c and 1c respectively. PAGE 98 98 APPENDIX B MATLAB/COMSOL SCRIPTS Solve.m function [fem,map] = solve(i nput,inputNeg,angle,anoIndex) % script designs and solves a forward model fo r given input and specified anomaly location % input & inputNeg : pair of electrodes used for current injection % angle : angle of the axis where anomaly li es on. zero corresponds to the positive y axis. % anoIndex : 1,2,3,4,5 (r = 0, 0.2, 0.4, 0.6, 0.8) % fem : solution (voltage field) structure flclear fem clear vrsn vrsn.name = 'COMSOL 3.4'; vrsn.ext = ''; vrsn.major = 0; vrsn.build = 248; vrsn.rcs = '$Name: $'; vrsn.date = '$Date: 2007/10/10 16:07:51 $'; fem.version = vrsn; % Geometry g1=circ2('1','base','center','pos',{'0','0'},'rot','0'); g2=square2('.1','base','corner','pos',{'.05','.05'},'rot','0'); g2=move(g2,[1,0]); g3=geomcomp({g1,g2},'ns',{'C1','SQ 1'},'sf','SQ1C1','edge','none'); g4=circ2('1','base','center','pos',{'0','0'},'rot','0'); [g5]=geomcopy({g3}); [g6]=geomcopy({g5}); g6=move(g6,[0,0]); g6=rotate(g6,0.7853981633974483,[0,0]); [g7,g8]=geomcopy({g3,g6}); [g9,g10]=geomcopy({g7,g8}); g9=move(g9,[0,0]); PAGE 99 99 g10=move(g10,[0,0]); g9=rotate(g9,1.5707963267948966,[0,0]); g10=rotate(g10,1.5707963267948966,[0,0]); [g11,g12]=geomcopy({g9,g10}); [g13,g14]=geomcopy({g11,g12}); g13=move(g13,[0,0]); g14=move(g14,[0,0]); g13=rotate(g13,1.5707963267948966,[0,0]); g14=rotate(g14,1.5707963267948966,[0,0]); [g15,g16]=geomcopy({g13,g14}); [g17,g18]=geomcopy({g15,g16}); g17=move(g17,[0,0]); g18=move(g18,[0,0]); g17=rotate(g17,1.5707963267948966,[0,0]); g18=rotate(g18,1.5707963267948966,[0,0]); g19=circ2('0.1','base','center','pos',{'0','0'},'rot','0'); [g20]=geomcopy({g19}); [g21]=geomcopy({g20}); g21=move(g21,[0,0.2]); [g22]=geomcopy({g21}); [g23]=geomcopy({g22}); g23=move(g23,[0,0.2]); [g24]=geomcopy({g23}); [g25]=geomcopy({g24}); g25=move(g25,[0,0.2]); [g26]=geomcopy({g25}); [g27]=geomcopy({g26}); g27=move(g27,[0,0.2]); g19=rotate(g19,angle*1.5707963267948966/90,[0,0]); g21=rotate(g21,angle*1.5707963267948966/90,[0,0]); g23=rotate(g23,angle*1.5707963267948966/90,[0,0]); g25=rotate(g25,angle*1.5707963267948966/90,[0,0]); PAGE 100 100 g27=rotate(g27,angle*1.5707963267948966/90,[0,0]); % Geometry objects clear s s.objs={g3,g4,g6,g9,g10,g15,g16,g17,g18,g19,g21,g23,g25,g27}; s.name={'CO1','C1','CO2','CO3','CO4','CO5','CO6','CO7','CO8','C2', ... 'C3','C4','C5','C6'}; s.tags={'g3','g4','g6','g9','g10','g15','g16','g17','g18','g19', ... 'g21','g23','g25','g27'}; fem.draw=struct('s',s); fem.geom=geomcsg(fem); % Initialize mesh fem.mesh=meshinit(fem, ... 'Report','off', ... 'hauto',5); % Application mode 1 clear appl appl.mode.class = 'ConductiveMediaDC'; appl.assignsuffix = '_dc'; clear bnd bnd.Jn = {0,0,1,1}; bnd.type = {'nJ0','cont','nJ','nJ'}; nbnd = flgeomnbs(fem.geom); bnd.ind = [1*ones(1,nbnd)]; bdd = find(bnd.ind == 1); arrowArrival = []; arrowDisplacement = []; for i = 1:length(bdd) [xy,dxy] = flgeomed(fem.geom,bdd(i),1); arrowArrival = [arrowArrival xy]; arrowDisplacement = [arrowDisplacement dxy]; PAGE 101 101 end arrowDeparture = arrowArrival (arrowDisplacement); leng = abs(complex(arrowDisplaceme nt(1,:),arrowDisplacement(2,:))); bnd.ind(bdd(find(abs(complex(arrowDeparture(1,:),arrowDeparture(2,:))) < 0.95))) = 2; temp1 = bdd(find(abs(complex(arrowArrival(1,:),arrowArrival(2,:))) > 0.95)); temp2 = bdd(find(abs(complex(arrowArrival(1,:),arrowArrival(2,:))) <= 1.05)); temp3 = bdd(find(abs(complex(arrowDeparture(1,:),arrowDeparture(2,:))) <= 1.05)); temp4 = bdd(find(leng <0.1)); temp = intersect(intersect(intersect(temp1,temp2),temp3),temp4); bnd.ind(temp) = 2; ebnd = [24 19 13 7 1 5 11 18]; bnd.ind(ebnd(input)) = 3; bnd.ind(ebnd(inputNeg)) = 4; appl.bnd = bnd; clear equ equ.sigma = {5.99e7,1,2}; equ.ind = [1,1,1,1,1,1 ,1,1,2,2,2,2,2,2]; if angle == 0 anoInd = [10,11,12,13,14]; else anoInd = [14,13,12,11,10]; end equ.ind(anoInd(anoIndex)) = 3; appl.equ = equ; fem.appl{1} = appl; fem.frame = {'ref'}; fem.border = 1; clear units; units.basesystem = 'SI'; fem.units = units; PAGE 102 102 % ODE Settings clear ode clear units; units.basesystem = 'SI'; ode.units = units; fem.ode=ode; % Multiphysics fem=multiphysics(fem); % Extend mesh fem.xmesh=meshextend(fem, ... 'Report','off'); % Solve problem fem.sol=femlin(fem, ... 'solcomp',{'V'}, ... 'outcomp',{'V'}, ... 'Report','off'); Solve_batch.m function [fem] = solve_batch() % script solves forward problems for all inputs prescribed % fem : solution (voltage field) structure nE = 8; for angle = [0:5:180] for kk = 1:5 for ii = 1:nE jj = mod(ii,8)+1; fem{ii,jj} = solve(ii,jj,angle,kk); cd forward eval(['save sol' num2str(angle) '_' num2str(kk) '.mat fem']); PAGE 103 103 cd .. end end end Measure.m function [v] = measure(fem,P) % script returns voltage values of specified locations for a given solution structure % fem : solution structure % P : points of measurements fem.xmesh = meshextend(fem); v = postinterp(fem, 'V', P); Measure_all_8.m function [Dv] = measure_all_8(fem) % script measures voltages from loaded fem structure in 8electrode case. % fem : solution structure % Dv : boundary voltage measurements nE = 8; clear i; Pcomplex = exp(i*2*pi/8.*[0:nE1]); x=real(Pcomplex);y=imag(Pcomplex);P = [x;y]; % P: electrode locations clear temp; index = []; for ii = 1:nE temp = sort(mod([1:nE3]+ii,nE)+1); index = [index; temp]; end Dv = zeros(1,nE*(nE3))'; for k = 1:nE PAGE 104 104 inPos = k; inNeg = mod(k,nE)+1; femIn = fem{inPos,inNeg}; [V] = measure(femIn,P); V = [V V(1)]; dVV = diff(V); outV = index(k,:); dV = dVV(outV); Dv([(nE3)*(k1)+1:(nE3)*k]) = dV'; end Check_sol.m function [M] = check_sol(fem,speed) % script plays movie file showing all solutions. % fem : solution structure % speed : play speed of the slideshow nE = 8; j = [1:nE]; i = [2:nE]; i = [i 1]; kk = 0; for k = 1:nE postplot(fem{j(k),i(k)},'tridata','V', ... 'flowdata',{'Jx_dc','Jy_dc'}, ... 'flowcolor',[1.0,0.0,0.0], ... 'flowlines',15, ... 'title','Streamline: Total current density [A/m^2]', ... 'axis',[2.08,2.08,1.2,1.2,1,1]) kk = kk+1; M(kk) = getframe; end movie(M,1,speed) Sensitivity.m % script is used to calcu lated sensitivity matrix PAGE 105 105 load sMesh.mat nsubd = 344; % area equ.init = {0}; equ.ind = ones(1,nsubd); appl.equ = equ; fem.appl{1} = appl; fem.frame = {'ref'}; fem.border = 1; fem.units = 'SI'; fem=multiphysics(fem); fem.xmesh=meshextend(fem); init = asseminit(fem); fem.sol = init; A = zeros(1,nsubd); for j = 1:nsubd A(j)=postint(fem,'1', 'dl',[j]); end save sMeshArea.mat A fem0 = fem; nsubd = flgeomnmr(fem0.geom); nbnd = flgeomnbs(fem0.geom); cd forward load solhomo.mat, cd .. nE = 8; for i = 1:nsubd pd(i) = posteval(fem0,'1','dl',i); vert{i} = pd(i).p; fac{i} = pd(i).t; % just in case end clear pd PAGE 106 106 numb = [1:nE]; for n = 1:nsubd disp('element'); disp(n); for k = 1:nE Pos = k; numbe = setdiff(numb,[Pos]); Neg = mod(Pos,nE)+1; f = fem{Pos,Neg}; f.xmesh = meshextend(f); Ex = postinterp(f,'Ex_dc',vert{n}); Exx{k,n} = Ex; Ey = postinterp(f,'Ey_dc',vert{n}); Eyy{k,n} = Ey; clear Ex, clear Ey end end save Exx.mat Exx; save Eyy.mat Eyy S = zeros(nE*(nE3),nsubd); clear temp; temp = ones(1,nE3); inputV = [temp 2*temp 3*temp 4*temp 5*temp 6*temp 7*temp 8*temp]; clear temp outputV = []; for i = 1:nE temp = sort(mod([1:nE3]+i,nE)+1); outputV = [outputV temp]; end for n = 1:nsubd for i = 1:nE*(nE3) input = inputV(i); output = outputV(i); S(i,n) = meshintegrate(vert{n},fac{n}, ... PAGE 107 107 dot([Exx{input,n}; Eyy{input,n}], ... [Exx{output,n}; Eyy{output,n}])); end end Sfull = S; save Sfull.mat Sfull Reconstruct.m % script is used to obtain TSVD image reconstruction. clear load Sfull8.mat; S = Sfull; load b_0.mat k = 16; [U,s,V] = svd(S); sinv = zeros(size(s)); ss = svds(S,k); for ii = 1:k sinv(ii,ii) = 1/ss(ii); end Sinv = V*sinv'*U'; x_TSVD = Sinv*b; PAGE 108 108 LIST OF REFERENCES Adler A, Berthiaum e Y, Guardo R and Amyot R 1995 Imaging of pulmonary edema with electrical impedance tomography Proc. 17th Annu. Int. Conf. IEEE EMBS, Montreal 17 557558 Adler A and Guardo R 1996 Electrical impedan ce tomography: Regularized imaging and contrast detection IEEE Trans. Med. Imag. 15 170179 Adler A, Guardo R and Berthiaume Y 1996 Impedan ce imaging of lung ventilation: Do we need to account for chest expansion? IEEE Trans. Biomed. Eng. 43(4) 414420 Adler A, Amyot R, Guardo R, Bates J H T a nd Berthiaume Y 1997 Monitoring changes in lung air and liquid volumes with electrical impedance tomography J. Appl. Physiol. 83 17621767 Adler A, Shinozuka N, Berthiaume Y, Guar do R and Bates J H T 1998 Electrical impedance tomography can monitor dynamic hyperinflation in dogs J. Appl. Physiol. 84 726732 Alessandrini G and Rondi L 1998 Stable determin ation of a crack in a planar inhomogeneous conductor SIAM J. Math. Anal. 30 326340 Avis N J and Barber D C 1992 Single step al gorithms for image reconstruction. II. Generalization of the Sheffield filtered back projection algorithm. [Electric impedance tomography] EIT/APT, IEE Colloquium on 5 13 Avis N J and Barber D C 1994 Image reconstruc tion using nonadjacent dr ive configurations Physiol. Meas. 15 A153A160 Barber D C and Brown B H 1984 Applied potential tomography J. Phys. E: Sci. Instrum. 17 723733 Barber D C and Seagar A D 1987 Fast reconstruction of resistance image Clin. Phys. Physiol. Meas. 8 A47A54 Barber D C and Brown B H 1988 Errors in recons truction of resistivity images using a linear reconstruction technique Clin. Phys. Physiol. Meas. 9 A101A104 Barber D C 1990 Quantification of impedance imaging Clin. Phys. Physiol. Meas. 11 A45A56 Bayford R H, Boone K G, Ha nquan Y and Holder D S 1996 Improvement of the positional accuracy of EIT images of the head using a Lagrange multiplier reconstruction algorithm with diametric excitation Physiol. Meas. 17 A49A57 Blott B H, Daniell G J and Meeson S 1998 Elect rical impedance tomography with compensation for electrode positioning variations Phys. Med. Biol. 43 17311739 PAGE 109 109 Blott B H, Cox S J, Daniell G J, Caton M J a nd Nicole D A 2000 High fidelity imaging and high performance computing in nonlinear EIT Physiol. Meas. 21 713 Boone K G and Holder D S 1995 Design considerati ons and performance of a prototype system for imaging neuronal depolarization in the brai n using direct current electrical resistance tomography Physiol. Meas. 16(3A) A87A98 Borsic A, Comina C, Foti S, Lancellotta R and Musso G 1995 Imaging heterogeneities with electrical impedance tomogr aphy: laboratory results Geotechnique 55(7) 539547 Boyle A, Lionheart W R B, GmezLaberge C and Adler A 2008 Evaluating deformation corrections in electri cal impedance tomography IXth Conf. on EIT, Dartmouth college, New Hampshire Breckon W R and Pidcock M K 1988 Data errors a nd reconstruction algorithms in electrical impedance tomography Clin. Phys. Physiol. Meas. 9 A105A109 Brown B H, Leathard A, Sint on A, McArdle F J, Smith R W M and Barber D C 1991 Blood flow imaging using electrical impedance tomography Ann. Intl. Conf. IEEE EMBS 13(1) 307308 Campbell J H, Harris N D, Zhang F, Brown B H and Morice A H 1994 Clinical applications of electrical impedance tomography in the mon itoring of changes in intrathoracic fluid volumes Physiol. Meas. 15 A217A222 Cheney M, Isaacson D, Newell J C, Simske S and Goble J 1990 NOSER: An algorithm for solving the inverse conductivity problem Intl.J.Imag Sys.Technol. 2 6675 Cherepenin V, Karpov A, Korjenevsky A, Ko rnienko V, Mazaletskaya A, Mazourov D and Meister D 2001 A 3D electrical impedance tomography (EIT) system for breast cancer detection Physiol. Meas. 22 918 Clay M T and Ferree T C 2002 Weighted regulariza tion in electrical impedance tomography with applications to acute cerebral stroke IEEE Trans. Med. Imag. 21 629637 CohenBacrie C, Goussard Y and Guardo R 1997 Regularized reconstruc tion in electrical impedance tomography using a vari ance uniformization constraint IEEE Trans. Med. Imag. 16 562571 Conway J 1987 Electrical impedance tomography for thermal monitoring of hyperthermia treatment: an assessment using in vitro and in vivo measurements Clin. Phys. Physiol. Meas. 8(A) 141146 Davalos R V, Otten D M, Mir L M and Rubi nsky B 2004 Electrical impedance tomography for imaging tissue electroporation IEEE Trans. Biomed. Eng. 51(5) 761767 Dehghani H, Soni N, Halter R, Hartov A a nd Paulsen K D 2005 Excitation patterns in threedimensional electrical impedance tomography Physiol. Meas. 26 S185S197 PAGE 110 110 Dickin F and Wang M 1996 Electrical resist ance tomography for process tomography Meas. Sci. Tech. 7 247260 Dijkstra A M, Brown B H, Leathard A D, Ha rris N D, Barber D C, Edbrooke D L 1993 Review clinical applications of el ectrical impedance tomography J. Med. Eng. Tech. 17(3) 8998 Dong G, Liu H, Bayford R H, Yerworth R, Gao S, Holder D and Yan W 2004 The spatial resolution improvement of EIT im ages by GVSMPFOCUSS algorithm Physiol. Meas. 25 209225 Dong G, Liu H, Bayford R H, Yerworth R, Schimpf P H and Yan W 2005 Spatial resolution improvement of 3D EIT images by th e shrinking sLORETAFOCUSS algorithm Physiol. Meas. 26 S199S208 Edic P M, Saulnier G J, Newell J C and Is aacson D 1995 A realtime electrical impedance tomography IEEE Trans. Biomed. Eng. 42 849859 Eyboglue B M, nert F A, Bays al U, Bibert Ihsan Keyf A, Yilmaz and Erdogan Y 1995 Application of electrical impe dance tomography in diagnosis of emphysema a clinical study Physiol. Meas. 16(3A) A191A211 Frerichs I, Hahn G and Hellige G 1999 T horacic electrical impedance tomographic measurements during volume controlled ventilationeffects of tidal volume and positive endexpiratory pressure IEEE Trans. Med. Imag. 18(9) 764773 Frerichs I 2000 Electrical impedance tomography (EIT) in applications related to lung and ventilation: a review of experi mental and clinical activities Physiol. Meas. 21 R1R21 Geselowitz D B 1971 An applica tion of electrocardiographic lead theory to impedance plethysmography IEEE Trans. Biomed. Eng. 18 38. Ghahary A and Webster J G 1989 Electrical safety for an electrical impedance tomography IEEE EMBS 11th annual intl. conf. 461462 Golub G H and Reinsch C 1970 Handbook series lin ear algebra: singular value decomposition and least squares solutions Numer. Math. 14 403420 GmezLaberge and Adler A 2006 Imaging of elect rode movement and conductivity change in electrical impedance tomography Proc. 19th Can. Conf. Electrical and Computer Eng. (Ottawa, Canada) 19 975978 Gorodnitsky I F and Rao B D 1997 Sparse signa l reconstruction from limited data using FOCUSS: A reweighted minimum norm algorithm IEEE Sig. Proc. 45(3) 600616 Halter R J, Hartov A, Heaney J A, Paulse n K D and Schned A R 2007 Electrical impedance spectroscopy of th e human prostate IEEE Trans. Biomed. Eng. 54(7) 13211327 Hansen P C 1987 The truncated SVD as a method for regularization BIT 27 543553 PAGE 111 111 Hansen P C 1992 Analysis of discrete ill posed problems by means of the Lcurve SIAM Review 34(4) 561580 Hansen P C and OLeary D P 1993 The use of the Lcurve in the regulariza tion of discrete illposed problems SIAM J. Sci. Comput. 14(6) 14871503 Hansen P C 1994 Regularization tools: A Matlab package for analys is and solution of discrete illposed problems Numerical Algorithms 6 135 Heikkinen L M, Vilhunen T, West R M and Vauhkonen M 2002 Simultaneous reconstruction of electrode contact impedances and internal el ectrical properties, Part II: applications Meas. Sci. Technol. 13 18551861 Heikkinen L M, Kourunen J, Savolainen T, Vauhkonen P J, Kaipio J P and Vauhkonen M 2006 Real time threedimensional electrical impeda nce tomography applied in multiphase flow imaging Meas. Sci. Technol. 17 20832087 Holder D S 2005 Electrical impedance tomography: Methods, histor y and applications Institute of Physics publishing Jianhua Li 1994 A method of reducing the er ror caused by boundary shape and electrode positions in electrical impedance tomography Physiol. Meas. 15 A169A174 Kao TJ, Isaacson D, Newell J C and Saulnier G J 2006 A 3D reconstruction algorithm for EIT using a handheld probe fo r breast cancer detection Physiol. Meas. 27 S1S11 Keshtkar A and Keshtkar A 2007 Measured and modelled electrical bioimpedance inside the human normal and malignant bladder epithelium International J. Biomed. Eng. Tech. 1(2) 127133 Kiber M A, Barber D C and Brown B H 1990 Estimation of object boundary shape from the voltage gradient measurements Proc. Mtg. EIT Copenhagen 5259 Kolehmainen V, Lassas M and Ola P 2005 The invers e conductivity problem with an imperfectly known boundary SIAM J. Appl. Math. 66(2) 365383 Kotre C J 1996 Subsurface electrical impedance im aging using orthogonal linear electrode arrays IEE Proc. Sci. Meas. Technol. 143 4146 Lionheart W R B, Lidgey F J, McLeod C N, Paulson K S, Pidcock M K and Shi Y 1997 Electrical impedance tomography for high speed chest imaging Physica. Medica. 13 Suppl. 1 Lionheart W R B 1998 Boundary shape and electrical impedance tomography Inverse Problems 14 139147 Lionheart W R B 2004 EIT reconstruction algorithms: Pitfalls, challenges and recent developments Physiol. Meas. 25 125142 PAGE 112 112 Ljaz U Z, Khambampati A K, Lee J S, Kim S and Kim K Y 2008 Nonstationary phase boundary estimation in electrical impedance to mography using unscented Kalman filter J. of Comput. Phys. 227(15) 70897112 Malmivuo J and Plonsey R 1995 Bioelectromagnetism: principles and applications of bioelectric and biomagnetic fields New york, Oxford university press. Mangnall Y F, Kerrigan D D, J ohnson A G and Read N W 1991 A pplied potential tomography: Noninvasive method for measuring gast ric emptying of a solid test meal Dig. Dis. Sci. 36(12) 16801684 McEwan A L and Holder D S 2007 Battery powered and wireless electrical impedance tomography spectroscopy imaging using Bluetooth IFMBE Proc. 11th Mediterranean Conf. on Medical and Biomedical E ngineering and Computing 2007 16 798801 Meeson S, Killingback L T and Blott B H 1995 The dependence of EIT images on the assumed initial conductivity di stribution: a study of pelvic imaging Phys. Med. Biol. 40 643657 Metherall P, Barber D C, Smallwood R H and Brown B H 1996 Threedi mensional electrical impedance tomography Nature 380(11) 509512 Molebny V V, Lionheart W R B, Vovk P P, Mykytenko Y H and Gouz V I 1996 Sensor position measurement for electroimpedance tomography Med. Biol. Eng. Comput. 34(1) Part 1 Morucci JP and Rigaud B 1996 Bi oelectrical impedance techni ques in medicine. Part III: impedance imaging. Third section: medical application Crit. Rev. in Biomed. Eng. 24(46)655677 Moskowitz M J, Ryan T P, Paulsen K D and M itchell S E 1995 Clinical implementation of electrical impedance tom ography with hyperthermia Intl. J. Hyperthermia 11(2) 141149 Mueller J L, Isaacson D and Newell J C 1999 A reconstruction algorithm for electrical impedance tomography data collected on rectangular electrode arrays IEEE Trans. Biomed. Eng. 46(11) 13791386 Murai T and Kagawa Y1985 Electrical impedan ce computed tomography based on a finite element model IEEE Trans. Biomed. Eng. 32 (3) 177184 Murphy D, Burton P, Coombs R, Tarassenko L and Rolfe P 1987 Impedance imaging in the newborn Clin. Phys. Physiol. Meas. 8A 131140 Murpy S C and York T A 2006 Electrical impedance tomography with nonst ationary electrodes Meas. Sci. Technol. 17 30423052 Nobes D C 1996 Troubled waters: Environmental applications of electri cal and electromagnetic methods Surveys in Geophysics 17 (4) 393454 PAGE 113 113 Nour S, Mangnall Y F, Dickson J A, Johns on A G and Pearse R G 1995 Applied potential tomography in the measurement of gastric emptying in infants J. Pediatr. Gastroenterol. Nutr. 20(1) 6572 Oh S, Tang T and Sadleir R J 2007 Quantitative analysis of shape change in Electrical Impedance Tomography (EIT) XIIIth conf. on Electrical bioimpedance and VIIIth conf. on Electrical Impedance tomography (ICEBI 07), IFMBE Proc. 17 424427 Oh S and Sadleir R J 2007 Compensating spatial variability of Quantity Index in 2D electrical impedance tomography: Comsol Multiphysics study COMSOL Conf. Boston, Massachusetts Oh S, Tang T, Tucker A S and Sadleir R J 2009 Normalization of a sp atiallyvariant image reconstruction problem in electrical impedance tomogr aphy using system blurring properties Physiol. Meas. 30 115 Oh T I, Woo E J and Holder D 2007 Multifreque ncy EIT system with radially symmetric architecture: KHU Mark1 Physiol. Meas. 28 S183S196 Pomerantz M, Delgado F and Ei seman B 1970 Clinical evaluation of transthoracic electrical impedance as a guide to in trathoracic fluid volumes Ann. Surg. 171(5) 686694 Powell H M, Barber D C and Freeston I L 1987 Impedance imaging using linear electrode arrays Clin. Phys. Physiol. Meas. 8(A) 109118 Pretty I A 2006 Caries detection and diagnosis: n ovel technologies J. of Dentistry 34(10) 727739 Romsauerova A, McEwan A, Horesh L, Yerworth R, Bayford R H and Holder D S 2006 Multifrequency electrical impedance tomography (EIT) of the adult human head: initial findings in brain tumors, arteriovenous malformations and chronic stroke, development of an analysis method and calibration Physiol. Meas. 27 S147S161 Sadleir R J and Fox R A 1998 Quantification of blood volume by electrical impedance tomography using a tissueequivalent phantom Physiol. Meas. 19 501516 Sadleir R J and Fox R A 2001 Detection and quan tification of intraper itoneal fluid using electrical impedance tomography IEEE. Biomed. Eng. 48(4) 484491 Sadleir R J, Zhang S U, Tucker A S and Oh Sungho 2008 Imaging and quantification of anomaly volume using an eight electrode h emiarray EIT reconstruction method Physiol. Meas. 29 913927 Sadleir and Tang 2009 Electrode conf igurations for detection of in traventricular haemorrhage in the premature neonate Physiol. Meas. 30 6379 Sha L, Ward E R and Stroy B 2002 A review of dielectric properties of normal and malignant breast tissue Proc. IEEE SoutheastCon. 457462 PAGE 114 114 Shini M and Rubinsky B 2007 Multiple biopsy pr obe sampling enabled minimally invasive electrical impedance tomography Physiol. Meas. 29 109126 Soleimani M, GmezLaberge C and Adler A 2006 Imaging of conductivity changes and electrode movement in EIT Physiol. Meas. 27 S103S113 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 Physiol. Meas. 25 301314 Taktak A, Record P M Gadd R and Ro lfe P 1995 Neonatal lung imaging using EIT Proc. IX Int. Conf. on Electrical BioImpedance, Heidelberg, Germany Heidelberg: University of Heidelberg 458459 Tang M, Wang W, Wheeler J, McCormick M and Dong X 2002 The number of electrodes and basis functions in EIT image reconstruction Physiol. Meas. 23 129140 Tang T, Zhang S U and Sadleir R J 2006 A porta ble 8electrode EIT measurement system Proc. VIIth Conference on Biomedic al Applications of EIT 190193 Thomas D C, SiddallAllum J N, Sutherla nd I A and Beard R W 1994 Correction of the nonuniform spatial sensitivity of elec trical impedance tomography images Physiol. Meas. 15 A147A152 Tidswell T, Gibson A, Bayford R and Holder D 2001 Threedimensional electrical impedance tomography of human brain activity Neuroimage. 13 283294 Tossavainen OP, Vauhkonen M, Heikkinen L M and Savolainen T 2004 Estimating shapes and free surfaces with electrical impedance tomography Meas. Sci. Technol. 15 14021411 Tucker A S, Sadleir R J, Oh Sungho and Tang Te 2008 Portable eightelec trode EIT system for detection and quantificati on of abdominal hemorrhage Electrical Impedance Tomography Conference, Hanover, New Hampshire Vauhkonen M, Vadsz D, Karj alainen P A, Somersalo E and Kaipio J P 1998 Tikhonov regularization and prior information in electrical im pedance tomography IEEE Trans. Med. Imag. 17(2) 285293 Vilhunen T, Heikkinen L M, Sa volainen T, Vauhkonen P J, Lappalainen R, Kaipio J P and Vauhkonen M 2002 Detection of faults in re sistive coatings w ith an impedancetomographyrelated approach Meas. Sci. Technol. 13 865872 Vonk Noordegraaf A, Kunst P W, Janse A, Smulders R A, Heethaar R M, Postmus P E, Faes T J and deVries P M 1997 Validity a nd reproducibility of electrica l impedance tomography for measurement of calf blood fl ow in healthy subjects Med. Biol. Eng. Comput. 35(2) 107112 PAGE 115 115 Walker D C, Brown B H, Hose D R and Smallwood R H 2000 Modelling the electrical impedivity of normal and premalignant cervical tissue IEEE Electronics letters 36(19) 16031604 Wheeler J L, Wang W and Tang M 2002 A comparis on of methods for measurement of spatial resolution in twodimensional circular EIT images Physiol. Meas. 23 169176 Xu C, Dong X, Fu F, Shuai W, Liu X and Zhang C 2007 A nove l image monitoring software system of electrical impedance tomography for internal hemorrhage IFMBE Proc. World Congress on Medical Physics and Biomedical Engineering 2006 14 38823885 Zadehkoochak M, Blott B H, Hames T K and Geor ge R F 1991 Spectral expansion analysis in electrical impedance tomography J. Phys. D: Appl. Phys. 24 19111916 Zlochiver S, Arad M, Radai M M, BarakShin ar D, Krief H, Engelman T, BenYehuda R, Adunsky A and Abboud S 2007 A portable bioimpedance system for monitoring lung resistivity Med. Eng. Phys. 29 93100 PAGE 116 BIOGRAPHICAL SKETCH Sungho Oh was born in 1975 in W o nju, South Ko rea. He graduated from Daewon Foreign Language High school with English and Spanish ma jor. He entered Hanyang University with radio science and engineering major. From 1996 to 1998, he served the Republic of Korea (ROK) Army as a Korean Augmentee To the US Army (KATUSA) in Youngsan, Korea. He graduated with a bachelors degree in electrical en gineering in 2000. Aft er graduation, he worked in marketing department, Mobile Internet Business Unit, LOCUS Co rporation. The next year, he moved to the USA and began his gra duate study in electrica l engineering at the University of Florida. He earned his masters degree in 2004 with a thesis: A new quality measure in electrocardiogram In 2005, he began his PhD study in biomedical engineering at the University of Florida. He earned his PhD in biom edical engineering in 20 09 with a dissertation: Compensation of shape change artifacts and spatia llyvariant image reconstruction problem in Electrical Impedance Tomography 