# Microwave imaging for two-dimensional electrical property distribution profiling

код для вставкиСкачатьINFORMATION TO USERS This manuscript has been reproduced from the microfilm master. UMI films the text directly from the original or copy submitted. Thus, some thesis and dissertation copies are in ^pewriter face, while others may be from any type of computer printer. The quality of this reprodnction is dependoit upon the quality of the copy submitted. Broken or indistinct print, colored or poor quality illustrations and photographs, print bleedthrough, substandard marginc, and inq)rqper alignment can adversely afEect reproduction. In the unlikely event that the author did not send UMI a complete manuscript and there are missing pages, these will be noted. Also, if unauthorized copyright material had to be removed, a note wiD indicate the deletion. Oversize materials (e.g., maps, drawings, charts) are reproduced by sectioning the original, beginning at the upper left-hand corner and continuing from left to right in equal sections with small overlaps. Each original is also photographed in one exposure and is included in reduced form at the back of the book. Photogr^hs inchided in the original manuscript have been reproduced xerographically in this copy. Higher quality 6" x 9" black and white photographic prints are available for any photogr^hs or illustrations appearing in this copy for an additional charge. Contact UMI directly to order. A Bell & Howell information Company 300 North Zeeb Road. Ann Arbor. Ml 48106-1346 USA 313.'761-4700 800.521-0600 MICROWAVE IMAGING FOR 2-D ELECTRICAL PROPERTY DISTRIBUTION PROFILING A Thesis Submitted to the Faculty in partial fulfillment of the requirements for the degree of Doctor of Philosophy by PAUL MORGAN MEANEY, JR. Thayer School of Engineering Dartmouth College Hanover, New Hampshire MAY 1995 Examining Committee: Chairman: Member: Member: Member: /^ ^ Member: Dean of Graduate Studi s / /'/ Author / �95 OMI Number: 9539476 DMI Microform 9539476 Copyright 1995, by DMI Company. All rights reserved. This microform edition is protected against unauthorized copying under Title 17, United States Code. UMI 300 North Zeeb Road Aim Arbor, MI 48103 MICROWAVE IMAGING FOR 2 D ELECTRICAL PROPERTY DISTRIBUTION PROFILING A Thesis Submitted to the Faculty in partial fulfillment of the requirements for the degree of Doctor of Philosophy by PAUL MORGAN MEANEY, JR. Thayer School of Engineering Dartmouth College Hanover, New Hampshire MAY 1995 Examining Committee: Chairman Member: Member: ( Member:: j c Member: Dean of Graduate Stud s Author �95 ii Thayer School of Engineering Dartmouth College MICROWAVE IMAGING FOR 2-D ELECTRICAL PROPERTY DISTRIBUTION PROFILING PAUL MORGAN MEANEY, JR. Doctor of Philosophy Committee: Keith D. Paulsen (Chairman) Eric W. Hansen Robert K. Crane (University of Oklahoma) B. Stuart Trembly Alexander Hartov ABSTRACT A microwave imaging system with concomitant reconstruction algorithm has been implemented for 2-D electrical property distribution profiling in laboratory-scale phantom experiments. Achieving quantitative absolute images using this system is an important preparatory step towards achieving dynamic images for temperature field estimation during hyperthermia treatment for cancer. A Newton-Raphson method has been implemented with the hybrid element as its computational base for the image reconstruction. It also incorporates a dual mesh scheme for data set and computational problem size reduction without compromising image quality. Integration of this 2-D reconstruction algorithm with a microwave measurement system has also required mapping of the 3-D antenna radiation patterns into 2-D space. With this transformation, the measured electric field data has been shown to agree well with that of the computational model. Initial laboratory experiments show this system to be the first reported to successfully reconstruct quantitative 2-D images of both the real and reactive components of simple targets with large contrasts with respect to the background medium. Subsequent phantom experiments performed at several frequencies, investigated different sized targets with varying contrast levels and multiple configurations and have confirmed that this method can successfully reconstruct absolute images of the complex electrical property profile. A significant finding from these studies was that the receiving antenna configuration plays a vital role in the ability to reconstruct images, with data from a dual-tiered receiver antenna arrangement improving images obtained where only a single tier of antennas had been used. A second result showed an improvement in image quality with frequency. This was especially evident when the object being imaged was encased in a high contrast cylinder, in which case the electrical properties were accurately reconstructable only at higher frequencies. Lastly, statistical analysis of the positions of the reconstructed objects clearly indicated that the ability of the system to locate an object degrades in the presence of additional targets. Based on these experiences, it is concluded that a next generation imaging system, which would be able to achieve dynamic images, will have good potential as a clinical aid in imaging temperature changes during hyperthermia. iii Acknowledgements I would like to thank the main people without whose help this work would not be possible. In particular, I would like to thank Keith Paulsen for his work with me not only as an advisor but as a colleague. This work has been a constant trade-off of what can be accomplished in the numerical modeling and with the actual microwave hardware implementation. He has given me a lot of credit for my experience in the microwave industry, but then constantly stretches my imagination as to what we would really like to build. Through this interaction we have developed a strong synnergism that has been very effective. I am also grateful to Robert Crane, Alexander Hartov and Thomas Ryan. Through the various stages of this project, 1 have always been able to draw on their expertise. I would also like to thank these three along with Eric Hansen and Stuart Trembly; as members of my thesis committee, they have lent constructive suggestions which have helped strengthen this work. I would like to thank all the people in the Numerical Methods Lab for their support and for putting up with me. In particular, Jim Waugh and Chris Naimie have been valuable resources on questions regarding our wide range of software capabilities. I would like to thank Mike Moskowitz who was not only a valuable resource but also an inspiration for blazing a trail before me. I am indebted to the Machine Shop at Thayer School, Roger Howes, Robert Stromberg and Leonard Parker, for their support and suggestions. I have spent many hours down there working and getting a kind of needed therapy. I would also like to thank Benoit Cushman-Roisin for listening to my constant complaints about the financial aid issues and working with the DSE to fix it. As any researcher knows, nothing gets done before the funding issues are settled. My special thanks goes to Victoria Fullerton, without whose constant help and moral support this would not have happened. This work was supported by DHHS/NCI Grant # R01-CA55034 and by Paul and Margaret Meaney. iV TABLE OF CONTENTS Abstract ii Acknowledgments iii Table of Contents iv List of Tables vii List of Figures viii Chapter L Introduction Chapter IL Formulation of the Image Reconstruction Algorithm II. 1 Introduction 1 12 12 11.2 Reconstruction Algorithm 14 11.3 20 Results II.3.A Small Mesh II.3.A.i Discrete Regions Problem 21 II.S.A.ii Continuous Distribution Problem 27 II.3.B Large Mesh 29 II.3.B.i Discrete Regions Problem 29 II.3.B.ii Continuous Distribution Problem 30 II.3.C Simulated Measured Data 11.4 Discussion and Conclusions Chapter III. The Dual Mesh Scheme for Overall Problem Size Reduction 111.1 21 32 34 38 Introduction 38 111.2 Dual Mesh Method 40 111.3 Results 45 III.3.A Small Mesh 46 III.S.A.i Discrete Regions Problem 46 IIl.S.A.ii Continuous Distribution Problem 51 III.3.B Large Mesh 52 III.J.B.i Discrete Regions Problem 54 IlI.S.B.ii Continuous Distribution Problem 59 III.3.C Simulated Measured Data 111.4 Discussion and Conclusions 59 65 V Chapter IV. The Active Microwave Measurement System 70 IV. 1. Introduction 70 IV.2. Hardware Design 72 IV.2.A Microwave Circuit Design IV.2.A.i Modulation Scheme 72 IV.2.A.ii Transmitter/Receiver Design 76 IV.2.B Antenna Design and Arrangement IV.3 Calibration IV.4 72 77 80 IV.3.A Sources of Measurement Error 83 Results 84 IV.5 Discussion and Conclusions Chapter V. System Evaluation in Tissue-Equivalent Phantoms 92 96 V.l. Introduction 96 V.2. Receiver Antenna Configuration 98 V.3. Measurement Sites Outside the Target Region 101 V.4. 105 Results V.4.A Materials V.4.B Phantom Tests VI. 1 107 VA.B.i Filter Parameter Determination 107 VA.B.ii Use of two-tieredReceiver Sites 109 VA.B.iii Resolution with Frequency 113 VA.B.iv Single Object Image Reconstructions 127 V.4.B. V Reconstructed Image Error Analysis 127 V.5 Discussion and Conclusions Chapter VI. 105 Conclusions Overview VI.2 Future Directions 134 150 150 152 References 155 Appendix A. Forward Solution of the Electric Field Distribution using the Hybrid Element Method 162 Appendix B. Algorithmic Computational Time Enhancements 165 Appendix C. Dual Mesh Overhead 167 vi Appendix D. Noise Figure Analysis 170 Appendix E. Water-filled Waveguide Return Losses 172 Appendix F. Radiation Field Variation Outside of the 2-D Plane 174 Appendix G. Calibration of Variations between Transmitter Channels and Receiver Channels 177 Appendix H. Spatial Low Pass Filter of Reconstructed Property Distribution 180 vii LIST OF TABLES Table I.l. List of various tissue electrical properties and their variation with temperature 4 Table II.3.A.i.l. Sampling rates for discrete regions distribution on the small mesh 23 Table II.3.B.i.l. Sampling rates for discrete regions distribution on the large mesh 29 Table III.2.1. Basis function values at the nodes of the new sub triangles 44 Table III.3.A.i.l. Sampling rates for discrete regions distribution on the small mesh 46 Table III.3.B.i.1. Sampling rates for discrete regions distribution on the large mesh 54 Table IV.4.1. Comparison of the calculated and adjusted measured data for the seventeen measurement sites in the homogeneous medium 85 Table IV.4.2. Image reconstruction errors 91 Table V.4.A.1. Electrical properties of tissue equivalent materials 107 Table V.4.B.v.l. Average RMS errors of the background and object areas of the target region for single-object image reconstructions 130 Table V.4.B.V.2. Average RMS errors of the background and object areas of the target region for multi-object image reconstructions 131 Table V.4.B.V.3. Average diameters and position errors of fat/bone equivalent cylinders for single-object image reconstructions 132 Table V.4.B.V.4. Average diameters and position errors of fat/bone equivalent cylinders for multi-object image reconstructions 133 viii LIST OF FIGURES Figure II.3.A.i.l. Meshes used for small and large problems 22 Figure II.3.A.i.2. Reconstructed images for discrete regions distribution on the small problem 24 Figure II.S.A.i.S. RMS electric field error convergences vs iterations 25 Figure II.3.A.i.4. Final RMS electric field errors versus noise levels for several contrast levels 26 Figure II.3.A.i.5. RMS 27 errors for a range of |k^Ax^| values Figure II.3.A.ii.l. Reconstructed images for continuous distribution on the small problem 28 Figure II.3.B.i.l. Reconstructed images for discrete regions distribution on the large problem 30 Figure lU.B.ii.l. Reconstructed images for continuous distribution on the large problem 31 Figure II.3.C.1. Very fine mesh used to compute simulated measured data 32 Figure II.3.C.2. Reconstructed images for the small problem using the simulated measured data 33 Figure II.3.C.3(a). Comparison of reconstructed and exact values along a longitudinal transect through body of the discrete regions distributions 34 Figure II.3.C.3(b). Comparison of reconstructed and exact k^l? values along a azimuthal path around the body of the continuous distributions case 35 Figure III.2.1. Geometry of piecewise integrations in the dual meshing scheme 42 Figure III.3.A.1. Fine and coarse meshes for small mesh problem 47 Figure III.3.A.i.l. Images for discrete regions distribution on the small mesh problem 48 Figure III.3.A.i.2. RMS electric field errors per iteration 50 Figure III.3.A.i.3. Condition numbers vs. |k^Ax^| for the single and dual mesh methods on the small mesh problem 51 Figure III.3.A.ii. 1. Images for continuous distribution on the small mesh problem ... 53 Figure III.3.B.1. Fine and coarse meshes for the large mesh problem 55 ix Figure IIU.B.i.l. Images for discrete regions distribution on the large mesh problem 56 Figure III.3.B.i.2. RMS electric field errors per iteration 57 Figure III.3.B.ii.l. Images for continuous distribution on the large mesh problem .... 58 Figure III.3.C.1. Very fine mesh 60 Figure III.3.C.2. Reconstructed images for the discrete regions distribution on the small problem using electric field data computed on the very fine mesh 61 Figure III.3.C.3(a). Re( k^L^) along a longitudinal transect across the body with the discrete regions distribution for the fine-fine, fine-coarse, and coarse-coarse reconstructions 62 Figure III.3.C.3(b). Im( k^L^) along a longitudinal transect across the body with the discrete regions distribution for the fine-fine, fine-coarse, and coarse-coarse reconstructions 63 Figure III.3.C.4. Reconstructed images for the continuous distribution on the small problem using electric field data computed on the very fine mesh 64 Figure III.3.C.5(a). Re( k^L^) along a longitudinal transect across the body with the continuous distribution for the fine-fine, fine-coarse, and coarse-coarse reconstructions 65 Figure III.3.C.5(b). Im(k^L^) along a longitudinal transect across the body with the continuous distribution for the fine-fine, fine-coarse, and coarse-coarse reconstructions 66 Figure IV.2.A.i.l. Schematic diagram of the data acquisition system 74 Figure IV.2.A.i.2. Photograph of the data acquisition system 75 Figure IV.2.A.i.3. Photograph of the microwave circuit plate 75 Figure IV.2.B.1. Return Losses of the quarter wave monopole immersed in two different media 78 Figure IV.2.B.2. Photograph of the antenna plate with the transmit and receive antennas suspended from the plate 79 Figure IV.3.1. Diagram used for calculation of the phase center 81 Figure IV.4.1. Meshes used in reconstructions using measured data 87 Figure IV.4.2. Real and imaginary parts of the reconstructed images of the k^ distribution for the various phantoms 89 Figure IV.4.3. Electrical property profile along longitudinal transects through the target regionsfor the real and imaginary parts, respectively 90 X Figure V.2.1. Comparison of calculated field values on the perimeter of a homogeneous and inhomogeneous target region 99 Figure V.2.2. Diagram of electromagnetic waves as they are perturbed after passing t&ough and around arbitrary objects in the target region 100 Figure V.2.3. Diagram of electromagnetic waves as they are perturbed after passing through arbitrary objects in the target region 101 Figure V.3.1. Representation of the hybrid element based reconstruction algorithm with the incorporation of measurement sites outside of the finite element region 103 Figure V.4.1. Diagram of meshes used in the reconstruction algorithm 106 Figure V.4.B.i. 1. Reconstructed images at 500 MHz for the 4.3 and 2.5 cm fat/bone equivalent cylinders separated by 0.4 cm using various filter weightings 108 Figure V.4.B.ii. 1. Reconstructed images at 500 MHz for the 4.3 cm and 2.5 cm fat/bone equivalent cylinders at different separations using data from just the inner arc and from both arcs of measurement sites 110 & 111 Figure V.4.B.ii.2. Electrical property profile along longitudinal transects through the target region for images shown in Figure V.4.B.ii. 1 112 Figure V.4.B.iii.l. Reconstructed images for the 4.3 cm and 2.5 cm fat/bone equivalent cylinders at 300, 500 and 700 MHz for two separation distances 114 & 115 Figure V.4.B.iii.2. Electrical property profile along longitudinal transects through the target region for images shown in Figure V.4.B.iii.1 116 & 117 Figure V.4.B.iii.3. Reconstructed images for the 4.3 cm diameter fat/bone equivalent cylinder and the 3.8 cm diameter DI water cylinder at 300 and 700 MHz with two separation distances 119 & 120 Figure V.4.B.iii.4. Electrical property profile along longitudinal transects through the target region for images shown in Figure V.4.B.iii.3 121 & 122 Figure V.4.B.iii.5. Reconstructed images for the 4.3 cm diameter fat/bone equivalent cylinder and the 3.8 cm diameter agar gel cylinder at 300 and 700 MHz with two separation distances 123 & 124 Figure V.4.B.iii.6. Electrical property profile along longitudinal transects through the target region for images shown in Figure V.4.B.iii.5 125 & 126 Figure V.4.B.iv.l. Reconstructed images for the 1.6 cm fat/bone equivalent cylinder centered in the target region at 300, 500 and 700 MHz 136 Figure V.4.B.iv.2. Electrical property profile along longitudinal transects through the target region for images shown in Figure V.4.B.iv.1 137 Figure V.4.B.iv.3. Reconstructed images for the 2.5 cm fat/bone equivalent cylinder centered in the target region at 300, 500 and 700 MHz 138 xi Figure V.4.B.iv.4. Electrical property profile along longitudinal transects through the target region for images shown in Figure V.4.B.iv.3 139 Figure V.4.B.iv.5. Reconstructed images for the 4.3 cm fat/bone equivalent cylinder centered in the target region at 300,500 and 700 MHz 140 Figure V.4.B.iv.6. Electrical property profile along longitudinal transects through the target region for images shown in Figure V.4.B.iv.5 141 Figure V.4.B.iv.7. Reconstructed images for the 2.5 cm fat/bone equivalent cylinder offset by 2.5 cm in the target region at 300, 500 and 700 MHz 142 Figure V.4.B.iv.8. Electrical property profile along longitudinal transects through the target region for images shown in Figure V.4.B.iv.7 143 Figure V.4.B.iv.9. Reconstructed images for the 4.3 cm fat/bone equivalent cylinder offset by 2.5 cm in the target region at 300, 500 and 700 MHz 144 Figure V.4.B.iv.lO. Electrical property profile along longitudinal transects through the target region for images shown in Figure V.4.B.iv.9 145 Figure V.4.B.iv.l 1. Reconstructed images for the 2.5 cm fat/bone equivalent cylinder offset by 3.8 cm in the target region at 300, 500 and 700 MHz 146 Figure V.4.B.iv.l2. Electrical property profile along longitudinal transects through the target region for images shown in Figure V.4.B.iv.11 147 Figure V.4.B.iv.l3. Reconstructed images for the 4.3 cm fat/bone equivalent cylinder offset by 3.8 cm in the target region at 300, 500 and 700 MHz 148 Figure V.4.B.iv.l4. Electrical property profile along longitudinal transects through the target region for images shown in Figure V.4.B.iv.13 149 CHAPTER I. INTRODUCTION Hyperthermia is a cancer treatment which seeks to raise the temperature of the whole tumor to a therapeutic threshold with only a minimal increase in the temperature of the surrounding tissue. Hyperthermia treatments are generally performed using microwave or ultrasound radiation and have shown the most promise when performed as an adjuvant therapy to ionizing radiation (RT) or chemotherapy. As such, there is strong evidence from studies with spontaneous animal tumors [Dewhirst, et al., 1982] that properly applied hyperthermia can enhance the effects of RT by increasing local tumor control [Overgaard, et al., 1987; Valdagni, et al, 1988; Hynynen, et al., 1990]. Temperature field distributions in the body induced by hyperthermia can vary considerably over the tumor region primarily due to variations in blood flow and uneven heating. To date, unfortunately, temperatures can only be monitored at a finite number of points, thus making it difficult to determine the complete distribution in the presence of such large variations. This can result in the use of hyperthermia systems which are unable to control the deposition of heat within the tumor region. Subsequently, this has made it difficult to assess the efficacy of hyperthermia as an adjuvant treatment and has led to confusing and conflicting results regarding its benefits [Perez, et al., 1991]. It is clear that conclusions drawn about hyperthermia's success as an adjuvant cancer therapy will depend on the ability to improve the temperature field monitoring of the tumor region during treatment. Monitoring of temperature field distributions generally fall into two categories: (I) computing the temperature field distributions directly from a finite number of temperature measurements, and (2) computing the temperature field distributions indirectly by monitoring changes in temperature dependent tissue properties. The first method involves a heat transfer model which computes the temperature fields everywhere within a region given a finite number of temperature measurements and assumptions about the heat transfer characteristics of the 1 2 affected tissue [Clegg and Roemer, 1989; Oleson, et al., 1989; Edelstein-Keshet, et al., 1989; Strohbehn, et al., 1986]. An initial study by Edelstein-Keshet, et al. [1989] developed a simple mathematical representation which included the geometry, boundary conditions of the bioheat transfer equation and the specific absorption rate (SAR) patterns. They evaluated a set of numerical quantities, called descriptors, which included such things as: average temperature, fraction of tumor heated above a given temperature, mean edge temperature, and mean edge temperature gradient. These descriptors were calculated based on a simple mathematical representation that had been fitted to the data and studied using statistical tests to determine the model's effectiveness. Another study by Clegg and Roemer [1989] sought to combine measurement data at finite sampling points with a numerical model, where the blood perfusion is treated as an unknown. The perfusion is then estimated in an iterative process by minimizing the difference between predicted and measured temperature values at a finite number of sample positions, from which the steady-state temperature field distribution can be deduced. While it has been shown that this method can achieve 1 "C accuracy under ideal conditions, the performance degrades significantly in the presence of noise. As is discussed in Clegg and Roemer [1989], this method is largely dependent on estimating the blood perfusion within distinct tissue regions, which to date, is largely unknown and is the most critical parameter in temperature field estimation during hyperthermia treatment [Paulsen, et al., 1985; Roemer, 1991]. Various attempts have been made to accurately model the effects of blood perfusion on the temperature field distributions. These include initial simplified models by Pennes, et al. [1948] and Perl, et al. [1966] which assumed that most of the energy exchange occurred at the capillary level and the convective effects could be treated as a volumetric heat source. Subsequent models by Lagendijk, et al., [1984], Weinbaum, et al. [1985], and Roemer [1990] have shown that the more significant temperature variations can be due to the presence of larger blood vessels and countercurrent vessel pairs within the tumor region. Studies by Moros, et al. [1993] compared 3 the simple bioheat equation model of Pennes [1948] and that by Lagendijk [1987] to in vivo experimental temperature measurements under hyperthermic conditions and showed that the former matched the measured data to better than 1癈 at all measurement points while the latter matched only 79% of the data to better than 1癈. The second method of inferring temperature Helds is to deduce them indirectly from changes in temperature dependent material properties. Two such properties which can be monitored through the use of magnetic resonance (MR) imaging are (1) the molecular diffusion, Tl, and (2) the proton chemical shift [LeBihan, et al., 1989]. Because of the size of MR systems, integration of this type of imaging with hyperthermia treatment is quite difficult, but has been achieved with some success [Samulski, et al., 1992]. Experiments by MacFall, et al. [1995] using MRI during hyperthermic heating have achieved 0.55獵 accuracy with a 1.8 cm resolution for phantom and 0.9''C accuracy with a 0.9 cm resolution during in vivo experiments. Electrical impedance tomography (EIT) has also been used as a means of mapping temperature distributions based changes in the tissue electrical conductivity with temperature [Conway, 1987; Blad, et al., 1992; Moskowitz, et al., 1994]. By first reconstructing an image of the target region before the heating has begun and subtracting it from subsequent images obtained during heating, temperature fields can be inferred. In the work by Moskowitz, et al. [1995,1994] temperature accuracies of 1癈 in phantom and 1.4''C in vivo were achieved. EIT has found several applications beyond just temperature field monitoring, such as imaging of basic anatomy [Henderson and Webster, 1978] and in vivo biological activity [Barber, et al., 1992; McArdle, et al., 1992]. Reconstruction algorithms that have been employed have ranged from simple and filtered back-projection [Barber and Seager, 1987; Barber and Brown, 1990; Brown, et al., 1992] to iterative techniques employing variations of the Newton-Raphson method [Yorkey, 1986; Isaacson, 1986; Paulsen, et al., 1994]. Differential imaging has showed promise; however, after heating and cooling the tumor region, the baseline images are often not recovered [Moskowitz, 1994; McRae and Esrick, 4 1993]. This suggests that permanent damage may result in the treatment region, causing changes in the material electrical properties which increases the need for good absolute imaging capabilities. Active microwave imaging is a modality which appears to be ideally suited to imaging the complex electrical properties, electrical permittivity and conductivity, below the skin surface. It has been shown that these electrical properties vary with temperature in biological tissue over the microwave frequency range [Schwan, et al., 1953; Bottomley, et al., 1978; Burdette, et al., 1986]. Table I.l lists the electrical properties of various tissue types and their variation with temperature at 433 and 915 MHz. 433 MHz Material Er A/癱 O 915 MHz Act/ , /A/癈 Er AEF/ /Er/癈 (s/m) a ACT/ , /C/癈 (s/m) Fat 5 2.6% 0.26 1.8% 5 4.2% 0.35 1.1% Kidney 60 -0.2% 1.22 2.0% 55 -0.4% 1.41 1.3% Liver 47 -0.2% 0.89 1.8% 46 -0.4% 1.06 1.4% Muscle 57 -0.2% 1.12 1.3% 55 -0.2% 1.45 1.0% Table I.l. List of various tissue electrical properties and their variation with temperature. Thus, being able to monitor changes in the electrical properties during hyperthermia treatment could allow temperature field distributions to be inferred using differential imaging techniques such as those employed in EIT as long as good absolute images were available on which to base relative changes. Promising results have been achieved by Mallorqui, et al. [1992] and Jofre, et al. [1990] using an active microwave imaging system operating at 2.45 GHz to image dynamic changes in biological activity such as heating via hyperthermia and blood flow 5 changes employing these differential techniques. In monitoring temperature, accuracies of 0.5=0 have been achieved [Jofre, et al., 1990]; however, absolute images of phantoms have not been possible to date, primarily because of limitations in the reconstruction algorithm involved which utilizes diffraction tomography based on the Bom approximation. This approximation is only useful when the size of the targets are electrically small and their contrast with the background medium is within a couple of percent [Slaney, et al., 1986] which is a significant limitation in the context of biomedical imaging where contrasts are often greater than one order of magnitude. Improvements in the ability to achieve good absolute images will only enhance the quality of the temperature field estimates utilizing differential techniques with these images. A second technique to deduce temperature field distributions has been reported by Miyakawa [1993]. The approach reported employs chirp radar techniques to isolate the straight line transmission path while eliminating multipath signals such that conventional tomographic reconstruction algorithms used in x-ray CT imaging can be utilized. It has achieved 5癈 accuracy with 1 cm resolution, to date. More generalized microwave imaging of biological tissue has been the subject of several investigations over the past two decades [e.g. Mallorqui, et al., 1992; Broquetas, et al., 1991; Jofre, et al., 1990; Bolomey, 1989; Bolomey, et al., 1983, 1982; Larsen and Jacobi, 1986; Guo et al., 1989] for a broad range of applications including temperature field estimation. For example, investigators have sought to exploit the broad range of soft tissue electrical properties (as much as 20:1 - See Table I.l), whereas conventional techniques such as X-ray CT can only image densities on the order of a few percent [Larsen and Jacobi, 1986]. Reconstructing images of these properties could potentially provide more detailed information about tissue functionality than is presently available. Problems with this method generally arise because these standard radar imaging techniques require that the target be in the Fraunhofer region with respect to both the transmit and receiver antennas; this is not possible in this case due to the high signal attenuation of the media under study [Guo, et al., 1989]. 6 Initial work has been performed by Bolomey, et al. [Bolomey, 1989; Bolomey, et al., 1983, 1982] and is based on diffraction tomography. This system incorporates a coherent signal detection scheme which is required for recovering both the magnitude and phase of the electric fields. The sixty-four horn antennas arranged on the perimeter of a water tank have been used and made electronically selectable to operate in either the transmit or receive mode. For a given antenna transmitter, only the 32 antennas directly opposite it are selected for the receive mode, which yields a total of 2048 pieces of measurement data [Jofre, et al., 1990]. At an operating frequency of 3 GHz, the system has been able to achieve qualitative images with resolution on the order of half a wavelength [Broquetas, et al., 1991]. Larsen and Jacobi [1986] have reported the use of standard radar techniques to image a canine kidney. Generally, a signal transmitted from one antenna and received by another undergoes amplitude, phase and polarization changes due to the presence of dielectric objects in the propagation path. Imaging performed using the forward scattering mode is achieved in a similar manner to that used in backscattered techniques [Sinclair, 1951; Copeland, 1960; Huynen, 1970] which are quite common in conventional radar applications where the only data available is from reflected signals. For medical applications, forward scattered systems offer some advantages over backscattered systems in part because of the typically smooth radar cross-sections as a function of wavelength in comparison with the rapid oscillatory behavior of the associated backscattered radar cross-sections with wavelength. Qualitative images performed at 3.9 GHz have been achieved using just the amplitude or phase data in either the co-polarized or cross-polarized configuration which show features that can be related to known regional specializations within the kidney. Work presented by Quo, et al. [1987] has introduced some interesting and potentially useful concepts. Their image reconstruction method discretely samples the Green's function 7 which converts these discretized values to the Green's matrix for use in the calculation of the electric fields based on known material properties. Their work has suggested that in the presence of a lossy medium, this matrix is invertable and that there is a unique solution. This method attempts to reconstruct complete 3-D images which is a daunting task due to the size of the problem; in particular, since the primary assumptions of digitizing the Green's function have been that the electric fields do not vary significantly over each cell, the required cell sizes must be on the order of one tenth of a wavelength (Their initial work sampled the region every quarter wavelength). At 3 GHz, this would require cell sizes on the order of 3 mm which would lead to a densely-packed Green's matrix of size NxN, where N is the total number of cells. Solving this matrix in itself would be a computationally intensive problem; it would also require N pieces of measurement data which would be difficult to realize. The work has also suggested that resolution is not necessarily tied to the frequency of operation, but is linearly related to the amount of measurement data. Even with reduction in operating frequency, the approach still presents a large matrix to be solved and large amounts of data to be collected. As of yet, practical implementation of this concept has not been reported. Recent work in the area of inverse scattering has included two new image reconstruction algorithms by Joachimowicz , et al., [1991] and Caorsi, et al., [1990]. The former is an iterative scheme for calculating the complex electrical properties of the materials in 2-D and 3-D by updating the material properties and recomputing the electric fields based on the new values until the difference between the computed and measured electric fields is small. It has some unique features such as using a spatial Fourier transform technique to reduce the electric field computation time for each iteration. It has shown some success with simulations, but has several limitations including the need for a plane wave illumination and the fact that the forward solution algorithm, based on the Moment Method, has been shown to have distinct inadequacies in certain electric field calculations [Borup, et al., 1987]. 8 The second technique [Caorsi, et al., 1991, 1990] uses the same computational base as the previous algorithm, namely the Moment Method, but reconstructs the images of the electrical properties in a single step using the pseudo-inverse transformation [Golub, et al., 1989]. It seeks to overcome the inherent nonlinear nature of this inverse problem by increasing the amount of measurement data available while increasing the effectiveness of that data through unique antenna configurations. This is quite similar to the algorithm reported by Guo, et al. [1987] mentioned above, in that they both perform a discretization of the Green's function and both attempt to perform the full 3-D image reconstruction. To date, neither of these groups has reported on the implementation of these algorithms as part of a complete imaging system to test their effectiveness in a laboratory setting. Measurement systems for microwave imaging experiments have one major goal in common: they all must be coherent. Two different techniques are generally employed. The first is a heterodyning system that compares the transmitted signal with the original reference signal to detect the relative amplitude and phase changes which occurred during transmission [Jofre, et al., 1990]. The second employs a chirp signal which transmits a sweep of frequencies and compares it (after transmission) with the reference frequency sweep to extract the amplitude and time of transit which can then be converted into phase [Miyakawa, 1993]. This has the advantage of being able to distinguish between the primary signal and any multipath contribution. However, its applicability is limited in many biological investigations because it requires a very large bandwidth. The higher frequencies will often not be able to transit the entire target region due to attenuation in the medium. The goal of this thesis work has been to develop a 2-D image reconstruction algorithm, implement it along with a coherent microwave measurement system in a laboratory scale test set-up and assess the overall imaging capabilities of the combination. The initial implementation has been achieved using only 4 transmit and 4 receiver antennas which require 9 manual movement from measurement site to site. Because of this, the data acquisition time is quite long, rendering dynamic imaging of phenomena such as heating, impossible at this stage. However, this implementation has allowed us to experiment with various antennas, antenna configurations, and operating frequencies to determine design parameters for a future generation system. As mentioned earlier, simultaneous quantitative reconstruction of the real and imaginary parts of phantom electrical properties using microwave imaging has not been achieved to date. Thus, strides in improving absolute imaging capabilities will greatly enhance the ability to achieve dynamic images (using differential imaging techniques) whose quality is dependent on a quantitative image of the tissue electrical properties as a baseline for extracting the differential image. Extracting differential images based on inaccurate absolute images may yield sub-optimal dynamic images. While assessing this prototype's capabilities, we have achieved quantitative 2-D images of a variety of phantoms with a wide range of material contrasts. The work contained herein studies 2-D image reconstructions using data from microwave imaging simulations and laboratory studies. We develop a strategy for transforming the 3-D measured data into a 2-D plane such that 2-D images can be achieved. The work includes descriptions of the image reconstruction algorithms, design of the data acquisition system and comparison of results. Portions of this work have been published elsewhere [Meaney, et al., 1995; Paulsen, et al., (a); Meaney, et al., (a); Meaney, et al., (b)] but will be repeated here for completeness and continuity. The organization of this work proceeds as follows: Chapter II describes the image reconstruction algorithm based on the hybrid element method [Lynch, et al., 1986; Paulsen, et al., 1988] which has been used successfully in modeling electric field distributions in biological applications. It begins with the derivation of the reconstruction algorithm using the hybrid element method and discusses which of the hybrid element method implementations is 1 0 best suited for this appUcation. Numerous simulations are the performed to illustrate the effects of noise, problem size, various problem classifications, mesh sampling rate and material contrast ratios. For these comparisons, four criteria were used in evaluating the effectiveness of this algorithm, including: (1) visual inspection of the reconstructed images, (2) transects of the images to compare the reconstructed with the actual material values, (3) RMS electric field errors after each iteration, and (4) the Hessian matrix condition number. In general, the results of these simulations were promising for future implementation with actual measured data. Chapter III describes the implementation of the dual mesh scheme which was introduced as a method of decoupling the mesh used for the electric field computation from that used for the material property estimation as a way of reducing both the overall problem size and the amount of measured data required. Its implementation is described including a discussion of the computational overhead required for coordination both meshes. Tests similar to those in Chapter II are then performed to demonstrate the algorithm's strengths under various situations. Using the same criteria as used in Chapter II, images of comparable quality were achieved using the dual mesh scheme with respect to those obtained using the single mesh method while significantly reducing the problem size and amount of measured data required. It is also shown to have potential advantages in situations where a priori information is available since deployment of the reconstruction parameter nodes could then be made based on rational decision making. Nonuniformly distributing these nodes is easily accomodated by the technique. Chapter IV describes the design of the coherent heterodyne data acquisition system and all its related considerations including trade-offs required resulting from the nature of the reconstruction algorithm outlined in Chapters II and III, and various calibration procedures. Design discussions include component selection for broad bandwidth and large dynamic range operation, and the antenna design utilizing a highly conductive surrounding medium. A 11 calibration procedure is discussed which transforms the 3-D antenna radiating characteristics to 2-D. Results first compare the measured electric field data with that computed using the numerical algorithms discussed in Chapter II. Preliminary reconstructed images are compared for cylindrical objects of varying diameters having high contrast levels with the surrounding medium to demonstrate the strength of this integrated system - both the numerical algorithm and data acquisition system. Chapter V investigates various antenna configurations with regard to improving the reconstructed image quality without unduly increasing the amount of measured data. The implications of using such configurations and their effects on the reconstruction algorithm are also discussed. Various test cases were chosen to demonstrate the capabilities of this system including: (1) several single object (varying diameters) experiments having high contrast levels with the background at multiple locations; and (2) several double object experiments of varying contrast levels and multiple separation distances. These tests were performed over a range of frequencies to further examine the system's capabilities. The results were generally successful and showed improvement with increased frequency in the real part of the images while the imaginary part showed some degradation with increased frequency due in part to the inadequate spatial sampling rates in the fine mesh with respect to the signal exponential decay in the medium. Statistical analysis of the reconstructed object sizes and their positions was performed for both single- and multi-target cases in a fiirther effort to quantify this system. Chapter VI provides an overview of the results achieved in the preceding chapters and includes a discussion of ideas for the next generation system implementation, especially with respect to the goals of: (1) achieving dynamic images, (2) implementing a system which can be operated in a clinical setting, and (3) utilizing parallel computation techniques for realizing pseudo-real time images. CHAPTER II. FORMULATION OF IMAGE RECONSTRUCTION ALGORITHM H I Introduction As indicated in Chapter /, there have been two recent inverse scattering methods introduced which attempt to perform electrical property reconstructions for biological bodies [Joachimowicz, et al., 1991; Caorsi, et al., 1990]. These two models base their approach on a Moment Method formulation with pulse basis functions on a domain which is comprised of uniform square sub-areas [Ney, 1985; Harrington, 1982; Joachimowicz, et al., 1990]. The first method by Joachimowicz, et al. [1991] calculates the scattered fields with an initial guess for So (r) which represents the contrast between the wavenumber of the external medium and that of each subsection of the body. The differences in the calculated and measured scattered fields are then minimized in a Newton iterative scheme to determine the final values of s(r), from which the complex dielectric properties can be determined. Because the problem is ill-conditioned in determining the As(r)'s, a regularization process is employed to help stabilize the matrix decompositions that are involved [Tikonov, et al., 1977]. The calculation of the scattered fields is performed using the Moment Method with pulse basis functions over a domain of uniform, square sub-areas, where the discretized integral equation is expressed in a convolutional form. The whole equation is Fourier transformed into a domain represented in terms of spatial wavenumbers. There the convolution integral of the spatial functions is replaced by multiplication of the transformed functions which takes fewer calculations and uses much less storage [Chommeloux, et al., 1986]. For simple cases, this method has been shown to converge to an acceptable error level. The second method by Caorsi, et al. [1991,1990,1989] forms a scattering matrix also based on the Moment Method above. This method uses multiple excitations to yield more 12 13 measurements than degrees of freedom being estimated and performs a noniterative pseudoinverse transformation to yield a least squares solution for the electrical property distribution. It has been shown to yield good approximations for the dielectric properties of relatively small problems. While the approach used here is an iterative technique which updates the values of the electrical properties similarly to [Joachimowicz, et al., 1991], it deviates substantially from both [Joachimowicz, et al., 1991] and [Caorsi, et al., 1990] in its computational base. The computational methodology used here is a hybrid coupling of the Finite Element and Boundary Element methods. In calculating the electric fields from the current estimates of the electrical properties, it seeks to take advantage of the strong points of the Finite Element method [Lynch, et al. 1985] in the regions where the electrical properties are inhomogeneous and not known. It also utilizes the strengths of the Boundary Element method [Yashiro, et al., 1985] in the regions where the medium is homogeneous, unbounded in nature and the electrical properties are known. Coupling of these two methods occurs only at the boundary of the biological body. This coupling creates a hybrid of the two methods, the Hybrid Element method (HE) [Lynch, ci al., 1986], which computes the electric fields inside an inhomogeneous dielectric region that is being radiated by microwave sources positioned outside the body which are embedded in an infinite homogeneous region. Because the Newton's iteration scheme is computationally demanding on a per iteration basis, the computational overhead of the forward solver is important which provides further rationale for the selection of the hybrid approach in this context. The computational costs in calculating the forward solution are modestly increased in comparison to using the FE alone because only the target region must be discretized which does increase the bandwidth but only to the extent that the small but dense BE system of equations is involved. Paulsen et al. [Paulsen, et al., 1992a] show that in general there are three convenient methods for the coupling of the FE and BE problems, but that when considering only the forward 14 solution of the scattering problem, the so-called Hybrid Element method #2 is advantageous in terms of storage space and the number of computations. Quite interestingly, in the inverse or reconstruction problem, a different Hybrid Element method (referred to as HE #1) has the advantage especially in computation time and overall complexity, as will be illustrated. Further, the potential for non-uniform meshing without compromising computational efficiency is powerful for problems where nonuniform meshing strategies are appropriate, for example, when some form of a priori information is available. This feature is not generally available to the schemes presented in [Joachimowicz, et al., 1991; Caorsi, et al., 1990]. The purpose of this chapter is to introduce a Newton iteration scheme for image reconstruction which is based on the Hybrid Element method and to characterize its performance on several simple imaging problems. Particular attention will be paid to the derivation of the inverse problem solution where the implementation using the Hybrid Element method is involved. Test cases have been chosen using simulated measured data to examine this technique under various conditions such as small and large problems and also conditions where noise is present. They have also been performed for two different electrical property distributions to demonstrate the flexibility of the method. II.2 Reconstruction Algorithm The reconstruction routine is an iterative technique based on Newton's method. It begins with an initial guess for vector k^, where kj is the complex-valued wavenumber in which both the dielectric constant and conductivity are imbedded for a set of sample points in the FE mesh. The reconstruction routine can be summarized as involving the following steps: 1) Computing the electric fields based on the current for each different excitation being employed using Hybrid Element method # 1 [Lynch, et al., 1986; Paulsen, et al., 1992a] (henceforth called the forward solution.) 15 2) Forming a difference electric field vector by subtracting the calculated electric fields from the measured electric fields at a finite number of locations for a finite number of distinct incident fields. 3) Constructing the Jacobian matrix, J, required in Newton's method similar to that used in Paulsen et al. [1994]. 4) Computing a perturbation vector, Ak^ using the method of normal equations [Golub, et al., 1989] where the difference field vector serves as the known right hand side. 5) Updating the values of based on Ak^ and repeating steps 1 thru 4 until a suitable stopping criteria is reached. Step 1 is described in detail in the references cited above and a brief description is included in Appendix A for completeness. Steps 2 thru 4 require more elaboration for the hybrid element implementation and will be described presently. We first approximate the true electric field, E, at a point as a function of the true electric property distribution through a two term Taylor series expanded around an approximate electrical property profile represented by an N dimensional vector, k^. 'approx ({^approx}) where is a 1 X )k J approx (ak^j j (11.2.1) N matrix of first derivatives of the approximate solution with respect to each member kf. k^, and k^pp^, are the true and approximate k^ vectors of length N and Ak^ is a difference vector between k^, and k^pp?,.' One can easily determine the electrical properties for each complex-valued member k? given the operating frequency since kf =co^nej + jcoiiCj where |X is the magnetic permeability, e is the electrical permittivity and a is the electrical conductivity. To make use of the approximation in (II.2.1) in order to update the k^ profile estimate, true values of the electric field are needed. In practice these are taken as field observations that can be recorded in the region around the ' Unless otherwise stated bold lettering is used to distinguish a column vector from a matrix in the text, whereas the symbol |} refers to a column vector and [] indicates a matrix when these quantities are used in the equations. 16 body after it has been illuminated by a known source. To generate a system of equations which can be solved for Ak^, at least N observations are needed. Due to series truncation a linear system for Ak^ can be constructed by considering N observation points each of which is approximated by the relationship given in (II.2.1). However, the true electric field varies non-linearly in k^; hence, using more observations than reconstruction parameters tends to improve the estimate for the update vector Ak^ The values of k^ are also restricted to a certain bound which limits the number of possible solutions that occur naturally due to the inherent non-uniqueness of such inverse problems [Habashy, et al., 1986]. To generate sufficient numbers of observations, multiple measurement locations can be used along with multiple excitation arrangements. In matrix form, the process can be written as [j j{ak2} = {�"}-{�=} (11.2.2) where J is the Jacobian matrix of dimensions O x N where O is the number of total observations which can be considered as the product of the number of excitation patterns times the number of observations per excitation. E? and are the measured and calculated electric field vectors of length O, respectively. This set of equations can be solved as a least squares problem using the method of normal equations [Golub, et al., 1989] by multiplying both sides of the equation by to yield: = (11.2.3) where J^J is the N x N Hessian matrix. To find the Jacobian matrix, J, we need to differentiate the equations from the forward solution for each excitation with respect to each kf parameter. As shown in Paulsen et al. [1992a] and Appendix A, the finite element region produces a matrix system which can be decomposed as: 17 ?A? > 0 I lEbl 0 0 " � .'^bl J E . L "0 QQ A,bl fb (11.2.4) and the boundary element region yields a similarly written matrix form: ? Dbsl fe. (II.2.5) 1 T"!U\ d Csb Cb.1 p "Cbb where the subscripts I, b, and s refer to all points interior to the FE surface, all points on the boundary shared by the FE and BE problems and all points on the radiating source(s), respectively. In equation (II.2.4), the sparse matrix A is partitioned into four submatrices of various dimensions, A?, A,y, A,,,, and A^y where, for example, submatrix A,,, has dimensions n, x n^ where n, , n^, and n, are the number of nodes in regions I, b, and s, respectively. Similar partitioning is used in column vectors E and F where F=^ This on notation also appears in the expression of matrix system (II.2.5) as well as later in the text. For convenience in some instances, we express the partitioned matrices in (II.2.4) and (II.2.5) in a nonpartitioned form in which case the subscripts are dropped in favor of a single alphabetical symbol. We first differentiate equations (II.2.4) and (II.2.5) with respect to the j^h member k^, of which results in #{e}+[a]^ = [b]^ dkf akf and, (11.2.6) (II.2.7) since matrices B, C and D do not vary with respect to the k?'s, because the material properties of the external medium are assumed to be known. As a result the matrices in (11.2.7) are identical to those in (II.2.5), hence, analogously to that in Paulsen et al., [1992a] and Appendix A, one can write 18 :: G bb G bs % 'sb dk] dE. dkf (II.2.8) where [G] = [c-'J [D]. In this formulation we have assumed that the boundary conditions on the source will 3E not vary with the kj's in the body and hence, = 0. Therefore, ^fb akf aeb akj gbb j (II.2.9) Rewriting equation (II.2.6) in partitioned matrix form and substituting equation (II.2.9) yields, after some simplification: d A lb/ aa,i 'bl dkj n aey abb ~ bhhg bb^ bb,j akj ^dkf (II.2.10) d A bb/ 'ak: "bj which allows computation of the elements of J. There are several important points to note. First, the matrix on the left hand side in (II.2.10) is the identical system involved in the forward solution (See Paulsen et al., [1992a] and Appendix A). Second, at each iteration the forward solution is calculated using the current values of the electrical properties to get E, and Ey, so the right side of (II.2.10) becomes a known column vector provided conditions are such that n x H or 0n can oe computed. Third, if the boundary is specified on the source, we would switch F, and af. Eg in equation (II.2.5) and, after inverting the matrix, set ?= 0 which yields a set of akj relationships that have identical form to (II.2.10) but with an altered submatrix, G^y. Fourth, as noted in Paulsen et al., [1992a], HE #2, where the FE system is inverted and inserted into 19 the BE relations, has some computational advantage in the forward problem. Here, however, if HE #2 were to be used, equation (II.2.6) would be rearranged so that A could be inverted to solve for y in terms of the right hand side which would then be substituted into equation (II.2.7) to yield a relationship for (II.2.6), the inversion of A to solve for Because of the presence of E in equation Y would require the computation of both AI^, and AI(,b, where AI=A ', which involves one LU decomposition and n,+n(, back substitutions. Comparing this with the number of back substitutions needed in the HE #2 forward solution, where only Al^b is used, results in an increase of n, back substitutions that would become necessary at each iteration. HE #1, on the other hand, needs no additional back substitutions to produce the left-hand matrix in (II.2.10) whether it is used in the forward calculation or in the construction of J. Thus, HE #1 offers significant savings both in computation time and memory when the HE is used in this reconstruction algorithm. (Additional advantages of HE #1 with respect to computation time enhancements are discussed in Appendix B.) The only remaining item needed to compute akj is akj , which can be evaluated on a matrix entry-by-matrix entry basis. It is well known that the Galerkin weighted residual form of the Helmholtz equation yields matrix elements (ii.2.11) where ((?)) denotes domain integration over the problem space and (|) are locally varying Lagrangian basis functions centered at each node in the FE mesh. Note that here we have chosen to expand in the same basis as the field solution in which case N, the number of reconstruction parameters, is identical to the number of nodes in the FE mesh. This is 20 convenient, but not necessary and perhaps not even desirable in all cases as we show elsewhere [Meaney, et al.]. Taking the derivative with respect to k? is quite straight forward and yields: (II.2.12) The Hessian matrix typically is ill-conditioned and must be modified to give stable decompositions [Hua, et al., 1990]. To handle this problem, we employ an algorithm developed by Marquardt [Dennis, et al., 1983; Marquardt, 1963]. In this case, we scale the Hessian matrix and the vector on the right hand side of (II.2.3). Then a factor, X., is added to the diagonal of the Hessian matrix and (II.2.3) becomes: (II.2.13) Equation (II.2.13) is solved at each iteration to get the new values of Ak^. Note that as X is made arbitrarily small, equation (II.2.13) becomes very similar to (II.2.3). The whole process essentially amounts to weighting the main diagonal so that the solution of the matrix is more stable. A similar process is used by Joachimowicz et al. [1991] that also adds a regularization factor to the main diagonal based on a Tikhonov method [Tikonov, et al., 1977], II.3 Results The experimental evaluation of the algorithm described in the previous section is conducted on two types of imaging problems carried out on two levels of mesh discretization. The two types of imaging problems are representative of those of interest in biomedical applications and involve (1) a domain with discrete regions of homogeneous electrical properties and (2) a domain with a continuously varying electrical property profile. The two mesh discretizations involve (1) a relatively small number of nodal degrees of freedom and 21 (2) a larger number of mesh nodes. In Sections II.3.A and B below, the exact electric field data, with which the computed electric fields at each iteration are compared, are calculated on the same mesh as that used in the reconstruction. In Section 11.3.C the method is evaluated similarly to that in Sections II. 3.A and B except that the exact electrical field data has been computed on a mesh with eight times as many nodes as that used in the small problem. In all tests, the radiators which provide the incident fields were excited one at a time with data taken at all of the measurement sites for each incident excitation. In both the smaller and larger problems, the excitation on each radiator was a simple cosine distribution, zero at the ends and unity at the center, similar to that expected from the opening of a waveguide radiator. The images are compared using three different criteria: (1) qualitatively by viewing the images, (2) comparing the relative RMS errors for the electric fields. The relative RMS errors for and (3) the relative RMS errors for and E are defined as: rms(e)= ?i (II.3.1) The condition numbers for the Hessian matrices are also compared for some cases where the condition number has been calculated for the matrix J^J, given in equation (II.2.3) based on the ratio of the largest to smallest eigenvalues of that matrix. All calculations were performed on an IBM RISC 6000 workstation. Finally the stopping criterion was an electricfield RMS error less than 10'5 or a non-declining steady-state level in this quantity. n.3.A Small Mesh II.3.A.i Discrete Region Problem: The first problem to be considered uses a mesh with 117 nodes and 197 elements in the finite element part of the domain. Also included are five applicators containing seven nodes each (Figure II.3.A.i.1(a)). Experience suggests that the product of the number of 22 Regioni Region 2 Radiators Region 1. ^ ^Region 2 ^ 1^ 1^ 1^ li*1^ l!91^ ^ 1^ 1^^ lil^li!l|i?^l^l?l^li<l?>i^i^c^i^^'4'1 i? 1^ ^ i? ^ij^ 1^ 1^ lij ~ ifdsy 141^ 5< ^iJ^1^ BiIJ5iliJliJ^'i I5i 2 ^ 1^ (i1^ 1^ li ijrii uj ^ 1^ f lj< i? i?i fc? i? ^ fc i 5is ^ ^ij^ ijjij ^ 1^1^ ^ k'i fc? ?> ?> t ,1^^ 1^ 1^ ^1^ 1^ 1^ ^1^ 1^ 1^ if 1^^ 1^ i?1^^ � v<? i? ^1^ 1^ i? ^1^ i? ik!< lii 1^ ^ p^ ik! ci 1^ 1^ w. Radiators (b) Figure lU.A.i.l. Hybrid Element meshes for small and large problems. Small problem mesh (a) has 117 finite element nodes and 197 elements and five applicators each with seven nodes. Large problem mesh (b) has 481 finite element nodes and 900 elements and eight applicators. 23 observation sites times the number of excitations should be on the order of twice the degreesof-freedom on the mesh (e.g. [Yorkey, et al., 1987]). In these examples, exceeding this amount of observational data had relatively little effect on the reconstructed images. Accurate calculation of the forward solution requires a mesh size which is governed by the relationships: Re(kW)<0.4 and Im(kW)<0.04 (II.3.A.i.l) which for a characteristic mesh length. Ax, correspond to ten samples per wavelength and 3.5 samples per exponential decay length [Paulsen, et al., 1992b]. Tests were performed at a number of values of these dimensionless groups in order to characterize the algorithm over a broad range. The electrical properties of the two regions and frequencies considered were such that: Sampling Rate Sampling Rate re(k^ax^) f nodes/A, > im(k^ax^) Material I (Outer) 3.420* 10-5 to 0.4150 1073 to 9.753 2.173*10-5 to 0.02391 67.14 to 6.401 Material 2 (Inner) 1.110*10-5 to 0.1348 1882 to 17.11 7.250*10-5 to 0.007978 to 11.09 Region (nodes/Exp decay) 116.2 Table II.B.A.i.1. Sampling rates for discrete regions distribution on the small mesh. For this problem the overall diameter of the object being imaged ranged from 0.008908 to 0.9800 wavelengths and .1428 to 1.498 exponential decays where these numbers were based on a weighted average of the electrical properties of the body. For these tests an electrical contrast ratio of 3:1 was maintained with the outer region having the larger k^ value. Since k^Ax^ completely controls the numerical solution of the forward problem, reporting our results in these terms is the most informative since they can be related to any tissue property, frequency and mesh combination which produces these values. For example, a mesh with nodal spacing of 2 cm at a frequency of 180.5 MHz and 24 typical soft tissue properties of Cr = 70, CT = 0.0351 siemens/m, reconstructs equivalently to a tissue with er=36.5, c =0.1013 at a frequency of 1 GHz on a mesh with 5 mm resolution provided the total number of nodes is constant. The reconstructions shown in figure II.3.A.i.2 are representative of the images that are attainable. Based on the average values of for these images, the target was 0.5308 wavelengths and 0.352 exponential decays in diameter. The regularization employed was a blending of both Tikhonov and Marquardt approaches with regularization coefficients being REAL IMAGINARY 0.75 1.15 1.55 1.94 2.34 2.74 3.13 3.53 3.93 (a) (b) Figure II.3.A.i.2. Images for discrete regions distribution on the small problem - Real and imaginary parts, (a) exact image, (b) image using Hybrid Element method. 25 determined empirically.2 (The main diagonal of the Hessian matrix was augmented first by the Tikhonov regularization and then the Marquardt regularization was applied to the resultant matrix before it was solved using the method of normal equations.) The image using this method appears to converge exactly to better than 10'^ RMS electric-field error in 26 iterations (Figure II.3.A.i.3). The time per iteration using this method was 47.25 seconds. The final RMS electric field errors were also compared for known noise levels introduced in the measured data for cases of electrical contrast ratios of 3:1, 10:1, and 50:1 (Figure II.3.A.i.4). This gives a practical indication of what kind of signal to noise levels will be necessary to yield adequate images. Discrete Regions Continuous Distribution I 2 HI CO 1 -4 0 10 20 30 Iteration Number Figure II.3.A.i.3. RMS electric-field error convergence for the discrete regions distribution and continuous distribution on the small problem. This method was analyzed as a function of average jk^Ax^j because jk^Ax^j controls the performance of the imaging algorithm. To keep the analysis consistent as the ^ The Tikhonov coefficient, ( X , was 0.5 and the starting value of X for the Marquardt approach was 7.5 and was decreased by a factor of 10 at each iteration. 26 dimensionless wavenumber is varied, starting values were chosen to be a weighted average of the exact values in the two regions since Newton's method can be sensitive to starting conditions. Bounds were imposed on the values of at each iteration to restrict the number of possible solutions as a practical way of dealing with the inherent non-uniqueness. These bounds were 25% greater than the maximum and 25% below the minimum of the exact values of k^. Figure II.3.A.i.5 shows the final RMS errors for k^ vs. |k^Ax^| for both the small problem and the large problem (discussed in next section.) The condition numbers for the Hessian matrix before any regularization in the smaller problem were typically between 10 9 and 10'^� and between 10"12 and 10"14 in the larger problem. lu 2 0) u. u k- "5 0) lu (0 s E-Field Error (3:1) E-Field Error (10:1) E-Field Error (50:1) to 0) oc -160 -140 -120 -100 -80 Noise Level (dB) Figure II.3.A.i.4. Final RMS electric field errors versus noise levels in the measured data for three material contrast levels. Imposed noise is generated using a uniform probability density normalized to the radiator's power level. 27 10' ? ? Q. i lo-l m a 10-2 -Q? Small Problem Large Problem 10 -3 .0001 * 1 ? ? 1 ??i 1 ? � ? i ? 1.1 .001 i 1 i 1 1 ? ..i .01 2 2 k * Ax Figure II.3.A.i.5. RMS Error vs. |k^Ax^| for the small and large problems using the discrete regions distribution. II.J.A.ii Continuous Distribution Problem: The reconstructed images for the continuous distribution case used exactly the same mesh as the discrete regions problem discussed above. The continuous distribution for k^ is described by the following equation: where k^ is equal to the complex k^ for material 2 in the discrete regions problem which produces a maximum contrast ratio of roughly 7:1 between the largest and smallest values. Maintaining all else the same, such as observational sites, excitations, starting values and regularizations, guarantees that the Hessian matrix would be the same for both the discrete regions problem and the continuous variation problem during the first iteration. This implies 28 that their initial condition number would also be identical for each k^Ax^.^ The images shown in figure II.3.A.ii.l are representative of those attainable on this type of problem. These images were reconstructed for an average having the same value as that used in the discrete region problem. Similarly to the earlier results, the method converged to an RMS electric-field error of 10"^ but required only 11 iterations. REAL IMAGINARY Figure II.3.A.ii.l. Images for continuous distribution on the small problem - Real and imaginary parts, (a) exact image, (b) image using Hybrid Element reconstruction method. ^ At the first iteration, differences occur in the right hand side of equation II.2.3 due to differences in the observed data on the two distributions. 29 II.3.B Large Mesh II.S.B.i Discrete Regions Problem: This larger problem uses a mesh with 481 nodes and 900 elements in the finite element portion of the mesh. Also included are eight applicators with seven nodes each (Figure II.3.A.i.l(b)) which are external to the body. The electrical properties of the two regions and frequencies considered were such that; Sampling Rate Sampling Rate re(k^ax^) (nodes/X.1 im(k^ax^) Material 1 (Outer Reg.) 3.090*10-5 to 0.3744 1129 to 10.26 to 0.02158 70.64 to 6.736 Material 2 (inner Reg.) 1.006*10-5 to 0.1217 1980 to 18.00 6.545*10-5 to 0.007199 122.4 to 11.66 Region 1.961*10-4 (nodes/Exp decay) Table III.3.B.i.l. Sampling rates for discrete regions distribution on the large mesh. For this problem the overall diameter of the object being imaged ranged from 0.03538 to 3.893 wavelengths and 0.5668 to 5.945 exponential decays where these numbers were based on a weighted average of the electrical properties in the body. The image shown in figure II.3.B.i.l is representative of those attainable for this problem. The reconstruction employed a blending of both the Tikhonov and Marquardt regularizations with coefficients of 0.5 and 7.5, respectively as before. The reconstructed image is quite uniform over the two separate regions; however, the boundary is somewhat jagged even though the method appears to capture the steep gradient rather well in several locations. The convergence of the RMS electric-field errors to a value less than 10-^ is achieved after 12 iterations. The cost per iteration using this method on the mesh in Figure II.3.A.i.l(b) was 1917.77 seconds. 30 REAL IMAGINARY Figure II.3.B.i.l. Images for discrete regions distribution on the large problem - Real and imaginary parts, (a) exact image, (b) image using Hybrid Element reconstruction method. Il.S.B.ii Continuous Distribution Problem: The description of the continuous distribution compared with the discrete regions case for the large problem is analogous to that for the small problem. The exact same meshes, algorithms, regularizations, starting values and constraints are used. The continuous distribution is the same as that described in equation (II.3.A.ii.l). The images shown in Figure Il.S.B.ii.l are representative of the results that were achieved. The RMS electric-field error reached a steady-state level of 10'^ after only six iterations. The final errors are 31 similar in behavior to those noted in the smaller problem, and the condition numbers are again identical to those in the discrete regions case. REAL IMAGINARY Figure II.3.B.ii.l. Images for continuous distribution on the large problem - Real and imaginary parts, (a) exact image, (b) image using Hybrid Element reconstruction method. 32 II.3.C Simulated Measured Data In a more realistic test of the capabilities of the method, the results shown in this section are based on computation of the exact electric field data on a very fine mesh (Figure II.3.C.1) whereas reconstructions are performed on the mesh of the smaller problem (Figure II.3.A.i.l(a)). Calculation of the "measurement" data using the very fine mesh will not yield values exactly the same as those computed using the small mesh. This altering of the exact electric field data eliminates the possibility of error cancellation due to the same mesh being used in the forward solution and the reconstruction. Region 1 Region 2 Radiators siSiiiil Figure II.3.C.1. Very fine mesh with 940 finite element nodes and 1773 elements and five applicators each with seven nodes. Figure II.3.C.2(a) shows the results of this test on the discrete regions problem. The contrast ratio of the exact solution was maintained at 3:1. Clearly the method does not reconstruct the images exactly, but it does appear to detect a subregion which is distinct from the rest of the body in both the real and imaginary images. Figure II.3.C.3(a) plots vs. position along a longitudinal transect across the body for this test and the exact solution 33 where L is the object diameter. The data bears out the observations above, namely that the method generally is able to resolve the two distinct regions in both the real and imaginary parts of the image. Iteration counts were the same as those for the smaller problem. REAL IMAGINARY Figure II.3.C.2. Reconstructed images for the small problem. Exact electric-field data was computed on the very fine mesh, (a) discrete regions distribution, (b) continuous distribution. 34 A similar test was conducted for the continuous distribution problem (See figure II.3.A.ii.l) with these results shown in figure II.3.C.2(b). From these images it is difficult to tell visually whether the method was able to resolve the distribution. Figure II.3.C.3(b) shows the plots of versus position along an azimuthal cut through the body at a radius equal to 1/2 the outer radius for both the reconstructed and true distribution. For the real part, the reconstruction agrees reasonably well with the exact solution with only a modest amount of ripple. For the imaginary part, it follows the exact solution but with a higher level of ripple. 10 Re (Exact) Re (Image) Im (Exact) o- Im (Image) .1 -U.5 -0.3 -0.1 0.1 0.3 0.5 Relative X - Coordinate Figure II.3.C.3(a) Comparison of reconstructed and exact along a longitudinal transect through the body of the discrete regions distribution. L is the target radius. 35 10 Re (Exact) Re (Image) Im (Exact) Im (Image) .1 0 90 180 270 360 Theba Figure II.3.C.3(b) Comparison of reconstructed and exact along an azimuthal path around the body at a radius equal to 1/2 of the outer radius in the continuous distribution case. L is the target radius. II.4 Discussion and Conclusions A Newton's iterative scheme for the electromagnetic reconstruction of complex wavenumbers in inhomogeneous dielectric bodies has been proposed for the 2D TM polarization case. It uses the Hybrid Element method which constitutes a coupling of the Finite Element method and the Boundary Element method. The solution to the Helmholtz equation using this approach has been shown to be robust for the forward problems of both the 2D TM and TE polarization cases, as well as in the 3D [Paulsen, et al., 1988]. Thus, the Hybrid Element method creates an excellent foundation for such electromagnetic reconstructions. The simulations performed on the small problem show that this method is able to reconstruct images exactly for both the discrete regions and continuous distribution cases under the conditions of no noise and perfect match between model and data. This is important because it provides a benchmark for comparing the performances on larger and noisier problems. The algorithm was able to reconstruct images for the continuous distribution much faster than for the discrete regions distribution. In general it would appear that it is easier for this routine to image smoother functions than those with steep gradients due in part to the inherent smoothing resulting from the Tikhonov regularization. The simulations performed on the larger mesh point to potential problems with image reconstruction based on larger matrix systems in general. The technique was not able to achieve an exact image reconstruction on the large domain for the discrete regions distribution as it was for the continuous distribution. Although the images reconstructed on the discrete regions distribution were not exact, they were able to differentiate between the two regions quite well and had only minor difficulties near the steep gradients at the material interface. One significant difference between the small and large problems seems to be the condition numbers. In general the condition numbers were several orders of magnitude greater for the large problem versus those for the small problem. This is a clear indicator of a more difficult inverse problem. This also seems to amplify the assertion that the continuous distribution is easier to reconstruct than the discrete regions distribution. Application of this algorithm in Section II.3.A with noise introduced in the discrete regions problem for various material contrast ratios gives some insight into the signal to noise ratios that will be required to obtain acceptable images. It also demonstrates the performance of this algorithm on problems of very high contrast which will be encountered in biomedical applications. A more difficult test was performed in Section II.3.C, where it was shown that the algorithm was capable of reconstructing good images in less than ideal conditions. Reconstructed images for both the discrete regions and continuous distributions 37 were able to distinguish the distinct material features with only modest amounts of ripple and slight boundary artifacts. Even though the images for the continuous distribution are somewhat difficult to interpret, quantification by way of an azimuthal cut through the image clearly shows that the reconstruction is achieving the desired trends. The hybrid element method has been used successfully in accurately modeling the electric fields and provides a solid foundation for this type of image reconstruction. By itself, this algorithm still has limitations because of the computational costs and the amount of measured data required. This is overcome by combining it with the dual mesh scheme of Chapter III which directly addresses these issues. Subsequently, this allows for its use in laboratory scale image reconstructions which are discussed in Chapters IV and V. CHAPTER III. THE DUAL MESH SCHEME FOR OVERALL PROBLEM SIZE REDUCTION III.l Introduction As discussed in Chapter II and [Meaney, et al., 1995], a reconstruction algorithm for microwave imaging has been developed which uses finite element methods as its computational base. In this type of approach, imaging considerations such as resolution, computational costs, amount of measured data required, etc., are dictated by the sampling rate used to discretize the target region and the same sampling has been used both to compute the fields as well as to express the reconstructed electrical property profile. While the computation of the physical quantity of interest, in this case the electric fields, may be restricted by the constraints of the governing equation and thus require a finely discretized mesh, the reconstructed electrical properties may be quite uniform and hence representable with far fewer degrees-of-freedom. To take advantage of the disparate sampling rates that arise in the image reconstruction problem, a general method is introduced which exploits two separate meshes, one uniformly dense for field solution and a second less dense mesh for parameter recovery. Accurate calculation of electric fields requires a mesh discretization which can capture the characteristics of the physical wave propagation. Analysis and experience dictate that the meshing guidelines, Re(k^Ax^j<0.4 and Im^k^Ax^j< 0.04, which for a characteristic mesh length, Ax, correspond to ten samples per wavelength and 3.5 samples per exponential decay length are necessary in this regard [Paulsen, et al., 1992; Lynch, et al., 1985]. Experience also shows that geometrical factors are involved which require finer discretization of the body anatomy in order to characterize the rapid field variations that can occur. While wavelength considerations and body geometry effects demand mesh detail, the number of electrically distinct subregions is usually modest suggesting that the electrical property distribution should be representable through far fewer degrees-of-freedom. Hence, 38 39 a separate, more coarse mesh can be used for reconstruction purposes, while the field solution is calculated at each iteration using a finer mesh meeting at least the above criterion. This dual mesh scheme can significantly reduce the problem size during the reconstruction while still maintaining an accurate representation of the dielectric properties; hence it has important ramifications for overall computational costs and the amount of measured data that may be needed for image reconstruction. It also allows one to deploy the degrees-of-freedom associated with the reconstructed parameters independently of those of the field itself, and therefore, according to any a priori information on the electrical property distribution that may be available. It is interesting to note that others have recognized the advantages of having separate meshes for field solution and property parameter recovery. For example. Woo et al. [1994] have grouped adjacent elements into a single, larger element for reconstruction purposes in an attempt to reduce the number of reconstruction parameters in their finite element based electrical impedance tomography reconstruction algorithms. While they have been successful with this approach, the element grouping has been heuristically based and has required the reconstruction parameter mesh to consist of a subset of samples coincident with those in the field mesh. Hence, the two meshes have not been truly decoupled which has limited the generality of their approach in practical problems. In contrast, a systematic methodology is introduced for incorporating two meshes during finite element based image reconstruction. This dual mesh scheme allows the use of two arbitrarily arranged but overlaying meshes and thereby achieves a complete decoupling of the meshing issues associated with the field and reconstruction property discretizations. The power of this approach is demonstrated for microwave imaging simulations in the following sections. Integration of the dual mesh scheme with the HE reconstruction algorithm {Chapter II and Meaney, et al., 1995) is discussed especially in the context of the additional computational overhead. Image reconstructions using simulated measured data for both discrete and continuous property distributions for small and large levels of discretization are performed, in 40 each case providing direct comparison between dual and single mesh reconstructions. These illustrate that the general character of the reconstructed images are not materially different between the two methods; hence, the findings of the detailed study of the single mesh method with respect to noise and regularization strategies presented in Chapter II and [Meaney, et al., 1995] should not be substantively different for the dual mesh scheme. While only ideal simulation results are shown in this chapter, image reconstructions based on laboratory data have successfully exploited the dual mesh concept as demonstrated in Chapters IV and V and [Meaney, et al., (a)]; in fact, such reconstructions would not have even been possible without the dual mesh scheme due to the level of discretization required at these high frequencies (300, 500, and 700 MHz) and the limited amount of measured data available. III.2. Dual Mesh Method In this section the dual mesh method is described as it pertains to the hybrid element (HE) microwave imaging technique described in Chapter II and [Meaney, et al., 1995]. Implementation of the dual mesh scheme impacts two areas of the HE reconstruction algorithm: (i) the forward solution at each iteration for the electric fields which must account for the more coarsely defined electrical property distribution and (ii) calculation of the Jacobian matrix which is used to update the electrical property estimates during the nonlinear iterative process. The impact of the existence of the two meshes on the solution process for the field distribution, itself, is best understood in terms of the Galerkin weighted residual statement [Meaney, et al., 1995;Lynch, et al., 1986]: (III.2.1) where for simplicity the field has been assumed to be polarized perpendicularly to the analysis plane (i.e. the z direction); = co^|Lie* is the complex valued wavenumber and the 41 (|)j are the linear Lagrangian basis functions defined over the fine mesh, referred to here as the electric mesh. The inner products, ((?)), are also performed over the elements of this electric mesh. In the dual mesh implementation, is expanded in terms of basis functions defined over the more coarse mesh, referred to as the k^-mesh, such that: k'(x,y) = XkX(x,y) (in.2.2) n=l where Yn is a Lagrangian basis function over the k^-mesh, k^ is the k^ value at the n^h node of the k^-mesh, and M is its number of nodes. To calculate the electric fields using this formulation, the k^ values are first expressed at the nodal locations that exist within the electric mesh. To illustrate this point, Figure 111.2.1 shows a triangular element from the finer e l e c t r i c m e s h w h i c h i s e m b e d d e d i n s i d e a t r i a n g u l a r e l e m e n t f r o m a m o r e c o a r s e k ^ mesh. To calculate the values of k^ at nodes m, n, and p on the electric mesh, a weighted sum of the k^ values and the basis functions, , evaluated at the nodes /y, /2, and I3 of the k^-mesh is used; for example, k^ at node m becomes: M k^(Xm.ym)= (III.2.3) n=I While in general can have any pre-defined spatial dependency, we have chosen a linear Lagrangian interpolant as our expansion function with which to represent the k^ distribution. As long as the dielectric properties do not vary too rapidly over the k^ mesh this is a suitable representation for k^. The second impact of the dual mesh scheme occurs during the construction of the Jacobian matrix, [J], which is used to update electrical property estimates. The terms of [J] are composed of partial derivatives of the electric field at the observation sites with respect to the values of k^ at each node within the k^-mesh, where the subscript o refers to the 9kj 42 Element Electric Element (a) i3 Element A Electric Element E k Element B i4 (b) Figure III.2.1 Geometry of piecewise integrations in the dual meshing scheme, (a) electric element totally within a single element, and (b) electric element spans two elements. 43 piece of measurement data and j signifies the node number in the k^-mesh. As described in detail in Chapter II and [Meaney, et al., 1995], this requires calculating the partial derivatives with respect to the matrix equation given by (III.2.1) where the electric field has been expanded as an N term sum of basis functions, <t)j, each multiplied by an unknown coefficient (to be determined) where N is the number of nodes in the electric mesh. The individual elements of this matrix can be written as {Chapter II, Meaney, et al., 1995; Lynch, et al., 1986) ay = {-V <l)i. V (t.j) + which when differentiated with respect to the 0j ^ (III.2.4) member of the sum over the k^-mesh bases produces ^ = (4'i(|>jYf) (III.2.5) where the subscript � refers to a node on the more coarse k^-mesh, is the basis function for node f in this mesh, and the inner product is still performed over the elements in the electric mesh. Evaluating this inner product can be quite involved due to the fact that (j); and are defined on overlapping meshes which do not contain coincident node locations. To illustrate this we present two examples of cases which commonly arise in the dual mesh approach for arbitrarily overlapping discretizations. (Appendix C discusses the overhead required for this approach in more detail.) The simplest situation occurs when a given electrical element lies completely within a k^ element (Figure 111.2.1(a)). The integration over the electric element can be performed in a very straightforward manner because all three basis functions on the k^ element are non� zero over the entire electric element. Complications arise when the electric element spans more than one k^ element (Figure 111.2.1(b)). Here, the integrations must be carried out piecewise within each k^ element because and (in Fig. 111.2.1(b)) are non-zero only 44 in the shaded regions within their respective elements. (There are other cases where the integrations become even more complicated than these, such as when the electric element spans three or more elements, but the overall strategy highlighted here still applies.) When calculating (III.2.5) for during integration over the electric element E shown in Figure 111.2.1(b), it is important to note that Yy is non-zero at electric node m, is zero at nodes a and b, and zero over the entire quadrilateral abnp. Hence, integration over the quadrilateral will contribute nothing to (III.2.5) for k^. The situation is more complicated for integration over the sub-triangle and the values of the electric basis functions at the new nodes a and b are summarized below: Electric Element, E, Basis Functions Value at Node m <t'? 1 0 <t>p 0 <t>m Value at Node b <t>m (xfc,yfc) Value at Node a 4)m(^a'yo) 0 0 <t )p(x�,yj Table III.2.1. Basis function values at the nodes of the new sub triangles. Once these values are determined, the integration over the triangle mba is quite straight forward. Integration over the quadrilateral abnp, when calculating contributions to (III.2.5) for k^ is handled in much the same manner except that we first divide the quadrilateral into two sub-triangles, then calculate its new basis functions, and perform the integrations over both sub-triangles as described above. As can clearly be seen, these piecewise integrations require extra computations such as determination of the locations of the vertices a and h, calculation of the sub-triangle areas, and the evaluation of various basis functions at these coordinate locations. The most efficient way of handling this is through a pre-processing program which calculates the coordinates of the vertices along with the associated basis functions and sub-triangle areas, and stores this information for later use (See Appendix B). This increases overall memory requirements, but substantially reduces the computation time at each iteration. 45 III.3. Results The experimental evaluation of the algorithm described in the previous section is conducted on two types of imaging problems carried out on two levels of mesh discretization. The two types of imaging problems are representative of those of interest in biomedical applications and involve (i) a domain with discrete regions of homogeneous electrical properties and (ii) a domain with a continuously varying electrical property profile. The two mesh discretizations involve (1) a relatively small number of nodal degrees of freedom and (2) a larger number of mesh nodes. Note that within each of these overall mesh sizes, both a fine and a coarse mesh were utilized in order to evaluate the dual mesh approach versus the single mesh strategy. In Sections III.3.A and III.3.B the exact electric field data, with which the estimated electric fields at each iteration are compared, are computed on the fine mesh. Because this is somewhat of an idealized case. Section III.3.C evaluates these methods using exact electric field data which has been computed on a mesh with eight times the resolution of the fine mesh. This will be a more realistic and difficult test of these algorithms. We note that while we focus primarily on ideal cases with no noise, the results below indicate that the overall image reconstruction characteristics of the HE algorithm do not materially change when the dual mesh method is introduced. As a result the analysis of noisy data presented in Chapter II and [Meaney, et al., 1995] should apply to the dual mesh reconstruction contained herein as well. The images will be compared using three different criteria: (1) qualitatively by viewing the images, (2) quantitatively comparing property profiles along transects through the domain and (3) examining the relative RMS errors for the electric fields. The relative RMS errors for the electric fields are defined as: 1 N T^m cc 2x1/2 (I1I.3.1) , ^ i=i / 46 The condition numbers for the Hessian matrices (i.e. are also compared for some cases where the condition number is taken as the ratio of the largest and smallest singular values of that matrix [Golub, et al., 1989]. All calculations were performed on an IBM RISC 6000 workstation. Finally the stopping criterion for the iterative reconstruction was taken to be an electric-field RMS error of less than 10"^. III.3.A Small Mesh The small mesh problem uses an electric mesh with 117 nodes and 197 elements in the finite element part of the domain and includes five applicators containing seven nodes each as shown in Figure III.3.A.l(a). The k^-mesh has only 65 nodes and 116 elements in this case as illustrated in Figure III.3.A.l(b). III.S.A.i Discrete Region Problem: The discrete regions problem consists of an offset circular target at 3 o'clock which is embedded within a larger cylindrical background region. The dielectric properties of the two regions and frequencies considered were such that: Sampling Rate Sampling Rate Region Re(k^Ax^) Re(k^L^) (nodes/A.) Im(k^Ax^) Im(k^L^) (nodes/Exp decay) Material I (Outer Region) 3.420* 10-5 to 0.4150 0.003518 to 42.57 9.753 to 1073 2.173*10-5 to 0.02391 0.02230 to 2.453 6.401 to 67.14 Material 2 (Inner Region) 1.110*10-5 to 0.1348 0.001143 to 13.83 17.11 to 1882 7.250*10-5 to 0.007978 0.007440 to 0.8184 11.09 to 116.2 Table III.3.A.i. 1. Sampling rates for discrete regions distribution on the small mesh. In all cases an electrical contrast ratio of 3:1 has been maintained with the outer region having the larger k^ value. Since k^Ax^ completely controls the numerical solution of the forward problem, reporting our results in these terms is the most informative since they can be related to any tissue property, frequency and mesh combination which produces these values. For example, a mesh scale (i.e. Ax) of 2 cm at a frequency of 180.5 MHz and typical 47 (b) Figure II.3.A.1. Fine and coarse meshes for small mesh problem: (a) fine mesh has 117 interior nodes and 197 elements with five applicators each with 7 nodes; (b) coarse mesh has 65 interior nodes and 116 elements. Shaded region is material 2. 48 REAL IMAGINARY Figure III.3.A.i.l. Images for discrete regions distribution on the small mesh problem - real and imaginary parts, (a) exact image, (b) images using single mesh method and (c) images using the dual mesh method. 49 soft tissue properties of er = 70, a = 0.0351 s/m, reconstructs equivalently to a tissue with er =36.5, 0=0.1013 at a frequency of 1 GHz on a mesh with 5 mm resolution provided the total number of nodes is constant. Assuming an object size of 20 cm for the 180.5 MHz case, jk^L^I changes from 0.01602 to 0.001001 for a 5 cm object in this example where L is a measure of the object size (in this case its diameter). The images shown in Figure III.3.A.i.l are representative of the reconstructions attainable using both single and the dual mesh methods, respectively. These images have an average Re(k^L^) = 2.75 and an average Im(k^L^) = 0.583. The regularization employed was a blending of both Tikhonov, et al. [1977] and Marquadt [1963] approaches with regularization coefficients of 0.5 and 7.5, respectively. These coefficients were determined empirically as in [Meaney, et al., 1995]. The image using the single mesh method appears to converge exactly while the dual mesh method has some difficulty at the interface between the two materials. This is also bourne out in the convergence of the RMS electric-field errors with iteration number as shown in Figure III.3.A.i.2. The single mesh method converges to better than 10"^ RMS electric-field error in 26 iterations while the dual mesh method only reaches an error of 1.2544*10"^ after 50 iterations. The time per iteration using these methods was 47.25 seconds for the single mesh method and 27.91 seconds for the dual mesh method, respectively. The condition numbers for those two methods are compared as a function of average jk^Ax^l in Figure III.3.A.i.3. To keep the comparisons as consistent as possible as the dimensionless wavenumber varied, starting values were chosen to be a weighted average of the exact values in the two regions. This is important because the Newton's method 50 E-field Error Convergence Small Problem - Discrete Regions Single Mesh Method Dual Mesh Method 10 - 6 _i_ 10 20 30 40 50 Iteration Number Figure III.3.A.i.2. RMS electric field error behavior for discrete regions distribution on the small mesh problem using both the single and dual mesh methods. approach can be very sensitive to starting values. Both show a plateau below a |k^Ax^j of 0.0034 and a peak at |k^Ax^| = 0.073 after which the condition number decreases gradually. The lower condition numbers at the higher |k^Ax^| would imply that a lower RMS electric field error might be achieved; however, the values of |k^Ax^| which are greater than 0.4 are above the range where the forward solution can be expected to be accurate. This figure confirms the fact that the sensitivity to measurement error between the dual and single mesh methods should not be substantially different in overall character. 51 Condition Numbers Small Problem 10 11 I E 10 1 0 . 10- ?o c o u 10' Single Mesh Method Dual Mesh Method LUlJ 10' .0001 I .001 .01 .1 k**2 * delta x**2 Figure III.3.A.i.3. Condition numbers vs. |k^Ax^| for the single and dual mesh methods on the small mesh problem. III.3.A.a Continuous Distribution Problem: The continuous distribution problem consists of the same outer circular region where the profile is described by the following equation: .2, k2(p,6)= where (III.3.A.ii.l) is equal to the complex k^ for material 2 in the discrete region problem. The constant in the denominator was chosen so that the starting value for a given frequency was exactly the same as that for the discrete region problem. Keeping everything else the same, such as observational sites, excitations and regularizations, this guaranteed that the Hessian 52 matrix would be the same for both the discrete region and the continuous variation problems prior to completion of the first iteration. This implies that their condition number would also be the same at each frequency; hence the results in Figure III.3.A.i.3 also apply to this case. (At the first iteration, the differences come in due to differences in the observed data on the two distributions.) The images shown in Figure III.3.A.ii.l are representative of the reconstructions attainable using both the single and dual mesh methods with the same fine and coarse meshes shown in Figure III.3.A.1 as compared with the exact solution. These images have jk^L^j of 2.811 where the is an average of |k^| over the exact distribution given in equation (III.3.A.ii.l). Except for the observed data, these images were obtained using the exact same algorithms, regularizations, starting values and constraints used in the discrete regions problem discussed above. Similarly to the earlier results, the single mesh method converged to an RMS electric-field error of 10"^ but in only 11 iterations in this case. The dual mesh method asymptotically reached an RMS electric-field error of 7.2193*10"^ after 50 iterations. I1I.3.B. Large Mesh The large mesh problem uses an electric mesh with 481 nodes and 900 elements in the finite element portion of the domain and includes eight applicators with seven nodes each which are external to the body as shown in Figure III.3.B.l(a). The k^-mesh has 119 nodes and 220 elements in this case and is shown in Figure III.3.B.l(b). The mesh was designed to illustrate a situation where a priori information about a region of interest might be known. Areas of homogeneity are generally sparsely populated with nodes while those known to be near a boundary have a more dense sampling. This type of meshing scheme was chosen in part because of the experience with the dual mesh method on the small problem in the previous section where the results indicated a difficulty in identifying the boundary between 53 REAL IMAGINARY Figure III.3.A.ii. 1. Images for continuous distribution on the small mesh problem - real and imaginary parts: (a) exact image, (b) images using single mesh method and (c) the images using the dual mesh method. 54 two regions. The specialized discretization takes advantage of the flexibility of the dual mesh approach by allowing a denser sampling of certain areas of the domain. IIU.B.i Discrete Regions Problem: The dielectric properties of the two regions and frequencies considered for the discrete regions problem were such that: Sampling Rate Sampling Rate Refk^Ax^t Re(k^L^) Region (node.s/X,1 10.26 Material 1 (Outer Region) 3.090*10-5 to 0.3744 0.05911 to 715.3 to 112.9 Material 2 (Inner Region) 1.006*10-5 10 0.1217 0.01921 to 232.5 to 1980 18.00 Imrk^Ax^> (nodes/Expdecay) 1.961*10-4 to 0.02158 0.3742 to 41.22 6.736 to 70.64 6.545*10-5 to 0.007199 0.1250 to 13.75 to 122.4 11.66 Table III.3.B.i. 1. Sampling rates for discrete regions distribution on the large mesh. The image set shown in Figure III.3.B.i.l(c) is representative of those attainable using the dual mesh method for this problem class. It employed a blending of both the Tikhonov and Marquadt regularizations with regularization coefficients of 0.5 and 7,5 respectively. This image is quite uniform over the two separate regions and shows only modest distortion near the boundary of the two materials. The convergence of the RMS electric-field errors with iteration count is displayed in Figure III.3.B.i.2, which shows that the final error level is achieved after approximately 22 iterations. The computational time per iteration using this scheme was 219.07 seconds. In contrast, the image set shown in Figure III.3.B.i. 1(b) is representative of those attainable using only the fine mesh (i.e. Figure III.3.B.l(a)) on this same problem. The reconstruction employed the same Tikhonov and Marquadt regularizations with the same coefficients. Qualitatively, these images are not as good as those achieved using the dual mesh strategy. The RMS electric-field errors as displayed in Figure III.3.B.i.2 show that the final error level is an order of magnitude greater than that achieved with the dual mesh method. The difference is most likely attributable to the single mesh method placing a nearly 55 '� Ik-"* wT 'iiTi�U � I?f li? � f >? f I? i? � � I? U S>^^k?k?籭?i:*�i?k:<i?� i^T> 1^ ^^ ISliJ ^ 1^ 1^ 1^1^ "I f *Tii' ?^rj ^ 1^ 1^1^ ^ (a) (b) Figure III.3.B.1. Fine and coarse meshes for large mesh problem: (a) fine mesh has 481 interior nodes and 900 elements with eight applicators each with 7 nodes; (b) coarse mesh has 119 interior nodes and 220 elements. Shaded region is material 2. 56 REAL IMAGINARY Figure III.3.B.i.l. Images for discrete regions distribution on the large mesh problem - real and imaginary parts: (a) exact image, (b) images using single mesh method and (c) images using the dual mesh method. 57 E-field Error Convergence Large Problem - Discrete Regions Single Mesh Method Dual Mesh Method 10' 2w lU ?o g lU (/) S cc 0 10 20 30 40 50 iteration Number Figure III.3.B.i.2. RMS electric-field error behavior for discrete regions distribution on the large mesh problem using both the single and dual mesh methods. equal emphasis on all nodes, while the dual mesh method can emphasize certain areas in the reconstruction domain by more densely meshing those areas. It should also be noted that for the single mesh method, the iteration time was 1917.77 seconds; hence, for this larger problem, the dual mesh method also offers significant computational time savings. The condition numbers for the large mesh problem were quite similar in nature to those for the small mesh problem with the exception that they were roughly two orders of magnitude greater. This might account for the fact that the images for the larger problem were not as good as those for the smaller problem overall. 58 REAL IMAGINARY Figure III.3.B.ii.l. Images for continuous distribution on the large mesh problem - real and imaginary parts: (a) exact image, (b) images using single mesh method and (c) images using the dual mesh method. 59 III.S.B.ii Continuous Distribution Problem: The description of the continuous distribution problem relative to the discrete regions case for the large mesh situation is analogous to that for the small mesh scenario. The exact same meshes, algorithms, regularizations, starting values and constraints have been used. The continuous distribution is the same as that described in equation (III.3.A.ii.l). The image set shown in Figure III.S.B.ii.1(c) is representative of those attainable using the dual mesh method. These images are based on average values of Re(k^L^) = 134.3 and Im(k^L^) = 17.0. The RMS electric-field error behavior reveals that an error relatively close to its final value can be reached after only six iterations. The condition numbers are again identical to those shown for the discrete regions problem (e.g. Figure III.3.A.i.3). Comparison images for the exact solution and the single mesh approach are shown in Figures III.S.B.ii.1(a) and (b), respectively. III.3.C Simulated Measured Data In a more realistic test of the capabilities of both the single and dual mesh methods, the reconstructions appearing in this section utilize exact electric field data computed on the very fine discretization shown in Figure III.3.C.1. Since the calculation of the electric fields in the reconstruction process can never identically match the simulated measurement data in this case even if the exact electrical properties were supplied, this situation can be viewed equivalently to adding noise to the measurement data. We note that on average a 5% difference existed between the observed fields calculated on the two meshes when both contained the exact electrical properties (with this difference reaching as high as 10% for some sites); hence, the images shown below are representative of the impact that additive noise would have at this level. Combinations of the fine (Figure III.3.A.1(a)) and coarse meshes (Figure III.3.A.l(b)) from the small mesh problem are then used to perform the actual image reconstructions. Based on combinations of these meshes, three cases are studied which are referred to as: (I) 60 the fine-fine test, (2) the fine-coarse test, and (3) the coarse-coarse test. Case 1 (fine-fine) uses the single mesh method where both the electric fields computed and the values reconstructed at each iteration are obtained on the fine mesh (Figure III.3.A. 1(a)). Case 2 Figure III.3.C.1. Very fine mesh having 940 interior nodes and 1773 elements with five applicators each with 7 nodes. Shaded region is material 2. (fine-coarse) uses the dual mesh method where the electric fields are computed on the fine mesh and the values are reconstructed on the coarse mesh (Figure III.3.A.l(b)). Case 3 (coarse-coarse) is also a single mesh method except that it uses the coarse mesh (Figure III.3.A.l(b)) for both calculations as opposed to the fine mesh used in Case 1. Figure III.3.C.2 shows the results of these three tests for the discrete regions problem. The contrast ratio between the two regions was maintained at roughly 3:1. Clearly, none of the methods reconstructs the images exactly, but the first two, the fine-fine test and the finecoarse test, appear to (at least in the real component of the image) detect a subregion which is distinct from the rest of the body. Figures III.3.C.3(a) and (b) show plots of position along a longitudinal transect across the body for all three cases as well as for the vs. 61 IMAGINARY REAL (a) (b) (0 Figure III.3.C.2. Reconstructed images for the discrete regions distribution on the small mesh problem where the exact electric field data was computed on the very fine mesh: (a) fine-fine test, (b) fine-coarse test, and (c) coarse-coarse test. 62 k**2 * L**2 vs. Position Longitudinal Cut 8.00 7.00 I� I UJ QC 6.00 5.00 4.00 3.00 Exact fine-fine 2.00 fine-coarse coarse-coarse 1.00 -0.5 -0.3 -0.1 0.1 0.3 0.5 Figure III.3.C.3(a). Re(k L ) along a the longitudinal transect across the body with the discrete regions distribution on the small mesh problem for the fine-fine, fine-coarse, and coarse-coarse reconstructions, compared with the exact profile. exact profile. These figures bear out the observations above, namely that the first two methods generally are able to resolve the two distinct regions in the real part of the image, and in the imaginary component they seem to detect the distinct regions but the results are much noisier. In both cases the coarse-coarse test does not do an adequate job of resolving the two regions. With regard to the reconstruction times, the fine-fine, fine-coarse and coarse-coarse tests required 45.8, 27.67 and 11.67 seconds per iteration, respectively. From this data, the fine-coarse test offers a good compromise between fast reconstruction time and resolvable image profiles. In situations where there is a limit on the amount of observational data available, the comparison of the fine-coarse and coarse-coarse tests becomes very important. Both use exactly the same amount of measured information, and it is clear that 63 the former does recover the profile while the latter does not. This can have significant consequences on the overall performance of an imaging system. 1.20 Exact fine-fine fine-coarse coarse-coarse -0.5 -0.3 -0.1 0.1 0.3 Relative X - coordinate Figure III.3.C.3(b). Im(k^L ) along a the longitudinal transect across the body with the discrete regions distribution on the small mesh problem for the fine-fine, fine-coarse, and coarse-coarse reconstructions, compared with the exact profile. Similar tests have been performed on the continuous distribution problem with these results shown in Figure III.3.C.4. By qualitatively viewing these images it is difficult to determine whether one test was able to resolve the image profile significantly better than the others. Figures III.3.C.5(a) and (b) show quantitative plots of versus position for an azimuthal cut through the body at a fixed radius equal to one-half the outer radius of the domain. For the real component, the fine-fine test agrees well with the exact profile with only a modest amount of ripple. The fine-coarse test also agrees quite well with the exact distribution except near 90 degrees. The coarse-coarse test follows the general trends in the continuous profile but has several large fluctuations from the exact distribution. For the 64 REAL IMAGINARY (a) (b) (c) Figure III.3.C.4. Reconstructed images for the continuous distribution on the small mesh problem where the exact electric field data was computed on the very fine mesh: (a) fine-fine test, (b) fine-coarse test, and (c) coarse-coarse test. 65 imaginary component, both the fine-fine and fine-coarse tests follow the exact distribution but with a slightly higher level of ripple than in the reconstructed real part. Again, the coarse-coarse test has more significant fluctuations from the exact profile than do the other two tests. k**2 * L**2 vs. Position Longitudinal Cut 8.00 7.00 � * 6.00 5.00 i 4.00 CC 3.00 Exact fine-fine 2.00 fine-coarse coarse-coarse 1.00 -0.5 -0.3 -0.1 0.1 0.3 0.5 Figure III.3.C.5(a). Re(k L ) along an the azimuthal cut around the body at a fixed radius equal one-half of the outer radius with the continuous distribution on the small mesh problem for the fine-fine, fine-coarse, and coarse-coarse reconstructions, compared with the exact profile. III.4 Discussion and Conclusions In Chapter II and [Meaney, et al., 1995] a Newton's iterative scheme for the electromagnetic reconstruction of complex wavenumbers for inhomogeneous dielectric bodies was described for the 2D TM polarization case using the hybrid element method. Here, this particular implementation of the dual mesh scheme is well suited for integration with the hybrid element method which constitutes a coupling of the finite element and the boundary element methods. The solution to the Helmholtz equation using this approach has 66 1.20 100 Exact fine-fine fine-coarse coarse-coarse 0.5 Relative X - coordinate Figure III.3.C.5(b). Im(k L ) along an the azimuthal cut around the body at a fixed radius equal one-half of the outer radius with the continuous distribution on the small mesh problem for the fine-fine, fine-coarse, and coarse-coarse reconstructions, compared with the exact profile. been shown to be robust for the forward solutions of both the 2D TM and TE polarization cases, as well as in the 3D [Paulsen, et al., 1988]. Thus, the hybrid element method creates an excellent foundation for such electromagnetic reconstructions. In the formulation of this integrated reconstruction algorithm, the dual mesh scheme has been developed to decouple the sampling of reconstructed wavenumbers from the degrees of freedom used to compute the electric field distributions. The dual mesh method of reconstruction offers substantial flexibility by enabling dense discretization of critical areas of interest and coarse discretization in other areas while still maintaining accurate computation of the electric field at each iteration. This can be critical for high frequency imaging especially in the case where the amount of observational data from which image 67 reconstruction is to be accomplished is limited. It is also important under conditions where a priori information may be available to guide discretization decisions and can be used to improve the computational costs of reconstruction without significant compromise to image quality. The simulations performed on the small mesh problem show that the single mesh method is able to reconstruct images exactly for both the discrete regions and continuous distribution problems. This is important because it provides a benchmark for comparing the performances of the dual mesh method on the small mesh problem and subsequently the single and dual mesh methods on a large mesh problem. For the dual mesh method on the small mesh problem, the reconstructed images are not as good as those of the single mesh method. In these cases the condition numbers for the small mesh problem using both the single and dual mesh methods are quite similar; therefore, this is probably not the most likely explanation for why the dual mesh images are not as good as their single mesh counterparts. For the discrete regions image, differences in image quality may be attributable to the fact that the coarse mesh is uniformly sampled over the whole body without enough fine discretization near the interface between the two regions. The results show that the dual mesh method is ideal for problems which have distinct regions where k^ is known to have sharp contrasts in which case areas of uniformity can be coarsely discretized whereas areas of rapid change can be finely discretized. As expected the dual mesh method has more difficulty with the continuous distribution problem where k^ is continually varying over the whole body since accurate representation of this distribution requires a more uniform, finely discretized mesh. The simulations performed on the large mesh point to potential problems with larger matrix systems in general. Neither the single nor the dual mesh method was able to achieve an exact image reconstruction for the large domain. One significant difference between the small and large problems seems to be the condition numbers. In both the single and dual 68 mesh methods, the condition numbers were several orders of magnitude greater for the large problem versus those for the small problem. This is a clear indicator of a more difficult inverse problem. In this case the coarse mesh was specifically designed to be more densely discretized near the boundary of the two regions largely in response to the difficulty this method had in distinguishing the interface in the small mesh problem. Relatively good images were obtained using the dual mesh method. These were by no means exact, but they did differentiate well between the two regions with relatively constant values within each area. The single mesh method was not as successful. It still managed to distinguish the two regions, but the image profile was hardly constant over the inner material. The single mesh method was also the beneficiary of a condition number several orders of magnitude less than that for the dual mesh method. This suggests that the dual mesh method may have the effect of weighting the error function (in a least squares context) more in areas densely meshed as opposed to the single mesh method which applies more equal weight to all regions. Application of the single and dual mesh methods on the fine-fine, fine-coarse, and coarse-coarse test cases shows that they are both capable of reconstructing the desired image in this less than ideal, noisy case. Both reconstructed images from the fine-fine and finecoarse tests were able to distinguish the distinct material features with only modest amounts of ripple. The images from the coarse-coarse test did not do an adequate job of distinguishing these material features. This shows that simply using the single mesh method on the coarse mesh instead of the fine mesh will seriously degrade the image quality. However, using the dual mesh method instead of the single mesh method with the fine mesh will not seriously deteriorate the image quality. The fine-coarse test also provides a substantial computational time savings over the fine-fine test, albeit not as much as the coarse-coarse combination, while still maintaining good image quality. Therefore, the fine- 69 coarse approach provides a good compromise between reducing computation time and maintaining image quality. The reconstruction algorithm using the hybrid element method from Chapter II and the dual mesh method described in this chapter has demonstrated the ability to reconstruct images for a wide variety of simulated problems, especially for larger and noisier cases. In fact, experience to date indicates that the dual mesh method is essential for obtaining reconstructions from laboratory data as described later in Chapters IV and V where one is interested in high frequency illumination of targets with embedded heterogeneities which are multiple wavelengths in size [Meaney, et al., (a)]. Under these conditions, the need for accurate field computation drives the image reconstruction problem to excess without invoking an alternative discretization for the reconstruction parameter property profile; thus making realistic lab-scale imaging impractical without dual meshing. The ability to have arbitrarily overlapping discretizations offers a convenient and highly flexible framework for pursuing image reconstruction based on a practical setting such as presented in laboratory scale tests demonstrated in Chapters IV and V and ultimately in the clinical environment. CHAPTER IV. THE ACTIVE MICROWAVE MEASUREMENT SYSTEM IV. 1. Introduction The first step necessary for the eventual utilization of microwave based thermal imaging is to reconstruct good static images using microwave technology, and when this has been achieved, proceed to the problem of dynamic imaging which will yield data on changes in the temperature field distribution. Attempts have been made to image biological tissues using active microwave imaging with varying degrees of success. One system developed by Jofre et al. [1990] and Broquetas, et al. [1991] performs microwave imaging at 2.45 GHz. It has shown some promise, especially in the area of differential imaging; however, its major limitation appears to be in detecting the large electrical property contrasts that naturally occur in biological tissue and there have been some difficulties in resolving both the resistive and reactive components of the electrical property distributions as well. These problems appear to be due primarily to the use of the Bom approximation during image reconstruction - a method that has been shown to be limited in the presence of high contrast objects [Slaney, et al., 1984]. A related approach, which attempts to exploit conventional x-ray CT reconstruction algorithms by discriminating the straight line from multipath microwave signals (1-2 GHz) received through chirp-radar techniques, has also been reported [Miyakawa, 1993]. This work has demonstrated spatial resolutions on the order of 1 cm for static targets. Imaging of temperature changes in laboratory phantoms has also been achieved using difference methods, although quantitative measures of the actual temperature change have not been possible to date. Similar qualitative thermal imaging experiments have been reported in [Mallorqui, et al., 1992]. These works indicate that there is limited experience with the integration of a hardware data acquisition system with a reconstruction algorithm in order to produce actual images using microwave excitations in the biomedical context. Further, the ability to simultaneously recover absolute and quantitative values for both the resistive and reactive components of the imaged 70 71 electrical property distribution has remained elusive. The efforts reported in this chapter seek to advance the medical microwave knowledge base on both of these important fronts. From past research, plus our own simulations, it is clear that there are three key ingredients to a successful microwave imaging system: 1) the electric fields must be measured accurately and reliably by the measurement system: 2) the numerical algorithm for calculating the electric fields based on the electrical property distribution and electromagnetic excitation must be accurate; and 3) there must be some technique readily available which reduces the amount of measured data required without compromising image quality. This last item is important because it can seriously impact the complexity and cost of the measurement system. Items 2 and 3 have already been discussed in Chapters II and III, respectively. Item 1 is successfully addressed in this chapter, with the end result of the combination of all three being good reconstructions of the electrical property distribution in a simple experimental geometry where simultaneous recovery of both resistive and reactive components of the electrical property profile have been accomplished within an absolute imaging framework. The coherent measurement system (Section IV.2) retains all the phase and amplitude information of the received signals with the variations due to differences in the transmit/receive channels calibrated out at the data processing stage [Skolnik, 1980]. After the data acquisition is completed, the information requires modification (Section IV.3) because it is recorded by a system which radiates in three dimensions, while the image reconstruction algorithm is strictly two-dimensional at this stage. Once the measurements have been processed, the data are then used by the reconstruction algorithm to obtain profiles of the electrical properties. The image reconstruction algorithm using the hybrid element method is discussed in detail in Chapter II and [Meaney, et al., 1995] with its integration with a dual mesh scheme described in Chapter III and [Paulsen, et al., (a)]. The image reconstructions presented here are the first experimentally-driven confirmations of the practical merits of this reconstruction algorithm and its dual mesh concept which have been examined theoretically elsewhere. 72 In the following sections, the design and implementation of the microwave measurement system, including discussions of the various trade-offs concerning sensitivity and resolution and how these affected the design choices, is presented. The calibration procedure is also discussed in the context of using measured data from a system whose sources radiate into a 3-D space in a 2-D reconstruction algorithm. The results demonstrate the success achieved in using such a calibration procedure by comparing the data from the forward solution for the electrical fields in a known homogeneous medium with the measured data from the actual system. With this accomplished, the reconstruction algorithm is applied to several phantoms to show that objects on the order of a half wavelength can be quantitatively imaged (both resistive and reactive components) while objects of smaller size result in qualitatively correct images. IV.2 Hardware Design A 4 transmit channel/ 4 receiver channel measurement system has been implemented for use over the broad frequency range of 300 to 1100 MHz. This span allows for the determination of an optimal operating frequency based on trade-offs of resolution size and signal penetration depth and also presents the opportunity of having multi-spectral excitations which could have important biological implications such as in discrimination between blood flow and temperature related electrical property changes. The design utilizes low cost, wellcharacterized microwave components to reduce system integration problems. Its modularity allows for testing multiple antennas in various positions which is essential in determining an optimal configuration. 1V.2.A Microwave Circuit Design IV.2.A.i Modulation Scheme In designing the microwave imaging system, a major requirement is that it be coherent such that both the real and imaginary parts of the electrical properties can be extracted [Skolnik, 73 et al., 1980]. This design is a single side-band heterodyne system (Figure IV.2.A.i.l) for which the EG&G Model 5210 Lock-In Amplifier (EG&G Princeton Applied Research, Princeton, NJ, USA) is well suited. The reference continuous wave (CW) signal is modulated with an intermediate frequency (IF) using a single-sideband, carrier-suppressed upconverter. The received modulated signal is then compared with the original reference CW signal using a single mixer and the resulting IF is fed into the lock-in amplifier which extracts the I & Q components and also performs the A/D conversion. This arrangement is preferable because the I & Q discrimination can be accomplished in the KHz frequency range which is much more accurate and less expensive than attempting the phase and amplitude measurements directly from the microwave signal. The system has been constructed in a modular fashion such that all the microwave components reside on a microwave plate and the antennas are suspended from an antenna plate in such a way that the two plates are stacked on top of each other with appropriate spacers above the experimental imaging tank (Figure IV.2.A.i.2). This economical use of the available space makes it possible to keep the transmission lines between the components of the two levels as short as possible, thereby reducing potential phase shifts due to continual mechanical stresses during measurements. It allows ease of assembly or disassembly of the system in the event that components require alteration or servicing. To better illustrate the modularity of this system. Figure IV.2.A.i.3 shows the microwave circuit plate which can be easily removed to make changes in the antenna configuration or to perform any system servicing that might be necessary. A photograph of the isolated antenna plate with the various antennas suspended from it is shown in Figure IV.2.B.2, and is discussed in Section IV.2.B. 74 ator 1 L. - SD JUXUOi L .. P a s 3 � ?a I u 1 a Q 5 K19 SD � J M X WOS in Q(Jc/5 ill Figure IV.2.A.i.l. Schematic Diagram of the Data Acquisition System. Shows the microwave components of the coherent, heterodyne design in reference to the sahne tank, EG&G Lock-In Amplifier, and Computer. 75 Figure IV.2.A.i.2 Photograph of the Data Acquisition System. Shows the modular components of the coherent, heterodyne design in reference to the antenna plate and the saline tank. Figure IV.2.A.i.3 Photograph of microwave circuit plate. Shows modulator circuitry and transmitter multiplexer to the left and the receiver circuitry and multiplexer to the right. 76 IV.2.A.a Transmitter /Receiver Design: The reference CW signal is divided into two parts, one going to the transmitter and the other going to the receiver for later comparison. The former is amplified and attenuated prior to modulation to increase the isolation between the transmitter and receiver components. The transmitter signal is then modulated by a single-sideband, carrier-suppressed modulator, Merrimac IGM-9B Series, I and Q Phase Demodulator (Merrimac, West Caldwell, NJ). This yields an output signal of -16 dBm with the sideband and carrier suppressed to -25 dBc or better over the frequency range 300 to 1100 MHz. Greater sideband and carrier suppressions are generally available over narrow bandwidths; however, we have chosen this modulator for our preliminary investigations because of its wide bandwidth. This allows us to study the frequency performance of our system and contributes to the selection of an optimal operating frequency as well as keeping the option of multi-spectral excitation which may be of interest in vivo. The phase noise contributions due to these suppressed signals are discussed in Section IV.3. A. The modulated signal is fed into a power amplifier and then a multiplexer (SP4T switch), M/A-COM SW-255 (M/A-COM Microelectronics Division Lowell, MA) which selects the transmitting channel. The multiplexer has isolation between 50 and 60 dB across the band. Extra isolation was required to eliminate all leakage signals caused by the antenna configuration and the high attenuation of the saline bath. This was accomplished by incorporating a M/ACOM SW-215, switch at the four outputs of the multiplexer which was necessary to achieve the linear dynamic range goal of greater than 115 dB. Each receiver channel is comprised of a cascade of broadband components (M/A-COM) which include; (1) a low noise amplifier ? AMC-174, (2) an isolating switch (SPST switch) ? SW-215, (3) a multiplexer ? SlV-255, (4) two 5 bit digital attenuators ? AT-104, (5) two amplifiers ? AMC-184, and (6) a mixer ? MDC-149. The linear dynamic range is bounded on the low end of the spectrum by the signal-to-noise ratio of the system (System noise figure 77 is discussed in Appendix D) and on the high end by the saturation of the component most likely to see power levels approaching its 1 dB compression point. In this design, that component is the mixer. Attenuating the signal in known increments using the digital attenuators allows the mixer to operate in the linear region while still maintaining accurate phase and amplitude data. Finally, the resulting IF is fed into the lock-in amplifier for I & Q discrimination and eventually A/D conversion. Control of all the switches, attenuators, microwave source and lock-in amplifier is performed by a Dell Dimension 486 computer. IV.2.B Antenna Design and Arrangement Water-filled waveguide apertures were chosen as the transmit antennas. They are straightforward to design and are quite well understood for simulation purposes which is important when trying to construct an accurate numerical model of the system. They yield a certain degree of directivity without increasing size, an important consideration in reducing the signal loss from the transmitter to receiver in a lossy medium. For the receive antennas, monopoles were chosen because they closely emulate a point electric field detector in the 2-D plane. They were constructed by exposing a quarter wavelength (in the medium) of the center conductor of a semi-rigid coaxial cable. In a low loss medium such as air or deionized (DI) water, these monopole antennas are notoriously poorly behaved, exciting currents along the outside of the coaxial cable [Krause, 1950]. Methods of improving the performance of these antennas generally take two forms: (1) incorporating a balun in an attempt to properly terminate unwanted currents, and (2) use of an absorptive material such as a ferrite ring around the outside of the cable to attenuate the currents. Both of these generally have limited bandwidths. Our approach is to capitalize on the high attenuation of the surrounding saline solution. Figure IV.2.B. 1 shows the characteristic impedances of the monopole antenna in DI water and a 0.9% saline solution. The plot for the DI water case is quite uneven, displaying characteristics which are generally attributable to the lack of any sort of balun. The plot for the saline solution case is dramatically different; it is quite well-behaved over a broad bandwidth, with a nominal return loss of 9 dB for the frequency range of 300 to 1100 MHz. Return loss measurements for the water-filled waveguide radiators in DI water and saline are shown in Appendix E. Characteristic Impedance IMonopole Antenna 0 -10 -20 DI Water Saline Solution (0.9%) -30 100 300 500 700 900 1100 Frequency (MHz) Figure IV.2.B.1 Return Losses of the quarter wave monopole immersed in two different media; (1) DI water and (2) 0.9% saline solution. Similar advantages are also exploited to minimize mutual coupling of the various antennas. When the return loss of an antenna is measured in the presence of another receiver antenna (location separations are a minimum of 2.8 cm), there is no detectable change. This is due primarily to the high attenuation of a radiated signal during its round trip from the one antenna to another and back again. The large tapered structures (see discussion below) used to improve the antenna positioning also do not contribute to the mutual coupling because the radiated signals must travel out of the plane of interest and then be reflected back into that plane, thus yielding an even longer transit distance with even greater attenuation. This 79 source/detector arrangement is fairly close to the TM electromagnetic mode where the electric field is predominantly in the vertical direction (z) for semi-infmitely long cylinders and hence is appropriate for 2D image reconstruction (This is not entirely true and full compensation for this is discussed in Section IV.3). For completeness, a discussion of the signal variations due to changes in the antennas vertical positioning is included in Appendix F. Figure IV.2.B.2. Photograph of the antenna plate with the transmit and receive antennas suspended from the plate. The antenna plate with antennas has been removed from the tank for illustration purposes. During data collection, it is important that the positions of the various antennas be known quite accurately in order to be represented with the numerical model. If a dimension is in error by 1.0 centimeter it can contribute roughly a 55� phase shift at 500 MHz. It is also important that the fixtures used for alignment not be so cumbersome as to perturb the electric field patterns, themselves. The design suspends the antennas from the antenna plate situated above the experimental tank (Figure IV.2.B.2). Plexiglas? rods are fixed at known locations 80 on top of the waveguides and suspended through holes in the antenna plate. A tapered metal sleeve is soldered to the coaxial cable of the monopole antenna for stability and is also mounted to the antenna plate. All mounting holes in the antenna plate were drilled by a CNC milling machine to ensure antenna position accuracy. IV.3 Calibration Since the electric fields at the apertures of the transmit and receive antennas are polarized in the z direction, simulation using a 2D TM mode in the reconstruction algorithm becomes possible. In fact, this is only true after some modifications of the raw measured data. (An additional calibration procedure, to calibrate out the amplitude and phase differences between transmitters and receivers, is discussed in Appendix G.) All antennas have three dimensional radiation patterns and hence have a 3-D "free space loss factor" (FSLF) which is (X./47iR)^ [Balanis, 1982]. Those in the numerical algorithm have a 2-D FSLF which is A,/47iR. TO be able to compare the data accurately, the 3-D FSLF needs to be subtracted from the measured data and a 2-D FSLF added in. This requires knowledge of the complex electrical permittivity and a reference point from which to measure R, the distance from the phase center. For certain horn antennas, there exist algorithms to calculate the phase center [Ohtera, et al., 1975]. For most applications this is not possible, and experimental methods must be used [Dyson, 1967]. For points in the far field of an antenna, and well within its main lobe, the phase of the electric field measured at a point some distance R from its phase center varies linearly as a function of R. Thus, given the location of the measurement point and the phase of the electric field at that point (Figure IV.3.1), a relationship involving the phase center, to within some error due to noise and measurement imprecision, can be written as: Pi=Pc+tpR+ei (IV.3.1) 81 Pc (0,0) Figure IV.3.1 Diagram used for calculation of phase center which is used to calculate \htfree space loss factor , the electrical properties of the homogeneous medium and phase and amplitude offsets. where P, is the phase at the measurement point, Pc the phase at the phase center, tp is the slope of the phase with distance after taking into account the cyclical discontinuities of the phase, and Ej is a random error term associated with each measurement. This equation can be manipulated to give: Pi=Pc+tpV(xi-Xc)^+y? +ei (IV.3.2) where (x,, y;) are the coordinates of the measurement point and (xc ,0) are the coordinates of the phase center. This equation has three unknowns and three knowns along with the error term and requires measurements from at least three locations to determine Xc, Pc, and Xp. An 82 improved approach takes measurements at many locations and performs a least square regression to determine these values [Mines, et al., 1990]. Since it is difficult to linearize equation (IV.3.2), a simple three dimensional search over the three unknowns can be performed until the squared error is minimized. Once the phase center is known, the (X/47[R)^ term can be subtracted from the amplitude of the measured data for points directly in front of the radiator. This leaves data whose attenuation is due solely to the conductivity of the medium. Using a linear regression analysis on this reduced data set, the slope of the amplitude loss, Ta, without the influence of the FSLF can be determined: Ai=Ac+XA(xi-Xc}fei (IV.3.3) where Ai is the amplitude at the measurement point and Ac the amplitude at the phase center. Knowledge of Tp and Xa makes it possible to accurately determine the relative dielectric constant and the conductivity of the homogeneous medium [Liao, 1980]. This is critical because the HE method used in image reconstruction requires that the material properties outside the FE region be known. It also provides a good starting value for use in the iterative reconstruction algorithm. Once the electrical properties are known for the homogeneous medium, the same regression analyses are run for data calculated using the HE method forward solution. This yields calculated electric field data for those sites corresponding to the measurement locations. In general there will be an arbitrary phase and amplitude offset between the measured and calculated data. These phase and amplitude offsets for all measurement sites are averaged to give overall offset values which will be used to modify the measured data when targets are present. This along with the subtraction of the 3-D FSLF and the addition of the 2-D FSLF must be performed to the measured data before the image reconstruction algorithm can be executed. 83 IV.3.A Sources of Measurement Errors Several mechanisms exist which can contribute to overall errors in the measurement system. The primary of these include: ( 1 ) Thermal Noise - The available thermal noise is kTB; k is Boltzman's constant, T is the temperature in degrees Kelvin, and B is the system bandwidth. For this system, it is roughly -135 dBm. At present this is not a major contribution to measurement error since the dynamic range is limited to 115 dB primarily due to leakage signals. ( 2 ) Leakage Signals - The leakage signals are primarily due to the level of isolation between the transmitter or receiver channels that are ON and those that are OFF. In this case the isolation between channels is nominally 50 dB (due to the multiplexers - Figure IV.2.A.i.1) with an additional 60 dB (due to the 4 SPSTs - Figure IV.2.A.i.l) added for increased isolation. At present these leakage signals are not performance-limiting since we do not yet utilize the full dynamic range of our system, but this will have to be monitored carefully in the future as we involve measurement sites at larger distances from the transmitter. ( 3 ) Phase Error Contributions due t o Suppressed Sideband Signals - In the worst case scenario, where the relative phases of the undesirable sidebands are in quadrature with the primary transmitted signal, the phase error would be approximately 10"25/20 radians, or 3.2 degrees for the case of -25 dB sidebands. This is only possible when the signals have traveled many wavelengths; however, in our system, the maximum distance that a signal travels is roughly six wavelengths. Given the existing sideband level of -25 dBc, this would produce only a 0.0002 degree phase error which would not be detectable with our current system configuration. ( 4 ) Antenna Positioning Inaccuracies - In the high dielectric, lossy medium of interest, positioning inaccuracies contribute phase errors of 5.5�/mm and amplitude errors of 0.8 dB/mm. Inaccuracies in the machining of the antenna mounting plate have manifested themselves in two forms: (1) lateral positioning errors on the order of 0.5 mm have been observed, and (2) rotational errors of the transmit antennas on the order of 1" have also 84 been noted. These can contribute up to 0.4 dB amplitude and 2.7� phase errors in the former case, and 0.5 dB amplitude and 1.5" phase errors in the latter situation. Efforts will be necessary to improve this in the next generation system once an optimal antenna arrangement has been determined. (5) Mechanical Stress of Cables - Flexible coaxial cables (EZ Form Cable Corp., New Haven, CT) form the connections between the antennas and the microwave component module with the connectors being continually engaged and disengaged from the antennas (both transmit and receiver) in various positions in the present manual mode of operation. Data from the manufacturer predicts that with three successive bends of 90� in a cable length (typical of that found in our present configuration) we should observe less than 1� in phase error at 500 MHz. The sum of these predicted measurement error contributions is on the order of 0.9 dB in amplitude and 5.2" in phase. This error estimate will be compared in the next section with the actual observed measurement data. IV.4 Results The results presented in this section were obtained with the microwave measurement system configured as shown in Figure IV.2.A.i.2 using an operating frequency of 500 MHz with an IF of 5 KHz. The experimental set-up consisted of four transmit and four receive antennas such that, at any one time, 16 measurements could be recorded. Eight possible combinations of these antenna configurations existed making the total amount of data that could be recorded for any one phantom experiment equal to 128 measurements. The receive antenna locations consisted of two sets of eight positions equidistant from the center of the phantom target (radii of 7.37 and 3.68 cm, respectively). The eight transmit antenna locations positioned the face of the antenna apertures 11.7 cm from the target center. Both transmit and receive antennas were spaced at equal angular intervals around the periphery of the phantom object. 85 The first experiment performed was to determine the accuracy with which the numerical computations can predict the measured fields when the electrical properties are known. This information provides validation of both the calibration procedures identified in the previous section and the model of the radiating sources used in the computational procedures. Since the image reconstruction algorithm makes repeated use of these calculations as it updates an initial estimate of the electrical property distribution towards the final profile, accurate calculation of the measured field distribution will be essential for successful image formation when the electrical properties are unknown. Measured data is used to calculate the electrical properties of the medium, the phase center and the phase and magnitude offsets between the measured and calculated data as described earlier. Using this information, the comparison between the calculated data and the measured data has been performed by subtracting and adding the appropriate FSLF's. Table IV.4.1 Node Number 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 Measured Measured Calculated Calculated Amplitude (dB) Phase (degrees) Amplitude (dB) Phase (degrees) -15.97 -170.53 -15.68 -169.68 -34.20 30.41 -34.06 33.59 -52.28 -14.49 -52.33 -17.04 -63.70 -155.67 -64.24 -156.22 -67.53 -68.90 -67.40 -81.48 -64.54 -148.06 -64.24 -156.72 -53.04 -15.37 -52.33 -17.04 -34.90 33.23 -34.06 33.59 27.96 -29.55 -29.84 29.36 109.07 -35.43 -34.73 107.34 -45.03 -109.40 -45.27 -97.59 -52.16 29.58 -52.63 29.71 77.73 -55.22 -55.39 76.88 -52.59 29.61 -52.63 29.71 -97.99 -45.28 -45.27 -97.59 -35.23 107.66 -35.43 107.34 -128.94 -42.75 -42.41 -127.74 Table IV.4.1 Comparison of the calculated and adjusted measured data for the seventeen measurement sites in a homogeneous medium. 86 shows the comparison of field amplitudes and phases at 17 locations in a homogeneous phantom having electrical properties; = 74.4, o = 1.77 s/m. The agreement is quite good at all comparison points. The average phase difference was 3.2 degrees and the average amplitude difference was 0.32 dB which compares quite favorably with the expected sources of measurement error which have been estimated to be (see Section IV.3.A) 0.9 dB and 5.2�. This data indicates that the calibration procedures and the numerical waveguide model allow the radiation of the physical 3D antennas within the experimental environment to be simulated using a 2D calculation which can form the basis of a 2D material property reconstruction process. The numerical model used to perform the calculations appearing in Table IV.4.] has eight transmit antennas and utilized the finite element mesh shown in Figure IV.4.1(a) which consists of 2342 nodes and 4538 triangular elements along with 560 boundary elements representing the external radiators. The measured data was obtained at the locations identified in Figure IV.4.1(b). The computational model is based on the HE scheme outlined in Meaney, et al. [1995]. Because of the lossiness of the background medium, the presence of the seven other transmitters in the numerical model had a negligible effect on the calculated electric fields for the case when a single transmitter was being modeled as demonstrated by the data in Table IV.4.1. The simultaneous presence of all radiating positions in the computational model affords certain efficiencies in the reconstruction algorithm as described in Meaney et al. [1995]. Note that the reconstruction parameters are obtained on the more sparsely sampled mesh of Figure IV.4.1(b) and have linear spatial variation inbetween the discrete nodal positions shown in this figure. While the dual mesh reconstruction algorithm which uses both meshes in Figure IV.4.] has been described theoretically elsewhere [Paulsen, et al., (a)], the results shown here represent the first experimental confirmation of this concept. The degree to which this scheme is successful experimentally is significant since reconstruction using a single mesh method would require a domain sampled as finely as that shown in Figure IV.4.1(a) (due to wavelength considerations) and would not be possible with the present hardware system given 87 Homogeneous Region Radiator Target R e g i o n Measurement Sites Figure IV.4.1 Meshes used in the reconstructions, (a) is the electrical mesh used in the forward solut ion. It includes the finite element mesh and the boundary elements. (b) is the k mesh for reconstructing the electrical properties. Hiqli! ighted are the measurement sites. 88 the hmited amount of observational data that it can provide (i.e. 128 measurements per image reconstruction). As an initial series of experiments to demonstrate the imaging capability of the 4 channel prototype system, tests were performed on three cylindrical phantoms of diameters 1.55, 2.51, and 4.23 cm (denoted XI, X2 and X3). In all three cases the phantoms were positioned at the center of the target region and were constructed from a polyester resin with hardener, acetylene black and aluminum powder as prescribed by Guy [1971] for making equivalent phantoms for bone and fat tissue. As such, these phantom targets offered a very high contrast with the surrounding medium, their electrical properties being reported to be Er = 4.6 - 6.2 and c = 0.17 - 0.55 s/m [Guy, 1971]. The phantom targets were machined to lengths of 17.83 cm to simulate a semi-infinite cylinder consistent with the effort to simulate a 2-D problem. The reconstructed images are of the 2-D k^ distributions, where k is the complex wavenumber such that k^ = (i)^ne +(mis the radian frequency and is the magnetic permeability). In this way, the real part of the k^ image conveys information about the relative dielectric constant, e^, of the material and the imaginary part shows information about the conductivity, a. Figure IV.4.2 shows images of the real and imaginary parts of the reconstructed electrical profiles for all three phantoms. The image of X3 phantom is clearly the best in terms of both the real and imaginary components. It achieves good spatial resolution of the object diameter, and the k^ values reconstructed for the phantom agree quite well with that reported in Guy [1971]. Its diameter is roughly one-half wavelength which may be considered the lower limit of resolution for such an imaging system. The other two phantoms are substantially smaller than a half wavelength and their reconstructed images clearly degrade with decreasing size. Figure IV.4.3 shows plots of the electrical profile along longitudinal transects through the target region for the real and imaginary parts of k^, respectively. Comparisons are performed between the actual and reconstructed values for all three phantoms. Here, the real 89 Re(k2) (m-2) 2000 3000 4000 5000 6000 7000 8000 0000 Im(k2) (m"2) 0 1000 2000 3000 4000 5000 6000 7000 8000 9000 10000 XI Phantom (0.016 m) 0.148 m X2 Phantom (0.025 m) X3 Phantom (0.042 m) Figure IV.4.2 Real and imaginary parts of the reconstructed images of the k^ distribution for phantoms XI, X2, and X3. 90 RE(k ^) (m ^) IM(k ^) 10000 10000 7500 7500 5000 5000 2500 2500 (m ^) XI Phantom (1.55 cm) 10000 10000 7500 7500 5000 5000 2500 2500 X2 Phantom (2.51 cm) 10000 10000 7500 7500 5000 5000 2500 2500 0 -0.08 -0.04 0.04 X - position (m) 0.08 ?0.08 ?0.04 0.04 X - position (m) X3 Phantom (4.23 cm) Figure IV.4.3 Electrical property profile along longitudinal, transects through the target region for the real and imaginary parts of k , respectively. Comparisons are performed between the actual and reconsti uci.od values for all three phantoms. 0.08 91 part of the image for phantom X2 shows the correct tendencies, i.e. an object in the center is identifiable and of smaller diameter than the X3 phantom. However, the reconstruction in this case is less quantitative with respect to the calculated value of the phantom electrical properties. The corresponding imaginary part shows a similar trend with some significant fluctuations in the center of the phantom region. The images for the XI phantom are somewhat worse than those for the X2 phantom, especially in the real part where an erroneous local peak is reconstructed. It should also be noted that the values of k^ match quite well to the surrounding material at the outer edges of the target region in all cases. As a final quantitative measure of the quality of the images of targets of decreasing size shown in Figure IV.4.2, we provide Table IV.4.2 where we report the averaged difference in reconstructed k^ values normalized by their exact value for both the target and the background regions. The sampling of the reconstructed k^ distribution for purposes of this comparison RMS error of Re(k2) RMS error of Im(k2) Phantom background target background target XI 0.1133 10.24 0.1060 119.3 X2 0.09832 5.706 0.1228 100.1 X3 0.06925 2.932 0.1533 57.38 Table IV.4.2. Image Reconstruction Errors was performed on the finely discretized mesh of Figure IV.4.1(a). A clear quantitative improvement (factor of 2-3) is seen in the reconstructed images of the real part of the k^ distribution for both the target and background regions as the object size increases from less than one-fifth to near one-half wavelength. The same is true for the imaginary portion of the k2 distribution in the target region, although the background does not exhibit this trend. It should be noted that while the data reported for the background data in Table IV.4.2 can be viewed as a percent error (i.e. the background has been reconstructed quantitatively to 92 approximately 10% accuracy), the values for the target region are skewed by the normalization to the very small (relative to the background) exact value in the target region and the fact that smoothing of the jump discontinuity in the exact profile by the imaging algorithm (as illustrated in Figure IV.4.3) has occurred. The image reconstruction algorithm used to produce the data in Figures IV.4.2 and IV.4.3 is essentially the same as that reported elsewhere, [Paulsen, et al., (a)], where a dual mesh scheme is used to minimize the number of reconstruction parameters (i.e. Figure IV.4.1(b) shows the reconstruction parameter mesh whereas Figure IV.4.1(a) shows the electric field mesh.), but with several recent enhancements. A weighted least-squares minimization [Kallman, et al., 1992] involving pre- and post-multiplying of the Jacobian matrix with diagonal weights has been used to improve matrix conditioning which is particularly important in this case since the measured electric field strengths at different spatial positions can differ by several orders of magnitude. In addition, low pass filtering has been applied during the iterative reconstruction process by spatially averaging the updated reconstruction parameters at each iteration, which reduces the effects of noise and measurement imprecision (Appendix H). This prevents the algorithm from quickly seeking local minima associated with physically unrealistic electrical property distributions. Reconstructions were performed on the IBM RISC 6000 Model 590 work station and took 5.4 minutes per iteration and in all cases required five or fewer iterations to achieve adequate results. The stopping criterion used was when the RMS electric field difference between computed and observed values for each iteration reached a steady-state behavior. IV.5 Discussion and Conclusions A prototype microwave imaging system has been designed and constructed to produce multiple TM-type electric field illuminations of a phantom and subsequent measurements of the resultant field distribution around the periphery of a larger target region containing the object 93 interest. The present 4 channel device (4 transmit, 4 receive) is a broadband (300-1100 MHz) yet high quality data acquisition system for both amplitude and phase, which is based on heterodyning principles using an intermediate frequency in the KHz range. At this stage of development the system has been found to have a 115 dB dynamic range with suppression of carrier and sideband (lower) signals of better than -25 dBc. The construction of the imaging system has been modularized for compact and efficient testing and development of the hardware in an experimental setting involving a target region containing biologically-based phantom samples of varying diameter. The microwave circuitry has been confined to a single template which is attached closely to a support frame from which the transmit and receive antennas are suspended. Precision positioning of the transmitters and receivers (especially critical, being essentially point measurements) has been achieved by predrilling fixed locations in the support frame. At the present time signal transmission and reception is under computer control; however, repositioning of the sources and receivers remains a manual operation. Careful calibration of the device has been accomplished in two respects. First, phase shifts inherent in each transmit and receive channel have been eliminated by empirically determining the offset phase values needed to achieve a zero phase reference at the radiating antenna and the true phase at the tip of each receiver. In addition, loss factors associated with 3-D and 2-D radiation have been determined in the experimental system and used to calibrate the measured data prior to its input to the image reconstruction algorithm which is presently 2D. Data collected by the measurement system and compared with the numerical model of the physical set-up showed agreement between phases to within 3 degrees and amplitudes to within 0.3 dB when the electrical properties were known which is below estimates of the expected measurement errors in these quantities from our system. These results also demonstrated that the numerical model implemented (Chapters II and III and Meaney, et al., 94 [1995], Paulsen et al., (a)) provides a sound framework on which to base an image reconstruction process when the electrical properties are unknown. Several initial imaging experiments have been conducted on a bone/fat type of phantom centrally placed in the target area with different sized objects (diameters ranging from 1.5 to 4 cm) having been imaged successfully using the algorithms of Chapters II and III. The largest phantom which was greater than one-half wavelength (of the background medium) at 500 MHz produced quantitative images both in terms of reconstructing the electrical properties (simultaneous recovery of resistive and reactive components) of the target area and the phantom and determining the size and location of the phantom object. As the size of the phantom object decreases below one-half wavelength, the reconstructed information maintains a qualitative character, but becomes quantitatively less useful. The size and location of the next smallest phantom (2.5 cm diameter) is accurately determined, however, the values of its electrical properties are significantly in error, although the target region electrical properties are correctly recovered. An unexplained localized perturbation in the reconstructed profile of the imaginary component of the phantom electrical properties was observed towards one edge of the object in this case. For the smallest phantom, the location of the heterogeneity is correctly determined and general trends in its size and electrical composition relative to the surrounding target area medium are recoverable. However, more severe local anomalies are present in the reconstructed images, especially in the real component where an aberrant peak is observed in this case. Overall the images produced are encouraging, particularly in light of the fact that if the target meets intuitive size-related restrictions, quantitative and absolute images of both the resistive and reactive components of the electrical property distribution were reconstructed. The data suggests that somewhat higher frequencies may be desirable to increase spatial resolution. The present system has a 115 dB dynamic range which provides another 55 dB of capability 95 over the conditions presented in these experiments, hence, increases in frequency should be achievable. Having successfully completed these initial experiments, a variety of more complex situations are investigated in Chapter V. There, in the static imaging mode, off-axis and multiple targets are studied over a broad frequency range to investigate what improvements are achievable at higher frequencies. Eventually, data acquisition must be completely automated within the system before dynamic imaging experiments can begin. Here, system design trade� offs between adding more channels attached to fixed position radiators and receivers must be considered relative to using a small number of channels attached to automatically-moveable transmit/receive elements. The design of this next generation system is continuing on the basis of what has been accomplished with the results of these static experiments and those reported in Chapter V. CHAPTER V. SYSTEM EVALUATION IN TISSUE-EQUIVALENT PHANTOMS V.l. Introduction As indicated in Chapter I, microwave imaging applied to various medical problems can offer several benefits relative to traditional imaging modalities [Larsen, et al., 1986]. In particular, with regard to soft tissue, it provides a material property contrast of roughly 20:1 as compared to the contrast levels of only a few percent which are available with other modalities [Larsen, et al., 1986; Gregg, 1977]. Being able to image such contrast levels may offer more detailed information about the functionality of various tissues. Of particular interest is the area of hyperthermia, where these electrical properties can vary with tissue temperature, most notable being changes in conductivity with thermal variations [Schwan, et al., 1953; Bottomley, et al., 1978; Burdette, et al., 1986]. Several investigators have sought to exploit microwave imaging in this context [Jofre, et al., 1990; Miyakawa, 1993] and while fundamental groundwork has been laid by these studies, to date the temperature imaging experiments that have been reported have been crude [Jofre, et al., 1990; Miyakawa, 1993; Mallorqui,et al., 1992]. Much of the difficulty has resided in the fact that the static images produced by the integrated functioning of the hardware and software components of these systems have not been accurate enough to expect that quantitative temperature changes could be recoverable from the subtle changes in electrical properties that must be sensed. Based on the strong rationale for medical microwave imaging, and the overall promising, but qualitatively successful, preliminary work of others, we have been developing an imaging system with the goal of temperature estimation. As such, one of our primary foci has been the realization of an approach with quantitative imaging capabilities. Towards this end we have attempted to exploit certain intrinsic advantages in both our hardware and image reconstruction software. In this regard, a hardware system has been constructed, which is based on coherent measurement of field amplitude and phase using heterodyning principles, that operates over a 96 97 broad frequency range with high fidelity. This system and its ability to simultaneously reconstruct quantitative values of the resistive and reactive components of the electrical properties for simple single-target symmetric phantom configurations has been demonstrated in Chapter IV and [Meaney, et al., (a)]. As discussed in Chapters II and III, the key to this success has been the development of an iterative reconstruction algorithm which utilizes the hybrid element method for accurate field computation and a dual mesh scheme for reducing the problem size and amount of measurement data required. In the following sections, a thorough evaluation of this complete system [Meaney, et al., (a); Meaney, et al., 1995; Paulsen, et al., (a)] is performed through a series of more complex, multi-target, non-symmetric phantom experiments which further substantiate the merits of the overall approach. Several different frequencies and target contrast levels have been considered in these studies and the results reported herein confirm the expected advantages of higher frequency excitation relative to the theoretical studies of Chapters II and III and Meaney, et al. [1995] and Paulsen, et al. [(a)]. To achieve quantitative images under the more complex target conditions, a novel two-tiered receiver pattern has been identified whose implementation has required modifications to the original reconstruction algorithm as described herein. These findings suggest that the quality of the static target images obtained is sufficient to provide a framework with which to proceed towards the more difficult problem of dynamic imaging of temperature distributions with the next generation system. The following sections describe the new receiver antenna arrangement in detail, and derive the changes required in our reconstruction algorithm to allow for measurement sites outside the target region. The overall system performance is demonstrated by imaging various phantom configurations of different sized objects with several contrast ratios. Since this system has a useful frequency range of 300 to 1100 MHz, reconstructed images are also shown for similar phantoms at different frequencies. These studies allow for a detailed 98 characterization of system capability and offer insights into design enhancements appropriate for the next generation version. V.2. Receiver Antenna Configuration The reconstruction algorithm employed here is a total field formulation which differs significantly from other approaches [Guo, et al., 1989; Joachimowicz, et al., 1991; Caorsi, et al., 1991] where the total field is decomposed into its constituitive parts - the incident field and the scattered field - such that the scattered field is used in an integral formulation to recover the material property distribution. In addition to mathematical nuances which distinguish these two formulation strategies, measurement issues can also arise depending on the approach that is taken. For example, the measurement sites which are located on the near side of the target region with respect to the EM source in the highly conductive medium used in this work, the total field values are due almost entirely to the source incident fields with very little contribution resulting from the scattered fields of the objects in the target region. This feature is illustrated in Figure V.2.1 which compares computed electric fields (both magnitude and phase, respectively) at sites on a 14.7 cm. diameter target region perimeter (Figure V.2.2) for a homogeneous saline medium, and that same medium with a 4.28 cm bone/fat equivalent cylinder at its center irradiated by an electromagnetic source located at an angular position of 180�. It is clear from both the magnitude and phase plots that the maximum differences between these two sets of data occur between -90" and +90�. Within our reconstruction algorithm which uses a total field formulation, the measured data for points nearest the electromagnetic source produces almost no information about potential objects in the target region. In a scattered field formulation, subtraction of the incident field from the total field would yield relatively small scattered fields that might prove useful. However, in the presence of measurement and roundoff errors due to subtracting two large, nearly identical values, effective utilization of this data in the scattered field formulation could prove difficult. 99 Measured Field Data Comparison iVIagnitude vs. Angular Position m jj Homogeneous (dB) Fat/Bone Phantom (dB) o> TJ C O) (Q ?180 -90 0 90 180 Angular Position (degrees) (A) Measured Field Data Comparison Phase vs. Angular Position 300 % 200 - 2 ?D 0) W (0 JZ 0. Homogeneous (deg) Fat/Bone Phantom (deg) 100 - -100 - -200 -180 -90 0 90 180 Angular Position (degrees) (B) Figure V.2.1. Numerically simulated electric field values at measurement sites on the perimeter of the target region for 500 MHz plotted as a function of angular position under the following conditions: (A) a homogeneous saline medium; and (B) the same saline medium perturbed with a 4.3 cm diameter fat/bone equivalent cylinder at the center of the target region. 100 Measurement Sites Electromagnetic Source Incident Waves Figure V.2.2. Diagram of elctromagnetic waves as they are perturbed after passing through and around arbitrary objects in the target region with a conventional tomography measurement configuration. Based on these observations, we have chosen not to use the near-source portion of our potential measurement data set. Recognizing that this significantly reduces the amount of measured data available for image reconstruction, which can directly effect the number of reconstruction parameters that can be accurately determined, we have sought to compensate for the decrease by increasing the amount of measurement data from the far side (transmission side) of the target region. Specifically, we have utilized a combination of measurement sites on the target region perimeter along with a set of sites which are outside of this perimeter. Figure V.2.3 conceptually illustrates our receiver antenna positioning strategy relative to the more conventional tomographic measurement site configuration from which the results in Figure V.2.2 were obtained. Utilization of the two-tiered set of sites shown in Figure V.2.3 requires 101 some modification to the hybrid element reconstruction algorithm [Meaney, et al., 1995] which is described in the next section. Electromagnetic Source Incident Waves Measurement Sites Figure V.2.3. Diagram of elctromagnetic waves as they are perturbed after passing through objects in the target region with the two-tiered measurement configuration. V.3. Measurement Sites Outside the Target Region The hybrid element reconstruction algorithm utilizes measurement sites which coincide with nodes in the finite element mesh. Thus, when the electric field distribution is calculated at each iteration, the computed electric field values at the measurement sites are available without any additional effort. In principle the finite element discretization could be extended beyond the perimeter of the target region in order to encompass the new measurement sites which are exterior to the target region. Unfortunately, this would significantly increase the computational costs due to the unnecessary discretizion of that exterior region which is not only homogeneous 102 but also of known electrical properties. An alternative approach, is to allow the new measurement sites to reside in the boundary element region of the problem (i.e. the homogeneous region exterior to the finite element mesh) which can be accomplished without sacrificing computational accuracy or increasing the size of the finite element mesh and requires only minor additional computations at each iteration. The key issues in incorporating such measurements are involved in calculation of the electric fields at these points during each iteration and calculation of the derivatives required for construction of the Jacobian matrix. It is also imperative that these procedures not place an undue computational load on the overall algorithm. Calculation of the electric fields at these "exterior" points is governed by solution of the matrix system [Paulsen, et al., 1992b]: (V.3.1) where {E} and {F} are the discretized electric fields, E, and (nxH), in the exterior or boundary element region (external to the target area. See Figure V.3.1); and the subscripts b and s refer to those particular points which are located either along the boundary of the finite element region or the source, respectively. The electric fields at all the points in vector X = {x|,x2,...xn}^ (for example at points 10-18 in Figure V.2.3) can be calculated as an integral of boundary electric fields and fluxes by [Paulsen, et al., 1989]: where the subscript z refers to the fact that the electric field is only oriented in the z-direction due to the assumed TM illumination; G is the Green's function for the 2-D Helmholtz equation with singularity at Xj; and the subscript i refers to the i^^ member of X which is located in the exterior region along the perimeter of the boundary element region. Integral equation (V.3.2) can be written in matrix form as: (V.3.3) 103 Finite Element Region X Outside Measurement Sites n Electromagnetic Source / / Boundary Element Region ? ^ *3 Figure V.3.i. Representation of the hybrid element based reconstruction algorithm with the incorporation of measurement sites outside of the finite element region. where {E^ } is a vector of electric field values calculated at the points in X. It is important to note that the matrices [Q,] and [Q2] are solely functions of the material properties in the boundary element region; hence, their values do not change with iteration as the electrical properties within the interior finite element region are updated; hence, they can be calculated and stored by a pre-processing program. In the existing algorithm, {E,} is initialized to a set of values appropriate for the type of radiator used which have been determined empirically. {Ey} is also calculated at every iteration in the present algorithm; hence, the only quantities needed to compute {E*} are {F|,} and {Fg}. Equation (V.3.1) can be rewritten as: (V.3.4) 104 where [G] = [D]~'[C] which is already calculated as part of the forward solution as described in [Meaney, et al., 1995]. Thus, {Fy} and {F,} can be computed at nominal extra cost from a simple matrix-vector multiplication of a pre-existing matrix and a vector. {E*} can also be calculated from two matrix-vector multiplications as shown in (V.3.3). The calculation of the derivatives of the field solution which are required to construct the Jacobian matrix, [J], are handled in very much the same manner. differentiation of equation (V.3.3) with respect to the ak? =K faEwl Qis] ak/ For example, value at reconstruction node, j, yields; [<32. 0 raFb" ak/ Q2.] aF. (V.3.5) akfj Again it is important to note that the matrices [Q,] and [Q^] are not functions of the k^ values within the finite element region, and are, therefore, treated as multiplicative constants during this differentiation. Since {E,} has been dictated earlier, it, too, does not change with variations in jk^j. Similarly to the forward solution, iteration and it follows that j^rr ^akf is already computed at each can be found by differentiating (V.3.4) and recognizing that matrix, [G] is also not a function of jk^j which means that these quantities are obtained by multiplying [G] by known vector ?! [? Given that [akjj i a k j ji are precisely the terms needed in the Jacobian matrix which are obtained directly from (V.3.5) once aF ak? as ak? and are known, it is clear that both calculation of the electric fields and their derivatives can be obtained during image reconstruction by existing matrix/vector multiplies which pose no major computation load on the overall algorithm. 1 05 V.4. Results This section presents reconstructed images obtained with our imaging system operating at three different frequencies (300, 500, and 700 MHz) using the receive antenna configuration described in Section V.2. The meshes used for the forward solution calculation and the image reconstruction consisted of 2342 and 143 nodes, respectively, shown in Figure V.4.1. Eight transmitter locations, which were positioned equidistantly from the target region center in angular increments of 45 degrees, and 32 receiver locations formed by two rings of equi� angular positions, fixed at radii of 7.4 cm and 10.2 cm from the target region center, comprised the experimental system configuration. However, for each transmitter excitation, only 18 receiver positions have been used involving the two semi-circular arcs of the nine measurement sites directly opposite each transmitter. Since the present system only has four transmitter and four receiver channels, the antennas have to be manually moved to prescribed locations in order to obtain a complete data set for each phantom experiment. V.4.A Materials Several materials with various electrical properties have been used during these experiments to characterize the capabilities of this system. Their properties are listed in Table V.4.A. 1. Values for the saline solution were determined as a matter of course during system calibration [Meaney, et al., (a)]. All others were measured using the HP 85070B Dielectric Probe Kit in conjunction with the HP 8753C Network Analyzer. The agar and distilled (Dl) water phantoms were held in cylindrical tubes of Plexiglas? with inner diameter of 3.81 cm and outer diameter of 4.45 cm. The other two phantoms were solid cylinders constructed from a fat/bone equivalent material [Guy, 1971] and had diameters of 2.5 and 4.3 cm., respectively. 106 Target Region ii i (B) Figure V.4.1 Meshes used in the reconstructions, (a) is the electrical mesh used in the forward solution. It includes the finite element mesh and the boundary elements. (b) is the k mesh for reconstructing the electrical properties. Highlighted are the measurement sites. 107 Frequency (MHz) 500 300 700 Re(k2) Im(k^) Re(k^) Im(k^) Re(k2) Im(k^) (Er) (a - s/m) (Er) (a - s/m) (Gr) (a - s/m) Saline 2921.08 4163.24 8199.86 7316.14 16161.3 10669.8 (Background Medium) (73.891) (1.7576) (74.672) (1.8532) (75.088) (1.9305) Fat/Bone Cylinder 229.287 32.214 613.848 65.139 1198.84 103.354 (5.80) (0.0136) (5.59) (0.0165) (5.57) (0.0187) 2755.00 1870.80 7603.37 3252.63 14840.2 4816.21 (69.69) (0.7898) (69.24) (0.8239) (68.95) (0.8714) 2988.24 46.663 8298.48 224.237 16258.6 591.387 (75.59) (0.0197) (75.57) (0.0568) (75.54) (0.1070) 96.4587 19.9445 255.862 34.0304 514.403 51.4562 (2.44) (.00842) (2.33) (.00862) (2.39) (.00931) Material Agar (.3% NaCl) DI Water Plexiglas? Table V.4.A.I. Electrical properties of the tissue equivalent materials used in the reconstruction experiments V.4.B Phantom Tests V.4.B.i Filter Parameter Determination: As discussed in [Meaney, et al., (a)], the reconstruction algorithm uses a low pass filter to smooth the reconstruction parameter updates at each iteration. It acts to average the value of a given node with values of the surrounding nodes in a weighted manner such that the influence of the surrounding nodal values can be systematically controlled. This filtering can be represented by the following formula: N kj^new (x, y) = (1 - 8)kLoid y) + ^ (%. Y) (V.4.B.i.l) j=l where the summation is over the values of the N nodes surrounding node i. As 6 is increased, the resultant image should have a tendency to appear smoother. Figure V.4.B.i.l shows images reconstructed at 500 MHz for the target region consisting of the 4.3 cm and 2.5 cm 08 Re(k2) (m-2) Ini(k2) (m"2) 000 2000 3000 4000 5000 6000 7000 8000 9000 2000 3000 4000 5000 6000 7000 8000 9000 0000 0000 0=0.25 0=0.375 0=0.5 Figure V.4.B.i.l. Real and imaginary parts of distribution of the 4.cm and 2.5 cm fat/bone equivalent cylinders with separation distance of 0.4 cm. for varying filter coefficients (0). 109 fat/bone equivalent cylinders separated by 0.4 cm for three values of 0. (This was chosen as a benchmark case because it was one of the more difficult for the reconstruction algorithm.) Proceeding from 0=0.5 to 0=0.25, it is clear that the background becomes less smooth but the smaller object on the right becomes more visible. As 0 approaches 0.25, however, the images of the target objects also appear to be significantly degrading. A value of 0=0.375 appears to be fairly optimal in that the background is quite uniform, the larger object is quite visible and uniform, and smaller one is beginning to appear. This value has been used for all subsequent reconstructions. V.4.B.a. Use of two-tiered Receiver Sites: It is informative to compare reconstructions involving the use of only nine receiver sites from the inner arc with reconstructions based on 18 receiver sites from both arcs. These two measurement configurations yield totals of 72 and 144 pieces of data (taken at 500 MHz), respectively. This comparison not only yields some indication as to how useful the data from the outer arc really are, but also provides information on how effective the algorithm can be in the under-determined least squares problem since it will be reconstructing 143 separate parameters in both cases. Figure V.4.B.ii.l show the images of the 4.3 cm and 2.5 cm fat/bone equivalent cylinders separated by increasing distances for reconstructions using 72 and 144 pieces of measurement data, respectively. Longitudinal transects across these images of the k^ values are shown in Figure V.4.B.ii.2 and compared with the actual material distributions. In general, both reconstructions recover the larger object in the proper location, while the images using the 144 data points capture its material property values much more accurately. This may partly be due to the raw increase in data ; however, while the results of both sets of images are good, there is a clear benefit from the extra data coming from the two-tiered configuration. 1 10 Re(k2) (m-2) Im(k2) (m"2) 2000 3000 4000 5000 6000 7000 8000 9000 2000 3000 4000 5000 6000 7000 8000 9000 10000 10000 0.4 cm spacing 2000 3000 4000 5000 6000 7000 8000 9000 10000 2000 3000 5000 6000 7000 8000 9000 0000 1.7 cm spacmg 0 1000 2000 13000 4000 5000 6000 7000 8000 9000 10000 2000 4000 6000 7000 9000 10000 4.2 cm spacmg Figure V.4.B.ii.I(A). Real and imaginary parts of the distribution of the reconstructed images for the 4.cm and 2.5 cm fat/bone equivalent cylinders for varying separation distances at 500 MHz using only the inner bank of measurement sites. 1 1 Im(k2) (m"2) Re(k2) (m"2) 2000 3000 4000 5000 6000 7000 8000 0000 2000 3000 4000 5000 6000 7000 8000 9000 10000 4000 5000 8000 10000 0.4 cm spacmg 4000 7000 8000 9000 10000 1.7 cm spacmg 0 1000 12000 3000 4000 5000 6000 7000 8000 9000 10000 2000 3000 5000 6000 7000 8000 9000 10000 4.2 cm spacmg Figure V.4.B.ii.l(B). Real and imaginary parts of the distribution of the reconstructed images for the 4.cm and 2.5 cm fat/bone equivalent cylinders for varying separation distances at 500 MHz using both the inner and outer banks of measurement sites. 112 RE(k ^ ) (m ^) IM(k ^) 10000 10000 7500 7500 5000 5000 2500 2500 (m ^) 0.4 cm spacing 10000 10000 7500 7500 5000 5000 2500 2500 1.7 cm spacing 10000 10000 7500 7500 5000 5000 2500 2500 0 0 -0.08 ?0,04 X - position 0,04 0,08 ?0,08 (iti) ?0,04 0,04 X - position (m) -1 . 2 c m s p a c i n g Figure V.4.B.ii.2. Electrical property transects through the taiget parts of , respectively. profile along longitudinal, region for the real and imaginary Ilots correspond to reconstructed images shown in Figure V.4.P.ii.1(A) and (B), respectively. 0,08 1 13 VA.B.iii. Resolution with Frequency: While a variety of image reconstructions have been performed at selected frequencies, only those involving multiple objects will be shown explicitly. These generally pose a more difficult reconstruction problem and cases have been carefully chosen to illustrate the algorithm's ability to resolve object size, separation of objects and material contrast. The first test case to be demonstrated involves the 4.3 and 2.5 cm fat/bone equivalent cylinders separated by 0.4 and 1.7 cm, respectively. Figures V.4.B.iii.l.A and B show images for these cases obtained at 300, 500 and 700 MHz with longitudinal transects across the electrical property distributions displayed in Figures V.4.B.iii.2.A and B. For the 1.7 cm separation, the images distinguish the two objects and their sizes quite well, with the real part of the 700 MHz images showing the most uniformity in the background medium and the sharpest resolution at material interfaces when compared with the real part of the 300 MHz images. It is unclear whether there is any improvement with frequency in the imaginary part of the images. The images for the smaller separation (0.4 cm.) are quite informative. All three contain what appears to be an extension to the right of an otherwise normal large cylinder. Again, the 300 MHz case shows the greatest degree of blending with the background with the 700 MHz case exhibiting a much sharper contrast. Again, it appears that the real part of the background becomes more uniform with frequency, but no improvement in the imaginary part is discemable. It also appears that since there is less blending of the materials; the 700 MHz case is beginning to distinguish two separate objects. This is not totally unexpected, as in other frequency dependent imaging techniques, shorter wavelengths most often yield better resolution. 14 Re(k2) (m"2) Im(k2) (m"2) 600 2400 2800 3200 3600 4000 300 MHz 0 1000 2000 3000 4000 5000 6000 7000 8000 9000 10000 2000 4000 6000 8000 1 0000 12000 14000 16000 18000 20000 1 3000 2000 4000 5000 6000 7000 8000 9000 10000 500 MHz 2800 4200 5600 8400 9800 11200 12600 14000 700 MHz Figure V.4.B.iii.l(A). Real and imaginary parts of the k^ distribution of the reconstructed images for the 4.cm and 2.5 cm fat/bone equivalent cylinders with a separation distance of 0.4 cm at 300, 500 and 700 MHz using both the inner and outer banks of measurement sites. 1 15 Im(k^) (m"2) Re(k2) (ni"2) 0 400 800 1200 1600 2000 2400 2800 3200 3600 4000 2400 3000 3600 4200 4800 5400 300 MHz 2000 3000 2000 5000 6000 7000 8000 9000 5000 7000 8000 9000 0000 10000 500 MHz eoDo.o MOOD HODO.O 2800 D bODO.O 4200 D EODO.O lODOD.a XOO.D 7000.D 12DOO.O MOOD IIDOD.Q 16D0D.a 11200.0 18DOD.O 12600.0 700 MHz Figure V.4.B.iii.l(B). Real and imaginary parts of the k^ distribution of the reconstructed images for the 4.cm and 2.5 cm fat/bone equivalent cylinders with a separation distance of 1.7 cm at 300, 500 and 700 MHz using both the inner and outer banks of measurement sites. 116 RE(k ^) (m ^) IM(k ^ ) 4000 6000 3000 4500 2000 3000 1000 1500 (m "^ ) 0 300 M H z 10000 10000 7500 7500 5000 5000 2500 2500 0 500 MHz 20000 14000 15000 10500 10000 7000 5000 3500 0 ?0.08 ?0.04 0 0.04 0.08 ?0.04 -0.08 X - position (m) 0.04 0.08 X - position (m) 700 MHz F i g u r e V . 4 . B . i i i . 2 ( A ) . Electrical property p r o f i l e a l o n g l o n g i t u d i n a l , transects parts of through , the target respectively. region Plots images shown i n Figure V .-I . B . i 11 . l (A). for the real correspond to and imaginary reconstructed : 117 RE(k 2) (m "2, IM(k ^ ) 4000 6000 3000 4500 2000 3000 1000 1500 (m ^ ) 0 300 M H z 10000 10000 7500 7500 5000 5000 2500 2500 0 500 MHz 20000 14000 15000 10500 10000 7000 5000 3500 0 -0.08 0.04 0 0.04 0.08 0.08 X - position (m) ?0.04 0.04 0.08 X - position (m) 700 MHz F i g u r e V . 4 . B.ii i . 2 (B). Electrical property profile a l o n g l o n g i t u d i n a l , transects parts of through the tarqet region for the real and imaginary , respectively. Pl.jts correspond t o reconstructed : images shown in Figure V . 4 . P.. i i i . !( B). I18 The next set of comparisons involves a dual object phantom consisting of the 4.3 cm fat/bone equivalent cylinder and the 3.8 cm diameter DI water cylinder. These are particularly interesting because the DI water and saline should demonstrate little difference in their relative dielectric constants, e^, but their conductivity values, a, should be strongly distinct. Figures V.4.B.iii.3.A and B show the images for these two objects with separations of 1.0 and 2.3 cm, respectively, at 300 and 700 MHz. Longitudinal transects of these images are shown in Figures V.4.B.iii.4.A and B. In the 700 MHz case, the DI water object is correctly imaged; the imaginary part of the image (corresponding to a) clearly shows the object, which is essentially not visible in the real part of the image (corresponding to 6^). The situation is quite different at 3(X) MHz. The background of the real part is much less uniform than that at 700 MHz and a low object emerges erroneously. It is important to note here that the Plexiglas^^ wall thickness was roughly 0.3 cm, and may have had some influence in making the composite object appear as a low dielectric cylindrical target. Nonetheless, there is a clear advantage in the 700 MHz reconstructions. The last situation to be considered here is similar to the previous one in terms of the objects, sizes and material compositions used except that the DI water is replaced with an agar gel that is 0.3% NaCl. This target results in a dielectric constant slightly lower than that of the saline background, and a smaller conductivity value, but one not as dramatically different than that of the fat/bone equivalent cylinder. Figures V.4.B.iii.5.A and B show the results of these image reconstructions at 300 and 700 MHz with longitudinal image transects displayed in Figures V.4.B.iii.6.A and B. At both frequencies the algorithm does quite well in recognizing the fat/bone equivalent cylinder, and the conductivity of the agar object. However, the results of the dielectric constant for the agar object are quite similar to the DI water tests described above. At 700 MHz, the algorithm is able to fairly accurately recover the dielectric constant of agar, but the results at 300 MHz show that it is still detecting an object of significantly lower dielectric constant. 119 Im(k2) (m-2) Re(k2) (m-2) 600 2000 2400 3000 2400 2800 3200 3600 4000 4200 4800 5400 300 MHz 0 2000 4000 6000 8000 10000 12000 14000 16000 18000 20000 2800 4200 5600 7000 9800 200 2600 4000 700 MHz Figure V.4.B.iii.3(A). Real and imaginary parts of the k2 distribution of the reconstructed images for the 4.2 cm fat/bone equivalent and the DI water cylinders with a separation distance of 1.0 cm at 300 and 700 MHz using both the inner and outer banks of measurement sites. 120 Re(k2) (m-2) Im(k^) (m-2) 2000 2400 3600 4000 300 MHz 2000 4000 6000 8000 10000 12000 4000 6000 18000 20000 4200 5600 7000 8400 9800 11200 12600 14000 700 MHz Figure V.4.B.iii.3(B). Reai and imaginary parts of the distribution of the reconstructed images for the 4.2 cm fat/bone equivalent and the DI water cylinders with a separation distance of 2.3 cm at 300 and 700 MHz using both the inner and outer banks of measurement sites. 121 RE(k 2 ) (m 2 ) IM(k; ^ ) 4000 6000 3000 4500 2000 3000 1000 1500 (m ^ ) 300 MHz 20000 14000 15000 10500 10000 7000 5000 3500 0 -0.08 -0.04 0 0.04 0.08 -0.04 -0.08 X - position (ml 0 0.04 0.08 X - position (m) 700 MHz F i g u r e V . 4 . B . i i i. 4 ( A ) . E l e c t r i c a l p r o p e r t y transects through the target region profile along longitudinal, for the real parts of , respectively. Plots correspond images shown in Figure V.4.B.iii.3(A). to and imaginary reconstructed ! 122 RE (k ^ ) (m ^ ) IM(k ^ ) 4000 6000 3000 4500 2000 3000 1000 1500 (m ^ ) 0 300 M H z 20000 14000 15000 10500 10000 7000 5000 3500 0 ?0.08 -0.04 0.04 0.08 ?0.08 X - p o s i t i o n (in) -0.04 0.04 0.08 X - position (m) 700 MHz Figure V.4.B.iii.4(B) Electrical property profile along longitudinal, transects through the target region for the real parts of , respectively. Plots correspond images shown in Figure V.4.B.iii.3(B). to and imaginary reconstructed : 123 Re(k2) (m-2) Im(k2) (m-2) 0 400 1800 1200 1600 12000 2400 2800 3200 3600 14000 3600 5400 300 MHz 2000 4000 6000 8000 10000 2000 14000 2800 4200 5600 7000 9800 16000 11200 18000 20000 12600 4000 700 MHz Figure V.4.B.iii.5(A). Real and imaginary parts of the k2 distribution of the reconstructed images for the 4.2 cm fat/bone equivalent and the agar gel (0.3% NaCl) cylinders with a separation distance of 1.0 cm at 300 and 700 MHz using both the inner and outer banks of measurement sites. 1 24 Re(k2) (m-2) Im(k2) (m"2) 200 200 1600 2000 2400 2800 3200 3600 4000 2000 4000 6000 8000 10000 12000 4000 2400 3000 3600 4200 4600 6000 300 MHz 2800 4200 5600 8400 9800 16000 11200 18000 12600 20000 4000 700 MHz Figure V.4.B.iii.5(B). Real and imaginary parts of the k^ distribution of the reconstructed images for the 4.2 cm fat/bone equivalent and the agar gel (0.3% NaCl) cylinders with a separation distance of 2.3 cm at 300 and 700 MHz using both the inner and outer banks of measurement sites. 125 RE(k 2) (m ^ ) IM(k ^ ) 4000 6000 3000 4500 2000 3000 1000 1500 300 MHz 20000 14000 15000 10500 10000 7000 5000 3500 0 ?0.08 ?0.04 0 0,04 0.08 ?0.08 X - p o s i t i o n (rn) ?0.04 0 0.04 0.08 X - position (m) 7 0 0 MHz F i g u r e V . 4 . B . i ii. 6 ( A ) . Electrical p r o p e r t y p r o f i l e a l o n g l o n g i t u d i n a l , transects through the target region for the parts of , respectively. Plots correspond images shown in Figure V.4.B.iii.5(A). real to and imaginary reconstructed : 126 RE(k ^) (m ^) IM(k ^) 4000 6000 3000 4500 2000 3000 1000 1500 (m "^) 0 300 MHz 20000 14000 15000 10500 10000 7000 5000 3500 -0.08 ?0.04 0.04 0.08 ?0.08 X - position (m) ?0.04 0.04 0.08 X - position (m) 700 MHz F i g u r e V . 4 . B . i i i . 6( B ) . E l e c t r i c a l p r o p e r t y transects through the target region profile along longitudinal, for the parts of , respectively. Plots correspond images shown in Figure V.4.B.iii.5(B). real to and imaginary reconstructed : 127 VA.B.iv: Single Object Image Reconstructions: For completeness, at the end of this chapter are shown a number of reconstructed images of single objects. In all cases the targets consisted of varying diameter cylinders of the fat/bone equivalent material. Imaging was again performed at 300, 500 and 700 MHz for cylinders of diameters of 1.6, 2.5 and 4.3 cm. Tests were conducted with these objects placed at various positions in the target region to determine whether the position influenced the system's ability to reconstruct the objects. Longitudinal transects through the reconstructed images are also shown to indicate the quality of these images. For all the images shown in Figures V.4.B.iv.l thru V.4.B.iv.l4, similar trends to those demonstrated in the multi-target cases were observed, namely, that the quality of the image is better for the larger targets and that the real part of the image improves as the operating frequency is increased. V.4.B. V : Reconstructed Image Error Analysis: Along with viewing the reconstructed images and transects through the images as ways of assessing system performance, we have also conducted an error analysis to further quantify these results. For this examination we have chosen to analyze three quantities and their variation with frequency, object size and differences in the real and imaginary parts of the reconstructed property profiles. These quantities include: (1) RMS error of the values both in the object and the background subregions, (2) diameters of the reconstructed objects, and (3) differences in the location of the reconstructed and actual objects. A brief discussion of how these quantities were determined is necessary before the analysis is performed: (1) RMS errors - The first step in this procedure is to map the reconstructed image parameters of the coarse mesh onto the fine mesh. This will make the analysis less susceptible to problems arising from insufficient discretization in the target areas with steep material property gradients. Then the fine mesh elements are chosen for computing the background and object RMS errors separately. In the case of the background, this is achieved by selecting all fine mesh elements which are outside of the actual object area (with an extra buffer of 1 mm to eliminate elements in the steepest gradient area of the material interface). 128 In the case of the object itself, selection of all fine mesh elements which are inside of the actual object area, again with a 1 mm buffer, is made. The RMS error values for the real parts of the images are then calculated using the following formulas: BackgroundRc[kactuai (C(i))] RMS Error =. i y Rc[kfeconstn)cte<i (C(i))] (V.4.B. V .1) Re[kLa,(C(i))] Object- 14 RMS Error =. Rc[kactual(C(j))] RG[kreconstructed(('(j))] (V.4.B. V .2) Re[kLa,(C(j))] where N is the number of elements outside of the object and M is the number of elements inside, and C(^) refers to the center of element �. A similar procedure is performed for the imaginary components of the images by simple replacing the function Re() with Im(). (2) Diameter of reconstructed object - Because of the large contrast ratios between the fat/bone equivalent objects and the saline background medium, a full-width, half-height criterion is employed to determine the object diameter of the recovered object. Similarly to the calculation of the RMS error, the k^ values are mapped from the coarse mesh to the fine mesh. If there are multiple objects in the target region, a subregion is specified which fully encompasses only the object under consideration. Then a midpoint k^ value is determined by taking the average of the lowest k^ value in the selected subregion (In all of these cases the values of the real and imaginary parts of the objects are significantly less than that of the background medium) with the predetermined value of the background medium. Based on this cutoff, all the elements whose k^ values at their centers are less than the k^ midpoint value are selected and their individual areas are summed to produce the reconstructed object 129 size. Assuming a circular object shape, an effective diameter can be calculated from the following equation: (V.4.B. V .3) where Dgff is the effective diameter and A is the total area of the subregion selected. (3) Location of the reconstructed object - Similarly to the procedure for calculating the object diameters, the values are mapped from the coarse mesh to the fine mesh and a subregion of the image is selected in the case of multiple targets. Then the center of the object is determined in a manner similar to calculating the center of mass of an object. However, in this case there is the added complication that the values for the object are all less than those of the background. To overcome this, the center of the object is calculated for the inverse of the k2 distribution which also tends to reduce the effects of ripple in the background region which could skew the center position estimate. As an example, the center of the real part of an object in the x-direction is calculated using the following formula: ^Re(k^(C(i)))' XRE (center) = ?? (V.4.B. V .4) |j'Re(k^(C(i))) where (center) is the x-coordinate of the real part of the object, and similarly to earlier conventions, k^(C(i)) is the k^ value at the center of element /, A(i) is the area of element i and N is the total number of elements in the selected region. With minor alterations, this formula is used to calculate X,M (center), Y^g(center) and (center). In all of these tests, only the 4.28 and 2.54 cm diameter fat/bone equivalent phantoms have been analyzed to simplify the comparisons. Similar trends were noted with phantoms of varying contrast levels. Table V.4.B.v.l shows a summaiy of RMS errors for the single object experiments versus frequency for both the background and object regions. It is clear 130 4.28 cm diam. Object RMS ERRORS Frequency (MHz) REAL 300 500 700 IMAGINARY 300 500 700 Samples per wavelength 28.5 18.4 13.5 Samples per exponetid decay 8.74 7.70 7.18 Background Object .153 .095 .085 6.34 4.81 3.90 .132 .135 .121 52.8 57.8 69.7 2.54 cm diam Object RMS ERRORS Frequency (MHz) REAL 300 500 700 IMAGINARY 300 500 700 Samples per wavelength 28.5 18.4 13.5 Samples per exponetid decay 8.74 7.70 7.18 Background Object .133 .083 .062 9.49 5.38 3.07 .099 .103 .091 94.0 91.9 90.0 Table V.4.B.v.l. Average RMS enors for the Real and Imaginary parts of the reconstructed material electrical properties for tests performed at 300, 500 and 700 MHz. Error values were averaged for the 4.28 cm and 2.54 cm diameter fat/bone equivalent cylinders at three distances (0.0, 2.54, and 3.81 cm) from the center of the target region along the horizontal axis. that there is an improvement in the errors with frequency for the real part of the images for both sizes of objects in both the background and object regions. It is not clear, however, what trends are emerging for the imaginary part of the images, and in the case of the object region of the 4.28 cm diameter target, it appears that the error actually increases with frequency. Similar results were observed in the multi-target cases as shown in Table V.4.B.V.2. Again the real part of the images (background and object) show improvement with frequency. In these cases. 1 3I Left Object RMS ERRORS Frequency (MHz) REAL 300 500 700 IMAGINARY 300 500 700 Samples per wavelen^ 28.5 18.4 13.5 Samples per exponetied decay 8.74 7.70 7.18 Background Object .242 .163 .154 6.32 4.38 3.72 .201 .261 .275 49.1 59.0 63.2 Right Object RMS ERRORS Frequency (MHz) REAL 300 500 700 IMAGINARY 300 500 700 Samples per wavelength 28.5 18.4 13.5 Samples per exponetid decay 8.74 7.70 7.18 Background Object .242 .163 .154 11.1 3.91 3.99 .201 .261 .275 88.5 108.5 120.4 Table V.4.B. V .2. Average RMS errors for two cylindrical objects made out of fat/bone equivalent material with the actual diameters being 4.28 cm and 2.54 cm for the left and right objects, respectively. Averages were compiled for three test cases where the two objects were separated by 0.4, 1.7, and 4.2 cm, respectively. the imaginary part of the images show increased degradation with frequency. Some insight into this phenomenon can be achieved by observing the mesh sample rates of the fine mesh as shown in Table V.4.B.V.1 and 2. In all cases, the number of samples per wavelength is greater than 10 which is a benchmark established in Chapter II for accurately modeling the forward solution. However, as the frequency increases, the mesh sampling rate per exponential decay approaches 7 which is also a benchmark for accurate forward solution 132 modeling. Thus, because of the high conductivity of the medium, a higher level of discretization may be required for more accurate forward solution modeling in an effort to compensate for the observed increase in RMS error with frequency. Tables V.4.B.V .3 and 4 list the results of the diameter and center of object computations for the 4.28 and 2.54 cm fat/bone equivalent cylinders in both the single and multi-target experiments, respectively. In general the real part of the images more accurately captures the cylinder diameters, especially for the larger of the two objects. It does not appear 4.28 cm diam. Object Frequency (MHz) Average Diameter (cm) Average Distance from exact location (cm) 300 500 700 4.31 4.28 4.02 .127 .036 .078 300 500 700 4.09 3.63 3.16 .061 .068 .064 Average Diameter (cm) Average Distance from exact location (cm) 300 500 700 3.23 2.83 2.59 .079 .048 .035 300 500 700 3.42 2.95 2.37 .069 .040 .030 REAL IMAGINARY 2.54 cm diam. Object Frequency (MHz) REAL IMAGINARY Table V.4.B.v.3. Average diameters and position errors for the two fat/bone equivalent cylinders with the actual diameters being 4.28 cm and 2.54 cm, respectively, during singletarget imaging experiments. Values were averaged for the cylinders at three distances (0.0, 2.54, and 3.81 cm) from the center of the target region along the horizontal axis. 133 Left Object Frequency (MHz) Average Diameter (cm) Average Distance from exact location (cm) 300 500 700 4.04 4.00 3.85 0.572 0.517 0.398 300 500 700 2.96 3.14 2.94 0.642 0.648 0.721 Average Diameter (cm) Average Distance from exact location (cm) 300 500 700 2.96 3.14 2.94 0.751 0.828 0.830 300 500 700 3.27 2.81 2.38 0.857 0.973 0.973 REAL IMAGINARY Right Object Frequency (MHz) REAL IMAGINARY Table V.4.B. V .4. Average diameters and position errors for the two fat/bone equivalent cylinders with the actual diameters being 4.28 cm (left object) and 2.54 cm (right object), respectively, during multi-target imaging experiments. Averages were compiled for three test cases where the two objects were separated by 0.4, 1.7, and 4.2 cm, respectively. that the algorithm's ability to accurately capture the object diameters is significantly better in either the single or multi-target cases. However, the center of the object computations show a clear improvement when only one object is imaged at a time. In these cases, the positioning error is generally 1 mm or less for both object sizes and at all three frequencies. In contrast, the errors in the object positions of the multi-target cases were generally between 4 and 10 mm with roughly 90% of these errors occurring along the common axis of the two objects. Thus the accuracy of the object position is clearly affected by the presence of other objects. 1 34 Generally, the diameter estimates for the targets appear to be within 1 cm with typical errors being more on the order of a few millimeters. Similarly, the position estimates are correct to better than 1 cm and quite often within 5 mm or better. V.5. Discussion and Conclusion A prototype microwave imaging system has been developed which can achieve reconstructed images of the material property distributions for a wide range of material contrasts over a frequency range of 300 to 700 MHz. Particular attention has been paid to the amount of data necessary to reconstruct these images along with selecting that data which is most beneficial in such reconstructions. In the present configuration, the system was able to achieve relatively good results even for cases where the Jacobian was under-determined, with incremental improvement being achieved through the use of a two-tiered receiver system consisting of an additional arc of receiver sites outside the target region. Comparison of results at different frequencies demonstrated the capability of reconstructing relatively good images over the full 300-700 MHz frequency range with improvements occurring as the frequency is increased. These observations are based on the ability of the algorithm to achieve high contrast between two different materials, reconstruct a relatively uniform background medium and distinguish between two closely spaced objects. In the more challenging tests of the algorithm's capability, objects with an assortment of electrical properties were imaged. The results of these comparisons showed a clear advantage for the higher frequency excitations, where both the dielectric constants and conductivity were reconstructed quite well at 700 MHz, whereas similar tests at 300 MHz were only able to recover the conductivity of the objects. In practice, increasing frequency will have to be balanced with the increasing attenuation in the biological tissue. Since the microwave components of this system are capable of operating over a wide spectrum, 300-1100 MHz, we have some flexibility in this regard. An additional option which we have not explored to date is the combination of data from multiple frequencies in order to better determine target material properties. We also have 135 the flexibility of modifying the conductivity of the background medium. All of these will ultimately be key issues in determining the optimal size of the target region and the operating frequencies. This work also poses some interesting questions as to how much and what kind of measurement data is most useful for image reconstruction. For practical considerations the background medium has been manipulated and the reconstruction algorithm formulated in such a way as to render signal reflection data minimally useful. These imaging experiments are being conducted under conditions which represent a compromise between microwave systems that utilize reflection, refraction and transmission data and other imaging modalities such as CT which rely primarily on transmission data. These relationships will be further investigated with the implementation of the next generation system. 36 Re(k2) (m"2) Im(k2) (m'2) 1200 2000 3000 2400 2800 3600 4000 4200 4800 5400 300 MHz 2000 3000 4000 5000 6000 7000 8000 9000 10000 2000 4000 6000 8000 10000 12000 14000 16000 18000 20000 2000 3000 5000 6000 8000 9000 0000 500 MHz 400 2800 4200 8400 11200 12600 14000 700 MHz Figure V.4.B.iv. 1. Real and imagina^ parts of the distribution of the reconstructed images for the 1.6 cm diameter fat/bone equivalent cylinder centered in the target region at 300, 500 and 700 MHz using both the inner and outer banks of measurement sites. 137 RE(k 2) (m 2) IM(k ^ ) 4000 6000 3000 4500 2000 3000 1000 1500 (m ) 300 MHz 10000 10000 7500 7500 5000 5000 2500 2500 0 0 500 MHz 14000 10500 ?0.08 -0.04 0 0.04 0.08 -0.08 X - position (m) -0.04 0 0.04 X - position (m) 700 MHz F i g u r e V. 4 .B. i v . 2 . Electrical property p r o f i l e a l o n g l o n g i t u d i n a l transects parts of t h r o u g h the t.Trqet region , resper-t i'./>:? 1 y . for the real and imaginary P l o t s correspond to reconstructed images shown in Fiqnr'-.- V. 1.!. 0.08 138 Re(k2) (m"2) 1600 2000 2400 2800 3200 3600 4000 Im(k2) (m"2) 2400 3600 4200 4800 5400 300 MHz 2000 2000 3000 4000 5000 4000 5000 6000 6000 7000 8000 9000 10000 7000 8000 9000 0000 500 MHz 0 2000 4000 16000 8000 110000 2800 4200 2000 114000 16000 18000 20000 9800 11200 12600 14000 700 MHz Figure V.4.B.iv.3. Real and imaginary parts of the distribution of the reconstructed images for the 2.5 cm diameter fat/bone equivalent cylinder centered in the target region at 300, 500 and 700 MHz using both the inner and outer banks of measurement sites. 139 300 MHz 10000 10000 7500 7500 5000 5000 2500 2500 I 0 0 500 MHz 14000 10500 7000 3500 -0.08 -0.04 0 0.04 0,08 -0.08 X - position (m) -0.04 0 0.04 X - position (m) 700 MHz F i g u r e V. 4 .B. i v . 4 . Electrical property profile a l o n g l o n g i t u d i n a l , transects parts of through the target region for the real and imaginary , respectively. [lots correspond to reconstructed i m a g e s s h o w n i n F i q u i f- .-1 . P . i " . ?. 0.08 140 Re(k2) (m-2) Itn(k2) (m-2) 2400 2400 3600 2800 3600 4000 6000 300 MHz 3000 4000 5000 6000 7000 8000 9000 10000 0000 500 MHz 2000 4000 6000 8000 0000 2000 4000 6000 8000 20000 4200 8400 200 2600 14000 700 MHz Figure V.4.B.iv.5. Real and imagina^ parts of the distribution of the reconstructed images for the 4.3 cm diameter fat/bone equivalent cylinder centered in the target region at 300, 500 and 700 MHz using both the inner and outer banks of measurement sites. 141 RE(k 2 ) (m 2 IM(k ^) 4000 6000 3000 4500 2000 3000 1000 1500 (m ) 300 MHz 10000 10000 7500 7500 5000 5000 2500 2500 0 0 500 MHz 20000 14000 15000 10500 10000 5000 -0.08 -0.04 0 0.04 0.08 0.08 -0.04 0 0.04 X - position (m) X - position (ml 700 MHz F i g u r e V . 4 . B . i v . 6 . Electrical p r o p e r t y p r o f i l e a l o n g l o n g i t u d i n a l , transects through the taiqet parts of , respectively, images shown in Figuj e V . 4 .1-. region for the real and imaginary F'lots correspond to reconstructed 0.08 142 Re(k^) (m"2) Im(k2) (m'2) 800 2400 3000 3600 4200 4800 2400 2800 3200 3600 4000 6000 300 MHz 0 0 000 1000 2000 3000 4000 5000 6000 7000 8000 9000 1 0000 I 3000 4000 5000 "2000 6000 7000 8000 9000 10000 500 MHz 2000 4000 6000 8000 10000 2000 14000 0 J1400 12800 4200 5600 7000 8400 9800 16000 11200 18000 12600 20000 114000 700 MHz Figure V.4.B,iv.7. Real and imaginaiy parts of the k^ distribution of the reconstructed images for the 2.5 cm diameter fat/bone equivalent cylinder offset by 2.5 cm to the right of the target region center at 300, 500 and 700 MHz using both the inner and outer banks of measurement sites. 143 R E( k ^ ) (m ^ ) IM(k ^ ) 4000 6000 3000 4500 2000 3000 1000 1500 (m ^) 300 MHz 10000 10000 7500 7500 5000 5000 2500 2500 500 MHz 20000 14000 15000 10500 10000 7000 5000 3500 0 -0.08 -0.04 0 0.04 0.08 -0.08 X - position (m) -0.04 0 0.04 X - position (m) 700 MHz F i g u r e V .4.B. i v .8. Elefrtrical transects through property profile a l o n g l o n g i t u d i n a l , the target region parts of , respectively. i lots images shewn in Figure 1. :v . i. for the real and correspond to imaginary reconstructed 0.08 144 Re(k2) (m-2) Iin(k2) (m-2) 200 1600 2400 3000 3600 4200 4800 5400 2000 2400 2800 3200 3600 4000 6000 300 MHz 0 1000 1000 2000 3000 4000 6000 7000 8000 9000 0000 2000 4000 6000 8000 10000 12000 14000 16000 18000 20000 1 3000 - 24000 000 5000 ' 6000 7000 8000 9000 10000 500 MHz 0 ^1400 12800 4200 5600 7000 8400 19800 1200 112600 114000 700 MHz Figure V.4.B.iv.9. Real and imaginary parts of the k2 distribution of the reconstructed images for the 4.3 cm diameter fat/bone equivalent cylinder offset by 2.5 cm to the right of the target region center at 300, 500 and 700 MHz using both the inner and outer banks of measurement sites. 145 RE(k ^ ) (m ^) IM(k ^ ) 4000 6000 3000 4500 2000 3000 1000 1500 (m ) 300 MHz 10000 10000 7500 7500 5000 5000 2500 2500 500 MHz 20000 14000 15000 10500 10000 7000 5000 3500 -0.08 -0.04 0 0,04 0.08 -0.08 X - p o s i t i o n (rn) -0.04 0.04 X - position (m) 700 MHz Figure V.4.B.iv.lO. Electrical property profile along longitudinal, transects parts of through the target , respectively. region for the real and imaginary Plots correspond to reconstructed images shown in Figure V.4.B.iv.9. 0.08 146 Re(k2) (m"2) 2000 2400 3200 3600 4000 Im(k2) (m"2) 2400 3000 3600 4200 4800 5400 6000 300 MHz 2000 3000 4000 5000 4000 5000 6000 7000 8000 9000 10000 8000 9000 10000 500 MHz 0 2000 1400 4000 2800 6000 8000 4200 15600 7000 8400 9800 11200 10000 2000 14000 16000 18000 20000 12600 114000 700 MHz Figure V.4.B.iv.ll. Real and imaginary parts of the k^ distribution of the reconstructed images for the 2.5 cm diameter fat/bone equivalent cylinder offset by 3.8 cm to the right of the target region center at 300, 500 and 700 MHz using both the inner and outer banks of measurement sites. 147 RE(k ^) (m ^) IM(k ^) 4000 6000 3000 4500 2000 3000 1000 1500 (m ) 0 300 MHz 10000 10000 7500 7500 5000 5000 2500 2500 500 MHz 20000 14000 15000 10500 10000 7000 5000 3500 -0.08 ?0.04 0 0.04 0.08 -0.08 X - position (m) -0.04 0.04 X - position (m) 700 MHz F i g u r e V . 4 .B.iv. l 2 , Electrical property profile a l o n g l o n g i t u d i n a l , transects parts of through the target region for the real and imaginary , respectively. Plots correspond to reconstructed images shown In Figure , 1 . . i ?/. 1 I . 0.08 148 Re(k2) (m-2) 2000 2400 2800 3600 4000 Im(k^) (m'2) 1800 2400 3000 3600 4200 4800 6000 300 MHz 2000 3000 5000 6000 7000 8000 1 0000 3000 4000 5000 6000 7000 8000 9000 0000 500 MHz 0 2000 ; 1 400 4000 6000 |2800 10000 12000 14000 16000 4200 5600 7000 8400 9800 11200 18000 12600 20000 114000 700 MHz Figure V.4.B.iv.l3. Real and imaginary parts of the k^ distribution of the reconstructed images for the 2.5 cm diameter fat/bone equivalent cylinder offset by 3.8 cm to the right of the target region center at 300, 500 and 700 MHz using both the inner and outer banks of measurement sites. 149 IM(k ) (m 300 MHz 10000 10000 7500 7500 5000 5000 2500 2500 0 0 500 MHz 14000 10500 7000 3500 ?0.08 0.04 0 0.04 0 08 0.08 X - position (m) -0.04 0 0.04 X - position (m) 700 MHz Figure V.4.B.iv.14. Electrical property profile along longitudinal. transects region f o r t h e real and imaginary Plots correspond to reconstructed parts of through , the t arget respectively, images shown in F i g u r e Y. 1.1', i : 1 j. 0.08 CHAPTER VI. VI.1 CONCLUSION Overview The work presented in this thesis begins with a basic numerical algorithm, the Newton- Raphson iterative method, which is well suited for this class of non-linear problems. We have chosen the hybrid element method as its computational base because it yielded an accurate representation of the electric fields in a geometry which matched that needed for these imaging experiments. Issues associated with maintaining an accurate electric field representation while reducing the size of the reconstruction problem and the amount of measured data required were overcome through the introduction of a dual mesh scheme. A measurement system was developed to measure the electric fields in a saline tank for problems closely resembling a 2-D geometry. A calibration algorithm was realized to transform the measured 3-D data into 2-D for reconstruction purposes, and finally, tests were performed using a number of different phantoms to evaluate the capability of the system as a whole. In general this type of inverse problem is beset with difficulties due to its inherent nonuniqueness [Habashy et al., 1986]. The reconstruction technique developed here exhibits certain flexibilities which allow the implementation to utilize a priori information, place bounds on the intermediate electrical property distributions, apply a spatial filter on the material property distribution, and utilize various regularization procedures in an effort to stabilize the solution and limit the number of possible solutions in an effort to overcome that nonuniqueness. It also exhibits a certain parallelism, especially in the construction of the Jacobian matrix which can be exploited at a later date to further reduce the iteration time [Hua, et al., 1990]. The image reconstruction is performed over a small, defined region, but because of the unbounded nature of microwave radiation, accurate accounting of the fields radiating to infinity is required. As shown by Lynch, et al., [1986], this can be accomplished by implementing a 15 0 1 5 1 hybrid of the finite element method on the small defined region and the boundary element method for the region extending to infinity. In general, the number of reconstruction parameters is directly related to the amount of measured data that is available [Yorkey, et al., 1987; Guo, et al., 1989]. In microwave imaging, the amount of measured data simply becomes too great if the mesh used for the reconstruction parameters is the same as that used to represent the electric field distribution. The electric fields vary continuously over the entire region of interest due to their sinusoidal and exponential decay behavior. However, the material property distributions in biomedical applications are often nearly constant within various subregions, thus generally requiring far less discretization. A dual mesh scheme was introduced to decouple the mesh used for modeling the electric field distribution from that used for the material property mesh. This allowed for the fine discretization required in modeling the electric fields, and a much more coarse mesh for reconstructing the material property parameters. It also lends itself quite nicely to specialized coarse meshes that might be constructed based on a priori information. Simulations were performed comparing the dual mesh scheme with the single mesh scheme which show that reconstructed images were still quite good and achieved the desired effects of reducing the problem size and amount of measured data. A prototype microwave measurement system was developed to accurately detect the electric fields at various sites due to the radiation of water-filled waveguide radiators. It was designed around several criteria: (1) a large dynamic range given the high signal attenuation of biological tissue, (2) accurate detection of both phase and amplitude, (3) availability, cost and wideband performance of the microwave components, and (4) general modularity such that various antenna configurations could be tested without complete disassembly and reassembly. The system generally achieved these goals including a 1 IS dB dynamic range and well matched impedance characteristics over the frequency range of 300 to 1100 MHz. 152 Using measured data from the data acquisition system, a calibration procedure was implemented which transformed the measured data from the 3-D radiating pattern to the 2-D plane. Incorporation of the transformation was accomplished using a least squares technique to determine the phase center of the radiators, which was subsequently used to subtract out the requisite "free space loss factors." Agreement between the measured data modified by the calibration technique and that computed from the simulations was on the order of 0.3 dB in amplitude and 3" in phase. Such agreement is an essential component for using the NewtonRaphson method. Finally images were successfully reconstructed for a variety of targets using this system while operating over the range of 300 to 700 MHz. These showed general improvement as the frequency increased; in particular, detection of the real parts of the DI and agar gel phantoms were not accurately accomplished at 300 MHz but were recovered quite well at 700 MHz. This may be due in part to the ability of the system to detect the presence of plexiglas^^ cylinder around these objects at 700 MHz while not at 300 MHz, and hence, the exhibited smearing of the plexiglas? properties with those inside the cylinder. Otherwise, the targets were generally reconstructed quite well with some advantage going to the higher frequency in terms of background smoothness and delineation of material interfaces. VI.2 Future Directions These findings are the first reported implementation of a microwave imaging system for biomedical applications which has successfully reconstructed quantitative images of various biological equivalent tissue and thus, bodes well for the future implementation in a clinical setting. However, two major hurdles remain before clinical implementation is realized. First, dynamic images must be achieved. Such a system will either involve: (1) only few antennas 153 mounted on a motorized platform to achieve multiple views of the target region, or (2) a large number of fixed antennas which can be electronically selected. In various contexts, these both have advantages. The second hurdle is to implement a microwave imaging system without the targets and antennas being completely submerged in a water tank. To date, the results of the system described herein and that of other systems by Larsen, et al. [1984], Jofre, et al. [1990], and Miyakawa [1993] have all been in water tanks. There are several ways to simulate the desirable electrical properties of water and saline through the use of high dielectric ceramics. We plan to implement such an arrangement and continue on the lines of 2-D imaging which can be transported to a clinical setting for treatment monitoring. Future work will also need to be done on reducing the imaging time. The most significant contributions here will most likely involve utilizing the inherent parallelism of this algorithm. At present, hundreds of matrix back substitutions are performed using the same decomposed matrix. Computing these tasks in parallel will significantly reduce the computation time, but it will inevitably increase the system cost due to new levels of computational sophistication. Beyond that, the ideal goal would be to achieve 3-D images. This poses a substantially more difficult problems than dealt with here. Full characterization of the radiating fields in 3-D will be required. In addition, receiving antennas will need to detect all three spatial components of the electric fields, whereas, here we have only needed one component, assuming for practical purposes that our system is operating in the 2-D TM mode. It will also pose a much more substantial computational problem with all the associated issues of increased illconditioning. 154 At present, the envisioned 2-D imaging system would be a valuable asset in a clinical setting. In particular, using this for monitoring the temperature distributions for assessing the success of hyperthermic treatments and possible real-time feedback control will be necessary for the future success of hyperthermia. REFERENCES Balanis CA, Antenna Theory: Analysis and Design. New York, NY, Harper and Row Publishers, pp.64-67, 1982. Barber DC, Brown BH, Avis NJ, "Image reconstruction in electrical impedance tomography using filtered back-projection," Annual International Conference of the IEEE Engineering in Medicine and Biology Society, 1691-1692, Paris, 1992. Barber DC, Brown BH, "Reconstruction of impedance images using filtered back-projection," European Community Concerted Action in Electrical Impedance Tomography, pp. 1-7, Copenhagen, CAIT, 1990. Blad B, Perrson B, Lindstrom K, "Quantitative assessment of impedance tomography for temperature measurements in hyperthermia," International Journal of Hyperthermia, vol 8, pp. 33-43, 1992. Bolomey JC, Hawley MS, (Chapter 2) Non-invasive control of hyperthermia. Methods of Hyperthermia Control, Springer-Verlag, Berlin, Heidelberg, 1990. Bolomey JC, "Recent European developments in active microwave imaging for industrial, scientific, and medical applications," IEEE Transactions on Microwave Theory and Techniques, vol 37, pp. 2109-2117, 1989. Bolomey JC, Jofre L, Peronnet G, "On the possible use of microwave-active imaging for remote thermal sensing," IEEE Transactions on Microwave Theory and Techniques, vol 31, pp. 777-781, 1983. Bolomey JC, Izadnegahdar, Jofre L, Pichot CH, Peronnet G, Solaimani M, "Microwave diffraction tomography for biomedical applications," IEEE Transactions on Microwave Theory and Techniques, vol 30, pp. 1998-2000, 1982. Borup DT, Sullivan DM, Gandhi OP, "Comparison of the FFT conjugate gradient method and the finite-difference time-domain method for the 2-D absorption problem," IEEE Transactions on Microwave Theory and Techniques, vol 35, pp. 383-395, 1987. Bottomley PA, Andrew ER, "RF magnetic filed penetration, phase shift and power dissipation in biological tissue: implications for NMR imaging," Phys. Med. Biol, vol 23, pp. 630-643, 1978. Broquetas A, Romeu J, Rius JM, Elias-Fuste AR, Cardama A, Jofre L, "Cylindrical geometry: a further step in active microwave tomography," IEEE Transactions on Microwave Theory and Techniques, vol 39, pp. 836-844, 1991. Brown BH, Sinton AM, Barber DC, Leathard AD, McArdle FJ, "Simultaneous display of lung ventilation and perfusion on a real-time EIT system," Annual International Conference of the IEEE Engineering in Medicine and Biology, pp. 1710-1711, Paris, 1992. Burdette EC, Friedrich PG, Seaman RL, Larsen LB, "In situ permittivity of canine brain: regional variations and post-mortem changes," IEEE Transactions on Microwave Theory and Techniques, vol 34, pp. 38-50, 1986. 155 156 Caorsi S, Gragnani GL, Pastorino M, "A multiview microwave imaging system for twodimensional penetrable objects," IEEE Transactions on Microwave Theory and Techniques, vol 39, pp. 845-51, 1991. Caorsi S, Gragnani GL, Pastorino M, "Two dimensional microwave imaging by a numerical inverse scattering solution," IEEE Transactions on Microwave Theory and Techniques, vol 38, pp. 981-989, 1990. Chommeloux L, Pichot C, Bolomey JC, "Electromagnetic modeling for microwave imaging purposes," IEEE Transactions on Microwave Theory and Techniques, vol 34, pp. 1064-76, 1986. Clegg ST, Roemer RB, "Towards the estimation of three-dimensional temperature fields from noisy data during hyperthermia," International Journal of Hyperthermia, vol 5, pp. 918-924, 1989. Conway J, "Electrical impedance tomography for thermal monitoring of hyperthermic treatment: an assessment using in vitro and in vivo measurements," Clinical Physics and Physiological Measurements, vol 8 (Supplement A), pp. 141-146, 1987. Copeland JR, "Radar target classification by polarization properties," Proceedings of the IRE, vol 48, pp. 1290-1296, 1960. Curran WJ Jr., Goodman RL, "Hyperthermia 1991: a critical review," Radiation research - A twentieth-century perspective, Congress Proceedings, vol 2, pp. 883-888, 1991. Demoment G, "Image reconstruction and restoration: overview of common estimation structures and problems," IEEE Transactions on Acoustics, Speech, and Signal Processing, vol 37, pp. 2024-2036, 1989. Dennis JE, Schnabel R, Numerical Methods for Unconstrained Optimization and Non Linear Equations, Englewood Cliffs, NJ, Prentice-Hall, 1983. Dewhirst MW, Moon T, Carlin D, "Analysis of tumor volume and thermal dosimetric effects on tumor response to heat, radiation and heat plus radiation: results of phase III randomized clinical trial in pet animals," Physical Aspects of Hyperthermia. Nussbaum GH (ed.), (Medical Physics mono. no. 8), pp. 495-510, 1982. Dyson JD, "Determination of the phase center and phase patterns of antennas," Radio Antennas for Aircraft and Aerospace Vehicles, W.T. Blackband (ed.), AGARD Conference Proceedings, No. 15, Slough, England, Technivision Services, 1967. Edelstein-Keshet L, Dewhirst MW, Oleson JR, Samulski TV, "Characterization of tumor temperature distributions in hyperthermia based on assumed mathematical forms," International Journal of Hyperthermia, vol 5, pp. 757-777, 1989. Golub GH, van Loan CF, Matrix Computations. The Johns Hopkins University Press, Baltimore, MD, second edition, pp. 193-259, 1989. Gregg EC, "Radiation risks with diagnostic x-rays," Radiology, vol 123, pp. 447-453, 1977. Guo TC, Guo WW, "Three-dimensional dielectric imaging by inverse scattering with resolution unlimited by wavelength," Proceedings of the Conference on Electrical Insulation and Dielectric Phenomena, pp. 65-74, 1989. 157 Guo TC, Guo WW, "Physics of image formation by microwave scattering," SPIE Medical Imaging, vol 767, pp. 30-39, 1987. Guy AW, "Analysis of electromagnetic fields induced in biological tissues by thermographic studies and equivalent phantom models," IEEE Transactions on Microwave Theory and Techniques, vol 19, pp.205-213, 1971. Habashy TM, Chew WC, Chow EY, "Simultaneous reconstruction of permittivity and conductivity profiles in a radially inhomogeneous slab," Radio Science, vol 21, Number 4, pp. 635-645, 1986. Harrington RF, Field Computation by Moment Methods, Malabar, FL, Krieger Publishing Company, 1982. Henderson RP, Webster JG, "An impedance camera for spatially specific measurements of Thorax," IEEE Transactions on Biomedical Engineering, vol 25, pp. 250-254, 1978. Hines WW, Montgomery DC, Probability and Statistics in Engineering and Management Science. New York, NY, John Wiley and Sons, pp.455-558, 1990. Hua P, Woo EJ, "Reconstruction algorithms," Electrical Impedance Tomography, JG Webster (ed), Bristol, England, Adam Higler, Ch. 10, pp. 97-137, 1990. Huynen RJ, Phenomenological Theorv of Radar Targets. Drukkerij Bronder, 1970. Hynynen K, Shimm D, Anhalt D, Stea B, Sykes H, Cassady JR, Roemer RB, "Temperature distributions during clinical scanned, focused ultrasound hyperthermia treatments," International JoumcU of Hyperthermia, vol 6, pp. 891-908, 1990. Isaacson D, "Distinguishabilities of conductivities by electric current computed tomography," IEEE Transactions on Medical Imaging, vol 5, pp. 91-95, 1986. Joachimowicz N, Pichot C, Hugonin JR, "Inverse scattering: an iterative numerical method for electromagnetic imaging," IEEE Transactions on Antennas and Propagation, vol 39, pp. 174252, 1991. Joachimowicz N, Pichot C, "Comparison of three integral formulations for the 2-D TE scattering problem," IEEE Transactions on Microwave Theory and Techniques, vol 38, pp. 178-85, 1990. Jofre L, Hawley MS, Broquetas A, de los Reyes E, Ferrando M, Elias-Fuste AR, "Medical imaging with a microwave tomographic scanner," IEEE Transactions on Biomedical Engineering, vol 37, pp. 303-312, 1990. Kallman JS, Berryman JG, "Weighted least-squares criteria for electrical impedance tomography," IEEE Transactions on Medical Imaging, vol 11, pp.284-292, 1992. Krause JD, Antennas. McGraw-Hill, New York, NY, pp. 276-78, 1950. Lagendijk JJW, "Heat transfers in tissues," Phvsics and Technology of Hyperthermia. Field SB and Franconi C (eds.), Martinus Nijhoff, The Netherlands, pp. 517-552, 1987. 158 Lagendijk JJW, Schellekens M, Schipper J, van der Linden PM, "A three-dimensional description of heating patterns in vascularized tissues during hyperthermic treatment," Physics in Medicine and Biology, vol 29, pp. 495-507, 1984. Larsen LE, Jacobi JH, "Methods of active microwave imagery for dosimetric applications," Medical Applications of Microwave Imaging. Larsen LE and Jacobi JH (eds.) IEEE Press, pp. 118-137, 1986. Lebihan D, Delannoy J, Levin RL, "Temperature mapping with MR imaging of molecular diffusion: applications to hyperthermia," Radiology, vol 171, pp. 853-857, 1989. Liao SY, Microwave Devices and Circuits. Englewood Cliffs, NJ, Prentice-Hall, Inc., pp. 1863, 1980. Lynch DR, Paulsen KD, Strohbehn JW, "Hybrid element method for unbounded electromagnetic problems in hyperthermia," International Journal for Numerical Methods in Engineering, vol 23, pp. 1915-1937, 1986. Lynch DR, Paulsen KD, Strohbehn JW, "Finite element solution of Maxwell's equations for hyperthermic treatment planning," Journal of Computational Physics, vol 58(2), pp. 246-9, 1985. MacFall J, Prescott DM, Fullar E, Samulski TV, "Temperature dependence of canine brain tissue diffusion coefficient measured in vivo with magnetic resonance echo-planar imaging," International Journal of Hyperthermia, vol 11, pp. 73-86, 1995. Mallorqui JJ, Broquetas A, Jofre L, Cardama A, "Noninvasive active microwave thermometry with a microwave tomographic scanner in hyperthermia treatments," ACES Special Issue on Bioelectromagnetic Computations, AHJ Fleming and KH Joyner (Eds.), ACES Journal, vol 7, pp.121-127, 1992. Marquardt DW, "An algorithm for least-squares estimation of non-linear parameters," SIAM Journal of Applied Mathematics, vol 11, pp. 431-41, 1963. McArdle FJ, Brown B, Angel A, "Imaging cardiosynchronous impedance changes in the adult head," EC Concerted Action on EIT, London, CAIT, 1992. McRae DA, Esrick MA, "Changes in electrical impedance of skeletal muscle measured during hyperthermia," International Journal of Hyperthermia, vol 9, pp. 247-262, 1993. Meaney PM, Paulsen KD, Ryan TP, "Two-dimensional hybrid element image reconstruction for TM illumination," IEEE Transactions on Antennas and Propagation, vol 43, pp. 239-247, 1995. Meaney PM, Paulsen KD, Hartov A, Crane RK, "An active microwave imaging system for reconstruction of 2-D electrical property distributions," IEEE Transactions on Biomedical Engineering, (In press (a)). Meaney PM, Paulsen KD, Hartov A, Crane RK, "Multi-target microwave imaging for tissue assessment: initial evaluation in tissue equivalent phantoms," IEEE Transactions on Biomedical Engineering, (Submitted (b)). 159 Miyakawa M, "Tomographic measurement of temperature change in phantoms of the human body by chirp radar-type microwave computed tomography," Medical and Biological Engineering and Computing, vol 31, pp. S31-S36, 1993. Moros EG, Button AW, Roemer RB, Burton M, Hynynen K, "Experimental evaluation of two simple thermal models using hyperthermia in muscle in vivo," International Journal of Hyperthermia, vol 9, pp. 581-598, 1993. Moskowitz MJ, Ryan TP, Paulsen KD, Mitchell SE, "Clinical implementation of electrical impedance tomography with hyperthermia," International Journal of Hyperthermia, vol 11, pp. 141-149, 1995. Moskowitz MJ, "Electrical Impedance imaging for temperature field estimation during hyperthermia treatment of cancer," Ph.D. Dissertation, Thayer School of Engineering, Dartmouth College, pp. 151-184, 1994. Ney MM, "Method of moments as applied to electromagnetic problems," IEEE Transactions on Microwave Theory and Techniques, vol 33, pp. 972-980, 1985. Ohtera I, Ujiie H, "Nomographs for phase centers of connical corrugated and TE; i mode horns," IEEE Transactions on Antennas and Propagation, vol 23, pp. 858-9, 1975. Oleson JR, Dewhirst MW, Harrelson JM, Leopold KA, Samulski TV, Tso CY, "Tumor temperature distributions predict hyperthermia effect," International Journal of Radiation Oncology Biology and Physics, vol 16, pp. 559-570. Overgaard J, Overgaard M, "Hyperthermia as an adjuvant to radiotherapy in the treatment of malignant melanoma," International Journal of Hyperthermia, vol 3, pp. 483-501, 1987. Paulsen KD, Meaney PM, Moskowitz MJ, "A dual mesh scheme for finite element based reconstruction algorithms," IEEE Transactions on Medical Imaging, (In Press (a)). Paulsen KD, Moskowitz, MJ, Ryan TP, "Temperature field estimation using electrical impedance profiling: I. Reconstruction algorithm and simulated results," International Journal of Hyperthermia, vol 10, pp. 209-228, 1994. Paulsen KD, Lynch DR, Liu W, "Conjugate direction methods for helmholtz problems with complex-valued wavenumbers," International Journal for Numerical Methods in Engineering, vol 35, pp. 601-22, 1992a. Paulsen KD, Liu W, "Memory and operations count scaling for coupled finite-element and boundary-element systems of equations," International Journal for Numerical Methods in Engineering, vol 33, pp. 1289-1304, 1992b. Paulsen KD, Lynch DR, "Calculation of interior values by the boundary element method," Communications in Applied Numerical Methods, vol 5, pp. 7-14, 1989. Paulsen KD, Lynch DR, Strohbehn JW, "Three-dimensional finite, boundary, and hybrid element solutions of the Maxwell equations for lossy dielectric media," IEEE Transactions on Microwave Theory and Techniques, vol 36, pp. 682-693, 1988. 160 Paulsen KD, Strohbehn JW, Lynch DR, "Comparative theoretical performance of two types of regional hyperthermia systems," International Journal of Radiation Oncology Biology and Physics, vol 11, pp. 1659-1671, 1985. Pennes HH, "Analysis of tissue and arterial blood temperatures in the resting human forearm," Journal of Applied Physiology, vol 1, pp. 93-122, 1948. Perez CA, Pajak T, Emami B, Hombeck NB, Tupchong L, Rubin P, "Randomized phase III study comparing irradiation and hyperthermia with irradiation alone in superficial measurable tumors," American Journal of Clinical Oncology, vol 14, pp. 133-141, 1991. Perl W, Hirsch RL, "Local blood flow in kidney tissue by heat clearance measurement," Journal of Theoretical Biology, vol 10, pp. 180-251, 1966. Rius JM, Pichot C, Jofre L, Bolomey JC, Joachimowicz N, Broquetas A, Ferrando M, "Planar and cylindrical active microwave temperature imaging: numerical simulations," IEEE Transactions on Medical Imaging, vol II, pp. 457-469, 1992. Roemer RB, "Optimal power deposition in hyperthermia. I. The treatment goal: the ideal temperature distribution: the role of large blood vessels," International Journal of Hyperthermia, vol 7, pp. 317-341, 1991. Roemer RB, "Thermal dosimetry," Thermal Dosimetry and Treatment Planning. Gautherie M (ed.), Springer-Verlag, Berlin, pp. 119-208, 1990. Samulski TV, MacFall J, Zhang Y, Grant W, Charles C, "Non-invasive thermometry using magnetic resonance diffusion imaging: potential for application in hyperthermic oncology," International Journal of Hyperthermia, vol 8, pp. 819-830, 1992. Schwan HP, Li K, "Capacity and Conductivity of Body Tissues at Ultrahigh Frequencies," Proceedings of IRE, vol 41, pp. 1735-1740, 1953. Sinclair G, "The transmission and reception of elliptically polarized waves," Proceedings of IEEE, vol 39, pp. 535-540, 1951. Skolnik MI, Introduction to Radar Systems. New York, NY, McGraw-Hill Brook Company, 1980. Slaney M, Kak AC, Larsen LE, "Limitations of imaging with first-order diffraction tomography," IEEE Transactions on Microwave Theory and Techniques, vol 32, pp. 860-873, 1984. Strohbehn JW, Paulsen KD, Lynch DR, "Use of the finite element method in computerized thermal dosimetry," Physical techniques in clinical hyperthermia. Hand JW and James JR (ed.). Research Studies Press, Letchworth, England, pp. 383-451, 1986. Tikhonov AN, Arsenin VY, Solutions of Ill-Posed Problems, Washington, DC, V.H. Winston & Sons, 1977. Valdagni R, Amichetti M, Pani G, "Radical radiation alone versus radical radiation plus microwave hyperthermia for Ng (TNM-UICC) neck nodes: a prospective randomized clinical trial," International Journal of Radiation Oncology, Biology, Physics, vol 15, pp. 13-24, 1988. 161 Weinbaum S, Jiji LM, "A new simplified bioheat equation for the effect of blood flow on local average tissue temperature," Journal of Biomechanical Engineering, vol 107, pp. 131-139, 1985. Woo EJ, Hua P, Webster JG, Tompkins WJ, "Finite-element method in electrical impedance tomography," Medical & Biological Engineering and Computing, vol 32(5), pp. 530-536, 1994. Yashiro K, and Ohkawa S, "Boundary element method for electromagnetic scattering from cylinders," IEEE Transactions on Antennas and Propagation, vol 33(4), pp. 383-9, 1985. Yorkey TJ, Webster JG, and Tomkins WJ, "Comparing reconstruction algorithms for electrical impedance tomography," IEEE Transactions on Biomedical Engineering, vol 34, pp. 843-52, 1987. Yorkey T, "Comparing reconstruction methods for electrical impedance tomography," Ph.D. Dissertation, University of Wisconsin-Madison, 1986. APPENDIX A. Forward Solution of the Electric Field Distribution Using the Hybrid Element Method Starting with the system of equations (equations II.2.4-5), the finite element region > < ?1 o .^bl 1 ^11 CP produces a matrix relationship which can be decomposed as: J E . L "0 lEbJ 0 Bbb. (A.l) pj p [FsJ "Dbb 1 O DJ Csb u "Cbb D.b P 1 and the boundary element region yields a similarly written matrix form: fE. (/\.2) where the subscripts I, b, and s refer to all points interior to the FE surface, all points on the boundary shared by the FE and BE problems and all points on the radiating source(s), respectively. In equation (A.l), the sparse matrix A is partitioned into four submatrices of various dimensions, A?, A,b, A^, and Ayy where, for example, submatrix A,b has dimensions n, x n^ where n, , n^, and n, are the number of nodes in regions I, b, and s, respectively. Similar partitioning is used in column vectors E and F where F=^. This on notation also appears in the expression of matrix system (A.2). For convenience in some instances, we express the partitioned matrices in (A.l) and (A.2) in a nonpartitioned form in which case the subscripts are dropped in favor of a single alphabetical symbol. At this point, there are three methods of solving these two sets of simultaneous equations. For purposes of the appendix, only hybrid element methods 1 and 2 (HE #1 and HE #2) will be described here, and the reader is referred to Paulsen, et al. [1990] for a discussion of HE #3. A.l Hybrid Element Method #1 Using HE #1 technique, we first solve for equation (A.2): 62 > in terms of 4 E^ E. ^ by inverting [C] in 163 "Cbb Cbs" Csb Css. -1 "Dbb Dbs" mi "Gbb Gbs' El Dsb Dss. l E s J Gsb Gss. Es (A.3) solving this for {Fy} and multiplying by matrix [B] from equation (A.l) yields: (A.4) [B1{FJ = (BlGM,]{Eb} + [Bl[Ota]{E,} This expression can be substituted into equation (A.l) to give: E, Mb .^bl ^bb ~[B][Gbb]. (A.5) EbJ 1[B1[06.1{E,}J which will allow for the computation of the electric field at all the nodes within the finite element region and on its boundary given the radiation distribution on all exterior sources. A.2 Hybrid Element Method #2 The HE #2 approach is again based on equations (A.l) and (A.2), except that we first fEil solve for jg| i" equation (A.l) by inverting matrix [A]: AI? AI Albi AI| I:: (A.6) where [AI] = [A] '. Equation (A.6) provides a relationship between {Ej,} and {Fj,} alone: (A.7) {Eb} = [AIi,b][B]{Fb} which can be substituted into equation (A.2) to yield: u mi "Dbb 1" Q Csb Css [ F s J Dsb I J f[AIbbl[BKFb}l_ Hbb ' DJ Cbb (Ej r Hsb H bs E. (A.8) 164 where Hbb H(,s _ D(,b D|,s [AIbb][B] Hsb ^ss. ^ss. .^sb and [Ijs] is the identity matrix of rank s. [iss] Equation A.8 can be rearranged to produce: (A.9) Equation (A.9) can be solved for {Fy} which can then be substituted back into equation (A. 1) to recover the electric fields everywhere within the finite element mesh. When considering only the forward solution, it is important to note that HE #2 has slight computational advantages over HE #1 [Paulsen, et al., 1990]. In He #2, the initial step is to invert matrix [A] which involves decomposing a sparse system of equations and then back substituting various unit vectors to recover the necessary columns of [Al]. In this case, only [AIj,b] is needed, which implies only n^ back substitutions are required. For HE #1, the initial step is to invert matrix [C] which involves decomposing a full system of equations followed by n^ + n, back substitutions. Hence, in general, there is a slight computational advantage associated with HE #2, primarily due to its reduced number of back substitutions. APPENDIX B. Algorithmic Computational Time Enhancements There are numerous methods for reducing overall computation time of any computer program, from purchasing faster computers to using computer specific optimized functions to optimizing various features of the algorithm itself based on the problem at hand. This section will adopt the latter approach in optimizing two tasks: (1) reducing the cost per iteration resulting from calculating the boundary element contribution to the forward solution, and (2) reducing the cost per iteration associated with calculating various indices and quantities required for the coordination of the electric and B.l meshes. Boundary Element Contribution to Iteration Time As mentioned in Chapter II, hybrid element method #1 (HE #1) offers significant computational advantages over HE #2 when used as the basis of this reconstruction algorithm. As described here, it also offers significant advantages over HE #2 by further reducing the costs per iteration by computing matrices associated the boundary element contribution once in a preprocessor routine. In using any of the hybrid element methods, the matrices associated with the boundary element discretization are solely based on the geometry of the mesh in the boundary element region and the material properties of that region. Thus, their computation will not change with iteration. This means that these matrices can be computed once in a pre� processor routine and stored for later use at each iteration. Using HE #1, we can take this one step further: not only can we compute matrices [C] and [D] of equation (II.2.5) in a pre� processor routine, but we can also invert [C] and multiply its inverse with [D] in that same pre-processor and only store their product, [G] (equation II.2.8). This means that using HE #1, there are no full matrix inversions computed during each iteration. On the other hand, if HE #2 were to be used we would still have to perform a full inversion of matrix [A], because its entries are functions of the material properties in the finite element region which change with 165 166 iteration. The algorithm would also have to store both [C] and [D] which would be nontrivial since both are full complex matrices. B.2 Dual Mesh Computational Overhead and its Contribution to Iteration Time It is understood that when using this imaging system, that the target region, measurement sites and radiator positions are all pre determined. So, too, are the electric mesh and its associated mesh. This is important because, as described in Chapter III and Appendix C, introduction of the dual mesh method has significant computational ramifications related to the coordination of the two meshes. However, these relationships are strictly functions of the geometry of the two meshes and not of the various material properties. Thus, information can be computed prior to execution of the reconstruction algorithm which significantly reduces the computational costs associated with each iteration. APPENDIX C. Dual Mesh Overhead As described in Chapter III, the dual mesh method requires significant overhead in the form of various indices and vectors of basis function values in order to make its implementation efficient. The implications in terms of the computational costs of the method are felt primarily in computing (1) the forward solution, and (2) the Jacobian matrix used for determining the material property updates at each iteration. Figure C.l illustrates the problem for the forward solution: k element electric element Figure C. 1. Diagram of relationship between solution. 167 mesh and electric mesh during the forward 168 In calculating the matrix, [A] (described in equation 11.2.4), integrations of the active basis functions are performed over each electric element. To accomplish this, the material properties at nodes Ei, E2, and Eg must be calculated. In general, however, these values are only known at the nodes of the k^ mesh. To compute these material property values, several pieces of information need to be calculated beforehand. These are; (1) NIK(i) - This arrayt stores the k^ element number which contains the i^h electric node number. (2) NIN(i,l-3) - This array stores the three k^ node numbers associated with the k^ element which contains the i^^ electric node. (3) PHINK(i,l-3) - This array stores the three basis function values associated with the three k^ nodes in the NIN array for the i^h electric node. Given this information and the material property values at the k^ nodes, the inner products required in forming [A] can be performed. In calculating the Jacobian matrix, [J] (described in equation II.2.2), determining the appropriate matrix contributions is a bit more involved. Figure C.2 illustrates one configuration for this situation. Essentially, the elements of [J] are the derivatives of the elements of [A] with respect to the k^ values at the k^ nodes. These derivatives involve integrations which are performed over the electric elements. To accomplish this, the algorithm (1) loops over all of the k^ nodes, (2) finds all of the electric elements for which the basis function centered at the current k^ node is non-zero, (3) calculates the basis function value (whose center is located at the current k^ node) at each of the electric element nodes, and (4) computes the required integrations over the electric elements using this information. To accomplish these tasks efficiently, several more pieces of information need to be pre-computed and stored including: 169 (1) ISKE(j,l-NEAMAX) - This array stores the electric element numbers associated with the node. Since there is no set number of electric elements associated with each node, a maximum value is set to limit the array size. (2) lSKN(j,NNAMAX) - This array stores the electric node numbers associated with the k^ node number. A maximum number of electric nodes associated with each k^ node is set. (3) PHIK(j,I-NNAMAX) - This array stores the basis function value at the electric nodes given in ISKN due to the k^ basis function centered at the k^ node, (|>j. Given this information, the inner products required in forming [J] can be performed. k node Electric Element k elements associated with |*^k node Figure C.2. Diagram of relationship between k^ mesh and electric mesh during calculation of the Jacobian matrix, [J]. APPENDIX D. Noise Figure Analysis As mentioned in Chapter IV, one of the keys to achieving a large dynamic range is to improve the receiver signal-to-noise ratio. The most widely recognized figure of merit for this is the noise figure, Fn. By definition, the noise figure can be written: /'^out where S is the signal power, N is the noise power and the subscripts in and out refer to the input and outputs of the receiver channel, k is Boltzmann's constant, Tq is the receiver temperature in " Kelvin, Bn is the noise bandwidth and G is the gain of the receiver, where S / G= . Fn in its simplest terms is an indication of how much the signal-to-noise ratio / ^in has been degraded as a result of the signal having passed through the cascade of receiver components. In general, the noise figure of individual components is readily available, and from this data a composite noise figure can be determined for the whole receiver. Figure D. 1 shows the cascade of components through which a signal passes in one of the system channels. (Note that the gain and noise figure for the digital attenuators are for base attenuation levels.) The total noise figure for the system, \ o , can be written, + P k -_l . ^2zl , + llzl+_+ G, G,G2 G,G2...Gk_, (D.2) where k is the total number of components in the cascade. From equation D.2, it makes sense to have a low noise amplifier, LNA, as the first component in the cascade which will have the effect of keeping the first term, F i, small and will reduce all subsequent terms as long as the gain products are greater than 1. For the system designed the cascaded noise figure is, F^io^i =2.8 dB. This is quite good given the wide band, 3(X) - 1100 MHz, of operation. 170 171 LNA Switch Multiplexer Amplifier Digital Attenuator F. = 6dB 4 Fg= 6 dB In R, = 2.4 dB = +15 dB Amplifier F^=1 dB G, = -1 dB Digital Attenuator F3= 1 dB = -1 dB = +20 dB Amplifier Mixer Fg=4dB F^=6dB Fo=4dB F?=7dB C^ = +15dB (^ = -6 dB (^ = +15dB C^ = -6dB Figure D. 1 Diagram of receiver channel with associated noise figures and gains. % == -6?' dB APPENDIX E. Water-filled Waveguide Return Losss As was shown in Chapter IV, the measurement system utilizes the highly lossy nature of the 0.9% saline medium to improve the characteristic impedance of the monopole receiving antennas to such a degree that they are usable over a broad frequency range. In contrast, the water-filled transmitting waveguide generally has a broad usable bandwidth relative to that of the receiving monopole antenna. As seen in Figure E.1, the return loss is generally better than -15 dB from 400 to 900 MHz when the waveguide is immersed in DI water. Characteristic Impedance Water-filled Waveguide -10 CD ?o -20 C k- 3 -30 CC -40 DI Water Saline Solution (0.9%) -50 L- 100 300 500 700 900 1100 Frequency (MHz) Figure E.1. Characteristic impedance of water-filled waveguide antenna when immersed in DI water and a 0.9% saline solution, respectively. 1 72 173 Its characteristic impedance is further improved when immersed in saline solution (Figure E.1). in which case it has a usable bandwidth from 300 to 900 MHz, with return losses of better than -10 dB. The curve is also dramatically smoother over the full bandwidth, indicating that such things as standing waves and normally unterminated currents are significantly attenuated by the highly conductive medium. Thus, from the available data, this antenna's characteristics appear well suited for this application. APPENDIX F. Radiation Field Variation Outside of the 2-D Plane The calibration procedures described in Chapter IV assume that the main beams of both the transmitting and receiving antennas are centered in the target plane. Positioning either of these antennas above or below such a level could contribute significant amplitude and phase errors and subsequently distort the reconstructed images. Given the present antenna configuration, it is difficult to measure a true radiation pattern by taking measurements in an arc in the far field. However, it is possible to measure the variations in the received signals txansmitted from these antennas due to vertical position variations. Transmitted l\1agnitude vs. Position Water-filled Waveguide and Monopole Antennas m TJ -2 0) "D 5 -4 C O) (0 > 300 MHz (Waveguide) -6 500 MHz (Waveguide) 700 MHz (Waveguide) 0> oc 300 MHz (Monopole) 500 MHz (Monopole) 700 MHz (Monopole) -10 -3.5 -2.5 -1.5 -0.5 0.5 1.5 2.5 3.5 Position (cm) Figure F.l. Variations in the magnitude of the transmitted power for the water-filled waveguide and monopole antennas as a function of relative vertical position at 300, 500 and 700 MHz. 174 175 In this case, the antennas under test are either the water-filled waveguide or the monopole antenna, while the probe antenna is a specially modified monopole antenna whose vertical position can be changed. The antenna under test and the probe antenna are placed at set distances from each other and the transmitted signals are recorded for the various probe antenna vertical settings. Figure F.l shows the variations in the measured magnitude with probe antenna position for the waveguide and monopole positions. (The position on the abscissa is relative to the point where the center of the probe antenna is aligned with the center of the antenna under test.) Figure F.2 shows a similar plot for the phase variations for both types of antennas. Transmitted Phase vs. Position Water-filled Waveguide and Monopole Antennas -20 O) -40 300 MHz (Waveguide) 500 MHz (Waveguide) 700 MHz (Waveguide) -60 300 MHz (Monopole) 500 MHz (Monopole) 700 MHz (Monopole) -80 -3.5 -2.5 1.5 -0.5 0.5 1.5 2.5 3.5 Position (cm) Figure F.2. Variations in the phase of the transmitted power for the water-filled waveguide and monopole antennas as a function of relative vertical position at 300, 500 and 700 MHz. 176 In general, there appears to be a region near the zero relative position where the phase and amplitude are quite flat for the frequency range of 300 to 700 MHz. This suggests that we are not operating these antennas in a highly variable region of their respective antenna patterns. Thus, if care is taken to use constant height settings, vertical positioning should not contribute significantly to the overall measurement error. APPENDIX G. Calibration of Variations Between Transmitter Channels and Receiver Channels The first stage of the calibration process is to remove differences in the phase and amplitude of all four transmit and all four receiver antennas. We accomplish this by using a least squares technique. Figure G. 1 illustrates the antenna configuration for our calibration scheme 0 - - z '4 V \ o / \ I \ / ' / \ l / Figure G.l. Antenna arrangement for calibration of amplitude and phase offset differences between antennas. where T1-4 correspond to the four transmitters and R1.4 correspond to the four receivers. Lj is the distance from a given transmit antenna to its closest receiver antenna and L] is the 177 178 distance from that same transmit antenna to the two equidistant receiver antennas. The data acquisition system is set up such that each transmit and receive channel has an arbitrary amplitude and phase offset with respect to an amplitude of 0 dBm and a phase of 0 degrees associated with it. For transmitters 1-4 these are At(1-4) and Pt(1-4). respectively, and for the receivers are Ar(1-4) and Pr(1.4). The amplitude and phase offsets due to the signals travelling lengths, Li and L2, are Ai, A2 and Pi, P2, respectively. Using this notation, equations for the amplitude and phase can be written for the signals detected at receiver, i, due to a signal sent from transmitter, j, having traveled a length Lk: ^i.j - ^Ri + ^Tj + (G-1) Pi.j = ^Ri + ^Tj + Pk + Pg (G-2) where Ag and Pg are the amplitude and phase error terms associated with each measurement. Note that for each transmitter, we are only using three of the four receiver measurements sites possible. This is mainly due to the fact that measurements at the farthest receiver site are partially disturbed by the presence of the nearest receiver antenna. (From experience with this measurement system, it is known that the measurements at a selected receiver antenna will be slightly perturbed if another receiver antenna is in the direct line of propagation from the given transmitter antenna. Therefore, care must be taken to use only measurement data that is unaffected by such perturbations.) As a result, there are twelve measurements that can be used. This yields twelve equations and eight unknowns for both the amplitude and phase. (For both amplitude and phase, values are referenced to transmitter and receiver number 1.) For the amplitudes, the unknowns are: At2. At3, At4, Ar2, AR3, AR4, Ai, and A2 and similarly for the phases, they are: Pt2. PT3. PT4. PR2. PR3. PR4. Pi, and p2 amplitude relationships can be written in matrix form as: The 179 ^4,1 "0 0 0 0 0 1 0 r 0 0 0 0 0 0 1 0 0 0 0 1 At2 ^2,1 AT3 ^1.2 0 AT4 ^2,2 0 0 1 AR2 ^3,2 0 0 0 1 AR3 ^2,3 0 1 AR4 A 3.3 A, A 4,3 0 0 0 1 1 0 0 0 0 0 0 1 1 0 0 1 I 0 0 0 1 0 0 1 0 1 0 1 0 1 0 0 1 0 1 0 0 0 1 0 0 1 0 1 0 0 1 0 0 1 .0 0 1 ^1.1 0 0 1 0 0 1 1 0 A; 0 0 0 0 1. (G.3) ^3,4 ^4,4 i.4 where the matrix on the left is called the Jacobian, J, the vector on the left contains the desired unknowns, x, and the vector on the right contains the measured values, m. Multiplying both sides of the equation by yields: [fj]{x} = [jt]{m} Here, (g.4) is the Hessian matrix and can be solved easily using a Cholesky factorization. This yields the relative values of At2. At3, and At4 with respect to the amplitude offset from transmitter 1 and Ar2, Ar3, and Ar4 with respect to the amplitude offset from receiver 1. The exact same procedure is applied to the relative phase values. By using these quantities to modify the subsequent measured values, all data will be referenced to transmitter 1 and receiver 1. At this point, all of the data is processed by the transformation from 3-D to 2-D radiation via the method outlined in Chapter IV. The data is then ready for use in the reconstruction algorithm. APPENDIX H. Spatial Low Pass Filtering of Reconstructed Property Distribution As mentioned in Chapter IV, when using measured data with the reconstruction algorithm, a low pass filter must be applied to the updated material property distribution at each iteration to keep the solution from diverging. This is accomplished through a weighted averaging of the material property value, k^, at each node with the estimated values associated with its surrounding nodes. Figure H. 1 shows a portion of a mesh surrounding node A for illustration purposes: Figure H.1. Portion of k mesh to illustrate low pass filter technique. 180 181 The averaging for each node is accomplished using the following formula; k('?ew)(A) = (1-8) kfo,d)(A) + e ) (H. 1) i=l where N is the number of surrounding nodes (6 in this example) and 0 is the weighting factor. Equation (H.l) is then applied to the values for every node. Because each new k^ value is a linear combination of all the old k^ values, the procedure can be written as a system of equations such that: Fi where jk^^ew)} and k(old) 1,2 ?^(new) are vectors of the new and old k^-distribution values over the mesh. F] is a square sparse matrix with the associated weightings from equation (H.l). This matrix can then be applied multiple times such that: [F2]{k(old)} = Knew)} where [F2] = [F, (H.3) which has the effect of extending the spatial reach of the averaging at each node (For the images in Chapters IV and V, the filter was applied twice.). For our purposes, the value of 0 has been determined empirically. Because this matrix is only a function of 0 and the k^ mesh, it can be calculated prior to execution of the reconstruction algorithm and applied at each iteration simply by pre-multiplying the material property vector. As a result, the filtering routine adds minimally to the overall computational costs. r the inverse of the k2 distribution which also tends to reduce the effects of ripple in the background region which could skew the center position estimate. As an example, the center of the real part of an object in the x-direction is calculated using the following formula: ^Re(k^(C(i)))' XRE (center) = ?? (V.4.B. V .4) |j'Re(k^(C(i))) where (center) is the x-coordinate of the real part of the object, and similarly to earlier conventions, k^(C(i)) is the k^ value at the center of element /, A(i) is the area of element i and N is the total number of elements in the selected region. With minor alterations, this formula is used to calculate X,M (center), Y^g(center) and (center). In all of these tests, only the 4.28 and 2.54 cm diameter fat/bone equivalent phantoms have been analyzed to simplify the comparisons. Similar trends were noted with phantoms of varying contrast levels. Table V.4.B.v.l shows a summaiy of RMS errors for the single object experiments versus frequency for both the background and object regions. It is clear 130 4.28 cm diam. Object RMS ERRORS Frequency (MHz) REAL 300 500 700 IMAGINARY 300 500 700 Samples per wavelength 28.5 18.4 13.5 Samples per exponetid decay 8.74 7.70 7.18 Background Object .153 .095 .085 6.34 4.81 3.90 .132 .135 .121 52.8 57.8 69.7 2.54 cm diam Object RMS ERRORS Frequency (MHz) REAL 300 500 700 IMAGINARY 300 500 700 Samples per wavelength 28.5 18.4 13.5 Samples per exponetid decay 8.74 7.70 7.18 Background Object .133 .083 .062 9.49 5.38 3.07 .099 .103 .091 94.0 91.9 90.0 Table V.4.B.v.l. Average RMS enors for the Real and Imaginary parts of the reconstructed material electrical properties for tests performed at 300, 500 and 700 MHz. Error values were averaged for the 4.28 cm and 2.54 cm diameter fat/bone equivalent cylinders at three distances (0.0, 2.54, and 3.81 cm) from the center of the target region along the horizontal axis. that there is an improvement in the errors with frequency for the real part of the images for both sizes of objects in both the background and object regions. It is not clear, however, what trends are emerging for the imaginary part of the images, and in the case of the object region of the 4.28 cm diameter target, it appears that the error actually increases with frequency. Similar results were observed in the multi-target cases as shown in Table V.4.B.V.2. Again the real part of the images (background and object) show improvement with frequency. In these cases. 1 3I Left Object RMS ERRORS Frequency (MHz) REAL 300 500 700 IMAGINARY 300 500 700 Samples per wavelen^ 28.5 18.4 13.5 Samples per exponetied decay 8.74 7.70 7.18 Background Object .242 .163 .154 6.32 4.38 3.72 .201 .261 .275 49.1 59.0 63.2 Right Object RMS ERRORS Frequency (MHz) REAL 300 500 700 IMAGINARY 300 500 700 Samples per wavelength 28.5 18.4 13.5 Samples per exponetid decay 8.74 7.70 7.18 Background Object .242 .163 .154 11.1 3.91 3.99 .201 .261 .275 88.5 108.5 120.4 Table V.4.B. V .2. Average RMS errors for two cylindrical objects made out of fat/bone equivalent material with the actual diameters being 4.28 cm and 2.54 cm for the left and right objects, respectively. Averages were compiled for three test cases where the two objects were separated by 0.4, 1.7, and 4.2 cm, respectively. the imaginary part of the images show increased degradation with frequency. Some insight into this phenomenon can be achieved by observing the mesh sample rates of the fine mesh as shown in Table V.4.B.V.1 and 2. In all cases, the number of samples per wavelength is greater than 10 which is a benchmark established in Chapter II for accurately modeling the forward solution. However, as the frequency increases, the mesh sampling rate per exponential decay approaches 7 which is also a benchmark for accurate forward solution 132 modeling. Thus, because of the high conductivity of the medium, a higher level of discretization may be required for more accurate forward solution modeling in an effort to compensate for the observed increase in RMS error with frequency. Tables V.4.B.V .3 and 4 list the results of the diameter and center of object computations for the 4.28 and 2.54 cm fat/bone equivalent cylinders in both the single and multi-target experiments, respectively. In general the real part of the images more accurately captures the cylinder diameters, especially for the larger of the two objects. It does not appear 4.28 cm diam. Object Frequency (MHz) Average Diameter (cm) Average Distance from exact location (cm) 300 500 700 4.31 4.28 4.02 .127 .036 .078 300 500 700 4.09 3.63 3.16 .061 .068 .064 Average Diameter (cm) Average Distance from exact location (cm) 300 500 700 3.23 2.83 2.59 .079 .048 .035 300 500 700 3.42 2.95 2.37 .069 .040 .030 REAL IMAGINARY 2.54 cm diam. Object Frequency (MHz) REAL IMAGINARY Table V.4.B.v.3. Average diameters and position errors for the two fat/bone equivalent cylinders with the actual diameters being 4.28 cm and 2.54 cm, respectively, during singletarget imaging experiments. Values were averaged for the cylinders at three distances (0.0, 2.54, and 3.81 cm) from the center of the target region along the horizontal axis. 133 Left Object Frequency (MHz) Average Diameter (cm) Average Distance from exact location (cm) 300 500 700 4.04 4.00 3.85 0.572 0.517 0.398 300 500 700 2.96 3.14 2.94 0.642 0.648 0.721 Average Diameter (cm) Average Distance from exact location (cm) 300 500 700 2.96 3.14 2.94 0.751 0.828 0.830 300 500 700 3.27 2.81 2.38 0.857 0.973 0.973 REAL IMAGINARY Right Object Frequency (MHz) REAL IMAGINARY Table V.4.B. V .4. Average diameters and position errors for the two fat/bone equivalent cylinders with the actual diameters being 4.28 cm (left object) and 2.54 cm (right object), respectively, during multi-target imaging experiments. Averages were compiled for three test cases where the two objects were separated by 0.4, 1.7, and 4.2 cm, respectively. that the algorithm's ability to accurately capture the object diameters is significantly better in either the single or multi-target cases. However, the center of the object computations show a clear improvement when only one object is imaged at a time. In these cases, the positioning error is generally 1 mm or less for both object sizes and at all three frequencies. In contrast, the errors in the object positions of the multi-target cases were generally between 4 and 10 mm with roughly 90% of these errors occurring along the common axis of the two objects. Thus the accuracy of the object position is clearly affected by the presence of other objects. 1 34 Generally, the diameter estimates for the targets appear to be within 1 cm with typical errors being more on the order of a few millimeters. Similarly, the position estimates are correct to better than 1 cm and quite often within 5 mm or better. V.5. Discussion and Conclusion A prototype microwave imaging system has been developed which can achieve reconstructed images of the material property distributions for a wide range of material contrasts over a frequency range of 300 to 700 MHz. Particular attention has been paid to the amount of data necessary to reconstruct these images along with selecting that data which is most beneficial in such reconstructions. In the present configuration, the system was able to achieve relatively good results even for cases where the Jacobian was under-determined, with incremental improvement being achieved through the use of a two-tiered receiver system consisting of an additional arc of receiver sites outside the target region. Comparison of results at different frequencies demonstrated the capability of reconstructing relatively good images over the full 300-700 MHz frequency range with improvements occurring as the frequency is increased. These observations are based on the ability of the algorithm to achieve high contrast between two different materials, reconstruct a relatively uniform background medium and distinguish between two closely spaced objects. In the more challenging tests of the algorithm's capability, objects with an assortment of electrical properties were imaged. The results of these comparisons showed a clear advantage for the higher frequency excitations, where both the dielectric constants and conductivity were reconstructed quite well at 700 MHz, whereas similar tests at 300 MHz were only able to recover the conductivity of the objects. In practice, increasing frequency will have to be balanced with the increasing attenuation in the biological tissue. Since the microwave components of this system are capable of operating over a wide spectrum, 300-1100 MHz, we have some flexibility in this regard. An additional option which we have not explored to date is the combination of data from multiple frequencies in order to better determine target material properties. We also have 135 the flexibility of modifying the conductivity of the background medium. All of these will ultimately be key issues in determining the optimal size of the target region and the operating frequencies. This work also poses some interesting questions as to how much and what kind of measurement data is most useful for image reconstruction. For practical considerations the background medium has been manipulated and the reconstruction algorithm formulated in such a way as to render signal reflection data minimally useful. These imaging experiments are being conducted under conditions which represent a compromise between microwave systems that utilize reflection, refraction and transmission data and other imaging modalities such as CT which rely primarily on transmission data. These relationships will be further investigated with the implementation of the next generation system. 36 Re(k2) (m"2) Im(k2) (m'2) 1200 2000 3000 2400 2800 3600 4000 4200 4800 5400 300 MHz 2000 3000 4000 5000 6000 7000 8000 9000 10000 2000 4000 6000 8000 10000 12000 14000 16000 18000 20000 2000 3000 5000 6000 8000 9000 0000 500 MHz 400 2800 4200 8400 11200 12600 14000 700 MHz Figure V.4.B.iv. 1. Real and imagina^ parts of the distribution of the reconstructed images for the 1.6 cm diameter fat/bone equivalent cylinder centered in the target region at 300, 500 and 700 MHz using both the inner and outer banks of measurement sites. 137 RE(k 2) (m 2) IM(k ^ ) 4000 6000 3000 4500 2000 3000 1000 1500 (m ) 300 MHz 10000 10000 7500 7500 5000 5000 2500 2500 0 0 500 MHz 14000 10500 ?0.08 -0.04 0 0.04 0.08 -0.08 X - position (m) -0.04 0 0.04 X - position (m) 700 MHz F i g u r e V. 4 .B. i v . 2 . Electrical property p r o f i l e a l o n g l o n g i t u d i n a l transects parts of t h r o u g h the t.Trqet region , resper-t i'./>:? 1 y . for the real and imaginary P l o t s correspond to reconstructed images shown in Fiqnr'-.- V. 1.!. 0.08 138 Re(k2) (m"2) 1600 2000 2400 2800 3200 3600 4000 Im(k2) (m"2) 2400 3600 4200 4800 5400 300 MHz 2000 2000 3000 4000 5000 4000 5000 6000 6000 7000 8000 9000 10000 7000 8000 9000 0000 500 MHz 0 2000 4000 16000 8000 110000 2800 4200 2000 114000 16000 18000 20000 9800 11200 12600 14000 700 MHz Figure V.4.B.iv.3. Real and imaginary parts of the distribution of the reconstructed images for the 2.5 cm diameter fat/bone equivalent cylinder centered in the target region at 300, 500 and 700 MHz using both the inner and outer banks of measurement sites. 139 300 MHz 10000 10000 7500 7500 5000 5000 2500 2500 I 0 0 500 MHz 14000 10500 7000 3500 -0.08 -0.04 0 0.04 0,08 -0.08 X - position (m) -0.04 0 0.04 X - position (m) 700 MHz F i g u r e V. 4 .B. i v . 4 . Electrical property profile a l o n g l o n g i t u d i n a l , transects parts of through the target region for the real and imaginary , respectively. [lots correspond to reconstructed i m a g e s s h o w n i n F i q u i f- .-1 . P . i " . ?. 0.08 140 Re(k2) (m-2) Itn(k2) (m-2) 2400 2400 3600 2800 3600 4000 6000 300 MHz 3000 4000 5000 6000 7000 8000 9000 10000 0000 500 MHz 2000 4000 6000 8000 0000 2000 4000 6000 8000 20000 4200 8400 200 2600 14000 700 MHz Figure V.4.B.iv.5. Real and imagina^ parts of the distribution of the reconstructed images for the 4.3 cm diameter fat/bone equivalent cylinder centered in the target region at 300, 500 and 700 MHz using both the inner and outer banks of measurement sites. 141 RE(k 2 ) (m 2 IM(k ^) 4000 6000 3000 4500 2000 3000 1000 1500 (m ) 300 MHz 10000 10000 7500 7500 5000 5000 2500 2500 0 0 500 MHz 20000 14000 15000 10500 10000 5000 -0.08 -0.04 0 0.04 0.08 0.08 -0.04 0 0.04 X - position (m) X - position (ml 700 MHz F i g u r e V . 4 . B . i v . 6 . Electrical p r o p e r t y p r o f i l e a l o n g l o n g i t u d i n a l , transects through the taiqet parts of , respectively, images shown in Figuj e V . 4 .1-. region for the real and imaginary F'lots correspond to reconstructed 0.08 142 Re(k^) (m"2) Im(k2) (m'2) 800 2400 3000 3600 4200 4800 2400 2800 3200 3600 4000 6000 300 MHz 0 0 000 1000 2000 3000 4000 5000 6000 7000 8000 9000 1 0000 I 3000 4000 5000 "2000 6000 7000 8000 9000 10000 500 MHz 2000 4000 6000 8000 10000 2000 14000 0 J1400 12800 4200 5600 7000 8400 9800 16000 11200 18000 12600 20000 114000 700 MHz Figure V.4.B,iv.7. Real and imaginaiy parts of the k^ distribution of the reconstructed images for the 2.5 cm diameter fat/bone equivalent cylinder offset by 2.5 cm to the right of the target region center at 300, 500 and 700 MHz using both the inner and outer banks of measurement sites. 143 R E( k ^ ) (m ^ ) IM(k ^ ) 4000 6000 3000 4500 2000 3000 1000 1500 (m ^) 300 MHz 10000 10000 7500 7500 5000 5000 2500 2500 500 MHz 20000 14000 15000 10500 10000 7000 5000 3500 0 -0.08 -0.04 0 0.04 0.08 -0.08 X - position (m) -0.04 0 0.04 X - position (m) 700 MHz F i g u r e V .4.B. i v .8. Elefrtrical transects through property profile a l o n g l o n g i t u d i n a l , the target region parts of , respectively. i lots images shewn in Figure 1. :v . i. for the real and correspond to imaginary reconstructed 0.08 144 Re(k2) (m-2) Iin(k2) (m-2) 200 1600 2400 3000 3600 4200 4800 5400 2000 2400 2800 3200 3600 4000 6000 300 MHz 0 1000 1000 2000 3000 4000 6000 7000 8000 9000 0000 2000 4000 6000 8000 10000 12000 14000 16000 18000 20000 1 3000 - 24000 000 5000 ' 6000 7000 8000 9000 10000 500 MHz 0 ^1400 12800 4200 5600 7000 8400 19800 1200 112600 114000 700 MHz Figure V.4.B.iv.9. Real and imaginary parts of the k2 distribution of the reconstructed images for the 4.3 cm diameter fat/bone equivalent cylinder offset by 2.5 cm to the right of the target region center at 300, 500 and 700 MHz using both the inner and outer banks of measurement sites. 145 RE(k ^ ) (m ^) IM(k ^ ) 4000 6000 3000 4500 2000 3000 1000 1500 (m ) 300 MHz 10000 10000 7500 7500 5000 5000 2500 2500 500 MHz 20000 14000 15000 10500 10000 7000 5000 3500 -0.08 -0.04 0 0,04 0.08 -0.08 X - p o s i t i o n (rn) -0.04 0.04 X - position (m) 700 MHz Figure V.4.B.iv.lO. Electrical property profile along longitudinal, transects parts of through the target , respectively. region for the real and imaginary Plots correspond to reconstructed images shown in Figure V.4.B.iv.9. 0.08 146 Re(k2) (m"2) 2000 2400 3200 3600 4000 Im(k2) (m"2) 2400 3000 3600 4200 4800 5400 6000 300 MHz 2000 3000 4000 5000 4000 5000 6000 7000 8000 9000 10000 8000 9000 10000 500 MHz 0 2000 1400 4000 2800 6000 8000 4200 15600 7000 8400 9800 11200 10000 2000 14000 16000 18000 20000 12600 114000 700 MHz Figure V.4.B.iv.ll. Real and imaginary parts of the k^ distribution of the reconstructed images for the 2.5 cm diameter fat/bone equivalent cylinder offset by 3.8 cm to the right of the target region center at 300, 500 and 700 MHz using both the inner and outer banks of measurement sites. 147 RE(k ^) (m ^) IM(k ^) 4000 6000 3000 4500 2000 3000 1000 1500 (m ) 0 300 MHz 10000 10000 7500 7500 5000 5000 2500 2500 500 MHz 20000 14000 15000 10500 10000 7000 5000 3500 -0.08 ?0.04 0 0.04 0.08 -0.08 X - position (m) -0.04 0.04 X - position (m) 700 MHz F i g u r e V . 4 .B.iv. l 2 , Electrical property profile a l o n g l o n g i t u d i n a l , transects parts of through the target region for the real and imaginary , respectively. Plots correspond to reconstructed images shown In Figure , 1 . . i ?/. 1 I . 0.08 148 Re(k2) (m-2) 2000 2400 2800 3600 4000 Im(k^) (m'2) 1800 2400 3000 3600 4200 4800 6000 300 MHz 2000 3000 5000 6000 7000 8000 1 0000 3000 4000 5000 6000 7000 8000 9000 0000 500 MHz 0 2000 ; 1 400 4000 6000 |2800 10000 12000 14000 16000 4200 5600 7000 8400 9800 11200 18000 12600 20000 114000 700 MHz Figure V.4.B.iv.l3. Real and imaginary parts of the k^ distribution of the reconstructed images for the 2.5 cm diameter fat/bone equivalent cylinder offset by 3.8 cm to the right of the target region center at 300, 500 and 700 MHz using both the inner and outer banks of measurement sites. 149 IM(k ) (m 300 MHz 10000 10000 7500 7500 5000 5000 2500 2500 0 0 500 MHz 14000 10500 7000 3500 ?0.08 0.04 0 0.04 0 08 0.08 X - position (m) -0.04 0 0.04 X - position (m) 700 MHz Figure V.4.B.iv.14. Electrical property profile along longitudinal. transects region f o r t h e real and imaginary Plots correspond to reconstructed parts of through , the t arget respectively, images shown in F i g u r e Y. 1.1', i : 1 j. 0.08 CHAPTER VI. VI.1 CONCLUSION Overview The work presented in this thesis begins with a basic numerical algorithm, the Newton- Raphson iterative method, which is well suited for this class of non-linear problems. We have chosen the hybrid element method as its computational base because it yielded an accurate representation of the electric fields in a geometry which matched that needed for these imaging experiments. Issues associated with maintaining an accurate electric field representation while reducing the size of the reconstruction problem and the amount of measured data required were overcome through the introduction of a dual mesh scheme. A measurement system was developed to measure the electric fields in a saline tank for problems closely resembling a 2-D geometry. A calibration algorithm was realized to transform the measured 3-D data into 2-D for reconstruction purposes, and finally, tests were performed using a number of different phantoms to evaluate the capability of the system as a whole. In general this type of inverse problem is beset with difficulties due to its inherent nonuniqueness [Habashy et al., 1986]. The reconstruction technique developed here exhibits certain flexibilities which allow the implementation to utilize a priori information, place bounds on the intermediate electrical property distributions, apply a spatial filter on the material property distribution, and utilize various regularization procedures in an effort to stabilize the solution and limit the number of possible solutions in an effort to overcome that nonuniqueness. It also exhibits a certain parallelism, especially in the construction of the Jacobian matrix which can be exploited at a later date to further reduce the iteration time [Hua, et al., 1990]. The image reconstruction is performed over a small, defined region, but because of the unbounded nature of microwave radiation, accurate accounting of the fields radiating to infinity is required. As shown by Lynch, et al., [1986], this can be accomplished by implementing a 15 0 1 5 1 hybrid of the finite element method on the small defined region and the boundary element method for the region extending to infinity. In general, the number of reconstruction parameters is directly related to the amount of measured data that is available [Yorkey, et al., 1987; Guo, et al., 1989]. In microwave imaging, the amount of measured data simply becomes too great if the mesh used for the reconstruction parameters is the same as that used to represent the electric field distribution. The electric fields vary continuously over the entire region of interest due to their sinusoidal and exponential decay behavior. However, the material property distributions in biomedical applications are often nearly constant within various subregions, thus generally requiring far less discretization. A dual mesh scheme was introduced to decouple the mesh used for modeling the electric field distribution from that used for the material property mesh. This allowed for the fine discretization required in modeling the electric fields, and a much more coarse mesh for reconstructing the material property parameters. It also lends itself quite nicely to specialized coarse meshes that might be constructed based on a priori information. Simulations were performed comparing the dual mesh scheme with the single mesh scheme which show that reconstructed images were still quite good and achieved the desired effects of reducing the problem size and amount of measured data. A prototype microwave measurement system was developed to accurately detect the electric fields at various sites due to the radiation of water-filled waveguide radiators. It was designed around several criteria: (1) a large dynamic range given the high signal attenuation of biological tissue, (2) accurate detection of both phase and amplitude, (3) availability, cost and wideband performance of the microwave components, and (4) general modularity such that various antenna configurations could be tested without complete disassembly and reassembly. The system generally achieved these goals including a 1 IS dB dynamic range and well matched impedance characteristics over the frequency range of 300 to 1100 MHz. 152 Using measured data from the data acquisition system, a calibration procedure was implemented which transformed the measured data from the 3-D radiating pattern to the 2-D plane. Incorporation of the transformation was accomplished using a least squares technique to determine the phase center of the radiators, which was subsequently used to subtract out the requisite "free space loss factors." Agreement between the measured data modified by the calibration technique and that computed from the simulations was on the order of 0.3 dB in amplitude and 3" in phase. Such agreement is an essential component for using the NewtonRaphson method. Finally images were successfully reconstructed for a variety of targets using this system while operating over the range of 300 to 700 MHz. These showed general improvement as the frequency increased; in particular, detection of the real parts of the DI and agar gel phantoms were not accurately accomplished at 300 MHz but were recovered quite well at 700 MHz. This may be due in part to the ability of the system to detect the presence of plexiglas^^ cylinder around these objects at 700 MHz while not at 300 MHz, and hence, the exhibited smearing of the plexiglas? properties with those inside the cylinder. Otherwise, the targets were generally reconstructed quite well with some advantage going to the higher frequency in terms of background smoothness and delineation of material interfaces. VI.2 Future Directions These findings are the first reported implementation of a microwave imaging system for biomedical applications which has successfully reconstructed quantitative images of various biological equivalent tissue and thus, bodes well for the future implementation in a clinical setting. However, two major hurdles remain before clinical implementation is realized. First, dynamic images must be achieved. Such a system will either involve: (1) only few antennas 153 mounted on a motorized platform to achieve multiple views of the target region, or (2) a large number of fixed antennas which can be electronically selected. In various contexts, these both have advantages. The second hurdle is to implement a microwave imaging system without the targets and antennas being completely submerged in a water tank. To date, the results of the system described herein and that of other systems by Larsen, et al. [1984], Jofre, et al. [1990], and Miyakawa [1993] have all been in water tanks. There are several ways to simulate the desirable electrical properties of water and saline through the use of high dielectric ceramics. We plan to implement such an arrangement and continue on the lines of 2-D imaging which can be transported to a clinical setting for treatment monitoring. Future work will also need to be done on reducing the imaging time. The most significant contributions here will most likely involve utilizing the inherent parallelism of this algorithm. At present, hundreds of matrix back substitutions are performed using the same decomposed matrix. Computing these tasks in parallel will significantly reduce the computation time, but it will inevitably increase the system cost due to new levels of computational sophistication. Beyond that, the ideal goal would be to achieve 3-D images. This poses a substantially more difficult problems than dealt with here. Full characterization of the radiating fields in 3-D will be required. In addition, receiving antennas will need to detect all three spatial components of the electric fields, whereas, here we have only needed one component, assuming for practical purposes that our system is operating in the 2-D TM mode. It will also pose a much more substantial computational problem with all the associated issues of increased illconditioning. 154 At present, the envisioned 2-D imaging system would be a valuable asset in a clinical setting. In particular, using this for monitoring the temperature distributions for assessing the success of hyperthermic treatments and possible real-time feedback control will be necessary for the future success of hyperthermia. REFERENCES Balanis CA, Antenna Theory: Analysis and Design. New York, NY, Harper and Row Publishers, pp.64-67, 1982. Barber DC, Brown BH, Avis NJ, "Image reconstruction in electrical impedance tomography using filtered back-projection," Annual International Conference of the IEEE Engineering in Medicine and Biology Society, 1691-1692, Paris, 1992. Barber DC, Brown BH, "Reconstruction of impedance images using filtered back-projection," European Community Concerted Action in Electrical Impedance Tomography, pp. 1-7, Copenhagen, CAIT, 1990. Blad B, Perrson B, Lindstrom K, "Quantitative assessment of impedance tomography for temperature measurements in hyperthermia," International Journal of Hyperthermia, vol 8, pp. 33-43, 1992. Bolomey JC, Hawley MS, (Chapter 2) Non-invasive control of hyperthermia. Methods of Hyperthermia Control, Springer-Verlag, Berlin, Heidelberg, 1990. Bolomey JC, "Recent European developments in active microwave imaging for industrial, scientific, and medical applications," IEEE Transactions on Microwave Theory and Techniques, vol 37, pp. 2109-2117, 1989. Bolomey JC, Jofre L, Peronnet G, "On the possible use of microwave-active imaging for remote thermal sensing," IEEE Transactions on Microwave Theory and Techniques, vol 31, pp. 777-781, 1983. Bolomey JC, Izadnegahdar, Jofre L, Pichot CH, Peronnet G, Solaimani M, "Microwave diffraction tomography for biomedical applications," IEEE Transactions on Microwave Theory and Techniques, vol 30, pp. 1998-2000, 1982. Borup DT, Sullivan DM, Gandhi OP, "Comparison of the FFT conjugate gradient method and the finite-difference time-domain method for the 2-D absorption problem," IEEE Transactions on Microwave Theory and Techniques, vol 35, pp. 383-395, 1987. Bottomley PA, Andrew ER, "RF magnetic filed penetration, phase shift and power dissipation in biological tissue: implications for NMR imaging," Phys. Med. Biol, vol 23, pp. 630-643, 1978. Broquetas A, Romeu J, Rius JM, Elias-Fuste AR, Cardama A, Jofre L, "Cylindrical geometry: a further step in active microwave tomography," IEEE Transactions on Microwave Theory and Techniques, vol 39, pp. 836-844, 1991. Brown BH, Sinton AM, Barber DC, Leathard AD, McArdle FJ, "Simultaneous display of lung ventilation and perfusion on a real-time EIT system," Annual International Conference of the IEEE Engineering in Medicine and Biology, pp. 1710-1711, Paris, 1992. Burdette EC, Friedrich PG, Seaman RL, Larsen LB, "In situ permittivity of canine brain: regional variations and post-mortem changes," IEEE Transactions on Microwave Theory and Techniques, vol 34, pp. 38-50, 1986. 155 156 Caorsi S, Gragnani GL, Pastorino M, "A multiview microwave imaging system for twodimensional penetrable objects," IEEE Transactions on Microwave Theory and Techniques, vol 39, pp. 845-51, 1991. Caorsi S, Gragnani GL, Pastorino M, "Two dimensional microwave imaging by a numerical inverse scattering solution," IEEE Transactions on Microwave Theory and Techniques, vol 38, pp. 981-989, 1990. Chommeloux L, Pichot C, Bolomey JC, "Electromagnetic modeling for microwave imaging purposes," IEEE Transactions on Microwave Theory and Techniques, vol 34, pp. 1064-76, 1986. Clegg ST, Roemer RB, "Towards the estimation of three-dimensional temperature fields from noisy data during hyperthermia," International Journal of Hyperthermia, vol 5, pp. 918-924, 1989. Conway J, "Electrical impedance tomography for thermal monitoring of hyperthermic treatment: an assessment using in vitro and in vivo measurements," Clinical Physics and Physiological Measurements, vol 8 (Supplement A), pp. 141-146, 1987. Copeland JR, "Radar target classification by polarization properties," Proceedings of the IRE, vol 48, pp. 1290-1296, 1960. Curran WJ Jr., Goodman RL, "Hyperthermia 1991: a critical review," Radiation research - A twentieth-century perspective, Congress Proceedings, vol 2, pp. 883-888, 1991. Demoment G, "Image reconstruction and restoration: overview of common estimation structures and problems," IEEE Transactions on Acoustics, Speech, and Signal Processing, vol 37, pp. 2024-2036, 1989. Dennis JE, Schnabel R, Numerical Methods for Unconstrained Optimization and Non Linear Equations, Englewood Cliffs, NJ, Prentice-Hall, 1983. Dewhirst MW, Moon T, Carlin D, "Analysis of tumor volume and thermal dosimetric effects on tumor response to heat, radiation and heat plus radiation: results of phase III randomized clinical trial in pet animals," Physical Aspects of Hyperthermia. Nussbaum GH (ed.), (Medical Physics mono. no. 8), pp. 495-510, 1982. Dyson JD, "Determination of the phase center and phase patterns of antennas," Radio Antennas for Aircraft and Aerospace Vehicles, W.T. Blackband (ed.), AGARD Conference Proceedings, No. 15, Slough, England, Technivision Services, 1967. Edelstein-Keshet L, Dewhirst MW, Oleson JR, Samulski TV, "Characterization of tumor temperature distributions in hyperthermia based on assumed mathematical forms," International Journal of Hyperthermia, vol 5, pp. 757-777, 1989. Golub GH, van Loan CF, Matrix Computations. The Johns Hopkins University Press, Baltimore, MD, second edition, pp. 193-259, 1989. Gregg EC, "Radiation risks with diagnostic x-rays," Radiology, vol 123, pp. 447-453, 1977. Guo TC, Guo WW, "Three-dimensional dielectric imaging by inverse scattering with resolution unlimited by wavelength," Proceedings of the Conference on Electrical Insulation and Dielectric Phenomena, pp. 65-74, 1989. 157 Guo TC, Guo WW, "Physics of image formation by microwave scattering," SPIE Medical Imaging, vol 767, pp. 30-39, 1987. Guy AW, "Analysis of electromagnetic fields induced in biological tissues by thermographic studies and equivalent phantom models," IEEE Transactions on Microwave Theory and Techniques, vol 19, pp.205-213, 1971. Habashy TM, Chew WC, Chow EY, "Simultaneous reconstruction of permittivity and conductivity profiles in a radially inhomogeneous slab," Radio Science, vol 21, Number 4, pp. 635-645, 1986. Harrington RF, Field Computation by Moment Methods, Malabar, FL, Krieger Publishing Company, 1982. Henderson RP, Webster JG, "An impedance camera for spatially specific measurements of Thorax," IEEE Transactions on Biomedical Engineering, vol 25, pp. 250-254, 1978. Hines WW, Montgomery DC, Probability and Statistics in Engineering and Management Science. New York, NY, John Wiley and Sons, pp.455-558, 1990. Hua P, Woo EJ, "Reconstruction algorithms," Electrical Impedance Tomography, JG Webster (ed), Bristol, England, Adam Higler, Ch. 10, pp. 97-137, 1990. Huynen RJ, Phenomenological Theorv of Radar Targets. Drukkerij Bronder, 1970. Hynynen K, Shimm D, Anhalt D, Stea B, Sykes H, Cassady JR, Roemer RB, "Temperature distributions during clinical scanned, focused ultrasound hyperthermia treatments," International JoumcU of Hyperthermia, vol 6, pp. 891-908, 1990. Isaacson D, "Distinguishabilities of conductivities by electric current computed tomography," IEEE Transactions on Medical Imaging, vol 5, pp. 91-95, 1986. Joachimowicz N, Pichot C, Hugonin JR, "Inverse scattering: an iterative numerical method for electromagnetic imaging," IEEE Transactions on Antennas and Propagation, vol 39, pp. 174252, 1991. Joachimowicz N, Pichot C, "Comparison of three integral formulations for the 2-D TE scattering problem," IEEE Transactions on Microwave Theory and Techniques, vol 38, pp. 178-85, 1990. Jofre L, Hawley MS, Broquetas A, de los Reyes E, Ferrando M, Elias-Fuste AR, "Medical imaging with a microwave tomographic scanner," IEEE Transactions on Biomedical Engineering, vol 37, pp. 303-312, 1990. Kallman JS, Berryman JG, "Weighted least-squares criteria for electrical impedance tomography," IEEE Transactions on Medical Imaging, vol 11, pp.284-292, 1992. Krause JD, Antennas. McGraw-Hill, New York, NY, pp. 276-78, 1950. Lagendijk JJW, "Heat transfers in tissues," Phvsics and Technology of Hyperthermia. Field SB and Franconi C (eds.), Martinus Nijhoff, The Netherlands, pp. 517-552, 1987. 158 Lagendijk JJW, Schellekens M, Schipper J, van der Linden PM, "A three-dimensional description of heating patterns in vascularized tissues during hyperthermic treatment," Physics in Medicine and Biology, vol 29, pp. 495-507, 1984. Larsen LE, Jacobi JH, "Methods of active microwave imagery for dosimetric applications," Medical Applications of Microwave Imaging. Larsen LE and Jacobi JH (eds.) IEEE Press, pp. 118-137, 1986. Lebihan D, Delannoy J, Levin RL, "Temperature mapping with MR imaging of molecular diffusion: applications to hyperthermia," Radiology, vol 171, pp. 853-857, 1989. Liao SY, Microwave Devices and Circuits. Englewood Cliffs, NJ, Prentice-Hall, Inc., pp. 1863, 1980. Lynch DR, Paulsen KD, Strohbehn JW, "Hybrid element method for unbounded electromagnetic problems in hyperthermia," International Journal for Numerical Methods in Engineering, vol 23, pp. 1915-1937, 1986. Lynch DR, Paulsen KD, Strohbehn JW, "Finite element solution of Maxwell's equations for hyperthermic treatment planning," Journal of Computational Physics, vol 58(2), pp. 246-9, 1985. MacFall J, Prescott DM, Fullar E, Samulski TV, "Temperature dependence of canine brain tissue diffusion coefficient measured in vivo with magnetic resonance echo-planar imaging," International Journal of Hyperthermia, vol 11, pp. 73-86, 1995. Mallorqui JJ, Broquetas A, Jofre L, Cardama A, "Noninvasive active microwave thermometry with a microwave tomographic scanner in hyperthermia treatments," ACES Special Issue on Bioelectromagnetic Computations, AHJ Fleming and KH Joyner (Eds.), ACES Journal, vol 7, pp.121-127, 1992. Marquardt DW, "An algorithm for least-squares estimation of non-linear parameters," SIAM Journal of Applied Mathematics, vol 11, pp. 431-41, 1963. McArdle FJ, Brown B, Angel A, "Imaging cardiosynchronous impedance changes in the adult head," EC Concerted Action on EIT, London, CAIT, 1992. McRae DA, Esrick MA, "Changes in electrical impedance of skeletal muscle measured during hyperthermia," International Journal of Hyperthermia, vol 9, pp. 247-262, 1993. Meaney PM, Paulsen KD, Ryan TP, "Two-dimensional hybrid element image reconstruction for TM illumination," IEEE Transactions on Antennas and Propagation, vol 43, pp. 239-247, 1995. Meaney PM, Paulsen KD, Hartov A, Crane RK, "An active microwave imaging system for reconstruction of 2-D electrical property distributions," IEEE Transactions on Biomedical Engineering, (In press (a)). Meaney PM, Paulsen KD, Hartov A, Crane RK, "Multi-target microwave imaging for tissue assessment: initial evaluation in tissue equivalent phantoms," IEEE Transactions on Biomedical Engineering, (Submitted (b)). 159 Miyakawa M, "Tomographic measurement of temperature change in phantoms of the human body by chirp radar-type microwave computed tomography," Medical and Biological Engineering and Computing, vol 31, pp. S31-S36, 1993. Moros EG, Button AW, Roemer RB, Burton M, Hynynen K, "Experimental evaluation of two simple thermal models using hyperthermia in muscle in vivo," International Journal of Hyperthermia, vol 9, pp. 581-598, 1993. Moskowitz MJ, Ryan TP, Paulsen KD, Mitchell SE, "Clinical implementation of electrical impedance tomography with hyperthermia," International Journal of Hyperthermia, vol 11, pp. 141-149, 1995. Moskowitz MJ, "Electrical Impedance imaging for temperature field estimation during hyperthermia treatment of cancer," Ph.D. Dissertation, Thayer School of Engineering, Dartmouth College, pp. 151-184, 1994. Ney MM, "Method of moments as applied to electromagnetic problems," IEEE Transactions on Microwave Theory and Techniques, vol 33, pp. 972-980, 1985. Ohtera I, Ujiie H, "Nomographs for phase centers of connical corrugated and TE; i mode horns," IEEE Transactions on Antennas and Propagation, vol 23, pp. 858-9, 1975. Oleson JR, Dewhirst MW, Harrelson JM, Leopold KA, Samulski TV, Tso CY, "Tumor temperature distributions predict hyperthermia effect," International Journal of Radiation Oncology Biology and Physics, vol 16, pp. 559-570. Overgaard J, Overgaard M, "Hyperthermia as an adjuvant to radiotherapy in the treatment of malignant melanoma," International Journal of Hyperthermia, vol 3, pp. 483-501, 1987. Paulsen KD, Meaney PM, Moskowitz MJ, "A dual mesh scheme for finite element based reconstruction algorithms," IEEE Transactions on Medical Imaging, (In Press (a)). Paulsen KD, Moskowitz, MJ, Ryan TP, "Temperature field estimation using electrical impedance profiling: I. Reconstruction algorithm and simulated results," International Journal of Hyperthermia, vol 10, pp. 209-228, 1994. Paulsen KD, Lynch DR, Liu W, "Conjugate direction methods for helmholtz problems with complex-valued wavenumbers," International Journal for Numerical Methods in Engineering, vol 35, pp. 601-22, 1992a. Paulsen KD, Liu W, "Memory and operations count scaling for coupled finite-element and boundary-element systems of equations," International Journal for Numerical Methods in Engineering, vol 33, pp. 1289-1304, 1992b. Paulsen KD, Lynch DR, "Calculation of interior values by the boundary element method," Communications in Applied Numerical Methods, vol 5, pp. 7-14, 1989. Paulsen KD, Lynch DR, Strohbehn JW, "Three-dimensional finite, boundary, and hybrid element solutions of the Maxwell equations for lossy dielectric media," IEEE Transactions on Microwave Theory and Techniques, vol 36, pp. 682-693, 1988. 160 Paulsen KD, Strohbehn JW, Lynch DR, "Comparative theoretical performance of two types of regional hyperthermia systems," International Journal of Radiation Oncology Biology and Physics, vol 11, pp. 1659-1671, 1985. Pennes HH, "Analysis of tissue and arterial blood temperatures in the resting human forearm," Journal of Applied Physiology, vol 1, pp. 93-122, 1948. Perez CA, Pajak T, Emami B, Hombeck NB, Tupchong L, Rubin P, "Randomized phase III study comparing irradiation and hyperthermia with irradiation alone in superficial measurable tumors," American Journal of Clinical Oncology, vol 14, pp. 133-141, 1991. Perl W, Hirsch RL, "Local blood flow in kidney tissue by heat clearance measurement," Journal of Theoretical Biology, vol 10, pp. 180-251, 1966. Rius JM, Pichot C, Jofre L, Bolomey JC, Joachimowicz N, Broquetas A, Ferrando M, "Planar and cylindrical active microwave temperature imaging: numerical simulations," IEEE Transactions on Medical Imaging, vol II, pp. 457-469, 1992. Roemer RB, "Optimal power deposition in hyperthermia. I. The treatment goal: the ideal temperature distribution: the role of large blood vessels," International Journal of Hyperthermia, vol 7, pp. 317-341, 1991. Roemer RB, "Thermal dosimetry," Thermal Dosimetry and Treatment Planning. Gautherie M (ed.), Springer-Verlag, Berlin, pp. 119-208, 1990. Samulski TV, MacFall J, Zhang Y, Grant W, Charles C, "Non-invasive thermometry using magnetic resonance diffusion imaging: potential for application in hyperthermic oncology," International Journal of Hyperthermia, vol 8, pp. 819-830, 1992. Schwan HP, Li K, "Capacity and Conductivity of Body Tissues at Ultrahigh Frequencies," Proceedings of IRE, vol 41, pp. 1735-1740, 1953. Sinclair G, "The transmission and reception of elliptically polarized waves," Proceedings of IEEE, vol 39, pp. 535-540, 1951. Skolnik MI, Introduction to Radar Systems. New York, NY, McGraw-Hill Brook Company, 1980. Slaney M, Kak AC, Larsen LE, "Limitations of imaging with first-order diffraction tomography," IEEE Transactions on Microwave Theory and Techniques, vol 32, pp. 860-873, 1984. Strohbehn JW, Paulsen KD, Lynch DR, "Use of the finite element method in computerized thermal dosimetry," Physical techniques in clinical hyperthermia. Hand JW and James JR (ed.). Research Studies Press, Letchworth, England, pp. 383-451, 1986. Tikhonov AN, Arsenin VY, Solutions of Ill-Posed Problems, Washington, DC, V.H. Winston & Sons, 1977. Valdagni R, Amichetti M, Pani G, "Radical radiation alone versus radical radiation plus microwave hyperthermia for Ng (TNM-UICC) neck nodes: a prospective randomized clinical trial," International Journal of Radiation Oncology, Biology, Physics, vol 15, pp. 13-24, 1988. 161 Weinbaum S, Jiji LM, "A new simplified bioheat equation for the effect of blood flow on local average tissue temperature," Journal of Biomechanical Engineering, vol 107, pp. 131-139, 1985. Woo EJ, Hua P, Webster JG, Tompkins WJ, "Finite-element method in electrical impedance tomography," Medical & Biological Engineering and Computing, vol 32(5), pp. 530-536, 1994. Yashiro K, and Ohkawa S, "Boundary element method for electromagnetic scattering from cylinders," IEEE Transactions on Antennas and Propagation, vol 33(4), pp. 383-9, 1985. Yorkey TJ, Webster JG, and Tomkins WJ, "Comparing reconstruction algorithms for electrical impedance tomography," IEEE Transactions on Biomedical Engineering, vol 34, pp. 843-52, 1987. Yorkey T, "Comparing reconstruction methods for electrical impedance tomography," Ph.D. Dissertation, University of Wisconsin-Madison, 1986. APPENDIX A. Forward Solution of the Electric Field Distribution Using the Hybrid Element Method Starting with the system of equations (equations II.2.4-5), the finite element region > < ?1 o .^bl 1 ^11 CP produces a matrix relationship which can be decomposed as: J E . L "0 lEbJ 0 Bbb. (A.l) pj p [FsJ "Dbb 1 O DJ Csb u "Cbb D.b P 1 and the boundary element region yields a similarly written matrix form: fE. (/\.2) where the subscripts I, b, and s refer to all points interior to the FE surface, all points on the boundary shared by the FE and BE problems and all points on the radiating source(s), respectively. In equation (A.l), the sparse matrix A is partitioned into four submatrices of various dimensions, A?, A,b, A^, and Ayy where, for example, submatrix A,b has dimensions n, x n^ where n, , n^, and n, are the number of nodes in regions I, b, and s, respectively. Similar partitioning is used in column vectors E and F where F=^. This on notation also appears in the expression of matrix system (A.2). For convenience in some instances, we express the partitioned matrices in (A.l) and (A.2) in a nonpartitioned form in which case the subscripts are dropped in favor of a single alphabetical symbol. At this point, there are three methods of solving these two sets of simultaneous equations. For purposes of the appendix, only hybrid element methods 1 and 2 (HE #1 and HE #2) will be described here, and the reader is referred to Paulsen, et al. [1990] for a discussion of HE #3. A.l Hybrid Element Method #1 Using HE #1 technique, we first solve for equation (A.2): 62 > in terms of 4 E^ E. ^ by inverting [C] in 163 "Cbb Cbs" Csb Css. -1 "Dbb Dbs" mi "Gbb Gbs' El Dsb Dss. l E s J Gsb Gss. Es (A.3) solving this for {Fy} and multiplying by matrix [B] from equation (A.l) yields: (A.4) [B1{FJ = (BlGM,]{Eb} + [Bl[Ota]{E,} This expression can be substituted into equation (A.l) to give: E, Mb .^bl ^bb ~[B][Gbb]. (A.5) EbJ 1[B1[06.1{E,}J which will allow for the computation of the electric field at all the nodes within the finite element region and on its boundary given the radiation distribution on all exterior sources. A.2 Hybrid Element Method #2 The HE #2 approach is again based on equations (A.l) and (A.2), except that we first fEil solve for jg| i" equation (A.l) by inverting matrix [A]: AI? AI Albi AI| I:: (A.6) where [AI] = [A] '. Equation (A.6) provides a relationship between {Ej,} and {Fj,} alone: (A.7) {Eb} = [AIi,b][B]{Fb} which can be substituted into equation (A.2) to yield: u mi "Dbb 1" Q Csb Css [ F s J Dsb I J f[AIbbl[BKFb}l_ Hbb ' DJ Cbb (Ej r Hsb H bs E. (A.8) 164 where Hbb H(,s _ D(,b D|,s [AIbb][B] Hsb ^ss. ^ss. .^sb and [Ijs] is the identity matrix of rank s. [iss] Equation A.8 can be rearranged to produce: (A.9) Equation (A.9) can be solved for {Fy} which can then be substituted back into equation (A. 1) to recover the electric fields everywhere within the finite element mesh. When considering only the forward solution, it is important to note that HE #2 has slight computational advantages over HE #1 [Paulsen, et al., 1990]. In He #2, the initial step is to invert matrix [A] which involves decomposing a sparse system of equations and then back substituting various unit vectors to recover the necessary columns of [Al]. In this case, only [AIj,b] is needed, which implies only n^ back substitutions are required. For HE #1, the initial step is to invert matrix [C] which involves decomposing a full system of equations followed by n^ + n, back substitutions. Hence, in general, there is a slight computational advantage associated with HE #2, primarily due to its reduced number of back substitutions. APPENDIX B. Algorithmic Computational Time Enhancements There are numerous methods for reducing overall computation time of any computer program, from purchasing faster computers to using computer specific optimized functions to optimizing various features of the algorithm itself based on the problem at hand. This section will adopt the latter approach in optimizing two tasks: (1) reducing the cost per iteration resulting from calculating the boundary element contribution to the forward solution, and (2) reducing the cost per iteration associated with calculating various indices and quantities required for the coordination of the electric and B.l meshes. Boundary Element Contribution to Iteration Time As mentioned in Chapter II, hybrid element method #1 (HE #1) offers significant computational advantages over HE #2 when used as the basis of this reconstruction algorithm. As described here, it also offers significant advantages over HE #2 by further reducing the costs per iteration by computing matrices associated the boundary element contribution once in a preprocessor routine. In using any of the hybrid element methods, the matrices associated with the boundary element discretization are solely based on the geometry of the mesh in the boundary element region and the material properties of that region. Thus, their computation will not change with iteration. This means that these matrices can be computed once in a pre� processor routine and stored for later use at each iteration. Using HE #1, we can take this one step further: not only can we compute matrices [C] and [D] of equation (II.2.5) in a pre� processor routine, but we can also invert [C] and multiply its inverse with [D] in that same pre-processor and only store their product, [G] (equation II.2.8). This means that using HE #1, there are no full matrix inversions computed during each iteration. On the other hand, if HE #2 were to be used we would still have to perform a full inversion of matrix [A], because its entries are functions of the material properties in the finite element region which change with 165 166 iteration. The algorithm would also have to store both [C] and [D] which would be nontrivial since both are full complex matrices. B.2 Dual Mesh Computational Overhead and its Contribution to Iteration Time It is understood that when using this imaging system, that the target region, measurement sites and radiator positions are all pre determined. So, too, are the electric mesh and its associated mesh. This is important because, as described in Chapter III and Appendix C, introduction of the dual mesh method has significant computational ramifications related to the coordination of the two meshes. However, these relationships are strictly functions of the geometry of the two meshes and not of the various material properties. Thus, information can be computed prior to execution of the reconstruction algorithm which significantly reduces the computational costs associated with each iteration. APPENDIX C. Dual Mesh Overhead As described in Chapter III, the dual mesh method requires significant overhead in the form of various indices and vectors of basis function values in order to make its implementation efficient. The implications in terms of the computational costs of the method are felt primarily in computing (1) the forward solution, and (2) the Jacobian matrix used for determining the material property updates at each iteration. Figure C.l illustrates the problem for the forward solution: k element electric element Figure C. 1. Diagram of relationship between solution. 167 mesh and electric mesh during the forward 168 In calculating the matrix, [A] (described in equation 11.2.4), integrations of the active basis functions are performed over each electric element. To accomplish this, the material properties at nodes Ei, E2, and Eg must be calculated. In general, however, these values are only known at the nodes of the k^ mesh. To compute these material property values, several pieces of information need to be calculated beforehand. These are; (1) NIK(i) - This arrayt stores the k^ element number which contains the i^h electric node number. (2) NIN(i,l-3) - This array stores the three k^ node numbers associated with the k^ element which contains the i^^ electric node. (3) PHINK(i,l-3) - This array stores the three basis function values associated with the three k^ nodes in the NIN array for the i^h electric node. Given this information and the material property values at the k^ nodes, the inner products required in forming [A] can be performed. In calculating the Jacobian matrix, [J] (described in equation II.2.2), determining the appropriate matrix contributions is a bit more involved. Figure C.2 illustrates one configuration for this situation. Essentially, the elements of [J] are the derivatives of the elements of [A] with respect to the k^ values at the k^ nodes. These derivatives involve integrations which are performed over the electric elements. To accomplish this, the algorithm (1) loops over all of the k^ nodes, (2) finds all of the electric elements for which the basis function centered at the current k^ node is non-zero, (3) calculates the basis function value (whose center is located at the current k^ node) at each of the electric element nodes, and (4) computes the required integrations over the electric elements using this information. To accomplish these tasks efficiently, several more pieces of information need to be pre-computed and stored including: 169 (1) ISKE(j,l-NEAMAX) - This array stores the electric element numbers associated with the node. Since there is no set number of electric elements associated with each node, a maximum value is set to limit the array size. (2) lSKN(j,NNAMAX) - This array stores the electric node numbers associated with the k^ node number. A maximum number of electric nodes associated with each k^ node is set. (3) PHIK(j,I-NNAMAX) - This array stores the basis function value at the electric nodes given in ISKN due to the k^ basis function centered at the k^ node, (|>j. Given this information, the inner products required in forming [J] can be performed. k node Electric Element k elements associated with |*^k node Figure C.2. Diagram of relationship between k^ mesh and electric mesh during calculation of the Jacobian matrix, [J]. APPENDIX D. Noise Figure Analysis As mentioned in Chapter IV, one of the keys to achieving a large dynamic range is to improve the receiver signal-to-noise ratio. The most widely recognized figure of merit for this is the noise figure, Fn. By definition, the noise figure can be written: /'^out where S is the signal power, N is the noise power and the subscripts in and out refer to the input and outputs of the receiver channel, k is Boltzmann's constant, Tq is the receiver temperature in " Kelvin, Bn is the noise bandwidth and G is the gain of the receiver, where S / G= . Fn in its simplest terms is an indication of how much the signal-to-noise ratio / ^in has been degraded as a result of the signal having passed through the cascade of receiver components. In general, the noise figure of individual components is readily available, and from this data a composite noise figure can be determined for the whole receiver. Figure D. 1 shows the cascade of components through which a signal passes in one of the system channels. (Note that the gain and noise figure for the digital attenuators are for base attenuation levels.) The total noise figure for the system, \ o , can be written, + P k -_l . ^2zl , + llzl+_+ G, G,G2 G,G2...Gk_, (D.2) where k is the total number of components in the cascade. From equation D.2, it makes sense to have a low noise amplifier, LNA, as the first component in the cascade which will have the effect of keeping the first term, F i, small and will reduce all subsequent terms as long as the gain products are greater than 1. For the system designed the cascaded noise figure is, F^io^i =2.8 dB. This is quite good given the wide band, 3(X) - 1100 MHz, of operation. 170 171 LNA Switch Multiplexer Amplifier Digital Attenuator F. = 6dB 4 Fg= 6 dB In R, = 2.4 dB = +15 dB Amplifier F^=1 dB G, = -1 dB Digital Attenuator F3= 1 dB = -1 dB = +20 dB Amplifier Mixer Fg=4dB F^=6dB Fo=4dB F?=7dB C^ = +15dB (^ = -6 dB (^ = +15dB C^ = -6dB Figure D. 1 Diagram of receiver channel with associated noise figures and gains. % == -6?' dB APPENDIX E. Water-filled Waveguide Return Losss As was shown in Chapter IV, the measurement system utilizes the highly lossy nature of the 0.9% saline medium to improve the characteristic impedance of the monopole receiving antennas to such a degree that they are usable over a broad frequency range. In contrast, the water-filled transmitting waveguide generally has a broad usable bandwidth relative to that of the receiving monopole antenna. As seen in Figure E.1, the return loss is generally better than -15 dB from 400 to 900 MHz when the waveguide is immersed in DI water. Characteristic Impedance Water-filled Waveguide -10 CD ?o -20 C k- 3 -30 CC -40 DI Water Saline Solution (0.9%) -50 L- 100 300 500 700 900 1100 Frequency (MHz) Figure E.1. Characteristic impedance of water-filled waveguide antenna when immersed in DI water and a 0.9% saline solution, respectively. 1 72 173 Its characteristic impedance is further improved when immersed in saline solution (Figure E.1). in which case it has a usable bandwidth from 300 to 900 MHz, with return losses of better than -10 dB. The curve is also dramatically smoother over the full bandwidth, indicating that such things as standing waves and normally unterminated currents are significantly attenuated by the highly conductive medium. Thus, from the available data, this antenna's characteristics appear well suited for this application. APPENDIX F. Radiation Field Variation Outside of the 2-D Plane The calibration procedures described in Chapter IV assume that the main beams of both the transmitting and receiving antennas are centered in the target plane. Positioning either of these antennas above or below such a level could contribute significant amplitude and phase errors and subsequently distort the reconstructed images. Given the present antenna configuration, it is difficult to measure a true radiation pattern by taking measurements in an arc in the far field. However, it is possible to measure the variations in the received signals txansmitted from these antennas due to vertical position variations. Transmitted l\1agnitude vs. Position Water-filled Waveguide and Monopole Antennas m TJ -2 0) "D 5 -4 C O) (0 > 300 MHz (Waveguide) -6 500 MHz (Waveguide) 700 MHz (Waveguide) 0> oc 300 MHz (Monopole) 500 MHz (Monopole) 700 MHz (Monopole) -10 -3.5 -2.5 -1.5 -0.5 0.5 1.5 2.5 3.5 Position (cm) Figure F.l. Variations in the magnitude of the transmitted power for the water-filled waveguide and monopole antennas as a function of relative vertical position at 300, 500 and 700 MHz. 174 175 In this case, the antennas under test are either the water-filled waveguide or the monopole antenna, while the probe antenna is a specially modified monopole antenna whose vertical position can be changed. The antenna under test and the probe antenna are placed at set distances from each other and the transmitted signals are recorded for the various probe antenna vertical settings. Figure F.l shows the variations in the measured magnitude with probe antenna position for the waveguide and monopole positions. (The position on the abscissa is relative to the point where the center of the probe antenna is aligned with the center of the antenna under test.) Figure F.2 shows a similar plot for the phase variations for both types of antennas. Transmitted Phase vs. Position Water-filled Waveguide and Monopole Antennas -20 O) -40 300 MHz (Waveguide) 500 MHz (Waveguide) 700 MHz (Waveguide) -60 300 MHz (Monopole) 500 MHz (Monopole) 700 MHz (Monopole) -80 -3.5 -2.5 1.5 -0.5 0.5 1.5 2.5 3.5 Position (cm) Figure F.2. Variations in the phase of the transmitted power for the water-filled waveguide and monopole antennas as a function of relative vertical position at 300, 500 and 700 MHz. 176 In general, there appears to be a region near the zero relative position where the phase and amplitude are quite flat for the frequency range of 300 to 700 MHz. This suggests that we are not operating these antennas in a highly variable region of their respective antenna patterns. Thus, if care is taken to use constant height settings, vertical positioning should not contribute significantly to the overall measurement error. APPENDIX G. Calibration of Variations Between Transmitter Channels and Receiver Channels The first stage of the calibration process is to remove differences in the phase and amplitude of all four transmit and all four receiver antennas. We accomplish this by using a least squares technique. Figure G. 1 illustrates the antenna configuration for our calibration scheme 0 - - z '4 V \ o / \ I \ / ' / \ l / Figure G.l. Antenna arrangement for calibration of amplitude and phase offset differences between antennas. where T1-4 correspond to the four transmitters and R1.4 correspond to the four receivers. Lj is the distance from a given transmit antenna to its closest receiver antenna and L] is the 177 178 distance from that same transmit antenna to the two equidistant receiver antennas. The data acquisition system is set up such that each transmit and receive channel has an arbitrary amplitude and phase offset with respect to an amplitude of 0 dBm and a phase of 0 degrees associated with it. For transmitters 1-4 these are At(1-4) and Pt(1-4). respectively, and for the receivers are Ar(1-4) and Pr(1.4). The amplitude and phase offsets due to the signals travelling lengths, Li and L2, are Ai, A2 and Pi, P2, respectively. Using this notation, equations for the amplitude and phase can be written for the signals detected at receiver, i, due to a signal sent from transmitter, j, having traveled a length Lk: ^i.j - ^Ri + ^Tj + (G-1) Pi.j = ^Ri + ^Tj + Pk + Pg (G-2) where Ag and Pg are the amplitude and phase error terms associated with each measurement. Note that for each transmitter, we are only using three of the four receiver measurements sites possible. This is mainly due to the fact that measurements at the farthest receiver site are partially disturbed by the presence of the nearest receiver antenna. (From experience with this measurement system, it is known that the measurements at a selected receiver antenna will be slightly perturbed if

1/--страниц