# Fundamental, formational, and practical challenges in microwave inverse scattering for medical imaging of the breast

код для вставкиСкачатьNOTE TO USERS This reproduction is the best copy available. UMI* F U N D A M E N T A L , FORMATIONAL, A N D P R A C T I C A L CHALLENGES IN MICROWAVE I N V E R S E S C A T T E R I N G FOR MEDICAL I M A G I N G OF THE BREAST by Jacob Dan Shea A dissertation submitted in partial fulfillment of the requirements for the degree of Doctor of Philosophy (Electrical Engineering) at the UNIVERSITY OF WISCONSIN-MADISON 2010 UMI Number: 3448781 All rights reserved INFORMATION TO ALL USERS The quality of this reproduction is dependent upon the quality of the copy submitted. In the unlikely event that the author did not send a complete manuscript and there are missing pages, these will be noted. Also, if material had to be removed, a note will indicate the deletion. UMI Dissertation Publishing UMI 3448781 Copyright 2011 by ProQuest LLC. All rights reserved. This edition of the work is protected against unauthorized copying under Title 17, United States Code. uest A ProQuest LLC 789 East Eisenhower Parkway P.O. Box 1346 Ann Arbor, Ml 48106-1346 ® FUNDAMENTAL, FORMATIONAL, AND PRACTICAL CHALLENGES IN MICROWAVE INVERSE SCATTERING FOR MEDICAL IMAGING OF THE BREAST submitted to the Graduate School of the University of Wisconsin-Madison in partial fulfillment of the requirements for the degree of Doctor of Philosophy By Jacob D. Shea Date of final oral examination: December 14, 2010 Month and year degree to be awarded: December 2010 The dissertation is approved by the following members of the Final Oral Committee: Susan C. Hagness, Professor of Electrical & Computer Engineering Barry D. Van Veen, Professor of Electrical & Computer Engineering Nader Behdad,, Professor of Electrical & Computer Engineering Stephen J. Wright, Professor of Computer Science Frederick Kelcz, Associate Professor of Radiology i To Baby Buster. ii ACKNOWLEDGMENTS With deepest gratitude I thank my advisors, Susan and Barry, for patiently and sagely guiding me through this academic adventure. Their intellect, enthusiasm, and warmth were equally indispensable. I give my utmost love and thanks to my parents, Dan and Elyse, who deserve full recognition for every success they have enabled me to achieve. My research was supported by the National Institutes of Health under grant R01 CA112398 awarded by the National Cancer Institute, a Graduate Fellowship from the University of Wisconsin - Madison, and a Doctoral Research Award from the IEEE Antennas and Propagation Society. iii TABLE OF CONTENTS Page LIST OF TABLES vi LIST OF F I G U R E S vii ABSTRACT xv 1 Introduction 1 2 Background 6 2.1 2.2 3 Microwave Medical Imaging Microwave Breast Imaging 2.2.1 Physical Justification 2.2.2 Prior Work 2.2.3 Inverse Scattering Microwave imaging of normal breast tissue distributions 3.1 3.2 3.3 3.4 Introduction Numerical Breast Phantoms 3.2.1 Test Phantoms 3.2.2 Scattered-Field Data Acquisition Imaging Method 3.3.1 Electromagnetic Inverse Scattering and the Distorted Born Iterative Method 3.3.2 Forward Solution 3.3.3 Frequency-Dependent Model of Unknown Dielectric Properties . . . . 3.3.4 Linear System of Scattering Equations 3.3.5 Inverse Solution 3.3.6 Assumptions and Initial Conditions Results and Discussion 3.4.1 Simple Phantom 3.4.2 Realistic Phantoms 7 9 10 13 19 24 25 29 30 32 35 35 36 38 39 41 44 46 48 54 iv Page 3.4.3 3.5 4 Other Performance Considerations Conclusion Contrast-enhanced tumor detection using microwave imaging 4.1 4.2 4.3 5 Introduction Microwave Imaging of Small Tumors Testbeds and Data Generation 4.3.1 Heterogeneous Numerical Breast Phantoms 4.3.2 Homogeneous Numerical Breast Phantoms 4.3.3 Data Acquisition 4.4 Imaging Methodology 4.4.1 Electromagnetic Inverse Scattering 4.4.2 Distorted Born Iterative Method 4.4.3 Differential Imaging 4.5 Results 4.5.1 Imaging a Tumor in a Homogeneous Phantom Interior 4.5.2 Imaging a Tumor in a Heterogeneous Phantom Interior 4.5.3 Differential Imaging of a Tumor in a Heterogeneous Phantom Interior 4.6 Discussion and Conclusion Analysis of Microwave Scattering in Breast Imaging 5.1 5.2 5.3 5.4 5.5 5.6 5.7 6 Introduction Background Methods Models Results 5.5.1 Single-Frequency Formulation 5.5.2 Multiple-Frequency Formulation 5.5.3 System Design 5.5.4 Error Evaluation Discussion Conclusion Future Work 6.1 6.2 Summary of research Directions for future research 6.2.1 System Design via Scattering Information Analyses 60 66 68 69 71 73 74 77 77 79 79 82 84 84 86 88 92 96 99 100 103 105 108 Ill 112 114 117 119 124 126 127 127 130 131 V Page 6.3 6.2.2 Custom Spatial Basis for Pre-Regularization 6.2.3 Generalization of Radar and Inverse Scattering 6.2.4 Data Acquisition via Beamforming 6.2.5 Generalization of Sourcing Functions Conclusion LIST OF R E F E R E N C E S 131 133 134 135 137 139 vi LIST OF TABLES Table 2.1 2.2 3.1 3.2 4.1 4.2 4.3 5.1 Page Cole-Cole parameters of the median dielectric properties of normal and malignant breast tissue samples excised during cancer surgeries [1] 12 Feature comparison of inverse scattering contributions to microwave breast imaging research 16 Debye parameters (infinite and delta relative permittivity, e ^ and Ae, and static conductivity, crs) of the materials modeled in the numerical phantoms (valid from 0.5 to 3.5 GHz) 31 Debye parameters (infinite and delta relative permittivity, el00 and Ac. and static conductivity, a s ) of the homogeneous initial profile assumed by the DBIM algorithm for each realistic phantom 46 Debye parameters of the tissue models used in the construction of the numerical phantoms (valid from 0.5 to 3.5 GHz). The relaxation time constant parameter is t=15.0 ps for all tissues 76 Peak contrast of estimated complex permittivity at 2.5 GHz of a malignant spherical inclusion in a homogeneous breast interior. For reference, the electrical size of the inclusion diameter in each background tissue is given in wavelengths. . . 89 Changes in dielectric properties and location errors at 2.5 GHz in the differential imaging of a 1.0-cm-diameter malignant spherical inclusion in realistic heterogeneous phantoms using microbubble and carbon nanotube contrast agents. . . . 92 Debye parameters (static and infinite relative permittivity, e s and e ^ , and static conductivity, ers) of the tissues modeled in the numerical breast phantoms. . . . 109 vii LIST OF FIGURES Figure 2.1 2.2 3.1 3.2 3.3 3.4 3.5 Cole-Cole models of the frequency dependence of the median (a) relative permittivity and (b) effective conductivity (S/m) of breast tissues excised during cancer surgeries. Variability bars show the 25 th — 75th percentile ranges Page 11 (a) Propagation loss of electromagnetic waves in breast tissue, in dB/cm. Typical breast dimensions are on the order of 5 cm deep and 8 cm in diameter, (b) Penetration depth of electromagnetic waves in breast tissue. Penetration depth is defined as the propagation distance in a medium before conductive and/or dispersive losses reduce the wave amplitude by a factor of 1/e (-8.686 dB). . . . 13 Frequency dependence of the single-pole Debye models for representative tissues used in the numerical breast phantoms, (a) relative permittivity, e r , and (b) effective conductivity, crefF 31 3-D numerical breast phantoms and dipole antenna arrays, (a) Class 1, (b) Class 2, (c) Class 3, (d) Class 4 33 At a sample DBIM iteration, the norm of the CGLS solution is plotted vs. the norm of the CGLS residual at each iteration of the algorithm. The knee of the 'L' curve (indicated by the ' Q ' marker) locates the iteration at which to terminate the CGLS algorithm 44 3-D, inverse crime reconstruction of the simple phantom with spherical inclusion, shown in coronal (x — y, top row), sagittal (y — z, middle row), and axial (x — z, bottom row) cross-sections of the reconstructed Debye parameters, (a)-(c) e ^ , (d)-(f) Ae, (g)-(i) a s (S/m) 50 3-D reconstruction of the simple phantom with spherical inclusion, shown in coronal (x — y, top row), sagittal (y — z, middle row), and axial (x — z, bottom row) cross-sections of the reconstructed Debye parameters, (a)-(c) (d)-(f) Ae, (g)-(i) (S/m) 51 viii Figure 3.6 3.7 3.8 3.9 Page Orthogonal transections of the inverse crime reconstruction of the simple phantom with spherical inclusion. The dashed line is the transection of the phantom, and the solid line is the transection of the reconstruction. The transections are taken through the center of the true position of the inclusion, x — 6.6 cm, y = 6.0 cm, z = 4.8 cm 52 Orthogonal transections of the reconstruction of the simple phantom with spherical inclusion. The dashed line is the transection of the phantom, and the solid line is the transection of the reconstruction. The transections are taken through the center of the true position of the inclusion, x = 6.6 cm, y = 6.0 cm, z = 4.8 cm 53 Convergence of the DBIM algorithm for the four heterogeneous phantoms. Lefthand y-axes: normalized convergence of the residual scattering norm. Right-hand y-axes: similarity metrics between the reconstruction at each iteration and the actual phantom profile. Here the cos(4>e) metric is computed for the concatenation of all three Debye parameters, relative to the exact phantom profile 55 Class 1 phantom: Exact (a)-(c) and reconstructed (d)-(f) profile of Ae shown in coronal (x — y, top row), sagittal (y — z, middle row), and axial (x — z, bottom row) cross-sections. Axes in cm 56 3.10 Class 2 phantom: Exact (a)-(c) and reconstructed (d)-(f) profile of Ae shown in coronal (x — y, top row), sagittal (y — z, middle row), and axial (x — z, bottom row) cross-sections. Axes in cm 57 3.11 Class 3 phantom: Exact (a)-(c) and reconstructed (d)-(f) profile of Ae shown in coronal (x — y, top row), sagittal (y — z, middle row), and axial (x — z, bottom row) cross-sections. Axes in cm 58 3.12 Class 4 phantom: Exact (a)-(c) and reconstructed (d)-(f) profile of Ae shown in coronal (x — y, top row), sagittal (y — z, middle row), and axial (x — z, bottom row) cross-sections. Axes in cm 59 3.13 Coronal, sagittal, and axial cross-sections of Class 2 phantom reconstructions: (a)-(c) scalar field approximation, (d)-(f) inverse crime with scalar field approximation. Cross-sections are taken from the reconstructed profiles of Ae, and are shown in coronal (x — y, top row), sagittal (y — z, middle row), and axial (x — z, bottom row) cross-sections. Axes in cm 61 ix Figure Page 3.14 The similarity metric, cos(</>e), between the actual and reconstructed Debye parameter profiles of each of the four realistic phantoms. The metrics for the £00, Ae, and a s profiles are denoted by ' x ' , 'A', and 'o', respectively. The left-hand side markers for each phantom are the results for unbounded, unconstrained reconstructions. The right-hand side markers for each phantom, enclosed by ' • ' markers, are the results for bounded, constrained reconstructions. Here the cos(0 e ) metrics are the results for the individual vectorizations of each Debye parameter 62 3.15 Effect of bounding and constraint techniques on reconstruction of static conductivity. Coronal cross-sections of static conductivity, a s , in (S/m) of the Class 2 reconstruction: (a) exact profile, (b) constrained reconstruction, (c) unconstrained reconstruction 62 3.16 Coronal, sagittal, and axial cross-sections of Class 2 phantom reconstructions: (a)-(c) scalar field approximation, (d)-(f) vector field formulation. Cross-sections are taken from the reconstructed profiles of Ae, and are shown in coronal (x — y, top row), sagittal (y—z, middle row), and axial (x—z, bottom row) cross-sections. Axes in cm 63 3.17 Reconstruction error vs. signal to noise ratio, (a) e 2 error metric, referenced to noiseless reconstruction, (b) similarity metrics: cos(0 e ) and cos(^ r ) are referenced to exact profile and noiseless reconstruction, respectively (both are computed for the concatenation of all three Debye parameters) 65 3.18 Comparison of Coronal cross-sections of the Ae reconstruction of the Class 2 phantom for (a) noiseless data, (b) 20 dB SNR, (c) 10 dB SNR, (d) 0 dB SNR. 65 4.1 Scattering cross-section of a homogeneous dielectric sphere of radius a as a function of diameter. The diameter of the sphere is given in units of wavelengths in the background medium, A6. The sphere has the properties of malignant tissue and is placed in a lossless version of three representative backgrounds: adipose, transitional, and normal fibroglandular tissues. The contrast in relative permittivity between the sphere and the background at 2.5 GHz for these three cases are 11.6, 2.55, and 1.15, respectively. 73 X Figure 4.2 4.3 4.4 Page (a) Relative permittivity, erei, and (b) effective conductivity, creff (S/m), of the ranges of adipose, transitional, and normal fibroglandular tissues used in phantom construction, (c) Relative permittivity, erei, and (d) effective conductivity, creff (S/m), of the malignant tissue models with and without the assumed effect of carbon nanotube or microbubble contrast agents 75 Diagram of the 3-D forward solution domain for the Class 2 numerical phantom, showing the dipole array, the downsampled skin layer, the measurement region V, and the imaging region V of the interior breast volume 78 Comparison of the exact (top row) and reconstructed (bottom row) coronal crosssections (in color) of relative permittivity for a single spherical inclusion of malignant tissue properties within a homogeneous breast volume of transitional tissue properties. The inclusion diameters are (a)-(b) 6 mm, (c)-(d) 14 mm, and (e)-(f) 22 mm. Spatial smoothing of the inclusion boundary and increasing underestimation of the inclusion properties with decreasing diameter are evident in the reconstructions. (Axes in cm.) 87 4.5 The normalized error in the estimation of the dielectric properties of the spherical malignant inclusion versus diameter in a homogeneous breast interior of (a) fibroglandular tissue, (b) transitional tissue, and (c) adipose tissue. The estimate is defined as the average of the reconstructed properties over the actual inclusion volume. The error is calculated as (u* — uest) / (u* — fx"11'), where u is either relative permittivity, erei, or effective conductivity, aeff, at 2.5 GHz. . . . 88 4.6 Reconstruction of the Class 1 phantom with a single 1.0-cm-diameter spherical malignant inclusion. Coronal (top row) and sagittal (bottom row) cross-sections (in color) of the exact and reconstructed volumes are shown at 2.5 GHz for (a)(d) relative permittivity and (e)-(h) effective conductivity in (S/m). (Axes in cm.) 90 Reconstruction of the Class 2 phantom with a single 1.0-cm-diameter spherical malignant inclusion. Coronal (top row) and sagittal (bottom row) cross-sections (in color) of the exact and reconstructed volumes are shown at 2.5 GHz for (a)(d) relative permittivity and (e)-(h) effective conductivity in (S/m). (Axes in cm.) 90 4.7 xi Figure 4.8 Page Reconstruction of the Class 3 phantom with a single 1.0-cm-diameter spherical malignant inclusion. Coronal (top row) and sagittal (bottom row) cross-sections (in color) of the exact and reconstructed volumes are shown at 2.5 GHz for (a)(d) relative permittivity and (e)-(h) effective conductivity in (S/m). (Axes in cm.) 91 Reconstruction of the Class 4 phantom with a single 1.0-cm-diameter spherical malignant inclusion. Coronal (top row) and sagittal (bottom row) cross-sections (in color) of the exact and reconstructed volumes are shown at 2.5 GHz for (a)(d) relative permittivity and (e)-(h) effective conductivity in (S/m). (Axes in cm.) 91 4.10 Differential images of the Class 1 phantom. Coronal (top row) and sagittal (bottom row) cross-sections of the complex permittivity at 2.5 GHz are taken through the peaks of the changes in (a)-(b) relative permittivity and (c)-(d) effective conductivity (S/m) due to microbubble contrast agent, and (e)-(f) relative permittivity and (g)-(h) effective conductivity (S/m) due to carbon nanotube contrast agent. The ' x ' markers show the actual location of the tumor as projected onto each cross-sectional plane. (Axes in cm.) 93 4.11 Differential images of the Class 2 phantom. Coronal (top row) and sagittal (bottom row) cross-sections of the complex permittivity at 2.5 GHz are taken through the peaks of the changes in (a)-(b) relative permittivity and (c)-(d) effective conductivity (S/m) due to microbubble contrast agent, and (e)-(f) relative permittivity and (g)-(h) effective conductivity (S/m) due to carbon nanotube contrast agent. The ' x ' markers show the actual location of the tumor as projected onto each cross-sectional plane. (Axes in cm.) 93 4.12 Differential images of the Class 3 phantom. Coronal (top row) and sagittal (bottom row) cross-sections of the complex permittivity at 2.5 GHz are taken through the peaks of the changes in (a)-(b) relative permittivity and (c)-(d) effective conductivity (S/m) due to microbubble contrast agent, and (e)-(f) relative permittivity and (g)-(h) effective conductivity (S/m) due to carbon nanotube contrast agent. The ' x ' markers show the actual location of the tumor as projected onto each cross-sectional plane. (Axes in cm.) 94 4.9 xii Figure Page 4.13 Differential images of the Class 4 phantom. Coronal (top row) and sagittal (bottom row) cross-sections of the complex permittivity at 2.5 GHz are taken through the peaks of the changes in (a)-(b) relative permittivity and (c)-(d) effective conductivity (S/m) due to microbubble contrast agent, and (e)-(f) relative permittivity and (g)-(h) effective conductivity (S/m) due to carbon nanotube contrast agent. The ' x ' markers show the actual location of the tumor as projected onto each cross-sectional plane. (Axes in cm.) 94 4.14 Peak pertubation of relative permittivity in microbubble difference images reconstructed from noise-corrupted measurement data from each of the four heterogeneous numerical breast phantoms. The plots indicate the maximum absolute perturbation (over all voxels in V) due to a given noise level, relative to the noiseless difference image 96 5.1 Cross sections of the static permittivity parameter, e s , of a) the scattered fibroglandular (Class 2) phantom, and b) the heterogeneously dense (Class 3) phantom. Axes in cm 108 (a) Class 2 and (b) Class 3 breast phantoms with array locations shown by vertically polarized (z-directed) 4 mm sources 110 5.3 Measured gain of the FDTD sources Ill 5.4 Fidelity of TSVD estimates of complex permittivity using single-frequency data from 820 measurement channels of the Class 3 phantom versus a) index of singular value truncation, and b) SNR. The SNR results are normalized to the gain of the source at each frequency. The required measurement sensitivity for an operating point on one of the SNR curves is determined by subtracting the gain, in dB, of the transmit and receive antennas 113 Coronal cross-sections of static permittivity, e s , a) the actual Class 3 breast phantom, and the optimal TSVD projections of the scattering data at b) 0.5, c) 1.0, d) 1.5, e) 2.0, and f) 2.5 GHz. Axes in cm 114 Fidelity versus SNR of the TSVD estimates of relative permittivity, e r , and effective conductivity, creff, using single-frequency data at 0.5 and 2.5 GHz. The SNR results are normalized to the transmit and recieve gain at each frequency. 115 5.2 5.5 5.6 Figure 5.7 Coronal cross sections of exact phantoms and TSVD estimates of static permittivity, e s , using multiple-frequency data at 0.5, 1.0, 1.5, 2.0, and 2.5 GHz. a) exact Class 2 phantom, b) Class 2 estimate, c) exact Class 3 phantom, d) Class 3 estimate. Axes in cm Page 116 5.8 a) Fidelity versus SNR of TSVD estimated Debye parameters for the Class 3 phantom using multiple-frequency data at 0.5, 1.0, 1.5, 2.0, and 2.5 GHz over 820 measurement channels, b) Comparison of the fidelity versus SNR of the static permittivity parameter by frequency selection. The fidelity of the relative permittivity for the 2.5 GHz single-frequency formulation is included for reference. 117 5.9 Fidelity of the static permittivity parameter versus SNR for vertical and horizontal polarizations. The TSVD estimates are computed from multiple-frequency data at 0.5, 1.0, 1.5, 2.0, and 2.5 GHz over 820 measurement channels 118 5.10 Fidelity of the parametric Debye estimate versus SNR for "through" channels and "reflect" channels using multiple-frequency Debye formulation for data at 1.5, 2.0, and 2.5 GHz 120 5.11 Coronal cross sections of multiple-frequency (1.5, 2.0, and 2.5 GHz) optimal TSVD image estimates for a) through channels and b) reflect channels. Axes in cm 120 5.12 Comparison of the fidelity versus truncation index using the exact vector field formulation and the approximate scalar field formulation. Class 3 phantom for a) single-frequency data at 2.5 GHz. b) multiple-frequency data at 1.5, 2.0, and 2.5 GHz 121 5.13 Coronal cross-sections of the estimated Class 3 phantom for multiple-frequency data at 1.5, 2.0, and 2.5 GHz. a) Background includes exact skin region and homogeneous interior of the volume-averaged tissue properties, b) Background is the homogeneous immersion medium. Axes in cm 122 5.14 Fidelity of the e s estimate versus SNR for estimation for the Class 3 phantom with and without the exact skin model included in the background media. Multiplefrequency data at 1.5, 2.0, and 2.5 GHz are used in the estimation 122 5.15 Fidelity of the e s and a s estimates versus truncation index with and without the Born field approximation for multiple-frequency data at 1.5, 2.0, and 2.5 GHz. . 123 xiv Figure 6.1 Reconstructions of static permittivity, e s , the Class 2 (ID 010204) numerical breast phantom for (a) the 2 mm voxel basis of over 65,000 basis functions, and (b) the custom eigenvector basis of fewer than 2,300 basis functions Page 133 XV ABSTRACT Microwave frequency electromagnetic illumination is widely investigated for application to medical imaging. In particular, the breast is an excellent candidate for microwave medical imaging research given the application's technical viability and epidemiological relevance. From a technical perspective, waves scattered from the fibroglandular tissue structures are measurable with practical hardware, due to the accessibility of the breast to multi-view illumination and the relatively low loss of the adipose tissue which comprises much of the interior. From a public health perspective, since breast cancer has the highest morbidity and second highest mortality rates among women, any improvements over the current screening and detection methodologies may have a significant impact. Microwave imaging has the potential to address several shortcomings in contemporary medical imaging technology: the illumination is low-power and non-ionizing, the hardware comparatively portable and inexpensive, and the microwave frequency dielectric contrast between malignant and normal tissue provides a physical basis for tumor detection. Investigations by the research community into the clinical utility of microwave images of the breast have been largely inconclusive. Challenges have been encountered including insufficient image resolution, uncertain sensitivity and specificity of tumor detection, and stringent requirements on measurement and modeling accuracies. Our research objective is to further the understanding of the fundamental limitations and potential imaging performance of microwave breast imaging. State-of-the-art numerical breast phantoms, realistic in both anatomical structure and physical properties, are used to ensure faithful simulations of clinical data. Phantoms with varied density of healthy fibroglandular tissue are studied to reveal modest imaging performance with potential for density assessment. Phantoms with xvi contrast-enhanced tumors are imaged differentially to show the measurable perturbation of the scattering data by exogenous agents. A method of analysis of the scattering data is implemented to estimate an upper-bound on image quality relative to various design and error considerations. The totality of this research suggests that optimal design of measurement systems and further development of imaging methods may lead to clinically relevant microwave images of the breast. 1 Chapter 1 Introduction The incidence rate of breast cancer in the United States is the highest of any cancer among women and is twice as high as the incidence rate of the second most prevalent cancer, lung cancer [2], A woman has a 12% lifetime risk of developing breast cancer and a 2.8% of the disease becoming terminal [3]. The mortality rate of breast cancer is second only to lung cancer, and there is a strong correlation of the mortality rate to tumor size at detection [4], suggesting the value of regular screening for early detection of breast cancer. The World Health Organization has forecast that while cancer mortality rates are higher in developed countries, globalization will result in increased cancer rates in developing countries [5]. Breast cancer is thus a major contemporary health concern not only because of its prevalence, but also because of low availability of regular screening in underdeveloped and underprivileged parts of the world. Although many imaging modalities have been proposed for breast cancer screening, only x-ray imaging and magnetic resonance imaging are currently included in screening guidelines for early detection [6], with ultrasound imaging used as a secondary diagnostic tool. X-ray mammography has been the standard screening technology for decades and its use in regular screening has been shown to improve survival outcomes [7]. Mammography is capable of imaging small, non-palpable tumors and calcifications, but weaknesses such as the ambiguity in a two-dimensional (2-D) projection image, particularly of radiographically dense breasts, contribute to false positive error rates around 10% (low specificity) and false negative error rates as high as 30% (low sensitivity) [8]. In addition, both the discomfort of 2 breast compression and the risk of the radiation exposure make it an unappealing medical procedure from a patient perspective and may have been factors for some of the nearly 30% of women over age 50 who have not had a mammogram in at least two years [9]. Ultrasound imaging is perhaps more amenable to the patient than mammography and is used as an ancillary imaging procedure after a suspicious finding by mammogram or clinical breast exam. However, there is low acoustic contrast between malignant and normal tissues, detection performance is operator dependent, and there has been no evidence to support broad use of ultrasonography in general population screening [10]. Magnetic resonance imaging (MRI) produces very detailed three-dimensional (3-D) images of the breast with high sensitivity, but low specificity, to malignancies. Factors such as high cost, long imaging time, and uneven availability make it better suited to screening high risk candidates than screening the general population. Among emerging technologies in medical imaging, the application of microwave imaging to breast cancer screening has been widely investigated as an alternative or complementary modality, given the potential for low relative cost, 3-D images, a comfortable and quick procedure using low-power, non-ionizing electromagnetic waves, and the contrast available in dielectric properties of tissues at microwave frequencies. Microwave imaging techniques include both radar-inspired methods and tomographic methods [11,12]. After significant recent effort by several groups on both approaches, tomographic inverse scattering methods appear to hold the most promise for microwave breast imaging. The increased emphasis on inverse scattering techniques in the microwave breast imaging application follows from a few notable considerations. First, a new understanding of the microwave properties of breast tissues suggests that the viable technique must be able to image both healthy and malignant glandular tissue. Second, spurious image artifacts caused by an assumption of linear scattering suggests that the successful imaging method must account for the nonlinearity of microwave scattering from the high contrast between the tissue properties within and surrounding the breast. 3 Radar-based imaging methods are useful for imaging the strongest scattering sites in an unknown distribution of tissue. In such techniques an array of microwave-frequency antennas is used to illuminate the breast, and the reflections and transmissions received by the array are focused to the locations of strongly scattering tissue structures in the interior breast volume [13-15], Radar-based methods of microwave imaging are similar to synthetic aperture radar techniques, and will be referred to in this work as tissue-penetrating radar (TPR) in acknowledgment of the similarity to the well-established application of ground-penetrating radar (GPR). In both applications, transmit and receive antennas are positioned above a surface and the sub-surface scatterers are located by beamforming techniques that spatially focus the received signals to the scattering sources. Unlike the half-space propagation model for the earth's surface in GPR, in T P R there is no analytical Green's function expression for the skin surface enclosing the breast volume to be imaged. Instead, simple models of the complicated wave interactions in the skin and breast tissues are used, and the scattering is usually approximated as linear. The resulting error manifests as poor localization and focusing accuracy. Furthermore, fibrous and glandular tissue structures are larger and more extensive than most malignant tissue formations and have similar dielectric property contrast in a background of adipose tissue. This results in a low signal-to-clutter ratio when using radar to detect and locate small tumor signals in the presence of interfering clutter from nearby fibroglandular structures. Inverse scattering methods employ the same basic measurement approach as in TPR, using a transmit-receive antenna array surrounding the breast, but differ significantly in image formation by providing an estimate of the complete spatial profile of dielectric properties [16,17]. Since the full dielectric profile yielded by the estimation completely defines a propagation environment, the error introduced by the simplified propagation models of radar methods can be mitigated. Inverse scattering methods are thus able to account for the nonlinearity of the multiple scattering interactions that contribute to the measured data. Furthermore, rather than assuming tumors to be the primary scattering sites, the properties and the distributions of all constituent breast tissues are estimated. The resulting 3-D 4 images could potentially provide a clinician with dielectric property estimates of tissue formations, tissue density analysis for risk assessment, and detection of tumors by either visual assessment, quantitative analysis, or the introduction of exogenous contrast agents. In this work, microwave breast imaging is investigated using techniques of inverse scattering. The research builds directly upon prior work by the University of Wisconsin Computational Electromagnetics laboratory (UWCEM) in the area of microwave breast imaging. Lazebnik et al. [1,18] published extensive data on the microwave dielectric properties of ex-vivo breast tissue samples, significantly revising the common understanding of the absolute and relative properties of the constituent breast tissues. Zastrow et al. [19] developed anatomically realistic numerical breast phantoms derived from magnetic resonance images (MRI) of breasts from each classification of mammographic breast density [20]. The accurate information on dielectric properties from the ex-vivo tissue study and the realistic 3-D tissue profiles from the MRIs result in the first-of-their-kind, high fidelity, numerical testbeds capable of faithfully simulating clinical data. Winters et al. [21], imaged the numerical breast phantoms with a simple tumor model added to each using an inverse scattering algorithm to reconstruct the tissue profiles. The images were the first fully 3-D solutions of both the direct and inverse problems for realistic breast targets. However, these images (as well as other prior results of microwave breast imaging research) have presented results that arguably neither confirm nor deny the viability of microwaves for breast imaging. Both numerical [21] and clinical [22] investigations have shown images with such low spatial resolution that the performance of the systems in correctly detecting and locating malignant tissues is ambiguous at best. Low spatial resolution alone does not disqualify microwave medical imaging, so long as the technology is of some clinical utility in diagnosis or risk assessment. However, evaluation of the clinical utility of microwave imaging requires that researchers fully understand the best case performance expectations and carefully and conclusively corroborate the images with known benchmarks. 5 Thus, the viability of microwave imaging in medicine has yet to be determined. This dissertation addresses the need to further evaluate the fundamental, formational, and practical challenges in microwave breast imaging with the goal of furthering our understanding of practical performance expectations. Fundamental issues include the critical evaluation of microwave tomography in application to breast imaging with respect to the underlying physics of microwave illumination and measurement. Formational issues are those that influence the way the problem is posed and solved, the measurement data of interest, and the use of a priori information specific to breast imaging. Practical issues include the evaluation of likely sources of measurement and model error, the development and implementation of imaging algorithms for experimental and clinical data, and the design of the test fixture and microwave measurements. The history of the use of microwaves in medical imaging and the theoretical scattering relationship underlying electromagnetic imaging are detailed in Chapter 2. Numerical breast phantoms containing only normal tissue are imaged in Chapter 3 to first explore the ability to reconstruct and distinguish adipose from fibroglandular tissues before considering the issue of tumor detection. A technique for differentially imaging the effect of exogenous contrast agents is described in Chapter 4, and the preliminary results show the potential for contrastenhanced tumor detection. In Chapter 5, a general method of analysis for inverse problems is applied to the microwave breast imaging problem in order to evaluate the information content of the scattering data for various measurements and formulations. A summary of the progress embodied by this research is given in Chapter 6, followed by details of a number of further research opportunities inspired by this work. 6 Chapter 2 Background Microwave-frequency radiation, broadly defined as the range from 300 MHz to 300 GHz, has a long history in the field of medicine. Although time-harmonic currents at sub-microwave frequencies were applied with therapeutic intention as early as the 1840's, technological limitations dictated that investigation of the potential therapeutic or diagnostic effects of microwave frequencies didn't begin in earnest until the 1930's [23]. The earliest application was electromagnetic diathermy, the therapeutic heating of deep tissue, which continues to find use in therapeutic applications [24], In contemporary medicine, microwaves can be exploited for a range of biophysical phenomena, including interactions at the cellular level, conversion of electromagnetic to thermal energy in tissues and vice-versa, and field shaping by wave superposition. In microwave ablation, cellular death is induced in malignant tissue by a sufficient rise in temperature due to the microwave heating. A localized field for spatially-selective heating may be introduced either non-invasively by the focusing of an external antenna array [25] or invasively by a percutaneous microwave applicator tipped with a single antenna [26]. Similar to focused ablation, microwave hyperthermia is a noninvasive technique which focuses an external array to heat a region of malignant tissue, thereby improving the response to radiation or chemotherapy [27]. Several other examples of therapeutic applications of microwaves are detailed in [28]. The reciprocal to the process of focusing waves is the scattering of waves from tissue structures. Similar to conventional radar, the scattering response can be used to localize the scattering sources 7 formed by contrasts in the dielectric properties of different tissue types. Microwaves have thus been studied in a number of medical imaging techniques. Microwave medical imaging, specifically microwave breast tomography, is the area in which the proposed work is focused. In this chapter, a brief history of the use of microwaves in medical imaging is given, techniques and results of existing microwave breast imaging research are summarized, and the method of electromagnetic inverse scattering, the technique of interest to the proposed work, is detailed. 2.1 Microwave Medical Imaging The potential for the use of microwave illumination in creating medical images was generally viewed with pessimism until the later decades of the 20th century. Unlike X-rays, microwaves undergo significant reflection at the skin, have poor tissue penetration depth at the frequencies of desirable spatial resolution, do not travel in a straight-line collimated beam, and exhibit multiple reflections between objects illuminated by the wide beam of a source antenna. In the late 1970's, Larsen and Jacobi [29-32] addressed several of these issues to produce some of the first successful active microwave images of biological tissues. A mechanically scanned transmitter and receiver antenna pair was placed in a water bath to improve both the resolution of, and power transmitted into, a canine kidney immersed in the water between the antennas. The problem of the wide microwave beam width causing scattered power to be received from angles other than line-of-sight was tackled by Larsen and Jacobi using two distinct methods. The first method was diffraction tomography. It resembles the application of the Fourier slice theorem to X-ray computed tomography in its use of the Fourier diffraction projection theorem to transform a set of transmission measurements to the spatial frequency domain [29,31,33,34]. Multiple angles of illumination populated the Fourier representation of the object, which was then inverse transformed to obtain the image. The second method was referred to as time-delay spectroscopy [30,32,35], and can be considered a radar technique in that it directly considers the differences in multipath propagation time to spatially distinguish scattering contributions. Both methods reconstructed 8 impressively detailed images of the interior kidney structures. Their results represent some of the earliest work in the application of both microwave tomography and microwave radar to medical imaging, two approaches to breast imaging that have received significant attention from researchers over the last decade. Although computational, algorithmic, and measurement techniques varied widely, much of the subsequent work on microwave medical imaging over the last 30 years could be classified by whether or not the technique addresses the challenge that Larsen and Jacobi did not propose a solution for in their initial work: that electromagnetic scattering is nonlinear. Specifically, scattering interactions occur not only between the antennas and individual features of the tissues in the imaging region, but also between the features themselves. The difficulty in the nonlinearity of the problem was manifest in the limitations of imaging methods that either 1) quantitatively reconstructed the object by solving for the dielectric properties while incurring significant errors due to a linear scattering approximation, or 2) qualitatively reconstructed the object by solving a linear formulation of the problem that recovers only the spatial profile of equivalent currents (the product of the total electric field and the dielectric properties contrast) [36]. Early work on qualitative microwave medical imaging included ongoing development of Larsen and Jacobi's work by Guo and Guo [37,38] and the microwave camera system of Bolomey, Pichot, Jofre, and others [39-41], Quantitative imaging work included numerical investigations of thoracic imaging by Ghodgaonkar et al. [42] and the use of a cylindrical scanner to image experimental phantoms of tissue structures by Jofre et al. [43], Given the limitations of the earlier methods, by the early 1990's iterative techniques for nonlinear optimization were widely studied for application to quantitative microwave medical imaging. Iterative methods continue to be a focal point of microwave medical imaging research, and the research aims proposed herein fall into that category. In short, iterative methods make successive refinements to an estimate of an unknown object, typically using a linear scattering approximation for each refinement, such that there is less error in the field approximation at each step and the estimate converges to a solution close to 9 the exact object. The Born iterative method (BIM) and distorted Born iterative method (DBIM) [16] are two ubiquitous examples of iterative methods in microwave imaging. DBIM is an electromagnetics-specific application of Gauss-Newton (GN) or Newton-Kantorovich (NK) methods of nonlinear optimization [44]. Contrast-source inversion (CSI) [45] is another popular method, which operates by alternately estimating the equivalent source profile and the dielectric profile. The publication record for nonlinear inversion methods for electromagnetic imaging is extensive (e.g. [46-51]), but, although the potential for biomedical applications surely was not lost on the authors or the community, the techniques have only been widely applied to microwave medical imaging in the last two decades. In 1991, Joachimowicz et al. [52] proposed using the NK method for microwave imaging along with a priori information about the tissue or organ to be imaged. They applied the method to simple numerically modeled arrangements of pixels having the nominal properties of water, fat, muscle, or bone. Meaney et al. [53] (1996) reconstructed experimental phantoms of bone and fat tissue-equivalent materials. In the very same journal issue, Semenov et al. [54] asserted a lack of interest in microwave medical imaging in light of the established techniques of X-ray imaging and NMR tomography, but highlighted the potential for the correlation of physiological conditions to changes in microwave-frequency dielectric properties. They placed a perfused and beating mongrel heart in their experimental imaging system and demonstrated microwave images of various phases of the muscle contraction cycle. The same group published reconstructions of a simple numerical model of a heart in 2-D [55] and later in 3-D [56,57], and then compared the performance of the Newton method and CSI in reconstructing microwave images of a ham shank [58]. Abubakar et al. applied their CSI method to experimental models of the arm and neck in 2002 [59] and to a numerical model of the thorax in 2008 [60]. 2.2 Microwave Breast Imaging Techniques of microwave breast imaging can be classified into passive methods, active methods, or a hybrid of passive and active methods [61]. Passive methods measure existing 10 physiological conditions, such as the increased temperature of malignant tissue relative to the surrounding breast tissue. Barrett and Myers [62] proposed microwave thermography in 1975 to exploit this phenomena. In this technique, external radiometers are used to measure microwave radiation from the surface of the breast due to thermal conditions of the interior breast tissue. Active methods, such as inverse scattering, illuminate the breast with electromagnetic waves and form an image of the breast interior based on the measurements of the waves scattered from interior tissue structures. Hybrid methods combine the two approaches; for example, microwave illumination from an active method can enhance the temperature differential that is measured passively. 2.2.1 Physical Justification The physiological basis of active microwave imaging of biological tissues is the scattering of electromagnetic waves from interfaces between tissues with different dielectric properties. The dielectric properties of tissues is primarily related to their water content. The amount of free water in constitutive breast tissues can vary from 40% to 80% by volume [28], and thus the high permittivity and high conductivity of water creates a large contrast between the dielectric properties of the low water content adipose (fatty) tissues and the high water content stromal (fibrous) and glandular tissues. The latter two are often collectively referred to as fibroglandular tissue. Cancer typically forms in glandular tissues and may exhibit higher perfusion (vascularization and blood flow), potentially resulting in higher dielectric properties than healthy glandular tissue. Prior to the recent publication of dielectric measurements of freshly excised malignant and healthy breast tissues [1,18], the contrast between malignant and normal glandular tissues was widely accepted to be significant, with assumed contrast ratios of at least 2:1 and as high as 10:1 [28,61]. The thorough, large-scale study of excised tissues by the University of Wisconsin and University of Calgary credibly showed 10:1 contrast in complex permittivity between the median properties of adipose and glandular tissues, but only 1.1:1 contrast between the median properties of malignant and 11 normal glandular tissues. A recent clinical study by Dartmouth College suggested a similarly small mean contrast in permittivity of about 1.2:1, but a larger contrast in conductivity as high as 2:1. Given this new understanding, one might speculate that the discrepancies in the previously published data may have arisen from inaccurate characterization of the proportion of adipose tissue in the tissue samples. (a) • (b) Figure 2.1 Cole-Cole models of the frequency dependence of the median (a) relative permittivity and (b) effective conductivity (S/m) of breast tissues excised during cancer surgeries. Variability bars show the 2hth — 75th percentile ranges. The frequency dependence of biological tissues in the microwave frequency range is shaped by a number of relaxation phenomena, including the relaxations of cellular membranes and protein structures around 300 MHz, water bound to cells around 3 GHz, and free water around 25 GHz [28]. The single-pole Cole-Cole model is sufficient to capture the aggregate frequency dependence of the complex permittivity of tissues from 500 MHz to 20 GHz [63]. The Cole-Cole model expresses complex permittivity, e, as a function of angular frequency, uj, and the parameters of the model. i(uj,£00,A£,T,a,as) = £00 + 1 + {JUT) , + (2.1) JUEQ The relative permittivity is e^ at frequencies far above the relaxation pole at fT « 1/27it. At frequencies well below the relaxation pole the relative permittivity is + As and the 12 conductivity is <rs. In the exponent of the relaxation pole, a is a tuning parameter for the transition in the relaxation region between decreasing permittivity and increasing dispersive losses. The dielectric properties of adipose and fibroglandular tissues as reported by [1] are shown over frequency in Fig. 2.1. The corresponding Cole-Cole parameters are given in Table 2.1. Table 2.1 Cole-Cole parameters of the median dielectric properties of normal and malignant breast tissue samples excised during cancer surgeries [1], Tissue Type % Adipose As (S/m) r (ps) a Normal, Fibroglandular (group 1) 0-30 7.24 46.0 0.81 10.3 0.049 Transitional (group 2) 31-84 6.08 19.3 0.30 11.5 0.057 Adipose (group 3) 85-100 3.58 3.34 0.053 15.2 0.052 6.75 50.1 0.79 10.5 0.051 Malignant Since the contrast between normal and malignant glandular tissues is roughly constant over lower microwave frequencies, the selection of operating frequencies for breast imaging systems is generally guided by a trade off between the spatial resolution and the penetration depth of microwaves in the breast volume. Imaging resolution is nominally limited to half of the wavelength of the highest spatial frequency of the illumination within the breast volume. Penetration decreases with frequency, as illustrated by the propagation loss and 1/e depth of adipose and fibroglandular tissues in Fig. 2.2. Thus, the upper end of the operating frequency range is ultimately limited by the dynamic range of the measurement system, since the high frequency waves scattered from deep within the breast volume must exceed the noise floor of the measurement after attenuation by the round trip propagation loss. Note also that practical considerations specific to a particular imaging technique may further restrict the upper frequency limit. Lin [64] compared the wavelength and attenuation of fat and muscle over frequency and concluded that 2-8 GHz was an optimal range for microwave imaging. The FCC has designated the ultra wideband (UWB) of 3.1-10.6 GHz for use in medical applications, although some existing approaches use frequencies as low as 100 MHz. Most 13 radar-based methods have used the UWB band, while inverse scattering methods often use the UHF band (0.3-3 GHz). (a) (b) Figure 2.2 (a) Propagation loss of electromagnetic waves in breast tissue, in dB/cm. Typical breast dimensions are on the order of 5 cm deep and 8 cm in diameter, (b) Penetration depth of electromagnetic waves in breast tissue. Penetration depth is defined as the propagation distance in a medium before conductive and/or dispersive losses reduce the wave amplitude by a factor of 1/e (-8.686 dB). 2.2.2 Prior Work Active methods that have been investigated for microwave imaging of the breast include thermo acoustic tomography (TAT), tissue-penetrating radar (TPR), and microwave tomography (MT) based on spectral or inverse scattering techniques. In TAT techniques, thermal heating of tissues by externally applied microwaves enhance ultrasonic imaging of the breast [65-67]. While TAT need not be considered further in this context, the background of TPR and inverse scattering methods will both be covered in more detail, given their significant similarities. Radar-based methods of microwave medical imaging, including space-time beamforming [14], synthetic aperture radar [15], and time-reversal [68] techniques, are collectively referred to by this author as tissue-penetrating radar (TPR) in acknowledgment of the similarity to the well-established application of ground-penetrating radar (GPR). In the GPR class of 14 radar techniques, subsurface features of the earth, or layered media in general, are illuminated by transmit and receive antennas above the surface of the earth. Strong scatterers such as land mines or transitions from soil to bedrock may be revealed through the use of the halfspace Green's function and beamforming techniques to spatially focus the received signals to the source of their scattering. Analogously in TPR, transmit and receive antennas are positioned above a skin surface and the underlying tissues are illuminated so that the return signal may be processed to locate strongly scattering tissue formations. In contrast to GPR, in microwave TPR there is no analytical Green's function expression for the skin surface or heterogeneous breast interior. T P R methods were attractive under the earlier assumption that there is high contrast (2:1 or higher) between malignant and normal breast tissues, and thus prompted a significant amount of recent research on their application to breast imaging. In 1999, Hagness and Taflove [69] used a 3-D FDTD analysis analogous to a GPR half-space configuration to show the feasibility of using UWB confocal imaging for detecting early-stage tumors as small as 2 mm in diameter with a 5:1 contrast in permittivity. Fear and Stuchly [11] expanded the T P R concept in 2000, with an investigation of confocal imaging of small tumors (5:1 contrast) in a cylindrical numerical breast phantom measured monostatically by a circular antenna array. They concluded the technique was feasible, but noted a significant challenge presented by the large interfering scattering response of the skin layer. In 2003, Davis [70] and Bond [71] et al. at the University of Wisconsin advanced a space-time beamforming technique (MIST) for UWB confocal imaging and proposed an algorithm for cancellation of the interference from the skin layer. The method produced a clear image of 2 mm tumors (5:1 contrast in permittivity and 12:1 in conductivity). Around the same time, at the University of Calgary a similar UWB method was developed and dubbed tissue-sensing adaptive radar (TSAR) in observance of the similarity to traditional adaptive radar methods. A 2-D TSAR technique was successfully applied to experimental data by Sill and Fear [15] in 2005 for a 1 cm tumor model (10:1 contrast in permittivity and 35:1 in conductivity) in a simple cylindrical breast model. A hemispherical 3-D experimental and clinical measurement system was reported in 15 2009 by Klemm et al. [72] at the University of Bristol. The system uses UWB multistatic array measurements and the multistatic adaptive beamforming method proposed by Xie et al. [73], and imaged 4-6 mm tumors (5:1 contrast in permittivity and 30:1 in conductivity) in a simple hemispherical breast phantom. Although multistatic array measurements significantly improved the beamforming performance, the authors acknowledge that radar-based methods face a major challenge in imaging the low contrast of tumors buried in normal fibroglandular tissue. The newly refined understanding of the comparative dielectric properties of the constituent breast tissues indicates not only that contrasts as low as 1.1:1 must be resolved in order to detect a tumor in dense fibroglandular tissue, but that the low contrast of the tumor must be resolved amid the 10:1 contrast of a distribution of fibroglandular tissue against a background of adipose tissue. Therefore, the near-term emphasis in microwave breast imaging must include the potential of the modality to image complete fibroglandular tissue distributions. The problem of tumor imaging thus becomes somewhat subordinate to the ability of microwave imaging techniques to image realistic fibroglandular distributions with sufficient spatial and contrast resolution. The diversity of the density of fibroglandular distributions of the breast (see BIRADS breast density classifications [20]) from patient to patient further complicates the challenge of imaging realistic tissue distributions, particularly in comparison to the simple homogeneous adipose backgrounds of previously studied numerical and experimental phantoms. Despite the high contrast between glandular and adipose tissues, T P R methods are likely as ill-suited to normal tissue imaging as they are to tumor imaging because the extent and distribution of fibroglandular structures result in significantly more complicated scattering responses than those of small, isolated, high contrast tumors. Many of the considerations discussed thus far point to microwave tomography based on inverse scattering as the technique with the most potential for microwave breast imaging. Indeed, various approaches to the application of inverse scattering to breast imaging have long been established in the publication record. Inverse scattering contributions in the 16 T a b l e 2.2 Feature comparison of inverse scattering contributions to microwave breast imaging research. Frequency Tumor Inverse Author Year Data Dims. (GHz) #Ant. Contrast Method [8] Meaney" 2000 Clin 3D/2D 0.3-1.0 16 n/a n/a [74] Bulyshev b 2001 Num 3D 2, 3.5, 5 32 2.3:1 Gradient [75] Fang a 2004 Num 3D/2D 0.9 16 1.7:1 GN [76] Fang a 2004 Num, E x p 3D/2D 0.3-0.9 16 8:1, 4:1 LM [77] Zhang c 2004 Num 3D 0.8 60 2-3:1 CSI [22] Meaney a 2007 Clin 3D 0.5-1.7 16 n/a GN [78] Rubaek d ' a 2007 N., E., C. 3D/2D 1.1 16 4:1, 1.7:1 GN [79] Johnson e 2008 Num 2D 1.4-2.6 16 2:1 FTBS [80] Gilmore-f' 9 2009 Num 2D 2-6 32 20:1 CSI [21] Winters^ 2009 Num 3D 1.2-2.7 40 1.5:1 GN [81] Rubaek d ' a 2009 Num 3D 1.3 32 2.7:1 GN [82] Irishina 1 2010 Num 2D 1-5 40 1.4:1 Lvl.Set [83] Shea' 1 2010 Num 3D 1-2.5 40 1.15:1 GN [84] Shea ft 2010 Num 3D 1-2.5 40 n/a GN "Dartmouth College, Hanover, NH 6 c Carolinas Medical Center, Charlotte, NC Duke University, Durham, NC ^Technical University of Denmark, Lyngby, Denmark e Nagasaki University, Nagasaki, Japan •^University of Manitoba, Winnipeg, MB s Schlumberger-Doll Research, Cambridge, MA ^University of Wisconsin, Madison, WI 'Universidad Carlos III de Madrid, Leganes, Spain c o n t e x t of b r e a s t imaging have p r e d o m i n a n t l y b e e n m a d e in t h e present decade. Table 2.2 c o m p a r e s some of t h e details of t h e i m a g i n g s y s t e m s a n d a l g o r i t h m s r e p o r t e d for t h i s 17 application. A preliminary clinical investigation was performed by Meaney et al. [8] based on their earlier reported transceiver system [85]. The resulting images were very low resolution, but showed agreement between coronal slices and breast shape, similarity in contralateral comparisons, and correlation to patients' radiographic breast density. Various improvements to their imaging algorithm were later contributed [86-88]. Bulyshev et al. [74] imaged a 3-D numerical model of a breast immersed in a fat-matching medium. Their model included a homogeneous hemisphere of fat, a single 6-mm-diameter spherical tumor, and a chest wall of layered slabs of fat, muscle and bone. The peak of the reconstruction was centered at the tumor, but was underestimated by about 80%. No skin layer was modeled, so the nearly perfect coupling from the fat-matched immersion was grossly unrealistic. Fang et al. [75] used a 2-D approximation to invert 3-D data in 2-D slices for a numerical breast model with a 3-cm-diameter tumor. Here a water immersion was used, but no skin region was modeled. In a separate publication, Fang et al. [76] incorporated multiple frequencies into the inverse solution using an a priori assumption of the frequency dependence of tissues. The reconstruction of a 3-cm-diameter tumor is shown for several combinations of frequencies, with the general result that more is better. Zhang and Liu [77] confirmed the ability of the CSI method to invert 3-D high contrast object domains and suggested the method for breast imaging, however their numerical models consisting of cubes or spheres were exceedingly simple. The results of a 150-patient clinical study at Dartmouth College were reported by Meaney et al. [22] and Poplack et al. [89]. Meaney examined the in vivo dielectric properties suggested by the resulting microwave images. The findings included high bilateral symmetry in the average microwave properties of the entire breast and a correlation between radiographic density and microwave dielectric properties. The images were registered against MRIs of the same patients, from which it is clear that while the breast is being imaged, the resolution is very low. Poplack compared microwave tomography to two other modalities (electrical impedance tomography and near-infrared spectral tomography) and reported the in vivo dielectric contrasts between various tissue types estimated by each method. X- ray images of patients with known lesions were correlated to what were deemed abnormal 18 areas of the microwave images, but the low microwave imaging resolution and mismatch between the compared imaging planes do not convincingly confirm the ability of microwaves to image tumors. Johnson et al. [79] used a time-domain inverse scattering technique to peform imaging of a 2-D MRI-derived breast phantom and reconstructed a 5-mm-diameter tumor of 2:1 contrast with impressive spatial and contrast resolution. The CSI method was applied to a 2-D MRI-derived breast model by Gilmore et al. [80] to reconstruct a large and high contrast tumor model with moderate resolution. Anatomically realistic, 3D, MRI-derived numerical breast phantoms were used in an inverse scattering study by Winters et al. [21] that used GPU computing and a novel spatial basis to solve the largescale problem and produce the first quantitative 3-D images for realistic phantoms of varying density. A color level-set inversion was employed by Irishina et al. [82] that used multiple stages to reconstruct first areas and then properties of each tissue type. Fang et al. [90] studied computational considerations of 3-D microwave medical tomography and proposed a technique using an initial guess of the field distribution for decreasing the computation time of an FDTD forward solution. Gilmore et al. used a wideband frequency selection scheme in [91] to avoid errors due to mutual coupling between antenna elements, and investigated experimental imaging resolution relative to object complexity in [92], Ali and Moghaddam used data from closely spaced frequency points to invert simple pixelated objects in 3D, showing super-resolution for isolated point objects with a limited number of unknowns and multiviews from all eight sides of a cube-shaped object region [93]. Shea et al. used anatomically realistic 3-D breast phantoms to investigate microwave imaging of healthy tissue distributions [84] and to differentially image malignant tumors placed among healthy glandular tissues [83]. Details of these investigations comprise Chapters 3 and Chapter 4, respectively. 19 2.2.3 Inverse Scattering Electromagnetic inverse scattering methods are based on the relationship between field scattering phenomena and the dielectric properties of the media under test. The mathematical relation is typically expressed as an integral equation that can be formed for each unique observation of the field. A set of observations forms a system of equations which must be inverted to solve for the unknown dielectric properties. In this section, the formation of the integral equation is summarized and an inverse problem is posed in the context of nonlinear optimization. 2.2.3.1 Electric field integral equation The electric field integral equation is the mathematical formulation of electromagnetic scattering phenomena typically used in inverse scattering methods. It can be derived from the Helmholtz wave equation. The homogeneous wave equation describes the electric field E at any position r in a source-free background defined by wavenumber k(,. V2E{r) + k2b{r)E(r) = 0 (2.2) The square of the wavenumber in a lossy, non-magnetic medium is complex and can be expressed in terms of its relative permittivity, e r , and effective conductivity, <xefr, at angular frequency u. k2(r) = u}2noeoeT(r) - ju^0(jeS(r) (2.3) The Green's function solution to the wave equation is the field produced by a infinitesimal point source placed at point r' in the background medium. In the following expression, the delta function has units of current density, and the source term for current density, J , in the inhomogeneous wave equation is given by jujfx0J. V2G\r, r') + k2Gb(r, r') = -jufi0 Su(r - r') The subscript u G {a;, y, z} denotes the polarization of the impulse source. (2.4) Combining the vector Green's function solutions for the three polarizations produces a dyadic Green's 20 function Gb, a 3 x 3 matrix which maps an infinitesimal vector source element at r' to a vector electric field at r . Since the superposition of fields is linear, the dyadic Green's function may be used to express the field due to a arbitrary source distribution J(r) in terms of an integral over a region V enclosing the distribution [94]. E(r) = —jujn0 j j J G V , r')J(r')dr' (2.5) v Let an unknown inhomogeneous object of finite extent and defined by wavenumber be placed in the background region. k(r) Interaction between the wave and the object will perturb the field solution to Et = El + Es, where El the original incident field and Es is defined as the scattered field. By the volume equivalence theorem [94], the object may be mathematically replaced by a volume distribution of equivalent current density, Jeq, related to the product between the total field Et and the wavenumber contrast between the object and the background. Jeq(r) = -l-(k2(r)-eb(r))E\r) CJ/i0 (2.6) By the linearity of field superposition, the Helmholtz equation must also be satisfied for the scattered field due to the equivalent source. W2Es(r) + kl(r)Es(r) = - (/c 2 (r) - k2b(r)) E\r) (2.7) Thus, the electric field integral equation relating dielectric contrast to scattered field is found using (2.5). Es(r) = H I G\r, r')E\r') ( j k \ r ' ) - k2b(r')) dr' (2.8) v The goal of inverse scattering is to recover the unknown wavenumber k(r) (or, equivalently, complex permittivity) at each point in the region of interest. The scattered field data, Es(r), can be measured by illuminating the imaging region and taking the difference between field observations at r with and without the unknown object present in the background. However, by definition the total field E 1 over V in the integrand of (2.8) is a function of the unknown wavenumber. Therein lies the nonlinearity of the electromagnetic inverse scattering problem. 21 The relation may be linearized by the Born approximation, which says that if the the magnitude or dimensions of the contrast are sufficiently small then the total field and the incident field are approximately equal. While this is far from valid in medical imaging applications, most of the imaging algorithms of record are based on first-order assumptions and thus begin with the Born approximation. Es(r) « J j J G (r, r')E\r') (k2(r') - k2(r')) dr' (2.9) v The scattered field of the left-hand side of (2.9) may be interpreted as the difference between the field measured at r in the presence of the unknown object and the field simulated at r in the presence of a numerical estimate of the unknown object. E ~ ( r ) - Ecalc(r) « J J J Gb(r, r ' ) ^ ( r ' ) ( k 2 ( r ' ) - k2b{r')) dr' (2.10) v Equation (2.10) is the mathematical formulation underlying most frequency-domain iterative inverse scattering methods. At each iteration, the previous estimate of k(r) becomes the background kb(r) in which to recompute Ecalc(r). A set of measurements forms a set of integral equations that are discretized into a system of linear equations, which may then be inverted to obtain an update, k2(r) — kb(r), to be added to the estimated object profile, kb(r). As the estimate converges to the properties and dimensions of the actual object, the error in the Born approximation converges towards zero. Note that the fields are typically measured by single-polarization antennas, in which case the left-hand side of (2.10) is a scalar field value, Emeas{r) 2.2.3.2 - Ecalc{r). Nonlinear optimization Many inverse scattering methods are specific applications of nonlinear optimization techniques. Thus, it will be instructive to express the problem in the common nomenclature of nonlinear optimization theory [95] in which a scalar-valued cost function / ( x ) is nonlinearly related to a length-K, complex vector x. min/(x) (2.11) 22 The traditional optimization approach to inverse scattering casts a set of equations of form (2.10) as a least-squares problem. M /(*) = g E l M * ) ! ! 2 (2.12) 3=1 In the context of electromagnetic scattering, the function r^(x) corresponds to (2.10), where x is the vector of unknown wavenumber contrasts over some discretization of V into K elements, and rj(x) maps the contrast to the residual scattering in the jth channel, E™eas — E^alc. Collecting the residual values of M measurements into a vector r, / ( x ) = | ||r(x)|| 2 . The Jacobian is the M x K matrix, J , formed by the gradients of the scattering residuals with respect to the unknown contrast vector. J = Vri(x) Vr 2 (x) . . . V r M ( x ) (2.13) An evaluation of the function / ( x ) carries a huge computational cost, as it requires N full, independent electromagnetic simulations of the imaging domain to obtain residual equations for M = N(N + l ) / 2 unique channels, where N is the number of source antennas. In the context of inverse scattering, this function evaluation is referred to as the direct problem or the forward solution. After a forward solution is obtained, the Jacobian can easily be computed from a discretization of the linearized scattering equation (2.10), and thereby the gradient of the function, V / ( x ) = J ^ r , is readily available. The number of rows in J is the number of unique measurements, and the number of columns is the number of discrete unknown regions in the imaging volume. The volume of the breast, the feature size of the interior tissue distribution, and the practical considerations of antenna and test fixture dimensions all tend to dictate that there will be significantly more unknowns than measurements for 3-D imaging, i.e. M K. The Hessian of / is the matrix H, which gives the quadratic dependence of the function on the contrast. In the least-squares formulation, the Hessian is given by, M 2 H H = V /(x) = J J + 5>,(x)V2r,(x) 3=1 (2.14) 23 In the nonlinear scattering relation (2.8), the Laplacian of the residual cannot be computed since the total field is an unknown function of the contrast. While higher-order formulations may be available, the linear scattering approximation dictates that V 2 r 7 (x) = 0 and the Hessian reduces to the first term of (2.14) which is readily available from the Jacobian. This approximation for the Hessian is valid near the solution if the rj(x) or V 2 r j l (x) are small in that neighborhood [95], both of which are true in iterative methods of electromagnetic inverse scattering. Of course, starting with an initial guess near enough to the solution to ensure convergence is an entirely different matter since both rj(x) and V 2 r^(x) are large for high-contrast objects. There are many well-known methods of optimization and inversion [95-97]. As evidenced by Table 2.2, Gauss-Newton (GN) is a common approach in this application. (Note that Gauss-Newton, or Newton-Kantorovich, are equivalent to the DBIM method that is often referred to in electromagnetic inverse scattering [44]). In the Gauss-Newton method, the Newton equation is simplified to JHJp = —JHr (2.15) where an estimate p of the solution is the search direction chosen to update the contrast vector, i.e. Xj + i = Xj + a p , where a is a step length that satisfies descent conditions. Equation (2.15) is ill-posed for 3-D microwave breast imaging, and thus some form of regularization is required in the inversion process. Traditional regularizations such as the Tikhonov method [98] suppress high spatial frequency components of the solution space and often result in excessive smoothing of the estimated dielectric profile. Alternative methods of regularization have been proposed [99-102] and have been shown to improve imaging results [103]. Chapter 3 Microwave imaging of normal breast tissue distributions Electromagnetic tomography of the breast at microwave frequencies is widely studied as an emerging modality for a range of clinical applications from risk assessment to cancer detection. Although early clinical studies have been performed in this area, understanding of the microwave properties of constituent breast tissues is still developing. Controlled studies of both numerical and experimental phantoms have largely been limited by computational and practical considerations, further complicating the evaluation of microwave breast imaging. The purpose of this study is to develop and test the performance of a three-dimensional (3D) microwave tomography algorithm for reconstructing the dielectric-properties distribution of normal fibroglandular and adipose tissue in the breast. We employ state-of-the-art 3D numerical breast phantoms that are realistic in both structural and dielectric properties. The phantoms include examples from each of the four classes of mammographic breast density. Since the properties of these phantoms are known exactly, this testbed serves as a rigorous benchmark for the imaging results. We apply the imaging algorithm to simulated array measurements of the numerical phantoms. A multiple-frequency, bound-constrained, vector field inverse scattering solution is implemented that enables practical computation of the large-scale three-dimensional problem. We exploit knowledge of the frequency-dependent Prom: [84] J. D. Shea, P. Kosmas, S. C. Hagness, and B. D. Van Veen, "Three-dimensional microwave imaging of realistic numerical breast phantoms via a multiple-frequency inverse scattering technique," Medical Physics, vol. 37, no. 8, pp. 4210-4226, Aug 2009. 25 characteristic of breast tissues at microwave frequencies to obtain a parametric reconstruction of the dispersive dielectric profile of the interior of the breast. Imaging is performed on a high-resolution voxel basis and the solution is bounded by a known range of dielectric properties of the constituent breast tissues. The imaging method is validated using a breast phantom with a single, high-contrast interior scattering target in an otherwise homogeneous interior. The method is then applied to simulated scattered field data from a set of realistic numerical breast phantoms of varied fibroglandular tissue density. Imaging results are presented for each realistic phantom and show robustness of the method relative to tissue density. In each case, the distribution of fibroglandular tissues is well represented in the resulting images. The resolution of the reconstructions at the frequencies employed is lower than the dimensions of smallest features of the normal tissue structures, resulting in a smearing of their reconstruction. The results of this study support the utility of microwave imaging to tissue density classification and further investigation of its use in other breast imaging applications. 3.1 Introduction Illumination of the human breast at microwave frequencies has been proposed for several medical imaging applications relevant to breast cancer. Microwave technology offers the potential for a low-cost, noninvasive modality in a non-ionizing range of the frequency spectrum. The dielectic properties contrast between the constituent breast tissues provide the physical basis for microwave imaging. To date, the primary motivation for developing microwave systems has been the need for improved detection of malignant breast tumors. However, the characterization of the relative density of tissue in a healthy breast is another compelling potential application for microwaves. Density characterization is a valuable component in an overall assessment of an individual's risk of breast cancer [104], Dense breast tissue results from a large percentage of epithelial and connective tissue, collectively referred to as fibroglandular tissue. The large microwave-frequency dielectric properties contrast between the dense fibroglandular tissue and less dense adipose tissue [18] offers a very tractable 26 mechanism for density characterization. This application of microwave imaging is less sensitive to resolution performance than early-stage tumor detection. Furthermore, the capability of microwave imaging to identify the spatial features and dielectric properties of the networks of normal fibroglandular tissue structures is an important baseline for the development of tumor imaging methods. High breast density is one of the strongest predictors of breast cancer risk [105] and is fairly common. Approximately half of all women aged 40-49 and a quarter of all women aged 70-79 have breasts that are at least 50% dense, as measured by mammography [106]. Women with high breast density have a four- to six-fold greater cancer risk than women with less dense tissue [107]. The relative risk of dense breast tissue is greater than many traditional risk factors [104] and the prevalence of this risk factor is greater than that of most others [108], An increase in breast density over time has been linked to an increase in breast cancer risk, while a decrease in density is associated with decreased risk [105]. These study findings point to the importance of breast density evaluation and monitoring and their potential role in individualized risk assessment. The most common clinical approach to evaluating breast density involves a qualitative visual interpretation of a mammographic image using the Breast Imaging Reporting and Data System [20] (BI-RADS) that was established to describe the effect of density on diagnostic accuracy [104]. Quantitative measurement has been made possible with computer-aided algorithms [109-111]. However, the accuracy of these techniques is fundamentally limited by the fact that the mammographic image is a 2-D projection of a 3-D phenomenon. Mammographic density is also sensitive to projection angle, level of compression, and X-ray intensity, which limit the use of mammography for monitoring changes in density [104], Breast density can be accurately measured from 3-D magnetic resonance images (MRI) using quantitative methods [112]. However, there are substantial concerns about high costs and limited access to high-quality MRI services. Furthermore, breast MRI is time-consuming, problematic for claustrophobic [113] or very obese patients, and prohibited for women with pacemakers. Thus, neither mammography nor MRI is ideally suited for evaluating and monitoring breast 27 density. Three-dimensional tomographic microwave images may therefore offer advantages in volumetric density measurement for risk assessment and the monitoring of density over time. In this study we investigate the performance of microwave tomography in imaging normal fibroglandular and adipose breast tissues and thereby evaluate the potential for breast density classification. Our investigation focuses on resolving, locating, and estimating the dielectric properties of healthy fibroglandular and adipose tissue distributions. We employ a method of 3-D microwave tomography for imaging the varied and complex dielectric spatial profiles of realistic breast tissue structures. The imaging method is applied to multi-frequency data synthesized from computational electromagnetics simulations of anatomically realistic 3-D numerical breast phantoms. The high-fidelity phantoms [19] contain MRI-derived distributions of normal adipose and fibroglandular tissues having dielectric properties based on large-scale dielectric spectroscopy studies [1,18]. Our imaging technique is based on a multiple-frequency inverse scattering method for reconstructing the 3-D profiles of the parameters of a dispersive dielectric model. In our inverse scattering approach the forward solutions are obtained from Maxwell's equations in the time-domain by the finite-difference time-domain method (FDTD). A set of integral Helmholtz equations describing the electromagnetic scattering are inverted to estimate the unknown dielectric profiles. We implement both the forward and inverse solutions on a grid of voxels of much smaller dimension than the wavelength of illumination so the resolution available to the implementation is not inhibited by the size of the voxels. The electromagnetic simulations and large-scale inverse solutions required by the imaging method are computationally costly on the voxel grid. Computation of the electromagnetic simulations is accelerated using graphics processing unit (GPU) hardware, while the computational challenge presented by the large-scale inverse problem is met with an efficient regularization and inversion scheme. In addition, we constrain the inversion algorithm based on a priori knowledge by restricting the inverse solution to he within the known physical range of dielectric properties of the constituent breast tissues. We have implemented a vector field formulation of the inverse scattering problem in order 28 to evaluate the validity of the scalar field approximation in this 3-D application. The effect of measurement noise on the imaging performance of the system is also studied. In addition to changing the focus of the investigation from tumor detection to the imaging of normal tissue profiles, this work is distinct in several respects from prior microwave breast imaging studies such as those involving synthetic-aperture-radar-based methods [13,15,71,114,115], time-reversal techniques [68], and tomographic approaches to inverse scattering [8,74-79]. Namely, the testbeds employed in this work are state-of-the-art 3-D numerical breast phantoms [19] that are realistic in both dielectric properties and tissue configurations, and for which the known actual properties distributions serve as a rigorous benchmark for the imaging results. Recent large-scale dielectric spectroscopy studies of freshly excised healthy and diseased breast tissue [1,18] suggest that much of the prior microwave breast imaging research assumed too large of a dielectric contrast between cancerous and healthy glandular tissues and too small of a contrast between normal glandular tissue and fat. Also, existing studies of 3-D microwave tomography of the breast have often focused on piecewise homogeneous object domains [74,77,81], which do not effectively capture the spatial complexity of the small features and continuously varying properties distribution of the constitutive tissues. Without the presence of a complex fibroglandular network adjacent to malignant tissue, the resulting phantoms overly simplify the tumor detection problem by creating unrealistically high dielectric contrast and low level of signal clutter from background scattering. This work also differs from our earlier work with realistic numerical breast phantoms [21] wherein the 3-D image was formed on a limited-resolution spatial basis and the phantoms included simple tumor models. Hence, the effect of the tumor could not be conclusively distinguished from the normal fibroglandular breast tissue in the reconstructed images. We proceed in the next section with a discussion of the test phantoms and the acquisition of simulated array measurements. Section 3.3 briefly summarizes the theoretical background for electromagnetic inverse scattering and gives the details of our imaging method. Imaging results for each test phantom are given in Section 3.4, along with evaluations of the scalar 29 field approximation, noise performance, and effects of modeling error. We make concluding remarks in Section 3.5. 3.2 Numerical Breast Phantoms Clinical studies provide the most faithful test domain for microwave breast imaging, but present a challenge in the interpretation and registration of the results because the actual dielectric properties of the breast interior are unknown. Although some initial clinical results of microwave breast imaging have been reported [22,89,115], clinical studies do not readily lend themselves to rigorous validations of imaging performance. Both numerical and experimental laboratory testbeds permit relatively straightforward performance evaluations since the dielectric properties distributions are known exactly or approximately. Experi- mental phantoms may be constructed from synthetic materials that accurately mimic tissue dielectric properties [116], but the diverse and complex structural distributions of the breast are difficult to mimic in an experimental phantom. Accordingly, many prior experimental laboratory studies have been restricted to models consisting of arrangements of homogeneous cylindrical [59,76,78] or spherical [15,56,114] targets which do not accurately mimic the complexity of the breast. Numerical phantoms offer more flexibility in capturing the structural realism of breast tissue. The dimensionality and complexity of numerical breast models have been limited in the past [13,74,77,78,80,81,117,118] by practical computing constraints and the lack of availability of 3-D anatomically realistic models. However, such limitations have been overcome with parallel computing strategies and the recent development of realistic MRI-derived numerical breast phantoms [19]. The 3-D numerical breast phantoms that serve as testbeds in this study are adapted from phantoms in the online repository developed and maintained by the University of Wisconsin Computational Electromagnetics Laboratory (UWCEM) [19,119]. The repository phantoms are derived from MRIs of healthy breasts, and thus convey the realistic shape and internal adipose (fatty) and fibroglandular tissue structure of the pendant breast. The range of dielectric properties assigned to each tissue type (adipose, transitional, or fibroglandular) 30 are derived from large-scale dielectric spectroscopy studies of freshly excised breast tissue specimens [1,18]. The repository phantoms are defined on a uniform 0.5 mm voxel grid and include a homogeneous skin layer and a chest wall. In this section we describe the modifications made to a subset of repository phantoms to adapt them for this study. We also describe the method of acquiring simulated scattered signals from the customized 3-D numerical breast phantoms. 3.2.1 Test Phantoms The dielectric properties assumed in the UWCEM numerical breast phantoms vary over the full range spanning the lower and upper bounds of all reported specimens in a study of ex-vivo tissues [18]. Cole-Cole curves reported [19] for those bounds as well as for the 25th, 50th, and 75th percentile values of adipose and fibroglandular tissues are used to generate a continuum of Cole-Cole curves that model a full range of dispersion properties of tissue in the breast. Over the frequency range of interest here, 0.5 to 3.5 GHz, a single-pole Debye model sufficiently models the frequency dependence of the complex permittivity of the tissues. e(u) e0 Ae as = £oo + t - — h1 + JUT JU£0 (3.1) The Debye parameters, e^,, Ae, and as, of the phantom profiles are chosen to fit the singlepole Debye model to the Cole-Cole model over the 0.5 to 3.5 GHz frequency range. The resulting Debye models for representative tissues are plotted in Fig. 3.1 and their parametric values are given in Table 3.1. The assumed dielectric properties of the skin layer [120] and the coupling medium (similar to vegetable oil) are also provided in Table 3.1. The relaxation time constant of the single pole, r, is set to 15 ps for all materials and is assumed to be spatially invariant to simplify our electromagnetic simulation code. The quality of the fit of Eqn. (3.1) to the Cole-Cole models is negligibly affected by this assumption. Additional modifications are made to the selected repository phantoms to reduce the computational burden of the imaging procedure. The 1.5-mm-thick homogeneous skin layer of each phantom is replaced with a 2.0-mm-thick homogeneous layer so that it can be better modeled on the 2.0 mm grid used in the imaging algorithm. In addition, we remove the 31 -g* 3.5 ^ ^ ^ Fibroglandular (minimax) — — Fibroglandular (75th pctile) - — - • Fibroglandular (25th pctile) <£L 3 Adipose (75th pctile) Adipose (25th pctile) >s I 0) 40(HO«O<3O«.O-8O-G-l)-O.O-9-O-8-O-S-O-6-O^OEK).i>oe>o-( CL < >1) 30 Adipose (maximin) 2.5 - o - - Dry Skin § 2 T3 §1.5 20 I •c 8 0.5 10 LU 3.5 1.5 2 2.5 Frequency (GHz) 0^ 0.5 (a) 1.5 2 2.5 Frequency (GHz) 3.5 (b) Figure 3.1 Frequency dependence of the single-pole Debye models for representative tissues used in the numerical breast phantoms, (a) relative permittivity, e r , and (b) effective conductivity, o ^ . Table 3.1 Debye parameters (infinite and delta relative permittivity, e<» and Ae, and static conductivity, a s ) of the materials modeled in the numerical phantoms (valid from 0.5 to 3.5 GHz). Material (percentile) ^oo Ae a s (S/m) Adipose tissue (maximin) 2.28 0.141 0.0023 Adipose tissue (25th) 2.74 1.33 0.0207 Adipose tissue (50th) 3.11 1.70 0.0367 Adipose tissue (75th) 4.09 3.54 0.0842 Fibroglandular tissue (25th) 16.8 19.9 0.461 Fibroglandular tissue (50th) 17.5 31.6 0.720 Fibroglandular tissue (75th) 18.6 35.6 0.817 Fibroglandular tissue (minimax) 29.1 38.1 1.38 Skin tissue 15.3 24.8 0.741 Coupling medium 2.6 0.0 0.0 32 slab regions of the chest wall model from the phantoms. Finally, the dimensions of the computational domain are adjusted to accomodate the phantom and a surrounding antenna array. The American College of Radiology (ACR) defines four classes of breast composition based on mammographic breast density: mostly fatty, scattered fibroglandular, heteroge- neously dense, and extremely dense [20]. We will henceforth refer to these as Class 1, Class 2, Class 3, and Class 4, respectively. We select one representative phantom of each of the four ACR classes from the repository. The repository ID numbers of the selected phantoms are 071904, 010204, 062204, and 012304. The exterior surface of the 3-D model of each phantom is depicted in Fig. 3.2. Orthogonal 2-D cross-sections of the Ae profile of the interior of each phantom are shown in Section 3.4 (see Figs. 3.9, 3.10, 3.11, and 3.12). The e ^ and o s Debye parameter profiles are similarly distributed, taking on the ranges of values given in Table 3.1. In addition to the four phantoms with realistic heterogeneous interiors, our testbeds include a simpler phantom with a single spherical inclusion in an otherwise homogeneous interior. This phantom is employed to validate the method and gain insight into the performance and limitations of the imaging system. The phantom is constructed by replacing the interior region of the Class 2 test phantom with a homogeneous medium having the mean properties of transitional tissue. Tissue is categorized as transitional when its dielectric properties fall between the 75t/l percentile for adipose tissue and the 25 t/l percentile for fibroglandular tissue (see Table 3.1). A 3-cm-diameter, homogeneous sphere having the 50th percentile properties of normal fibroglandular tissue is added to the interior of the phantom to provide a simple, high-contrast scattering target. Transections of this phantom are shown in Section 3.4 (see Figs. 3.6 and 3.7). 3.2.2 Scattered-Field Data Acquisition Each test phantom is surrounded by a cylindrical array of 40 electrically short dipole antennas, as shown in Fig. 3.2. Each dipole arm is 6 mm long and the excitation gap is 33 (a) (b) (c) (d) Figure 3.2 3-D numerical breast phantoms and dipoie antenna arrays, (a) Class 1, (b) Class 2, (c) Class 3, (d) Class 4. 2 mm. The cross-section of each antenna is 2 mm by 2 mm. The antennas are evenly distributed on five elliptical rings of eight antennas each, with adjacent rings rotated by 22.5°. The dimensions of the array are customized to fit each test phantom. The five rings are evenly spaced between the posterior and anterior coronal planes of each phantom. The dimensions of the major and minor axis of the array's elliptical cross-section are chosen such that the array conforms to the phantom with a minimum spacing of 1 cm between each antenna element and the skin surface. 34 We use the FDTD numerical method [121] to simulate array measurements of the numerical breast phantoms. In each simulation, one antenna in the array is sourced by exciting the electric field in the gap between the dipole arms.The source function is a modulated Gaussian pulse having a 3 dB bandwidth from 0.87 to 3.75 GHz.The 3-D computational domain is composed of 0.5 mm cubic grid cells and is terminated with uniaxial perfectly matched layer (UPML) boundary conditions for dispersive media [121]. The antenna measurements are observations of the co-polarized field component in the source gap. The time-domain scattered fields recorded at every antenna in the array are converted to phasors at the frequencies of interest by discrete Fourier transform. The four frequencies used in our reconstructions are 1.0, 1.5, 2.0, and 2.5 GHz. In addition, we simulate measurements on the simple phantom and Class 2 phantom after down-sampling them to a uniform 2.0 mm grid. The down-sampled grid matches the resolution of the grid used by the imaging algorithm described in Section 3.3. Imaging the down-sampled phantoms allows us to remove, and thereby identify, the effects of modeling error on the imaging results. Acquiring simulated measurements in the same computational domain as is used by the inverse scattering algorithm is often referred to as an "inverse crime." The difference in resolution between the 0.5 mm data-acquisition grid and the 2.0 mm imaging grid leads to antenna modeling mismatch and simulation artifacts such as increased numerical dispersion. We reduce the mismatch error using a calibration step [8,21]. Array measurements are made in a homogeneous immersion medium in both the data acquisition and imaging domains. Ratios of the array-only phasor measurements through each transmitreceive channel are applied as a complex correction factors to the acquired data. In this study, the calibration measurements of the data-acquisition grid and the imaging grid nominally agree within about 5%. 35 3.3 Imaging Method Multiple scattering interactions in a heterogenous target region result in nonlinearity of the electromagnetic scattering. Thus, an iterative method of nonlinear optimization is used to estimate the dielectric profile based on the measured scattered fields. Each iteration of the algorithm consists of a set of electromagnetic simulations and an inversion of a linear approximation to a system of scattering equations. The unknown dielectric properties are estimated simultaneously at multiple frequencies using a parametric model of the complex permittivity over frequency. This section describes electromagnetic scattering from a heterogeneous dielectric target, the nonlinear optimization method used to image the dielectric profile, the details of the implementation, and a discussion of the key assumptions made in our solution. 3.3.1 Electromagnetic Inverse Scattering and the Distorted Born Iterative Method Electromagnetic inverse scattering methods operate on a set of field measurements of a penetrable, unknown target region, V. A known source illuminates the region and the resulting field is measured at one or more observation points outside of V. An estimate of the unknown dielectric profile, e(r), in V is reconstructed based on a relationship between the fields scattered from the region and the dielectric profile within that region. For a measurement at location r at a given frequency co, this relationship can be expressed by an integral equation [16], E s (r) = E'(r) - E { (r) = u2n f G6(r|r')EV) [>(0 - eV)] dr ' (3-2) Jv where E s is the scattered electric field, E* is the total field, E® is the incident field in the presence of the background permittivity, e 6 (r), and G& is the dyadic Green's function for the background. The profile of the dielectric contrast over V is formed by the difference between the dielectric profile of the unknown region, e(r), and the background profile, c h (r). 36 A set of field measurements of a target space thus presents a system of equations of the form of Eqn. (3.2) in which the unknown quantity of interest is the dielectric contrast, e(r) — c 6 (r). However, there are several well-known complications in obtaining a solution to the system. The number of unknowns in V is often much greater than the number of measurements, yielding an under-determined system without a unique solution. The Green's function may not be available analytically when the background is not a canonical region such a homogeneous space or a half space. Furthermore, the total field within V is unknown and is a function of e(r), making the system nonlinear in the unknown contrast function. We employ the distorted Born iterative method [16] (DBIM) to obtain a solution to the nonlinear problem. At each iteration of the DBIM, the total field within the actual profile is approximated by the field in an estimated heterogeneous background profile. The fields within this background profile, E b , replace E 4 in Eqn. (3.2). The approach requires computation of the fields at the antennas and inside V for each iterate of the background profile. These field computations are collectively referred to as the forward solution. The DBIM is equivalent to a Gauss-Newton approach [44] to nonlinear least-squares optimization problems [95]. At each iteration, a system of scattering equations is constructed from the forward solution and the measurement data. The system is then inverted to find an approximate solution to the contrast between the current estimate of the background profile and the true phantom profile. This part of the algorithm determines the update to the background profile and will be referred to as the inverse solution. The DBIM algorithm alternates between forward and inverse solutions, updating the background profile, c 6 (r). at each iteration until convergence is reached in the minimization of the residual scattering. 3.3.2 Forward Solution We use the FDTD method in the forward solution to obtain 3-D, full-wave field solutions for an estimated background profile. FDTD is a useful method for the simultaneous aquisition of multi-frequency vector field data over the full bandwidth of interest. The sample spacing of the forward solution grid must be dense enough to limit numerical dispersion in 37 the FDTD simulation, but sparse enough to limit the computational cost of the simulations. We select a uniform 2.0 mm sample spacing for the forward FDTD solutions to balance this trade-off. Since multiple FDTD simulations of the numerical phantoms are still computationally costly on a 2.0 mm grid, we run the simulations on a GPU-based hardware accelerator to achieve a feasible imaging time. Alternative numerical techniques for reducing the computation time of the forward solution have been presented elsewhere [50,122,123]. An independent simulation is performed for each antenna in the array. The time-domain background field, E 6 . is recorded at all other antennas and at every voxel within the reconstruction region, V. The measurements are then converted to the phasor domain by discrete Fourier transform. The background fields are the fields associated with the current estimate of the unknown profile which is used as the background profile at each iteration of the DBIM. Although the fields observed in the FDTD forward solution are total field quantities, we refer to them in the inverse scattering context as background fields to distinguish them from the unknown total field E'* that existed in the actual profile during measurement. As there are no analytical Green's functions available for the heterogeneous background profiles estimated at each iteration of the DBIM, the Green's functions must be computed. For a unit source at a given antenna feed point r in V, the Green's function defines the field received at a particular antenna at rm. This quantity can be found using the principle of reciprocity [124] and the background field measurements already obtained by the forward solution. A particular advantage of computing the Green's function from the existing forward solution is that the exact radiative behavior of the antennas, including the mutual coupling, is built into the Green's functions and therefore no correction factors are necessary - with the exception of the calibration for antenna modeling error described in Section 3.2. The background Green's functions can be calculated from the source current and the incident field measurements, E*, based on a simple relationship [125]. Here, only the z-directed tensor elements of the background Green's functions are computed since the dipole antennas are z-polarized. For a ^-directed source current Iz of length L at the mth transmitting antenna 38 the effective Green's function tensor is calculated at each r 6 V as 0 0 0 0 0 0 Eix{ r\rm,uj) Ely(r\rm,u) (3.3) £*(r|rmi m ,w) The first two rows of the tensor are zeroed since no x- or y-directed sources are used. Equation (3.3) is used in the vector field formulation of the integral in Eqn. (3.2). When the scalar field approximation [21,57,74] is used, cross-polarization scattering effects are assumed to be negligible. Thus, Elx and Ely are set to zero in Eqn. (3.3) reducing the Green's function tensor to a single element. This approximation is employed to reduce the computational requirements of the inverse solution and reduce the memory required to store the forward solution. We note, however, that there is no scalar approximation employed in our full-wave 3-D FDTD computation of the forward solution itself; we use the approximation only in the construction of the inverse system. 3.3.3 Frequency-Dependent Model of Unknown Dielectric Properties The imaging algorithm seeks to estimate the complex permittivity profile of the unknown region over a set of discrete frequencies. Rather than reconstruct the profile independently at each frequency, it is advantageous to enforce an assumed frequency dependence of the tissue properties. This approach can reduce the number of degrees of freedom in the system, thereby restricting the size of the solution space and easing the computational burden of the imaging algorithm. In prior studies, a general dispersive formulation has been presented [76] along with sample solutions using linear and logarithmic two-parameter dispersion relations. A single-pole Debye model of dielectric dispersion has also been employed [21]. If the chosen model does not adequately capture the actual frequency dependence of the tissues, the model error will corrupt the resulting image. Except over narrowband ranges, an assumption of frequency independent dielectric properties is demonstrably inaccurate [1,18,63,120]. A linear dependence does not capture the dispersive effects at microwave 39 frequencies due to the water content in biological tissues. Instead, dielectric models which include the relaxation behavior explicitly, such as the Debye and Cole-Cole relations, are naturally adept at capturing this frequency dependence. We elect to use the single-pole Debye model as it employs one less parameter than the Cole-Cole model, it can be incorporated into FDTD in a straightforward manner, and it fits the dispersive behavior of breast tissues well over the frequency range of interest [126]. Furthermore, the single-pole Debye model uses the same number of parameters as would be needed to model both permittivity and conductivity with first-order models of frequency dependence. We assume the relaxation time constant parameter, r, to be invariant to reduce the number of unknowns and to keep the dispersion model linear with respect to the unknown parameters. The time constant is fixed in the inverse solution to the same value used by the data acquisition simulations and the forward solution. 3.3.4 Linear System of Scattering Equations The first-order Born approximation to Eqn. (3.2) is linear in the unknown contrast function, e(r) — £h(r). We will denote the contrast function over r G V more compactly as (S{s(r)j. To achieve a simultaneous solution of the multiple frequency system, the complex permittivity variables in Eqn. (3.2) are replaced by the frequency dependent Debye model of Eqn. (3.1) with known time constant r. The relationship between the field and the contrast function is then linear in the remaining three parameters of the Debye model: e^, Ae, and a s . The contrast functions of these three parameters, 5{e 00 (r)}, (5{Ae(r)}, and 5{crs(r)}, are the new unknowns over r e V. Since these parameters are all real-valued, Eqn. (3.2) is split into a pair of real and imaginary equations so that the solution space is limited to real values. For an array of N antennas there are N2F total frequency-domain measurements, where F is the number of discrete frequencies to be included in the solution. We discard redundant data from reciprocal channels to retain a single unique measurement for each transmit-receive antenna pair. We also discard monostatic data since these observations include the source 40 field. Each of the remaining MF measurements, where M = N(N — l)/2, yields a pair of real and imaginary vector equations. The resulting set of 2 M F equations are then discretized by a Riemann sum over the reconstruction region, V. We assume knowledge of the location and properties of the skin region and restrict V to the interior breast volume. The region is discretized using the same uniform 2 mm voxel basis used in the forward solution. The set of K voxels within V is vectorized so that the unknown contrast function for each Debye parameter forms a K-by-1 vector. These three vectors are collected into a single 3/T-by-l vector, x. For each equation, the remaining elements of the summand form a row of a 2MF-hy-3K matrix, A. The differences between the measured fields and the fields computed by the forward solution are collected into a vector of residual scattered fields, b. The resulting linear system, A x = b, is then structured as follows for each channel d between a transmitter at r m and a receiver at r n , 3R{Ef} 3{Bf} 3{Ef} tf(Ae) K{BS} where the 5R{} and (3.4) 6(<ra) operators denote the real and imaginary parts of the complex argument. In Eqn. (3.4), CpM K [6?(wi) ••• bjciui) = cp{cuF) [bf(uF) ••• bdK(tuF) El(r n |r m ,wi) - El(r n |r m ,Wi) E^ = ££(r„|r m ,o; F ) - £^(r„|r m , uF) 41 where bf(u>) = lo2/j,£0Gb(rn|rfc,o;)E1(rfc|rm,w), Ca(w) = 3.3.5 c^u) = 1, cA(u) = (1 + jut) 1 , and (ju£0)~\ Inverse Solution The linear model of the scattering system in Eqn. (3.4) is ill-posed. It is also highly under-determined on the selected voxel basis since 2 M F 3K. Our inversion strategy uses a conjugate gradient method to find a regularized, approximate solution to the system of normal equations, A 3.3.5.1 T A X = A T b . (3.5) Inexact N e w t o n Step As previously noted, the DBIM is equivalent to a Gauss-Newton method of optimization. An approximation, x, to the inverse of Eqn. (3.5) is then an inexact Newton step used to update the estimate of the background profile £b(r) over r e V. We use the conjugate gradient for least-squares (CGLS) algorithm to find an approximate inverse of (3.5). The CGLS minimizes the residual norm of the system, ||Ax — b|| 2 , over the Krylov subpace [97] span{ A T b, ( A T A ) A r b , . . . , ( A T A ) 3 ^ " 1 A T b }. The CGLS algorithm eases the memory requirements of the inverse solution considerably. Since only two matrix-vector products of A need to be computed per iteration, the complete normal matrix A T A never needs to be constructed and stored in memory. However, the matrix A on the 2.0 mm voxel basis often cannot be accomodated in memory either, in which case the matrix-vector products must be computed in blocks, with each block of A having to be recomputed twice per CGLS iteration. 3.3.5.2 Regularization In prior work [21], a Tikhonov regularization method was applied to a conjugate gradient inversion algorithm, using a series of trial solutions to determine the regularization parameter 42 [96]. This approach is reasonable for such reduced dimensionality approaches [21], but is far too computationally costly for large-scale linear systems. We avoid the need for any trial solutions by executing the CGLS algorithm without applying any regularization term to Eqn. (3.5) and terminating the algorithm prior to convergence [78]. The early termination method takes advantage of the self-regularizing properties of the CGLS [96]. This strategy is indispensible in light of the memory restrictions and resulting computational inefficiency faced by the CGLS algorithm. 3.3.5.3 Constraints Unconstrained methods of optimization will operate over the full solution space. However, at minimum we know that any nonphysical solution (e < eD and a < 0) is invalid. We can take further advantage of data on the bounds of breast tissue properties [1,18] and the observation that Ae > 0 in biological tissues. Applying these simple bound constraints to reduce the span of the solution space can aid in the solution of an under-determined, ill-posed system. In addition, we observe that the parameters £oo> Ae, and a s of the Debye models of breast tissues are roughly proportional in the sense that they are all at the low end of their range in adipose tissue and at the high end in fibroglandular tissue. We propose that bound constraints on these spatially correlated parameters can therefore be made interdependent. Specifically, we bound e ^ by its full range (see Table 3.1) but restrict the bounds on Ae and a s based on the largest and smallest ratios between the values of these two parameters and those of £oo over the adipose and fibroglandular entries in Table 3.1. At each DBIM iteration, current estimates of eoo are multiplied by these ratios to obtain new bounds for Ae and a s at each voxel of V. The tightened bounds on Ae and a s are intended to promote the spatial correlation of the Debye parameters that would be expected of accurately reconstructed profiles. Often bound constraints are enforced after completion of the optimization algorithm by simply projecting the inverse solution onto the allowed solution space. This approach can 43 undermine the optimality of the inverse solution in a non-sparse system where the elements of the solution have non-negligible interdependencies. A nonlinear transform to enforce the bounds within the contrast update has been proposed elsewhere [127]. In this work, we accomplish the constraint within the optimization routine using a projected-restarted method [128]. In this approach, the inverse solution is projected onto the allowed solution space after the initial CGLS algorithm terminates. The CGLS algorithm is then restarted using the projected solution as the initial guess. The new solution is then projected and CGLS algorithm is run again. The projecting and restarting is performed a predetermined number of times during each inverse solution and results in a bounded solution of improved optimality. 3.3.5.4 Termination Terminating the CGLS algorithm too early over-regularizes the solution and reduces the spatial resolution. Terminating too late under-regularizes the solution and allows the magnification of noise errors to corrupt the solution. We use an L-curve method to determine an appropriate termination condition for the algorithm. For ill-posed systems, a plot of the solution norm, ||xj|| 2 , versus the residual norm, ||Ax,- — b|| 2 , over the CGLS iterations j = 1, 2, 3... will depict an 'L' shape as the step length of the contrast update begins to grow faster than the scattering residual is reduced. Our criterion for choosing a solution is the iteration at which the plot has maximum curvature. This selection achieves a reduction in the residual while preventing small singular values from erroneously growing the solution. The method is illustrated by a sample L-curve in Fig. 3.3. We observe that for our numerical phantoms this heuristic generally selects the fourth or fifth CGLS iteration at every DBIM iteration. We fix the projected-restarted constraint method to three restarts of the CGLS algorithm. 44 E o "cc •g "w <D DC CO _i O o CGLS Solution Norm Figure 3.3 At a sample DBIM iteration, the norm of the CGLS solution is plotted vs. the norm of the CGLS residual at each iteration of the algorithm. The knee of the 'L' curve (indicated by the ' O ' marker) locates the iteration at which to terminate the CGLS algorithm. 3.3.5.5 Background Update The estimates of the contrast functions found by the inverse solution, ^ ( ^ ( r ) ) , (5(Ae(r)), and £(<rs(r)), are combined into an estimate of the complex permittivity contrast, £(e(r)), according to Eqn. (3.1). The contrast estimate is added to the K sample points of the ith background profile to create the new background for the next DBIM iteration, 4+i(r*) = 4 ( r f c ) + S(i(rk)), k = 1,..., K (3.6) The next iteration commences with the forward solution of the new background profile. 3.3.6 Assumptions and Initial Conditions Imaging performance is quite sensitive to the termination heuristics and regularization techniques employed by the optimization algorithm. Allowing the DBIM to run for too many iterations can allow noise and other errors to corrupt the image. Too few DBIM iterations restricts the ability of the method to overcome the nonlinearity of the problem. In this 45 work we determine convergence of the DBIM algorithm based on the norm of the residual scattering, 11112, computed after each forward solution. This quantity is a measure of the similarity between the fields scattered by the test object and the fields scattered by the reconstructed object. We define convergence to be the ith iteration at which ||bi_i|| 2 — ||bj||2 is less than 1% of ||bi|| 2 . The selection of this threshold is based on our observations of marginal improvement in image quality for additional iterations at the expense of increased influence of noise and model error. The initial guess provided to the DBIM routine is also crucial to the imaging performance of the DBIM. While convergence behavior can be relatively insensitive to initial conditions, the fidelity of the reconstruction can be quite sensitive to the initial guess. A reasonable choice for a homogeneous initial guess is to use the dielectric values of the often predominant adipose tissue. An initial reconstruction step can also be used [21,87] to find an improved homogeneous initial guess or to find a coarse heterogeneous initial guess. In the contrastsource method back-propagation of received fields can be used to obtain an initial estimate of the incident fields in the reconstruction region [60]. Such a profile of field magnitudes could be used to generate an initial profile of dielectric properties for the DBIM. In this work we have made use of a priori knowledge of the true average dielectric properties of the interior of the test phantoms, and use the volumetric average of the Debye parameter profiles as a homogeneous initial guess. The average Debye properties of each phantom are given in Table 3.2. Methods of estimating the average properties of the imaging volume have been reported [21,129], The scattering from the skin region and chest wall represent a primary source of interference in microwave breast imaging methods. There exist techniques for the cancellation of the skin response [71] and for imposing known boundary conditions at the chest wall [130]. Several approaches have been proposed for identifying the location and thickness of the skin layer [131,132]. The dielectric properties of the skin have been extensively documented [120] and there is potential for direct measurement in a clinical setting via dielectric spectroscopy [133]. These approaches to estimating the skin region may be superior to the 46 Table 3.2 Debye parameters (infinite and delta relative permittivity, e ^ and Ae, and static conductivity, crs) of the homogeneous initial profile assumed by the DBIM algorithm for each realistic phantom. Phantom £00 Ae a s (S/m) Class 1 3.62 2.66 0.0586 Class 2 3.97 3.12 0.0697 Class 3 7.42 9.58 0.219 Class 4 7.88 10.3 0.236 inclusion of the skin region in the reconstruction region of an inverse scattering technique, since its thickness and distinct surfaces are poorly reconstructed by the frequencies typically used in microwave tomography. We therefore assume in this work that the skin surface, thickness, and properties are known a priori, are included in the background profile, and are unchanged at each iteration of the DBIM. Despite this assumption, there remains a substantial difference between the fields scattered from the 0.5 mm resolution skin region in the data aquisition model and the downsampled 2.0 mm resolution skin region of the forward solution models. This error arises from both the stairstepping of the downsampled skin region and the additional numerical dispersion of the larger FDTD grid cells. 3.4 Results and Discussion We apply the DBIM algorithm to the data acquired from the set of numerical breast phantoms described in Section 3.2. We empirically validate the imaging method using a piece-wise homogeneous breast model with a single spherical inclusion. In addition, we investigate the imaging performance in realistic phantoms with respect to a number of important issues. Namely, we characterize the imaging artifacts arising from numerical and discretization errors in the forward model, examine the validity of the scalar field approximation, and evaluate the sensitivity of the imaging system to noise. 47 Each data set contains multistatic channel measurements from the 40-element array, excluding non-reciprocal and monostatic channels, at each of four frequencies (1.0, 1.5, 2.0, and 2.5 GHz). The frequency list is chosen to satisfy several trade-offs. The highest frequency is the primary design parameter in determining system resolution and must be balanced against the increase of the computation cost of the FDTD forward solution with frequency as well as the decrease in the stability of the inverse solution due to decreasing tissue penetration depth. Inclusion of lower frequencies help stabilize the inverse solution and increases the amount of information relative to measurement noise at the expense of increasing the size of the linear system. Finally, the bandwidth of the measurement must reflect practical antenna considerations, such as the trade-off between the antenna dimension and the upper and lower operating frequencies. In general, the selection of frequencies for practical imaging systems will be a function of considerations such as signal-to-noise ratio, immersion medium, array population and element size, and the conditioning of the inverse system. Reconstructions of dielectric profiles are often imaged as the complex permittivity at a particular frequency. Since we enforce the frequency dependence of the Debye model and solve for the contrast function parametrically, we elect to image the Debye parameters (e^, Ae, and a s ) directly. The resulting images provide complete multiple frequency information of the reconstruction. The diagnostic potential of dispersion profiles has previously been noted [76]. Presenting 3-D images of complex, heterogeneous objects in 2-D cross-sections is a challenge, particularly when the the reconstruction blurs detailed structures into adjacent grid cells. Simple and often-used 1-D transections of the reconstructions are even more limited in demonstrating the agreement between complicated 3-D profiles. Hence, we begin by imaging the simple case of a large, high-contrast, spherical scatterer placed within a homogeneous breast interior. The resulting reconstruction is presented in both 1-D and 2-D cuts to validate the imaging method. We then illustrate our 3-D reconstructions of realistic tissue 48 distributions using three orthogonal 2-D cross-sections through the parametric reconstructions. All 2-D cross-sections are selected to bisect the extent of the breast phantom along the Cartesian axes in the coronal, sagittal, and axial planes. We introduce two scalar metrics in support of quantitative comparisons between reconstructions. The first metric is an error measure based on the normalized root mean square. It is designed to represent the error relative to the actual profiles of all three Debye parameter estimates with a single metric. The proposed error metric is given as, e2 = 7 f ( Ae-Ae ^00 ^00 I1 Ae £00 2 + 2 CT a — crs (3.7) where K, the total number of voxels in V, is the length of each parameter vector. The hat notation, e.g. e, denotes an estimated quantity. The second metric is intended as measure of the qualitative similarity of two spatial profiles and is insensitive to the scaling of the profiles. The two 3-D profiles to be compared are vectorized and the metric is computed as the cosine of the angle between the vectors. Let 4> be the angle between two Debye parameter profiles represented by vectors pi and p 2 . Then the similarity measure is: COS ((f)) = (P1P2) (3.8) 1IP1II2 IIP2II2 When a reconstruction is referenced against an exact profile, the exact profile is obtained by down-sampling the phantom to the 2.0 mm grid so the dimensions of the phantom match the dimensions of the reconstruction. This case will be denoted by subscripted angle <f>e. The case of one reconstruction compared against another reconstruction will be denoted by subscripted angle <fir. The similarity measure of Eqn. (3.8) can be evaluated for a particular Debye parameter, e.g. p = A e , or each of the parameter profiles can be vectorized and concatenated as p = [e^, A e 1 .TIT to evaluate the similarity of the complete frequency- dependent complex permittivity profiles. 3.4.1 Simple Phantom The data used in the reconstruction of the simple phantom is noiseless and the DBIM algorithm employs the scalar-field approximation. We reconstruct images from two data 49 sets: the first is obtained from the phantom of 0.5 mm resolution and the second is obtained from the down-sampled phantom of 2.0 mm resolution. The latter is the inverse crime case, in which the forward model grid is identical to the data acquisition grid. The 3-D reconstructions of the 2.0 mm and 0.5 mm resolution cases are illustrated by the 2-D cross-sections of Fig. 3.4 and 3.5, and the similarity metrics with respect to the exact phantom are cos(0 e ) = 0.995 and cos(0 e ) = 0.967, respectively. Substantial imaging artifacts caused by the modeling error of the coarser resolution grid used in the forward solver are clearly observed in Fig. 3.5. Figures 3.6 and 3.7 show 1-D transections of the exact phantom and the corresponding reconstructions taken through the true center of the spherical inclusion. These plots confirm the accurate location and dimension of the reconstructed inclusion. The transections of Fig. 3.7 illustrate the reduced fidelity of the reconstruction caused by numerical and modeling mismatch, while those of Fig. 3.6 show the overestimation of the infinite permittivity, e ^ , as well as some artifacts in the static conductivity, <js, even when numerical and modeling errors are suppressed. 50 (c) (f) (i) Figure 3.4 3-D, inverse crime reconstruction of the simple phantom with spherical inclusion, shown in coronal (x ~ y, top row), sagittal (y — z, middle row), and axial (x — z, bottom row) cross-sections of the reconstructed Debye parameters, (a)-(c) e^, (d)-(f) Ae, (g)-(i) a, (S/m). 51 .30 Oi 2 1 4I 8 f 1ol 20 « * 10 12l 14l 0. 2 1 4 1 6 1 8 1 1ol (a) 0| 2 30 4l 6 1 8 20 » •# •» -m 12l 14l0 2 4 6 8 10 12 0 2 4 6 8 10 12 40 10 I0 0 (b) 0.6 1 0.4 10l 12l 14l 0.2 0 2 4 6 8 10 12 0 (g) 2 10 0.8 40 30 1 1 (d) 20 1 2 16 20 \ 6 8 8 10l 0 10l 0 0 (e) (h) 40 30 20 41 6 10 0 10l 0 Figure 3.5 3-D reconstruction of the simple phantom with spherical inclusion, shown in coronal (x — y, top row), sagittal (y — z, middle row), and axial (x — z, bottom row) cross-sections of the reconstructed Debye parameters, (a)-(c) e ^ , (d)-(f) Ae, (g)-(i) crs (S/m). 52 Infinite Permittivity 5 7 9 Delta Permittivity 5 11 x (cm) 5 7 9 4 6 z (cm) 9 11 5 x (cm) 11 5 7 8 2 9 4 6 z (cm) 7 9 11 x (cm) 3 5 y (cm) y (cm) 2 7 Static Conductivity 7 9 11 y (cm) 8 2 4 6 z (cm) Figure 3.6 Orthogonal transections of the inverse crime reconstruction of the simple phantom with spherical inclusion. The dashed line is the transection of the phantom, and the solid line is the transection of the reconstruction. The transections are taken through the center of the true position of the inclusion, x = 6.6 cm, y = 6.0 cm, z — 4.8 cm. 53 Infinite Permittivity Delta Permittivity Static Conductivity 30r 20 10 f A J ° 3 \f*A 5 7 9 11 3 x (cm) 3 5 7 9 x (cm) 11 3 y (cm) 30 20 10 0 5 7 9 11 -Av -AH z (cm) 7 9 11 x (cm) 3 y (cm) 40 30 20 10 0 z (cm) 5 5 7 9 11 y (cm) 0.9 0.6 0.3 0 z (cm) Figure 3.7 Orthogonal transections of the reconstruction of the simple phantom with spherical inclusion. The dashed line is the transection of the phantom, and the solid line is the transection of the reconstruction. The transections are taken through the center of the true position of the inclusion, x = 6.6 cm, y = 6.0 cm, z = 4.8 cm. 54 3.4.2 Realistic Phantoms Noiseless data from the simulation of the four MRI-derived phantoms on the 0.5 mm resolution grid are used to reconstruct each of the phantoms using the DBIM algorithm with the scalar field approximation. The convergence of the norm of the residual scattering fields is shown in Fig. 3.8 for each phantom. The convergence condition described in the previous section requires between four and six iterations, as shown the figure. In addition, Fig. 3.8 plots the similarity metric between the reconstructed and exact profiles, cos(0 e ), versus DBIM iteration. This metric is seen to increase convergently as the residual scattering decreases. A positive slope of the cos(<fie) curve at termination suggests that the reconstruction could be further improved with additional iterations. Visual agreement between the exact phantom profiles and the 3-D reconstructed profiles is evident in the cross-sectional images of Figs. 3.9, 3.10, 3.11, and 3.12. (Only Ae is presented in these figures, since the three reconstructed Debye parameter profiles are noted to be highly correlated in each phantom, as in the simple phantom images of Figs. 3.4 and 3.5.) These results demonstrate the ability of the imaging method to accurately locate areas of dense tissue in all four phantoms. There is a deficit of resolution available from the microwave illumination relative to the smallest dimensions of the fibroglandular features within the breast. The resolution is limited by the bandwidth of the data used in the imaging algorithm, effectively resulting in a spatially-averaged reconstruction of the actual distribution. The images show that reconstructions of dense areas of tissue are smeared and that the smallest features of the tissue structures are not reconstructed. The resolution deficiency also contributes to the underestimation of the absolute dielectric properties observed in some areas of fibroglandular tissues. The abrupt transition between tissue types in the phantom profiles is reconstructed as a more gradual change in contrast, and this behavior tends to broaden and underestimate the reconstruction of the fine details of fibroglandular tissue structures. In addition, the cylindrical array geometry does not fully illuminate the anterior and posterior of the breast volume, 55 further smearing the features of the reconstructed profile along the array axis. Finally, the effective resolution may also be limited by over-regularization of the inverse solution. Class 1 1 2 3 4 DBIM Iteration Class 2 2 3 4 DBIM Iteration 5 • Residual Norm DBIM Iteration DBIM Iteration Figure 3.8 Convergence of the DBIM algorithm for the four heterogeneous phantoms. Left-hand y-axes: normalized convergence of the residual scattering norm. Right-hand y-axes: similarity metrics between the reconstruction at each iteration and the actual phantom profile. Here the cos(<j>e) metric is computed for the concatenation of all three Debye parameters, relative to the exact phantom profile. Figure 3.9 Class 1 phantom: Exact (a)-(c) and reconstructed (d)-(f) profile of Ae shown coronal (x — y, top row), sagittal (y — z, middle row), and axial (x — z, bottom row) cross-sections. Axes in cm. 57 Figure 3.10 Class 2 phantom: Exact (a)-(c) and reconstructed (d)-(f) profile of Ae shown in coronal (x — y, top row), sagittal (y — z, middle row), and axial (x — z, bottom row) cross-sections. Axes in cm. 58 35 30 25 20 15 10 5 0 2 4 0 35 30 0. 2 1 1 61 4 8 I 10 12 l 0 °l ol Ciass 3 phantom: Exact (a)--(c) anc. reconstructed (a)-(f) profile of Ac shown t igure a in coronal (x — y, top row), sagittal (y — z, middle row), and axial (x — z, bottom row) cross-sections. Axes in cm. 59 Figure 3.12 Class 4 phantom: Exact (a)-(c) and reconstructed (d)-(f) profile of Ae shown in coronal (x — y, top row), sagittal (y — z, middle row), and axial (x — z, bottom row) cross-sections. Axes in cm. 60 3.4.3 Other Performance Considerations The effect of the forward model error on the realistic phantoms is investigated by reconstructing the inverse crime case using simulated measurements of the down-sampled Class 2 phantom. The reconstructed Ae profiles from the 0.5 mm and 2.0 mm phan- toms are compared in Fig. 3.13. The similarity measure between the two reconstructions is cos(4>r) = 0.970. Artifacts in the Class 2 reconstruction due to modeling error are minimal in comparison to the artifacts in the case of the simple phantom (see Fig. 3.5). The severity of the artifacts in the simple phantom and the comparative absence of such artifacts in the Class 2 phantom suggest that our fixed-iteration CGLS termination heuristic may under-regularize the reconstruction in the case of the simple scattering target. The benefit of the bounding and constraint method described in Section III is evaluated by comparing the similarity metric for reconstructions with and without the constraints. The similarity measures in both cases are shown for each Debye parameter in Fig. 3.14. These results show that for each of the four realistic phantoms the quality of the reconstruction is improved when the bounding and constraint techniques are employed. Figure 3.14 clearly indicates that the static conductivity profile is the least successfully reconstructed of the Debye parameters in this multi-frequency solution. This observation helps explain the benefit of the bounding method, which restricts the range of values of a s based on the estimated profile of a better reconstructed parameter, e ^ . Since the bounds on both Ae and <js are tightened considerably at each iteration based on eoo, the projected-restarted constraint technique aids in the optimization of the inverse solutions within those bounds. The improvement in the a s profile is illustrated in Fig. 3.15 by comparing the constrained and unconstrained reconstructions of the Class 2 phantom. The fidelity of the scalar field approximation in 3-D imaging is investigated by imaging the Class 2 phantom using the vector field formulation of Eqn. (3.2). The contribution of the additional information available to the vector formulation is revealed by comparing the reconstructions of the vector and scalar formulations. The Ae profiles reconstructed in each case are shown in Fig. 3.16. There is little difference between the solutions obtained by the 61 two formulations. The peak values of the dense areas of the reconstruction in the vector case are only a few percent higher than those in the scalar case. The similarity measure between the two reconstructions is cos(4>r) = 0.989, confirming the visual agreement. (c) (f) Figure 3.13 Coronal, sagittal, and axial cross-sections of Class 2 phantom reconstructions: (a)-(c) scalar field approximation, (d)-(f) inverse crime with scalar field approximation. Cross-sections are taken from the reconstructed profiles of Ae, and are shown in coronal (x — y, top row), sagittal (y — z, middle row), and axial (x — z, bottom row) cross-sections. Axes in cm. 62 1 0.9 X® X 0.8 A A O & CO 8 X® 0.7 A A O O 0.6 0.5 1 2 3 Phantom Class ID 4 Figure 3.14 The similarity metric, cos(</>e), between the actual and reconstructed Debye parameter profiles of each of the four realistic phantoms. The metrics for the £00, As, and a s profiles are denoted by ' x ' , 'A', and 'o', respectively. The left-hand side markers for each phantom are the results for unbounded, unconstrained reconstructions. The right-hand side markers for each phantom, enclosed by ' • ' markers, are the results for bounded, constrained reconstructions. Here the cos(<f)e) metrics are the results for the individual vectorizations of each Debye parameter. (a) (b) (c) Figure 3.15 Effect of bounding and constraint techniques on reconstruction of static conductivity. Coronal cross-sections of static conductivity, a s , in (S/m) of the Class 2 reconstruction: (a) exact profile, (b) constrained reconstruction, (c) unconstrained reconstruction. 63 Figure 3.16 Coronal, sagittal, and axial cross-sections of Class 2 phantom reconstructions: (a)-(c) scalar field approximation, (d)-(f) vector field formulation. Cross-sections are taken from the reconstructed profiles of Ae, and are shown in coronal (x — y, top row), sagittal (y — z, middle row), and axial (x — z, bottom row) cross-sections. Axes in cm. 64 Strictly speaking, we have not implemented a full vector formulation, since the polarization of the dipole antennas restricts the system to scalar sourcing and observation of ^-directed fields. We also note that the scalar field approximation is referred to in this work in the context of the inverse solution. There is no need for the scalar approximation in the forward solution or simulated data acquisition, since 3-D FDTD is a full-wave numerical method. The scalar approximation is employed in this study due to insufficient memory resources for applying the method to the vector formulation of the large Class 1 phantom. We evaluate imaging performance in the presence of noise by adding increasing levels of white Gaussian noise to the measurement data aquired from the Class 2 phantom. The profile reconstructed at each noise level is quantitatively compared to the profile reconstructed for the noiseless case. The error between the noiseless and noisy reconstructed profiles is plotted in Fig. 3.17 by the metrics e2, cos(0 e ), and cos(0 r ) versus signal to noise ratio (SNR). Figure 3.17(a) depicts the classical log-log relationship between estimation error and additive white Gaussian noise. The slope of this line depends in part on the regularization and termination techniques used in the imaging algorithm. Representative cross-sections of the reconstructions over a range of SNR are shown in Fig. 3.18. Corruption of the images due to noise is apparent, as the accurately located areas of fibroglandular density in the noiseless reconstruction begin to shift and disappear at SNR levels below about 20 dB. The decay of the similarity metrics with decreasing SNR, shown in Fig. 3.17(b), confirm this threshold. We note that the reference signal level of the SNR is defined here as the mean of the measured total fields at each frequency over all multi-static measurements of the phantom. The forward solution is executed on a 128-core GPU using Acceleware's FDTD API, while the inverse solution is executed on a 4-core CPU using MATLAB code. The 40 antenna simulations of the forward solution are run in serial on the GPU, each taking about 30 seconds. The inverse solution requires about 5 minutes on a memory-limited system, for a total time per iteration of about 25 minutes. Parallelization of the forward simulations, implementation of the inverse solution on the GPU, and advances in GPU performance is expected to reduce the time per DBIM iteration by about a factor of ten [134]. 65 20 30 SNR (dB) 20 30 SNR (dB) (a) (b) Figure 3.17 Reconstruction error vs. signal to noise ratio, (a) e 2 error metric, referenced to noiseless reconstruction, (b) similarity metrics: cos(</>e) and cos(<f>r) are referenced to exact profile and noiseless reconstruction, respectively (both are computed for the concatenation of all three Debye parameters). 8 10 12 (a) 2 4 6 8 10 12 (c) 8 10 12 (b) 8 10 12 (d) Figure 3.18 Comparison of Coronal cross-sections of the Ae reconstruction of the Class 2 phantom for (a) noiseless data, (b) 20 dB SNR, (c) 10 dB SNR, (d) 0 dB SNR. 66 3.5 Conclusion We have presented 3-D imaging results for microwave tomography of realistic numerical breast phantoms of varying fibroglandular tissue density. The inverse scattering method employs a constrained solution that is shown to be robust across the full range of ACR density classifications. The successful reconstruction of the fibroglandular tissue distributions suggest the utility of the method to density characterization and related risk assessments. As the ACR classifications are based on the percentage of fibroglandular tissue in a 2-D mammographic projection, 3-D density assessments may offer improved or complementary information. Microwave imaging thus has the potential to play a role in an individualized risk assessment which includes an estimate of cancer risk based on breast density characterization. Our implementation of a 3-D microwave inverse scattering method also serves as a reference point for more computationally efficient techniques. The use of the scalar field approximation in the complex 3-D scattering environment is supported by the strong similarity of the imaging results to those of the full vector formulation. Solution of the inverse problem on a high resolution voxel basis provides a benchmark for methods that seek to exploit the lower resolution of the microwaves by using lower-dimensionality spatial bases to reduce computational cost. Our solution of the large-scale inverse problem is obtained by the use of an efficient regularization and constrained inversion technique. Applying microwave inverse scattering techniques to data simulated from realistic numerical models offers valuable insight into the performance and challenges that can be expected in experimental laboratory and clinical investigations. Using simulated data allows for the selective treatment or idealization of practical issues whose effect on imaging performance is otherwise difficult to attribute. In this way we can analyze a complicated system in an incremental manner. Conversely, there are practical challenges which cannot be faithfully represented by a computational model and the associated simplifying assumptions. Though 67 several practical issues are demanding of further study, the presented imaging results underscore the feasibility of microwave inverse scattering techniques in application to 3-D breast imaging. Future work in system development and tumor detection is motivated by the performance of the microwave imaging system in reconstructing the dielectric properties and spatial profile of realistic distributions of normal breast tissue. The expectation of low contrast between malignant and normal glandular tissue increases the importance of spatial resolution and estimation accuracy to enable direct visual diagnosis from the morphology of an invasive malignant growth. The resolution of the images presented in this paper suggests that direct visual diagnosis is challenging without further resolution improvements, for example by adding higher frequency data or by adopting edge-preserving regularization techniques. Alternatively, different approaches to microwave tumor detection, such as the use of differential imaging with exogenous contrast agents to enhance the relative contrast of malignant tissue [83,135], may be considered. 68 Chapter 4 Contrast-enhanced tumor detection using microwave imaging The detection of early-stage tumors in the breast by microwave imaging is challenged by both the moderate endogenous dielectric contrast between healthy and malignant glandular tissues and the spatial resolution available from illumination at microwave frequencies. The high endogenous dielectric contrast between adipose and fibroglandular tissue structures increases the difficulty of tumor detection due to the high dynamic range of the contrast function to be imaged and the low level of signal scattered from a tumor relative to the clutter scattered by normal tissue structures. Microwave inverse scattering techniques, used to estimate the complete spatial profile of the dielectric properties within the breast, have the potential to reconstruct both normal and cancerous tissue structures. However, the ill-posedness of the associated inverse problem often limits the frequency of microwave illumination to the UHF band within which early-stage cancers have sub-wavelength dimensions. In this computational study, we examine the reconstruction of small, compact tumors in three-dimensional numerical breast phantoms by a multiple-frequency inverse scattering solution. Computer models are also employed to investigate the use of exogenous contrast agents for enhancing tumor detection. Simulated array measurements are acquired before and after the introduction of the assumed contrast effects for two specific agents currently From: [83] J. D. Shea, P. Kosmas, B. D. Van Veen, and S. C. Hagness, "Contrast-Enhanced Microwave Imaging of Breast Tumors: a computational study using 3D realistic numerical phantoms", Inverse Problems, vol. 26, no. 7, Sp. Iss. SI, July 2010. 69 under consideration for breast imaging: microbubbles and carbon nanotubes. Differential images of the applied contrast demonstrate the potential of the approach for detecting the preferential uptake of contrast agents by malignant tissues. 4.1 Introduction Consideration of the safety, cost, availability, sensitivity, and specificity of established breast cancer screening methods such as mammography and magnetic resonance imaging (MRI) motivates interest in alternative or complementary technologies. The last decade has seen expansive investigation of a number of techniques in which the interior of the breast is probed noninvasively with low-power, microwave-frequency electromagnetic waves. The ultimate goal of this research is a safe, inexpensive, and accurate approach to the detection, monitoring, and assessment of breast cancers. Tissue-penetrating radar techniques have been applied to the tumor detection problem to locate strong scattering targets within the breast [14,15]. Electromagnetic inverse scattering has been widely investigated in application to breast imaging as such methods offer a complete mapping of the dielectric profile of the breast [21,74-78,80,84,117]. Inverse scattering approaches have been applied to data from numerical and experimental models of the breast, and some preliminary clinical studies have been performed [22,89]. The three-dimensional (3-D) imaging results presented in [21,84] were obtained with a frequency-domain inverse scattering algorithm formulated using the distorted Born iterative method (DBIM) [16]. The algorithm reconstructs the complete 3-D profile of the dielectric properties of the interior of the breast using multistatic array data measured over a discrete set of frequencies. The imaging method was evaluated in [84] in application to realistic numerical breast phantoms of healthy tissues over a range of fibroglandular density classifications. The method was shown to produce representative reconstructions of the actual profile of complex permittivity, suggesting that it may be suitable for normal tissue density assessment. However, the low resolution of the images suggests that early-stage malignancies may not be identifiable by direct inspection of the images. 70 Endogenous contrast and spatial resolution are two important considerations in the detection and location of tumors by microwave imaging techniques. Recent studies of the histopathology and microwave-frequency dielectric properties of excised breast tissues suggested a much lower contrast between healthy and cancerous tissues than was previously understood [1]. The ex-vivo microwave dielectric properties of malignant glandular tissues were observed to be about ten times those of adipose tissue, but only one tenth higher than the properties of normal glandular tissues. Furthermore, the electrical dimension of earlystage cancers are at or below the nominal half-wavelength resolution limit of the UHF-band frequencies (0.3 - 3.0 GHz) typically employed in frequency-domain microwave inverse scattering. Although super-resolution has been observed and attributed to evanescent waves in near-field measurements or in multiple-scattering environments [136,137], microwave detection of early-stage malignancies is nevertheless challenged by the moderate endogenous dielectric contrast, the small scattering area of these malignancies, and the heterogeneous scattering environment of healthy glandular tissue in which tumors often form. These fundamental challenges can be addressed in part through the use of comparative information, such as a contralateral comparison of a patient's left and right breasts [89], or the monitoring of tissue properties and distribution over time. A simple proof-of-concept of the latter approach was illustrated for a numerical breast phantom with and without a malignant inclusion [21]. Several medical imaging modalities have made use of contrast agents to enhance imaging in the areas of concentration of the agent. For example, breast malignancies are imaged in MRI with the introduction of a gadolinium contrast agent [138]. The preferential uptake of a contrast agent by malignant tissue [139,140] could enhance the contrast of the malignant tissue relative to the healthy surrounding tissue and allow microwave detection via comparative imaging using pre- and post-contrast measurements. In this computational study we investigate the use of exogenous contrast agents in microwave imaging to detect tumors based on the change in their dielectric properties after the introduction of the agents. We assume an effective change in dielectric contrast as suggested 71 by preliminary studies of two potential contrast agents for use in microwave breast imaging: microbubbles [67] and carbon nanotubes [141,142], A homogeneous spherical inclusion is placed among healthy fibroglandular tissues in anatomically realistic numerical breast phantoms of varying fibroglandular density. The inclusion is given dielectric properties representative of either malignant tissue or malignant tissue under an assumed influence of a contrast agent. A Gauss-Newton nonlinear inverse scattering method is used to reconstruct 3-D estimates of the numerical phantom profiles in each case. Differential imaging then compares the estimated profiles to assess the tumor information available to this approach. The challenges inherent to detecting small, compact breast tumors with microwaves are discussed in section 4.2 to motivate study of contrast agents. In section 4.3 we review the numerical breast phantoms which serve as the computational test beds for the imaging investigation and we summarize the reported effects of the contrast agents at microwave frequencies. Section 4.4 details the inverse scattering algorithm used to reconstruct the dielectric profiles of the numerical phantoms. We present imaging results in section 4.5 for three groups of numerical phantoms, each having a single tumor: homogeneous interiors, realistic heterogeneous interiors, and realistic heterogeneous interiors with the tumor under the influence of a contrast agent. Further discussion and concluding remarks are offered in section 4.6. 4.2 Microwave Imaging of Small Tumors Before turning our attention to the imaging of malignant inclusions in numerical breast phantoms, we first evaluate a fundamental dynamic affecting tumor detection performance in microwave imaging systems: the amount of scattering information available from the tumor depends primarily on the size of the tumor and the contrast of the tumor relative to the surrounding tissue. To illustrate these considerations, we compute the scattering crosssection of a lossy dielectric sphere as the diameter and contrast are varied. The scattering cross-section, as defined in [143], is the area of the incident plane-wave power density equal to the total power scattered from the sphere over all 47T steradians. The performance of inverse 72 scattering techniques generally improves with additional information obtained by adding angles of illumination and multi-static observations to the system. Thus, the scattering cross section of an object is a measure of the relative amount of information available to an imaging system over a range of frequencies. In particular, it can illustrate the limitation of information from sub-resolution features and low contrast objects. The scattering crosssection given by equation (11-104) of [143] is evaluated and plotted in figure 4.1 versus electrical diameter for a sphere of malignant tissue in three different lossless background materials exhibiting relative permittivity values that span the range of normal breast tissue. The three background media are referred to as fibroglandular, transitional, and adipose; these representative normal breast tissue types are explained in further detail in section 4.3. In each case, figure 4.1 shows the sharp decline in scattering information in the region below the A/2 resolution limit, where the scattering response from the tumor will be overtaken by errors, noise, and scattering from larger tissue structures. These curves also show that the low contrast of a tumor buried in normal fibroglandular tissue presents a scattering cross-section that decays below the physical cross-section above the nominal A/2 resolution limit. Inverse scattering is a model-based approach to microwave imaging that operates by comparing the computed scattering from an estimated object profile to the measured scattering from the actual object. In the case of an electrically small object, for which a resolutionlimited reconstruction will overestimate the dimension of the object, figure 4.1 suggests that reconstructed contrast will be underestimated to limit the magnitude of the computed scattering response equal to that of the actual object. According to this reasoning, we expect a resolution-limited microwave breast imaging system to underestimate the properties of electrically small features of the tissue distribution, such as early-stage tumors. We illustrate this effect in section 4.5 by imaging a simple numerical breast phantom with a single spherical tumor of varying diameter placed in the three cases of homogeneous background tissue (normal fibroglandular, transitional, and adipose) considered in figure 4.1. 73 Adipose — Transitional Fibroglandular (D 0 0.5 1 1.5 b Diameter (d/^ ) 2 Figure 4.1 Scattering cross-section of a homogeneous dielectric sphere of radius a as a function of diameter. The diameter of the sphere is given in units of wavelengths in the background medium, Xb. The sphere has the properties of malignant tissue and is placed in a lossless version of three representative backgrounds: adipose, transitional, and normal fibroglandular tissues. The contrast in relative permittivity between the sphere and the background at 2.5 GHz for these three cases are 11.6, 2.55, and 1.15, respectively. 4.3 Testbeds and Data Generation Measurements for our study are simulated using anatomically realistic numerical models from the University of Wisconsin Computational Electromagnetics Laboratory's (UWCEM) breast phantom repository [19,119]. In this section we summarize the models selected from the repository and the modifications made to these numerical phantoms to create the test beds employed in this paper. We present the endogenous frequency-dependent dielectric properties of the breast tissue models as well as the assumed effect of the microbubble and carbon nanotube contrast agents on the dielectric properties of malignant tissue. 74 4.3.1 Heterogeneous Numerical Breast Phantoms We select one phantom from each of the four American College of Radiology classifications of fibroglandular breast tissue density [20]: ID number 071904 from the Class 1 group ("mostly fatty"), 010204 from the Class 2 group ("scattered heterogeneity"), 062204 from the Class 3 group ("heterogeneously dense"), and 012304 from the Class 4 group ("extremely dense"). The phantoms are defined on a 0.5-mm uniform Cartesian grid. The 1.5-mm-thick, homogeneous skin region of the repository phantoms is replaced by a 2.0-mm-thick region to match the cell dimension of our imaging grids. The dielectric properties of the constituent tissues of the MRI-derived phantoms are mapped as described by Zastow et al. [19]. We use a least-squares fit of the single-pole Debye relaxation model to the reported properties of exvivo breast tissues [18] over our frequency range of interest, 0.5 to 3.5 GHz. The single-pole Debye model for complex permittivity is given by, OeffM , Ae CrelM + — = too + 7 — JUt0 1 + JUT , <Ts h- JU€0 (4.1) The relaxation time constant, r, of the Debye model is fixed at 15 ps for all tissues. The values of the Debye model parameters (infinite permittivity, e^; delta permittivity, Ae; static permittivity, es = e^ + Ae; and static conductivity, a s ) for the breast tissue models in the phantoms are given in table 4.1. The table gives the range of parameters assigned to adipose, transitional, and fibroglandular tissue regions in the phantoms. The relative permittivity, erei, and effective conductivity, creff, of these models are plotted over frequency in figure 4.2 (a)-(b). The Debye parameters assigned to the homogeneous skin region [120] are also given in table 4.1. The phantom and array are immersed in a lossless medium with a relative permittivity similar to vegetable oil, erei = 2.6. In this study we investigate the image reconstruction of a compact malignant mass. Since each phantom in the repository is comprised only of normal tissues, we add a homogeneous spherical inclusion having the median properties of malignant glandular breast tissue, as characterized by Lazebnik et al. [1] and given in table 4.1. A single 1-cm-diameter inclusion is placed adjacent to or within an area of fibroglandular tissue in each of the four test 75 phantoms. A survey of the size of the tumors detected in a large collection of breast cancer studies of various modalities shows that the median diameter of detected tumors ranges between 10 to 20 mm [144]. Thus, a 1-cm diameter for the inclusion represents a clinically relevant dimension at which to evaluate tumor imaging performance. In addition, 1 cm is 4r 70 3.5- E 60 w >> Fibroglandular > 50 1 o^ <d 40 Q. <D ,2.5(Fibroglandular o Q •U 30 nJ a> * 3 1.5 b Transitional 20 1 10 0.5- Transitional Adipose 1 05 Aaipose" 1.5 2 2.5 Frequency (GHz) 0 5 3.5 3.5 1.5 2 2.5 Frequency (GHz) (a) (b) 4 3.5 60 • Endogenous malignant Malignant w / microbubbles 3 | y Malignant w/ nanotubes 70 / / ' 50 :> 2.5 E a> a. o o> 1.5 1 Malignant w/ nanotubes Endogenous malignant 0.5 Malignant w/ microbubbles 0 5 1 1.5 2 2.5 Frequency (GHz) (c) 3.5 5.5 1 1.5 2 2.5 Frequency (GHz) 3 3.5 (d) Figure 4.2 (a) Relative permittivity, erei, and (b) effective conductivity, <7eff (S/m), of the ranges of adipose, transitional, and normal fibroglandular tissues used in phantom construction, (c) Relative permittivity, erei, and (d) effective conductivity, creff (S/m), of the malignant tissue models with and without the assumed effect of carbon nanotube or microbubble contrast agents. 76 Table 4.1 Debye parameters of the tissue models used in the construction of the numerical phantoms (valid from 0.5 to 3.5 GHz). The relaxation time constant parameter is r=15.0 ps for all tissues. Tissue ^oo Ae a s (S/m) Adipose 2.42-7.63 2.28-4.09 0.141-3.54 0.0023-0.0842 Transitional 7.63-36.7 4.09-16.8 3.54-19.9 0.0842-0.461 Fibroglandular 36.7-67.2 16.8-29.1 19.9-38.1 0.461-1.38 endogenous 56.6 18.8 37.8 0.803 with //-bubbles 39.7 13.2 26.5 0.562 with nanotubes 69.3 14.8 54.5 1.47 40.1 15.3 24.8 0.74 Malignant, Skin on the order of one half wavelength - the nominal resolution limit - in fibroglandular tissue at the upper frequencies of the UHF band. To emulate an idealized preferential uptake of a contrast agent by malignant tissue, we vary the properties of the inclusion according to the contrast effects suggested by preliminary studies of two contrast agents. A computational study of air-filled microbubbles in malignant tissue [145] showed that the effective permittivity and conductivity were decreased by approximately 30% due to the presence of microbubbles at a concentration of 20% by volume - a concentration that is within the dosage of existing ultrasound procedures. We reduce the Debye parameters of the malignant inclusion accordingly to simulate this concentration of microbubbles; thus the frequency-averaged contrast of the Debye model given in table 4.1 for malignant tissue with microbubbles is -30% in both relative permittivity and effective conductivity relative to the endogenous properties. In the case of carbon nanotube contrast agents, we fit Debye models to the reported dielectric properties measurements for an experimental material mimicking the dielectric properties of malignant tissue with and without a 2 mg/mL concentration of carbon nanotubes [142], The resulting frequency-averaged contrast of the Debye model given in table 4.1 for malignant tissue with carbon nanotubes is 77 +22% in relative permittivity and +66% in effective conductivity relative to the endogenous properties. The effect of each contrast agent on the dielectric properties of the malignant tissue model is plotted over frequency in figure 4.2 (c)-(d). 4.3.2 Homogeneous Numerical Breast Phantoms An additional set of numerical phantoms is created by downsampling the Class 2 phantom from the 0.5-mm grid to a 2.0-mm grid. The interior of the downsampled phantom is replaced by a homogeneous material described by one of three Debye models representative of nominal properties of each tissue type: normal fibroglandular tissue (600=17.5, Ae=31.6, crs=0.720), transitional tissue (eoo=10.5, Ae=11.7, ers=0.273), or adipose (€00=3.11, Ae=1.70, crs=0.0367). A single homogeneous spherical inclusion having the endogenous properties of malignant tissue given in table 4.1 is placed in the homogeneous interior. 4.3.3 Data Acquisition A 40-element cylindrical array of 14-mm dipoles surrounds each numerical test phantom, as shown in figure 4.3. The antennas are distributed on five elliptical rings of eight antennas each, with adjacent rings rotated by 7r/8. The rings are evenly spaced between the posterior and anterior coronal planes of each phantom. The dimensions of each array are set to provide a minimum spacing of 1 cm between each antenna element and the skin surface. Array measurements of the phantoms are simulated using the finite-difference timedomain (FDTD) numerical method. An auxiliary differential equation approach [121] is used to simulate Debye relaxation in the dispersive dielectric materials, and a perfectly matched layer [121] is used to terminate the grid. The dipole antennas are driven with a modulated Gaussian pulse having a bandwidth covering 0.5 to 3.5 GHz. Each antenna is sourced individually in independent simulations and the time-domain received fields are recorded at every antenna and converted to phasors at the frequencies of interest via the discrete Fourier transform. 78 y (cm) Figure 4.3 Diagram of the 3-D forward solution domain for the Class 2 numerical phantom, showing the dipole array, the downsampled skin layer, the measurement region T>, and the imaging region V of the interior breast volume. Data is acquired from each of the four numerical phantoms from the repository for three cases of the dielectric properties of the spherical inclusion: the endogenous malignant properties, decreased properties due to microbubbles, and increased properties due to carbon nanotubes. Data is also acquired from the downsampled Class 2 phantom for each of three cases of homogeneous interior background properties. The diameter of the inclusion in the breast interior is varied from 2 mm to 22 mm in 4 mm increments. Simulations of the downsampled phantom are conducted on the same 2-mm grid used by the reconstruction algorithm in order to remove modeling error from the evaluation. These simulated measurement sets are used to evaluate resolution and estimation error of the imaging system. White Gaussian noise is added to the simulated data such that the signal-to-noise ratio (SNR) is 50 dB. The reference signal level is defined at each frequency as the mean of the measured total fields over all multi-static measurements of the numerical phantom - typically about 60 dB below the source level for the inefficient and unmatched antenna models used in this study. This level of noise is comparable to the noise levels specified in prior experimental and numerical microwave breast imaging studies [8,78,80]. 79 4.4 Imaging Methodology In section 4.4.1, we briefly review the well-known theoretical relationship between the field scattered from a domain of unknown objects and the dielectric profile of that domain, and formulate an approximate discrete linear system of equations using multi-static array observations at multiple frequencies. The linear system is constrained by restricting the frequency dependence of the solution to the space of single-pole Debye solutions. An iterative method for solving the nonlinear scattering problem is summarized in section 4.4.2, and a simple approach to creating a differential image from two sets of measurements is presented in section 4.4.3. 4.4.1 Electromagnetic Inverse Scattering The scattering of electromagnetic radiation from a non-magnetic, heterogeneous dielectric object domain V at an observation location described by coordinate vector rQbs can be written as an volume integral equation, E*(r obs ) = k20 [ G 6 (r 0 bs|r')E t (r') [e(r') - eb(r')] dr' Jv (4.2) where the scattered field, E s , is the difference between the total field observed in the object environment and the total field observed in the background environment. The scattered field results from the re-radiation of the total field, E*, in V from the dielectric contrast formed by the difference between the object profile, e(r), and the background dielectric profile, e6(r). The free-space wavenumber is ho. The scattering contribution at r obs due to the dielectric contrast at r € V is determined by the Green's function, G b (r 0 b s |r), of the background dielectric profile. The object and background profiles are complex permittivities and both may be heterogeneous. The total field, E4, in V depends on the multiple scattering interactions between the features of the dielectric profile. The scattered field, E s , is therefore nonlinearly related to the 80 contrast function due to the product E 4 (r) [e(r) — eb(r)] in the integrand of (4.2). The relation can be linearized using the Born approximation [16], in which E* in V is approximated by the background field Eb, that is, the total field in the background profile e6(r). The volume V enclosing the unknown object volume can be discretized by the set of constant basis functions of edge length A given by 4>k( r) =< 1 0 /orlkfc-rlL < f (4.3) forotherwise for all rk = A [(lk + \) ax + (mk + ±) a^ + (nk + a 2 ] e V, (/, m, n) G A/"3. Using this basis for voxels k = 1 , . . . , K within V, the unknown contrast function can be written as a length K vector of basis coefficients, v, where vk — ^ fv [e(r) — e 6 (r)] 4>k(x)dv'. The Born approximation of (4.2) can then be written as a discrete vector equation. Noting the vector form of the scattered field, E s (r 0 b S ) = [E® Ey E^]T, the linearized vector scattering equation for a source at r s r c is E s = H v , where H is the 3 x K matrix given by, H = fc02A3[G6(robs|r1)Eb(r1|rsrc) ••• G & (r obs |r K ) E ft (r K |r src ) ] . (4.4) A set of field observations due to a source at r src , as well as additional sets of field observations due to other source locations, can be used to create a larger system of equations in the unknown v, where the vector scattering equations for each of M multi-static observations are combined into a length 3M vector, e s , and a 3M x K matrix, B. In practice, fields are typically measured by single-polarization antennas. We therefore consider only the scalar field observations, Es, made by an array of ^-directed dipoles. In this case e s is only length M. In either the scalar or vector observation case, the matrices associated with each observation are collected as the block-rows of a larger matrix, B = [ H f H j • • • H ^ ] T , operating on the same unknown contrast function vector, v, such that e,s' = Bv. Observations at multiple frequencies can be added to the linear system of scattering equations by vertically concatenating the multiple-observation vectors, e s (to), and matrices, B(u;), at each of F frequencies. The resulting MF x K system is less ill-posed than 81 a single-frequency system [146] and can be approximately inverted to yield a simultaneous solution for a frequency-independent contrast function v. When the actual dielectric profile is frequency-dependent, the system can be transformed by a parametric model of the frequency dependence [76]. For a general parametric model of complex permittivity, e(oj) = J7 (lu,Pi,P2, • • •), the multiple-observation linear scattering system at each frequency is transformed to, Pi - Pi e » = P2 - P2 (4.5) We choose the Debye relaxation model of (4.1) for the parametric model of the contrast function to fit the dispersive properties of biological tissues. Note that the coefficients dJ-/dp must be constant so that the system remains linear in the unknown parameters. Since d T j d r is a function of r in the Debye model, we therefore assume a fixed value for the time constant. The remaining three degrees of freedom in the model are sufficient to accurately capture the frequency dependence of the tissues over the frequency range of interest. Further, noting that the Debye parameters are real-valued we split each complex equation into real and imaginary parts so that the solution space is limited to real values. A 2 M F x 3 K system is obtained using (4.5) and the three-parameter Debye model (4.1). 82 e 6 ool oo K e ool e ooK Aei - Ae 6 ! 9{e-(WF)}J [^{[-g-BM £B(Wf) §-B(uf) ]} (4.6) We will denote the system in (4.6) as y = Ax. It is a multi-static, multiple-frequency, linearized description of the scattering due to the contrast function, formulated in the Debye parameter space. The solution of this system, described below, yields an estimate of the contrast function of the unknown object. 4.4.2 Distorted Born Iterative Method The linear system y = A x of (4.6) approximates the nonlinear scattering relation of (4.2). The error in the Born approximation used to linearize (4.2) increases with the magnitude of the contrast function. The high contrasts between the constituent breast tissues of table 4.1 will result in a highly inaccurate solution based on the Born approximation alone. However, a nonlinear method of optimization can be used to overcome the limitation. The distorted Born iterative method (DBIM) [16], a form of Gauss-Newton optimization [44], is the approach we use here to estimate the contrast function from the observations of nonlinear scattering. This method iteratively refines the estimate, beginning with an initial guess for the background profile that may include any available a priori information about the object profile. A solution to (4.6) (the inverse solution) is found at each iteration i to obtain an update to the estimated profile. The updated estimate, q(r), is used as the background profile, eb(r), 83 in the subsequent iteration and the background fields and Green's functions are recomputed (the forward solution). The difference in the scattering from the background profile and the object profile decreases at each iteration as the estimate is improved. The DBIM algorithm is terminated when the norm of the residual scattering, ||y||, converges sufficiently. The full DBIM procedure is outlined in Algorithm 1. We use the FDTD method in the forward solution to efficiently compute the background fields and Green's functions at multiple frequencies. The FDTD simulations are performed as described in section 4.3 for the acquisition of measurements with the exception that the grid dimension in the forward model is 2.0 mm. The cells of the forward solution grid are the same as the voxels defined in (4.3), extended V k : £ D. where T> is a cuboid measurement region enclosing V and the array of N antennas. The forward solution domain is depicted by the 3-D diagram in figure 4.3. In addition, the background field due to each source, Eb(rfc|rsrc), is measured at each cell r^ £ V. The heterogeneous background Green's function is computed using the principle of reciprocity [125] and the field measured at e V due to a source current Ix, Iy, or Iz of length A at rsrc G V. G b (r obs |r fc ) = cu/^A i Li 4 Tz Eby(rk) |,si/>=0 o i i E z( r *)\i t ,i z = 0 (4-7) E»z(rfc)|Wy=0 The antenna array described in section 4.3 consists of ^-oriented dipoles sourced with I z . Thus, only the last row of the Green's function tensor of (4.7) is computed. The inverse solution is the inexact Newton step found using the conjugate gradient for least-squares (CGLS) inversion method [95], which finds a solution to the system of normal equations, A T y = A T Ax. The inverse solution approximates the contrast updates for each voxel rfc e V. The update step length is fixed at unity. The under-determined system in (4.6) is ill-posed and requires regularization. Regularization is accomplished by early termination of the CGLS algorithm after a fixed number of iterations [78,84], The termination condition is based on empirical observations of the trade-off between reduction of the norm of the 84 residual scattering, ||y||, and the divergent growth of the norm of the unknown contrast function, ||x||, that is indicative of the corruption of the solution by noise and other errors. 4.4.3 Differential Imaging A small change in contrast or a change to a localized area of the object profile may not be recognizable among the more sizable features of the reconstruction. The effect of a change in the object profile on the reconstructed profile can be visualized by the difference between the reconstructed profiles before and after the change. In this way, not only is the location of the change revealed but the ability of the system to measure information from that change is confirmed. Thus, by imaging the difference between the reconstructions of the phantoms before and after the introduction of a contrast agent, we evaluate both the scattering information available from a change in contrast of a compact malignant inclusion as well as the ability of the method to locate that change. First, we simulate measurements of a phantom containing a malignant lesion with endogenous dielectric properties. Next, we assume a change in dielectric properties of the inclusion due to the contrast agents as detailed in table 4.1 and simulate a second set of measurements. The "before" and "after" measurements are then individually reconstructed by the DBIM, and the difference between the two reconstructions is imaged. 4.5 Results We begin with the reconstruction of the downsampled homogeneous numerical phantom modified described in section 4.3.2. We then present images of the reconstruction of the four realistic heterogeneous numerical phantoms, each having a single inclusion of malignant tissue with endogenous properties. Finally, we present the differential images of all four realistic phantoms, using the reconstructions of each with and without the effect of the contrast agent applied to the properties of the inclusion. In all cases presented in this paper, the linear system is constructed with data acquired at 1.0, 1.5, 2.0, and 2.5 GHz from all multistatic array channels, excluding the monostatic Algorithm 1 Gauss-Newton for Distorted Born Iterative Method i 0 €q initial background profile repeat for m = 1 to N do FDTD: E 6 (r n |r m ) for n = 1,..., N ; E 6 (r fc |r m ) and G b (r fc |r m ) for k = 1,..., K E s (r„|r m ) E m e a s (r n |r m ) - E b (r n |r m ) for n = 1,..., N G 6 (r m |r fc ) 4r- G 6 (r fc |r m ) T for k = end for construct Aj, yj x, CGLS(A l ,y l ) e\ +1 ^ - e J + ^ X i ) i <r~ i + 1 until ||yi_i|| - ||y l || < 5 ||y 0 || \,...,K 86 channels. For each phantom, the initial background profile of the breast interior is a homogeneous interior having the average properties of the true phantom profile and includes a downsampled version of the skin region of the true phantom. The unknown imaging region V is defined to be the interior breast volume inside the skin layer. The dielectric properties of the immersion and skin regions are assumed known a priori. The DBIM is terminated after the 8th step and the CGLS algorithm is terminated at the 5th iteration at every step of the DBIM. These conditions ensure convergence at termination in all the cases shown such that the decrease in the residual scattering norm, ||i/||, is less than 1% of the residual scattering norm of the initial step. The CGLS termination condition is based on the results of an L-curve heuristic used for imaging these numerical phantoms in [84], All imaging results are visualized in orthogonal cross-sections of the profiles of relative permittivity and effective conductivity at 2.5 GHz, computed according to (4.1) from the reconstructed profiles of the Debye parameters. On the 2-mm voxel grid, the forward domains of the four realistic numerical phantoms contain between roughly 500,000 and 1,000,000 cells while the inverse domain contains between roughly 120,000 and 360,000 unknowns. The forward solution is executed on a 128-core GPU using Acceleware's FDTD API, while the inverse solution is executed on a 4-core CPU using MATLAB code. The 40 independent simulations of the forward solution are run in serial on the GPU, each taking 20 to 30 seconds. The inverse solution requires about 3 minutes on a memory-limited system, for a total time per step of the DBIM of around 20 minutes. Parallelization of the forward simulations, implementation of the inverse solution on the GPU, and advances in GPU performance are expected to reduce the time per DBIM iteration by about a factor of ten [134]. 4.5.1 Imaging a Tumor in a Homogeneous Phantom Interior The simplified breast model consisting of a skin region, a homogeneous background in the interior breast volume, and a single spherical inclusion is imaged over all combinations of three background tissue types (adipose, transitional, and fibroglandular tissue properties) 87 and six diameters (2, 6, 10, 14, 18, and 22 mm) of a homogeneous inclusion with malignant tissue properties. The effect of inclusion diameter on the estimation accuracy is shown in figure 4.4 by comparison of the coronal cross-sections of the reconstructed relative permittivity at 2.5 GHz of inclusion diameters 6, 14, and 22 mm in a transitional tissue background. The error in the estimate is normalized by comparing the estimated Debye curves to the actual curve relative to the initial guess curve, (e* — e) / (e* — e°), averaged over the frequencies used in the reconstruction. The real and imaginary parts of the error are computed separately and are plotted in figure 4.5 versus the diameter of the inclusion. The peak of the estimated contrast for each case of inclusion diameter and background tissue is compared to the actual contrast in table 4.2. These results show the improvement in the estimation accuracy of the imaging system as the inclusion diameter is increased through a range of sub-wavelength dimensions. The 0 2 4 6 8 (b) 10 12 14 0 2 ^ i 8 (d) 10 12 14 0 2 4 6 8 1*0 12 14 (f) Figure 4.4 Comparison of the exact (top row) and reconstructed (bottom row) coronal cross-sections (in color) of relative permittivity for a single spherical inclusion of malignant tissue properties within a homogeneous breast volume of transitional tissue properties. The inclusion diameters are (a)-(b) 6 mm, (c)-(d) 14 mm, and (e)-(f) 22 mm. Spatial smoothing of the inclusion boundary and increasing underestimation of the inclusion properties with decreasing diameter are evident in the reconstructions. (Axes in cm.) 88 \ ' 0.6 -o * . l/ue * * ^ efl - _est> , / Jrue Jnit, eft ^ ^ ell ' ^ 0.5 1 1 0.6 l U) L U V ; 0.4 •§0.4 N i : 0.2 Diameter ( d I X b ) (a) 1 o £ 0.8 c0 i X. "-A . 1 i i 0.8 0. .true •< a ,ll 0.2 i 0.2 o z _estv , - true Jnit> elf' ' eff - " e l l I 0.4 0.6 Diameter (d/X 6 ) (b) 0.S 0, • (<"" -E~V (E^-C'l ,_true •C.n _ 0.1 est, , .true inil efl ^ ^ ell " " e l l 0.2 0.3 Diameter (dA.b) (c) Figure 4.5 The normalized error in the estimation of the dielectric properties of the spherical malignant inclusion versus diameter in a homogeneous breast interior of (a) fibroglandular tissue, (b) transitional tissue, and (c) adipose tissue. The estimate is defined as the average of the reconstructed properties over the actual inclusion volume. The error is calculated as (u* — ii est ) / (u* — ii mit ), where u is either relative permittivity, erei, or effective conductivity, <reff, at 2.5 GHz. diameters are given in wavelengths in figure 4.5 and table 4.2 to demonstrate the dependence of estimation accuracy on the electrical dimension of the scatterer in a given background tissue. For example, despite the higher contrast of the inclusion in an adipose background, the longer wavelength (compared to that of other background tissues) results in increased estimation error. Thus, the sensitivity of tumor detection in the microwave imaging system will depend not only on the tumor contrast but also the wavelength in the tissue surrounding the tumor. 4.5.2 Imaging a Tumor in a Heterogeneous Phantom Interior In the Class 1 and 2 phantoms the inclusion is placed adjacent to areas of both adipose and fibroglandular tissue, while in the Class 3 and 4 phantoms the inclusion is mostly surrounded by fibroglandular tissues. Table 4.2 and figure 4.5 suggest that a 10 mm inclusion will not be well reconstructed in a background of primarily adipose tissue, and that while the same inclusion in a background of dense fibroglandular tissue may be more accurately reconstructed, the low contrast and the low imaging resolution may blur it into adjacent 89 Table 4.2 Peak contrast of estimated complex permittivity at 2.5 GHz of a malignant spherical inclusion in a homogeneous breast interior. For reference, the electrical size of the inclusion diameter in each background tissue is given in wavelengths. Exact Background Contrast Fibroglandular, Peak Estimated Contrast (inclusion diameter) (2 mm) (6 mm) (10 mm) (14 mm) (18 mm) (22 mm) (0.11A) (0.34A) (0.57A) (0.80A) (1.03A) (1.26A) max (e re i/4 el ) 1.15 1.00 1.05 1.10 1.17 1.18 1.19 max (aefr/cr'!ff) 1.13 1.00 1.03 1.07 1.12 1.14 1.15 (0.08A) (0.23A) (0.39A) (0.54A) (0.70A) (0.85A) Transitional, max (e r e i/4,) 2.55 1.02 1.27 1.77 2.57 2.98 2.88 max (<7eff /(Jgff) 2.99 1.02 1.28 1.78 2.47 2.73 2.94 (0.04A) (0.11A) (0.18A) (0.25A) (0.33A) (0.40A) Adipose, max (e re i/e 6 el ) 11.6 1.01 1.12 1.39 2.16 5.25 9.16 max (creff/ <Tgff) 22.1 1.01 1.17 1.57 2.71 7.15 16.3 features of the heterogeneous background. These predictions are confirmed in figures 4.64.9, where the reconstructions of the four realistic phantoms introduced in section 4.3 are compared to the actual phantom profiles. While the areas of dense fibroglandular heterogeneity are reconstructed, there is no clear visual evidence of the inclusion in any of the four reconstructions. These imaging results illustrate the expected challenges in directly imaging electrically small tumors, and motivate the use of contrast agents in microwave imaging. 90 Figure 4.6 Reconstruction of the Class 1 phantom with a single 1.0-cm-diameter spherical malignant inclusion. Coronal (top row) and sagittal (bottom row) cross-sections (in color) of the exact and reconstructed volumes are shown at 2.5 GHz for (a)-(d) relative permittivity and (e)-(h) effective conductivity in (S/m). (Axes in cm.) (b) (d) (f) (h) Figure 4.7 Reconstruction of the Class 2 phantom with a single 1.0-cm-diameter spherical malignant inclusion. Coronal (top row) and sagittal (bottom row) cross-sections (in color) of the exact and reconstructed volumes are shown at 2.5 GHz for (a)-(d) relative permittivity and (e)-(h) effective conductivity in (S/m). (Axes in cm.) 91 Figure 4.8 Reconstruction of the Class 3 phantom with a single 1.0-cm-diameter spherical malignant inclusion. Coronal (top row) and sagittal (bottom row) cross-sections (in color) of the exact and reconstructed volumes are shown at 2.5 GHz for (a)-(d) relative permittivity and (e)-(h) effective conductivity in (S/m). (Axes in cm.) Figure 4.9 Reconstruction of the Class 4 phantom with a single 1.0-cm-diameter spherical malignant inclusion. Coronal (top row) and sagittal (bottom row) cross-sections (in color) of the exact and reconstructed volumes are shown at 2.5 GHz for (a)-(d) relative permittivity and (e)-(h) effective conductivity in (S/m). (Axes in cm.) 92 4.5.3 Differential Imaging of a Tumor in a Heterogeneous Phantom Interior Differential images showing the effect of each contrast agent on the reconstruction of relative permittivity and effective conductivity in each of the four realistic phantoms are shown in figures 4.10 - 4.13. These images are produced by subtracting the profiles reconstructed with and without the effect of the contrast agents. The cross-sections are taken through the voxel having the peak estimated contrast and the values are clipped at zero since the direction of the change in contrast caused by each agent is known. The white ' x ' marker on each cross-section marks the true location of the center of the inclusion, as projected onto the cross-sectional planes of the estimated location. The estimated change in contrast of the inclusion and the error in the estimated location of the inclusion are given in table 4.3 for each agent in each phantom. Table 4.3 Changes in dielectric properties and location errors at 2.5 GHz in the differential imaging of a 1.0-cm-diameter malignant spherical inclusion in realistic heterogeneous phantoms using microbubble and carbon nanotube contrast agents. Phantom Exact Change (%) Est. Change (%) ||f-r*| | (mm) Class Agent ^rel OefT Crel Ceff Crel Ceff 1 //-bubble -30 -30 -1.8 -1.2 4.5 2.8 nanotube +22 +60 +0.9 +0.8 6.0 3.5 //-bubble -30 -30 -6.6 -5.6 8.2 7.2 nanotube +22 +60 +4.1 +4.8 12 11 //-bubble -30 -30 -9.1 -8.5 4.5 4.5 nanotube +22 +60 +4.6 +6.8 8.2 4.5 //-bubble -30 -30 -6.6 -4.5 8.5 7.5 nanotube +22 +60 +4.6 +3.9 8.9 6.3 2 3 4 93 -0.2 |-0.0 05 -0.4 -0.01 -0.6 -0.8 0 2 4 6 8 10 12 14 16 18 0 2 4 6 8 10 12 14 16 : -0.015 10 1-0.02 128 10 12 14 16 18 0 2 4 6 8 10 12 14 16 18 y. Figure 4.10 Differential images of the Class 1 phantom. Coronal (top row) and sagittal (bottom row) cross-sections of the complex permittivity at 2.5 GHz are taken through the peaks of the changes in (a)-(b) relative permittivity and (c)-(d) effective conductivity (S/m) due to microbubble contrast agent, and (e)-(f) relative permittivity and (g)-(h) effective conductivity (S/m) due to carbon nanotube contrast agent. The ' x ' markers show the actual location of the tumor as projected onto each cross-sectional plane. (Axes in cm.) (b) (d) (f) (h) Figure 4.11 Differential images of the Class 2 phantom. Coronal (top row) and sagittal (bottom row) cross-sections of the complex permittivity at 2.5 GHz are taken through the peaks of the changes in (a)-(b) relative permittivity and (c)-(d) effective conductivity (S/m) due to microbubble contrast agent, and (e)-(f) relative permittivity and (g)-(h) effective conductivity (S/m) due to carbon nanotube contrast agent. The ' x ' markers show the actual location of the tumor as projected onto each cross-sectional plane. (Axes in cm.) 94 82.5 •2 O r " 21 41 Or 2- 10.12 I0.1 4! I008 ll.5 ll •0.5 14L 0 2 4 6 B 10 0 12 2 4 6 8 10 12 80.06 10| |0.04 12! 10 02 lo 0 2 4 6 (e) (d) 10 12 10 12 lo (g) 0 (b) 8 2 4 6 6 (h) (0 Figure 4.12 Differential images of the Class 3 phantom. Coronal (top row) and sagittal (bottom row) cross-sections of the complex permittivity at 2.5 GHz are taken through the peaks of the changes in (a)-(b) relative permittivity and (c)-(d) effective conductivity (S/m) due to microbubble contrast agent, and (e)-(f) relative permittivity and (g)-(h) effective conductivity (S/m) due to carbon nanotube contrast agent. The ' x ' markers show the actual location of the tumor as projected onto each cross-sectional plane. (Axes in cm.) 0 2 4 6 8 10 12 14 16 I 1 I 0.06 -0.02 0.04 -0.04 -0.06 0 2 4 6 8 10 12 14 16 -0.08 0.02 0 2 4 6 0 10 12 14 16 2 4 6 8 (a) 10 12 14 16 (g) 1° 10.06 Jo.04 I:! 0 2 4 6 10.02 10 12 14 16 10 12 14 16 (b) (d) (f) 0 2 4 6 8 10 12 14 16 'o (h) Figure 4.13 Differential images of the Class 4 phantom. Coronal (top row) and sagittal (bottom row) cross-sections of the complex permittivity at 2.5 GHz are taken through the peaks of the changes in (a)-(b) relative permittivity and (c)-(d) effective conductivity (S/m) due to microbubble contrast agent, and (e)-(f) relative permittivity and (g)-(h) effective conductivity (S/m) due to carbon nanotube contrast agent. The ' x ' markers show the actual location of the tumor as projected onto each cross-sectional plane. (Axes in cm.) 95 In all cases, the differential images clearly show a compact change in the contrast of both relative permittivity and effective conductivity. The change is well resolved in all three dimensions, with some smearing along the coronal axis caused by the incomplete illumination of the breast by the cylindrical antenna array geometry. The error of the estimated location of the peak change in contrast is on the order of 10 mm or less and there is no discernible trend in the location error relative to the fibroglandular density of the phantom. The change in the dielectric properties of the inclusion is underestimated in all cases, consistent with the prediction of both the scattering analysis of section 4.2 and the reconstructions of the homogeneous phantoms with 1.0 cm inclusions. The estimated change in properties are about 10-20% of the actual change in Class 2, 3, and 4 phantoms, while the estimate was only a few percent of the actual change in the mostly-fatty Class 1 phantom. As noted in section 4.3, the data used in the reconstructions of the heterogeneous numerical phantoms include additive white Gaussian noise to create a 50 dB SNR. The effect of the noise on the difference images is evaluated over a range of noise levels. The perturbation of a difference image due to the noise is evaluated by subtracting the noiseless difference image from the noise-corrupted difference image. The peak noise perturbation for each numerical phantom is plotted in figure 4.14 for SNR from 10 to 50 dB. Comparison of these curves to the peak of the tumor spot in the difference images of figures 4.10a, 4.11a, 4.12a, and 4.13a, suggests accurate tumor identification in all cases studied for SNR greater than approximately 40 dB. We note that the image noise is higher in the reconstructions of the dense heterogeneous Class 3 and Class 4 tissue distributions since they reach higher dielectric properties over a more extensive fibroglandular distribution. 96 t« \\ \ \ -•©•-Class 1 • - * • - Class 2 - - • - - Class 3 - Class 4 • q* 75 10 yo o CO CD • » % Q. ' • V V, CD (/] °10 20 •w 30 SNR (dB) % 40 50 Figure 4.14 Peak pertubation of relative permittivity in microbubble difference images reconstructed from noise-corrupted measurement data from each of the four heterogeneous numerical breast phantoms. The plots indicate the maximum absolute perturbation (over all voxels in V) due to a given noise level, relative to the noiseless difference image. 4.6 Discussion and Conclusion The issue of the limited resolution of the imaging system has been carefully considered in this investigation. It is important to note, however, the role of the regularization technique in determining the detection sensitivity of the system. That the small tumors are not visible in the images of figures 4.6 - 4.9 is a consequence of both the loss of scattering information predicted by figure 4.1 as the tumor dimension decreases and the low-pass spatial filtering effect of the regularization technique employed in our algorithm. While there exist various edge-preserving approaches to regularization and global optimization techniques that can sharpen the resolution of the imaging system, there remains a fundamental issue of information deficiency that is independent of such techniques. Namely, as the dimension of a tumor decreases into the regime below one half-wavelength, the magnitude of the scattering decays exponentially until it decreases below the measurement sensitivity of the system. The problem of low level scattering information can exist even at larger dimensions due to the low dielectric contrast between malignant tissue and the surrounding healthy glandular tissue. 97 In this work we have addressed the deficiency through the enhancement of tumor contrast by exogenous contrast agents, and the use of a differential imaging scheme to capture their effect. There are practical issues to be considered in a real-world implementation of a differential imaging scheme. The assumptions of the contrast available from the agents and the selective uptake of the contrast agent in cancerous tissue are supported in the literature, but require further confirmation in a clinical setting. A differential image will reveal not only the change in the reconstruction due to a contrast agent but will also include artifacts due to any changes in the test fixture, location and position of the breast, noise, and any other measurement uncertainties. Our numerical models idealize many of these considerations, but provide an initial test bed in which to study the potential performance of a differential imaging scheme. Note that the assumed endogenous contrast between malignant and normal glandular tissue used in this initial test bed is based on ex vivo measurements and may differ slightly from the actual in vivo endogenous contrast. Our analysis and results have illustrated the dynamics of tumor detection with respect to contrast and resolution in microwave imaging, particularly in the lower microwave frequencies often employed for inverse scattering solutions. The differential images demonstrate a successful approach by which contrast agents can be used in conjunction with microwave imaging to detect tumors, including compact tumors of dimensions below the nominal resolution of the imaging system. The results of our investigation suggest a promising direction for further research on contrast-enhanced tumor detection by microwave imaging. Microwave breast imaging has been an active research area over the past two decades and has received considerable recent attention. As supported by the results and discussion in this paper, there are a number of fundamental and contemporary issues deserving of further consideration in ongoing research of the microwave modality. Resolution remains an important issue and its fundamental dependence on wavelength, contrast, and measurement sensitivity points to potential avenues towards improving microwave detection performance by improving the information content of the measurement data. In addition, there remains 98 significant opportunity for improvement of the performance through selection and development of imaging methods, including techniques of regularization, optimization, solution constraints, and modeling. Image validation is another important requirement in the development of microwave imaging. Computational studies are easily validated since the exact object is directly available for comparison; experimental studies offer practical measurement challenges with less realistic phantoms; clinical studies offer the best possible test case as well as perhaps the most challenging case for validation due to uncertainty of the true dielectric properties distribution. Lastly, there is a need for increased attention to the qualification of the potential clinical value of the microwave modality for medical imaging, including applications in screening, risk assessment, and tumor detection. 99 Chapter 5 Analysis of Microwave Scattering in Breast Imaging A wide variety of methods have been applied to the electromagnetic inverse scattering problem for medical imaging of the breast at microwave frequencies. While many techniques have been leveraged towards the development of a microwave imaging solution, they are all fundamentally dependent on the scattering data as a description from which to image the breast. Evaluating and optimizing the information contained in the scattering data is therefore instrumental in understanding and achieving optimal performance from any particular method of imaging. In this paper a method of analysis is employed for the evaluation of the information contained in the electromagnetic scattering from an unknown dielectric profile. The method obtains an estimation of imaging performance by mapping the data through the inverse of the scattering system onto the object space. The inverse is evaluated using a factorization of the linear system of equations resulting from the discrete approximation of a set of exact scattering relations. The truncated singular value decomposition (TSVD) is used to investigate the effective rank of the scattering system and to quantify and illustrate an upper bound on imaging performance when the solution space is limited to the span of the scattering matrix. The analysis is applied to three-dimensional, anatomically realistic numerical breast phantoms. The utility of the method is shown for various considerations in system design and problem formulation, providing design and development insights and best-case expectations of imaging performance. In particular, the analysis offers a method of decoupling the problem of selecting a set of data from the problem of forming an image from that data. 100 5.1 Introduction Research into the use of microwaves for medical imaging has primarily focused on breast imaging, as the location exterior to the high-loss thoracic wall and the low-loss adipose tissue content of the breast make the measurement of scattered waves feasible [61]. In addition, the dielectric contrast at microwave frequencies between the constituent tissue types of the breast provides an electromagnetic mechanism for sensing, locating, and identifying the tissue structures [18]. However, substantial technical challenges remain due in part to the propagation loss of the waves through the breast, the relatively low resolution of microwaves, and a recently revised understanding of only a modest contrast between normal and malignant glandular tissues [1]. Inverse scattering techniques provide an estimation of the complete dielectric property profile over the imaging region [16] and can therefore account for the scattering responses of both malignant and healthy fibroglandular tissue structures which have similarly high contrasts against a background of adipose tissue. However, inverse scattering systems for breast imaging are usually designed for frequencies in the UHF band (0.3 to 3 GHz), since the nonlinear inverse problem becomes increasingly difficult to solve at the higher frequencies of the UWB medical imaging band (3.1 to 10.6 GHz). The longer wavelengths of the UHF band scatter less efficiently from the millimeter-scale fine details of fibro-connective, ductal, and glandular tissue structures in the breast. The resulting loss of imaging resolution in UHFband inverse scattering has led to the use of techniques and algorithms capable of overcoming the nonlinearity of the scattering while improving image resolution and preserving the edges between different tissue types (see, e.g., [59,82,147-149]). While there is a considerable amount of contemporary research into imaging techniques that can be used to reconstruct the breast profile from the scattering data (see, e.g., [58, 60,74,81,84,85,103,123]), there is comparatively little evaluation of the quality of the scattering data. The information available from the scattering data is typically considered only indirectly in the evaluation of a final image, obscured by other uncertainties and errors. 101 However, there are two relatively distinct technical components to the microwave imaging problem. The first is the design of the measurements and the formulation of the inverse problem, and the second is the selection and development of an imaging method to efficiently and accurately image the breast from the data. It is natural, then, to decouple these two components to ensure that the success of an imaging algorithm is not hindered by poor data set selection or an error-prone problem formulation. Many prior contributions have considered the fundamental limitations of scattering data in the reconstruction of an unknown object or source profile. Slaney et al. showed the limitations of the Born and Rytov field approximations in microwave tomography [33]. Bucci and Franceschetti studied the information content of scattered electromagnetic fields, showing the degrees of freedom and band-limitation of scattered fields with finite measurement precision [150, 151] and investigating the maximum number of recoverable unknowns for radiative-region scattered field observations [152]. Hori presented an analysis method for geophysical inverse problems in which the Green's function is spectrally decomposed and the inverse problem is evaluated by a truncated set of eigenfunctions [153]. The use of singular-value decompositions (SVD) in the evaluation of ill-posed inverse problems is described in detail by Hansen [96]. Singular value analysis has been used previously in the study of microwave imaging. Fang et al. used the SVD to evaluate the effect of various design parameters on the invertibility of the microwave imaging problem [146]. The study was limited to a two-dimensional background region for which the singular spectrum of the scattering matrix was used to compare the degree-of-illposedness of the inverse problem as a parameter of interest was varied. Winters et al. used the SVD to create a custom spatial basis for the solution that was pre-regularized to remove basis functions that would not be resolved by the illumination [21]. In this paper we apply a truncated SVD (TSVD) method of analysis intended to decouple the obfuscating issues arising in the solution of nonlinear inverse problems from the evaluation of the quality of the scattering data with respect to system design and problem 102 formulation. The nonlinearity of the problem is avoided by making use of a priori knowledge of the exact fields over the object region. The problem is then linearized without field approximations. The exact fields are readily available from electromagnetic simulations when studying test objects in a purely computational environment. Our MRI-derived, three-dimensional numerical breast phantoms [19,119] therefore offer ideal test cases for the analysis. These high-fidelity numerical phantoms are accurate in dielectric properties and realistic in tissue distribution, allowing for a faithful simulation of loss, propagation, and scattering that is unmatched by 2-D test domains or studies of simple 3-D objects. We compute the vector fields and tensor Green's functions for the breast phantoms using finite-difference time domain (FDTD) simulations. The imaging region is illuminated by electrically short current sources for each orthogonal linear polarization. A system of scattering equations is constructed for inversion, but the idiosyncrasies of any particular choice of a wide variety of applicable imaging methods are deferred by employing the controlled regularization of the TSVD. The TSVD analysis provides a systematic investigation of the information contained in the scattering data and the performance relative to various system design considerations. Specifically, we investigate single-frequency formulations, multiple-frequency formulations, channel selection, noise performance, and field approximations. The analysis method provides a sense of an upper-bound on imaging performance to be used as a guide in design of the measurement system and formulation of the inverse problem, and as a benchmark for evaluating practical imaging performance. In the next section we provide a brief background for electromagnetic inverse scattering and TSVD analysis. In Section 5.3 the details of the analysis are presented in the context of the microwave breast imaging problem. The test phantoms are presented in Section 5.4, and in Section 5.5 the analysis is applied to the phantoms. A discussion of the method and the results follows in Section 5.6 and the investigation is concluded in Section 5.7. 103 5.2 Background We focus our attention in this investigation to breast imaging using methods of electromagnetic inverse scattering. In general, inverse scattering techniques reconstruct an estimate of the dielectric properties of an unknown object by inverting a system of equations describing electromagnetic scattering in the presence of an unknown object. In this section we briefly review the object-scattering relationship and the general form of the inverse problem. Though we define the problem in the frequency domain for three-dimensional space, the analysis in the time-domain or lower-dimensional space is analogous. Let an imaging target in a known background of complex wavenumber kb be described by an unknown spatial distribution, over a finite volume V, of the complex wavenumber k(r) where r is a position vector. Let Eb be the incident field over V in a known background, let E 1 be the total field in V in the presence of the object to be imaged, and define the scattered field as Es = El — Eb. The volume equivalency principle and Helmholtz wave equation in terms of the electric field may be used to formulate an equation for the field perturbation due to the unknown object in terms of the contrast between the unknown and background wavenumber profiles [94]. V2Es(r) + k2(r)Es(r) = - (k 2 {r) - k2b(r)) E\r) Given a tensor Green's function solution for the background profile, (5.1) a frequency domain integral equation is available to express the electromagnetic scattering from an object, due to a source at r s r c and observed at rQbs, in terms of the unknown wavenumber contrast. Es(rohs, rsrc) = J Gb(rohs, r')E\r\ 2 2 r s r c ) ( k ( r ' ) - k {r')) dr' (5.2) v The scattering relationship in (5.2) is a Fredholm integral equation of the first-kind [96], a class of equations often appearing in inverse problems. In the case of electromagnetic inverse scattering the equation is further characterized by its nonlinearity due to the dependency of the total field, E1, on the unknown object profile, k(r). 104 The scattering data, Es, are the difference between measurements of the total field and measurements or computations of the background field. Data can be collected over a set of spatial channels (selections of source and observation point pairs) and over a set of frequency channels (selections of source frequencies), resulting in a system of equations in the unknown wavenumber. Since closed-form, continuous, inverse operators are unavailable for the general scattering problem, the scattering equations are typically discretized over a finite spatial basis, yielding a system of linear equations. The system is composed of an M x N scattering matrix, A, an N x 1 vector, x, of the unknown contrast function, and an M x 1 a vector of scattered field data, b. In the basic formulation, M is the number of channel measurements and N is the number of spatial basis functions in the discretization. While the selection of field data, the choice of unknown parameters, and use of simplifying assumptions may vary, we will refer to any particular formulation of the linear problem in the generic form A x = b. Both linear and nonlinear inverse problems are often posed as a minimization of a cost function having the general form, min{|| A x — b|| + A p ( x ) } , (5.3) X where the residual error, A x —b, is minimized over a selected norm and a penalization term, p(x), is added (or multiplied, cf. [148]) in an amount controlled by a tuning parameter, A, to encourage or discourage certain solutions based on a prion assumptions about desirable solution characteristics. The residual term presents the phenomenological cost (i.e. how well the the scattering model and the estimate of the object predict the physical scattering phenomena of the true object), while the penalization term presents the presumptive cost (i.e. how well the estimate of the object agrees with known features of the true object). In this paper we are specifically interested in comparatively evaluating the phenomenological cost of the problem, such that any particular approach to the solution of (5.3) will be benefitted by superior data leading to higher quality images. We wish to evaluate the information available from the scattering data, b, and the ability of the formulation of A to estimate the unknown object from the data. To this end, we employ the singular value 105 decomposition of the scattering matrix. The SVD of the scattering matrix A is, A = U S VH (5.4) where £ is a diagonal matrix of singular values, and the left and right singular matrices U and V are unitary and form orthogonal bases for the data-space and object-space, respectively. The TSVD is a rank-limited, SVD-based inverse obtained by discarding a number of the least significant singular values and associated vectors. Written as a sum of left and right singular vectors, Uj and v.t. an estimate of the object is found from the data b as, (5.5) where in is the trucation index. The information available from the data is limited to the components supported by the span of the data-space basis formed by the left singular vectors, and the estimate of the object is limited to the components supported by the span of the object-space basis formed by the right singular vectors. Noise or errors obscuring the smaller components of the data are amplified by the inverse of the corresponding small singular values and mapped as spurious components of the object estimate. An optimal truncation index (with respect to noise and errors in A and b) can be found by maximizing the fidelity of the object estimate in (5.5) over the truncation index in G {1,..., N}. An objective evaluation of the formulation of A and the information available from b can thereby be compared for variations in system design, noise level, field approximations, among other issues arising in the design of an imaging system and the manner in which the problem is posed. 5.3 Methods The general background of the previous section is now extended to the application of interest. Three-dimensional microwave breast imaging presents a large-scale, nonlinear inverse problem characterized by an ill-posed, under-determined scattering matrix. The nonlinearity of (5.2) arises from the dependence of the total field in the unknown region on the unknown dielectric profile of the object to be imaged. An approximation of the total field is often 106 made to linearize the equations during solution the problem. However, for the purposes of this analysis we wish to evaluate the fundamental limits of microwave illumination in revealing information on the properties and structure of the unknown object. A computational electromagnetic solution is applied to the test domains (described in section 5.4) to acquire simulated channel measurements of the numerical breast phantoms. We further exploit the use of a full-wave field solver to compute the Green's functions of the background media and to record the total fields in the imaging region with the object present. The Green's function for the nth voxel of dimension A of the imaging region is computed using the principle of reciprocity and the background fields, Eb, observed at the nth voxel when the ith observation location is sourced for each orthogonal polarization. G b (r 2 ,r„)= 3 ujfx A Ix E x ( r n , n ) \ I y I z = 0 jx E y ( r n , r i ) \ I y I z = 0 jx Eb(rn, Tj) | I y J z = Q j-yEb(rn,ri)\Ix j-yEb(rn,ri)\Ix i- ^ ( r n , r t ) | ^ ^ 0 Iz=0 0 (5.6) The Green's functions may be computed for either homogeneous or heterogeneous background media. Heterogeneous backgrounds may include known boundary conditions, physical antenna models, or initial estimates of the boundary and properties of the unknown object. The total fields and Green's functions are used to construct a scattering matrix that is nearly exact, limited only by numerical precision, the accuracy of the electromagnetic solver, and the discretization of the integral equation. This exact formulation may be used as benchmark in evaluating the linearizing approximations and initial conditions used in practical formulations. To further idealize the benchmark, we include an exact skin model in the background media so that we can isolate the ability of the microwaves to image the finely detailed tissue structures of the breast interior [84]. The interior breast volume of the background is filled with a homogeneous material with properties equal to the volumeaverage of the breast phantom interior. The scattering system is also analyzed in Section 5.5 without the use of this a priori information to evaluate the cluttering effect of the skin layer's scattering response. 107 The cost of computing the SVD of the large-scale scattering matrix is prohibitive. The number of columns in A is equal to the number of unknowns incurred in the discretization of the imaging volume on a sufficiently fine grid such that it will not limit the image resolution achievable by the system. The number of rows is equal to the number of channel measurements, typically resulting in a highly underdetermined system of equations. Noting that this application often solves (5.3) using methods operating on the normal system of equations, we instead compute the eigen-decomposition of the symmetric matrix AAH. The dual of the normal matrix is chosen for decomposition (rather than the normal matrix AHA) since M <C N for typical array populations and volume discretizations of a 3-D breast region, and AHA. Furthermore, AAH share their M most significant eigenvalues and eigenvectors, and the rank of A cannot exceed M. A A h = WAWff (5.7) The eigenvalues, Aj, of this decomposition are equal to the square of the first M singular values of A, and the eigenvectors of W are equal to the vectors of the left singular matrix U of the SVD of A. The right singular vectors corresponding to the M most significant singular values can be expressed using (5.4) as V[x ... m] = A ^ W A " 1 / 2 . The TSVD object estimate of (5.5) may then computed from the eigen-decomposition for truncation index in < M. (5.8) The fidelity of the estimates provided by (5.8) must be quantitatively evaluated so that an optimal truncation index can be determined. Since the true object profile is readily available in the numerical test regime, the quality of the estimate is best judged in direct comparison to this exact solution, x*. Although such comparisons are often made according to an error norm ||x* — x||, we choose to quantify the fidelity by correlating the estimate and the exact solution using a normalized inner product. (5.9) 108 This is a scale-independent quantity that is equal to the cosine of the angle between the two vectors [84], This fidelity metric tends to demostrate more sensitivity to truncation index than an error norm, and in certain cases the error norm is observed to be non-decreasing over the truncation index. 5.4 Models Investigating microwave imaging in the numerical regime offers two significant advantages: the availability of the total field over the imaging volume, and the ability to create highly realistic breast phantoms with knowledge of the exact structure and electrical properties of the tissues. We use one Class 2 numerical breast phantom ('scattered fibroglandular', ID number 010204) and one Class 3 numerical breast phantom ('heterogeneously dense', ID number 062204) from the University of Wisconsin repository [119]. The x — y plane cross sections bisecting the phantoms on the coronal axis are shown in Fig. 5.1, illustrating the predomination of adipose tissue in the Class 2 phantom and the higher proportion of dense fibroglandular tissue in the Class 3 phantom. 60 40 20 I (a) (b) Figure 5.1 Cross sections of the static permittivity parameter, es, of a) the scattered fibroglandular (Class 2) phantom, and b) the heterogeneously dense (Class 3) phantom. Axes in cm. The dielectric properties of the breast tissue materials used in the construction of the numerical phantoms are based on published clinical data on the properties of ex-vivo tissue samples [18]. The properties of the skin layer are based on published data on dielectric 109 properties of various biological tissues [120]. The frequency dependent properties for each tissue type are fitted to a single-pole Debye model of the complex permittivity from 0.5 to 3.5 GHz. e c M = Coo + T 5 ^ + A - (5.10) The resulting model parameters are listed in Table 5.1. The time constant of the relaxation pole, r, is fixed to 15 ps for all materials. The numerical breast phantoms are mapped from MRI data as described in [19] where the Debye models of Table 5.1 are used as the bounds of the dielectic property ranges. The resulting phantoms are modified exactly as described in [84]. The test phantoms are downsampled to a grid of 2 mm x 2 mm x 2 mm voxels, in a lossless, non-dispersive background immersion with properties similar to vegetable oil (e r = 2.6). On the 2 mm grid, the Class 2 phantom contains about 65,000 interior voxels and the Class 3 phantom contains about 42,000 interior voxels. Table 5.1 Debye parameters (static and infinite relative permittivity, e s and e ^ , and static conductivity, ers) of the tissues modeled in the numerical breast phantoms. Material (percentile) £•00 a s (S/m) Adipose tissue (min) 2.42 2.28 0.0023 Adipose tissue (25th) 4.07 2.74 0.0207 Adipose tissue (50th) 4.81 3.11 0.0367 Adipose tissue (75th) 7.62 4.09 0.0842 Fibroglandular tissue (25th) 36.7 16.8 0.461 Fibroglandular tissue (50th) 49.1 17.5 0.720 Fibroglandular tissue (75th) 54.3 18.6 0.817 Fibroglandular tissue (max) 67.2 29.1 1.38 Skin tissue 40.1 15.3 0.741 Simulated field measurements are acquired using the finite-difference time domain (FDTD) method. Hertzian current sources allow sourcing and receiving of all three Cartesian field 110 components, either individually or simultaneously such that the vector sum achieves arbitrarilydirected source and receive polarizations. The z-directed vertical polarizations at the array locations are illustrated in Fig. 5.2 for each phantom. Horizontal polarizations at each source and observation location are simulated by simultaneously sourcing the x and y components of the current density vector with coefficients such that the source vector is parallel to tangent of the array surface. Since the FDTD field components are arranged on a staggered grid, the x, y, or z field components on both sides of the desired point are sourced (observed), so the active source (receiver) length is 4 mm. Forty source and observation points are distributed over a cylindrical surface sized to fit the breast phantoms with a minimum of 1 cm spacing to the skin surface. Locations are chosen every 45° on each of five coronal rings, with the rings equally spaced along the axial extent of the breast phantom. The second and fourth rings are rotated by 22.5° to increase the element spacing between adjacent rings. The distance between adjacent array locations is about 16 ± 2 mm for both phantoms. 10- 8 E E 5- N N 10 10 y (cm) x (cm) x (cm) (a) (b) Figure 5.2 (a) Class 2 and (b) Class 3 breast phantoms with array locations shown by vertically polarized (z-directed) 4 mm sources. An individual FDTD simulation is run for each source point and source polarization, with a UHF-band modulated Gaussian excitation pulse. Field measurements are recorded for each unique source-receive channel, and for each source configuration the vector electric field is recorded at every voxel in the imaging volume. The time domain data is transformed to the frequency domain at 0.5, 1.0, 1.5, 2.0, and 2.5 GHz, and the resulting phasor data Ill are normalized by the frequency response of the source pulse. The 4 mm source currents behave as short dipoles with an effective gain over frequency measured in the FDTD grid and plotted in Fig. 5.3. For single-frequency imaging formulations, the data at each frequency are gain-normalized to allow a fair comparison under the assumption that a single-frequency antenna design would have a fixed narrowband gain independent of operating frequency. In the multiple-frequency imaging formulations the data are not gain-normalized and remain weighted by the gain response toward higher frequencies where the source has higher radiation efficiency. Although data weighting is itself a good candidate for the TSVD inversion analysis, the gain response of Fig. 5.3 is characteristic of the electrically-small broadband antennas likely to be used in a densely populated array for multiple-frequency measurement. While the low gain of the short dipole sources results in channel measurement amplitudes averaging about -70 dB, the computational noise floor of the single-precision FDTD simulations is about -180 dB and the MATLAB single-precision noise floor of the post-processing and TSVD analysis is about -140 dB. Figure 5.3 Measured gain of the FDTD sources. 5.5 Results The analysis described in Section 5.3 is applied to field data collected from the Class 2 and Class 3 numerical breast phantoms defined in Section 5.4. We analyze both single-frequency and multiple-frequency formulations of (5.2), with the unknown contrast expressed in terms of complex permittivity, e c = fc2/cu2/i0, rather than the square of the wavenumber. In the 112 single frequency case, the unknown elements of x are complex permittivity contrasts of the form (e* — £c), where ec = e' — je". In the multiple-frequency formulation, the dispersive frequency dependence of the complex permittivity is parameterized [76] using the singlepole Debye model [21]. In addition, we investigate the use of the analysis in the selection of scattering measurements and the formulation of the inverse problem. The analysis is further employed to evaluate the effect of noise in the data, modeling errors, and field approximations. 5.5.1 Single-Frequency Formulation Single-frequency scattering systems at 0.5, 1.0, 1.5, 2.0, and 2.5 GHz are constructed for analysis using fields from the Class 3 phantom test bed. The systems are formulated with phasor data from vertically polarized channels and the scattering matrices are constructed using tensor Green's functions and vector total fields corresponding to the vertical source and receive polarizations. The system consists of 820 equations corresponding to 780 bistatic channel measurements and 40 monostatic channel measurements. The unknowns are the complex-valued permittivity contrasts over the N voxels in the imaging region, V. The imaging region is restricted to the volume interior to the skin layer. The fidelity is computed by applying Eqn. (5.9) to the N x 1 vectors of exact and estimated contrast in complex permittivity. Figure 5.4 charts the fidelity of the TSVD image as the truncation level and signal-tonoise ratio (SNR) are varied for each operating frequency, showing an improvement in the peak fidelity as frequency increases. In Fig. 5.4a, the truncation index corresponding to peak fidelity increases with frequency. The optimal truncation index provides a measure of the effective rank of the data and scattering matrix. The consideration of measurement sensitivity in the trade off between penetration and resolution is well illustrated by examining the fidelity of the TSVD estimate over a range of SNR. Figure 5.4b shows the effect of additive white noise in the data on the fidelity of the estimate for data sets at 0.5, 1.0, 1.5, 2.0, and 2.5 GHz. The plots demonstrate that the required measurement sensitivity increases as 113 frequency increases. In these results, note that the scattering matrices and the data have been normalized by the source gain, as explained in Section 5.4, such that the effective source and receiver gain is 0 dB at each frequency. The results may be interpreted for practical antennas by subtracting the transmit and receive gain from the SNR thresholds indicated by the knee of the curves in Fig. 5.4b. For example, the fidelity versus SNR curve for 2.5 GHz rolls off below about 70 dB. If the gain of the source and receive antennas are taken from Fig. 5.3 (about -25 dB) then the required measurement sensitivity is roughly 120 dB. (a) (b) Figure 5.4 Fidelity of TSVD estimates of complex permittivity using single-frequency data from 820 measurement channels of the Class 3 phantom versus a) index of singular value truncation, and b) SNR. The SNR results are normalized to the gain of the source at each frequency. The required measurement sensitivity for an operating point on one of the SNR curves is determined by subtracting the gain, in dB, of the transmit and receive antennas. Since the spatial frequency content of the singular vectors generally increases with singular value index [96], the increase in optimal truncation index with frequency demonstrated by Fig. 5.4a suggests that the resolution of the resulting images will improve as frequency increases. The trend in image resolution with frequency is visualized using the 3-D rendering of the TSVD estimate vector, x. The 3-D images at each frequency are represented in Fig. 5.5 by bisecting coronal cross-sections of the imaging volume. The corresponding cross-section of the exact phantom profile is shown for comparison in Fig. 5.5a. 114 (d) (e) (f) Figure 5.5 Coronal cross-sections of static permittivity, e s , a) the actual Class 3 breast phantom, and the optimal TSVD projections of the scattering data at b) 0.5, c) 1.0, d) 1.5, e) 2.0, and f) 2.5 GHz. Axes in cm. In addition to the comparison across frequencies, we investigate the information held by single-frequency scattering data with respect to the relative permittivity and effective conductivity parameters. The fidelity of each parameter is computed by applying (5.9) separately to the real and imaginary parts of the exact and estimated complex permittivity contrasts. Figure 5.6 clearly indicates that the effective conductivity is more poorly estimated by the data than the relative permittivity throughout the UHF frequency band. 5.5.2 Multiple-Frequency Formulation The TSVD analysis can be used to inform the choice of a particular formulation of the problem to be inverted. For example, simultaneous multiple-frequency solutions can be compared to single-frequency formulations and comparisons can be made of various choices of assumed frequency dependence of the unknown dielectric properties. Over the UHF band 115 E at 0.5 GHz - a „ at 0.5 GHz 0.8 erf • £ at 2.5 GHz a „ at 2.5 GHz eff >,0.6 a> •o 0 0 J. 20 40 / / 60 80 100 SNR (dB) Figure 5.6 Fidelity versus SNR of the TSVD estimates of relative permittivity, eT, and effective conductivity, creff, using single-frequency data at 0.5 and 2.5 GHz. The SNR results are normalized to the transmit and recieve gain at each frequency. the single-pole Debye model is a natural choice for constraining the frequency dependence of biological tissues, whose responses are dominated by their water content and the single-pole behavior of water in the band. With the time-constant of the pole fixed for all tissues, the unknowns are formulated as the three remaining free parameters of the Debye model: static permittivity, e s , infinite permittivity, e ^ , and static conductivity, a s . The Debye model is alternatively expressed with the static pemittivity parameter replace by a delta permittivity parameter, Ae = £s — £oo, which corresponds directly to the dispersive term in the Debye relation (5.10). We have included all three Debye permittivity parameters in this analysis for completeness, though since they collectively represent two degrees of freedom only Ae and Soo are estimated. A multiple-frequency scattering system is contructed using field data at 0.5, 1.0, 1.5, 2.0, and 2.5 GHz from the Class 3 phantom test bed. The systems are formulated with phasor data from vertically polarized channels and the scattering matrices are computed using tensor Green's functions and vector total fields corresponding to the vertical source and receive polarizations. The system consists of 8200 equations consisting of the real and imaginary equation pairs for each of 820 channels at each of the five frequencies. The 3N x 1 vector x is the concatenation of three JV x 1 vectors of the real-valued Debye parameters. 116 The construction of the linear system of equations in the Debye parameter unknowns is described in detail elsewhere [83]. Figure 5.7 Coronal cross sections of exact phantoms and TSVD estimates of static permittivity, £s, using multiple-frequency data at 0.5, 1.0, 1.5, 2.0, and 2.5 GHz. a) exact Class 2 phantom, b) Class 2 estimate, c) exact Class 3 phantom, d) Class 3 estimate. Axes in cm. As in the single-frequency cases above, the imaging region is restricted to the volume interior to the skin layer. Figure 5.7 shows coronal cross-sections of the estimated interior of the Class 2 and Class 3 phantoms. The fidelity is computed by applying (5.9) to the 3N x 1 vectors of exact and estimated Debye parameter contrasts. Figure 5.8a shows the fidelity versus SNR of each parameter of the optimal TSVD solution to the five-frequency Debye formulation of the inverse problem. Similar to the single-frequency results in Fig. 5.6, the static conductivity parameter of the multiple-frequency Debye formulation is shown to be poorly estimated by the scattering data compared to the estimation performance for the permittivity parameters. The effect of the selection of frequencies on the performance of the 117 (a) (b) Figure 5.8 a) Fidelity versus SNR of TSVD estimated Debye parameters for the Class 3 phantom using multiple-frequency data at 0.5, 1.0, 1.5, 2.0, and 2.5 GHz over 820 measurement channels, b) Comparison of the fidelity versus SNR of the static permittivity parameter by frequency selection. The fidelity of the relative permittivity for the 2.5 GHz single-frequency formulation is included for reference. Debye-parameterized multiple-frequency formulation is shown in Fig. 5.8b by comparison of the fidelity versus SNR curves for the static permittivity parameter. The fidelity versus SNR curve for relative permittivity in the 2.5 GHz single-frequency formulation without gain normalization is included for reference. The comparison shows only a slight degradation of the fidelity performance as well as an improvement in noise performance due to the inclusion of lower frequencies. 5.5.3 System Design The design of an imaging system is largely described by the selection of a set of measurement channels. A channel is defined by the locations of the source and receiver, the polarizations of the source and receiver, and the frequency of the source. The selection of frequency is investigated above; here we compare the information for varying selections of polarization and spatial channels. 118 The effect of the polarization on the information content of the scattering data can be studied by sourcing and observing the field components corresponding to the polarizations of interest. Any combination of co- and cross-polarization channels can be analysed in this way. For illustration, we compare vertical source and receive polarization to horizontal source and receive polarization. For the cylindrical array of source-receive locations described in Section 5.4, the Ez field component is defined as the vertical polarization. The horizontally polarized illumination is achieved by sourcing a both x- and y-directed current densities in a proportion such that vector sum of the current density is directed along the tangent of the cylindrical surface in the x — y plane. Similarly, the horizontally polarized observations are made by summing the appropriate proportions of the Ex and Ey fields corresponding to the observation location. For each case of polarization, the multiple-frequency, Deybe parameter, exact field formulation is constructed using data at 0.5, 1.0, 1.5, 2.0, and 2.5 GHz from all 820 spatial channels. The fidelity of the static permittivity estimate for each polarization versus SNR is shown in Fig. 5.9. The performance of the polarizations is similar, indicating that there is comparable imaging information available from either vertically or horizontally co-polarized scattering data. Figure 5.9 Fidelity of the static permittivity parameter versus SNR for vertical and horizontal polarizations. The TSVD estimates are computed from multiple-frequency data at 0.5, 1.0, 1.5, 2.0, and 2.5 GHz over 820 measurement channels. 119 The spatial selection of channels requires consideration of the array shape and the population density of elements in the array. In addition, it is of interest to consider whether certain channels contain better or different information than others. We investigate one example of channel selection by distinguishing through and reflect channels. In the cylindrical array configuration with the origin on the cylindrical axis, we define the through channels as the channels in which the receive location is equal to or greater than 90° from the source location, i.e. 14>src — (j)obs | > 7r/2, and the reflect channels are defined as the channels in which the receive location is equal to or less than 90° from the source location, i.e. \(f)src — (jfbs| < 7r/2. For the 40-element array described in Section 5.4, there are 556 through channels and 556 reflect channels. The fidelity versus truncation index of each channel set is charted in Fig. 5.10 for the full 3iV x 1 vector of Debye parameters. Although the peak fidelity is about equal each case, the difference in the optimal truncation indices suggest that through channel data favor lower spatial frequency solutions while reflect channel data favor higher spatial frequency solutions. This result is consistent with the Fourier-space treatment of microwave imaging in [33, 34] where it is shown for linear arrays of sources and receivers that the through data captures lower spatial frequencies and the reflect data captures higher spatial frequencies. Inspection of the cross-sections of the through and reflect image estimates in Fig. 5.11 suggests that the information content of each data set is complementary and that both should be used when possible. 5.5.4 Error Evaluation The preceding results have shown the use of the TSVD analysis for noise errors in the data. Others sources of errors in the data vector or in the scattering matrix can be analyzed in a similar manner. We illustrate the evaluation of such errors with a few representative examples. The performance of the scalar-field approximation is investigated, image degradation due to an imperfect estimation of the skin layer is visualized, and the Green's function approximation of the Born iterative method (BIM) is evaluated for the breast imaging application. 120 Figure 5.10 Fidelity of the parametric Debye estimate versus SNR for "through" channels and "reflect" channels using multiple-frequency Debye formulation for data at 1.5, 2.0, and 2.5 GHz. (a) (b) Figure 5.11 Coronal cross sections of multiple-frequency (1.5, 2.0, and 2.5 GHz) optimal TSVD image estimates for a) through channels and b) reflect channels. Axes in cm. The scalar field approximation is a computationally advantageous simplification of the inverse problem formulation which introduces an error in the scattering matrix. The validity of the scalar field approximation, in which the Green's functions and total fields calculated in the imaging region are assumed to be scalar, is compared to the exact formulation in which the tensor Green's functions and vector total fields are used. In our scalar example, only the z-component of the electric field and the zz-component of the Green's functions are used in the construction of the scattering matrix. For vertical source and observation polarizations, 820 data channels of the Class 3 phantom are used to evaluate the fidelity of the complex 121 permittivity estimate versus truncation index for the scalar and vector scattering matrix formulations. Figure 5.12 shows the limitation of the scalar-field-approximated scattering matrix in mapping the information available from single- and multiple-frequency data. (a) (b) Figure 5.12 Comparison of the fidelity versus truncation index using the exact vector field formulation and the approximate scalar field formulation. Class 3 phantom for a) single-frequency data at 2.5 GHz. b) multiple-frequency data at 1.5, 2.0, and 2.5 GHz. The results in Sections 5.5.1 to 5.5.3 used simplifying assumption that the skin region is exactly modeled in the background media, allowing us to isolate the performance of a microwave system in imaging the detailed heterogeneous tissue distributions of the breast interior. We can compare this to the performance of the system with no a priori information about the skin region by using a homogeneous background of the immersion medium. The imaging volume is defined to extend 6 mm away from the skin surface. The optimal TSVD image of multiple-frequency data at 1.5, 2.0, and 2.5 GHz is show in Fig. 5.13 with and without the background skin model. In the cross-section shown in Fig. 5.13b the skin layer is only partially reconstructed by the data. The resulting error is evident as clutter surrounds the fibroglandular regions that are seen to be correctly reconstructed by comparison to Fig. 5.13a. Figure 5.14 shows the decrease in image fidelity due to the imperfect imaging of the skin layer. 122 (a) (b) Figure 5.13 Coronal cross-sections of the estimated Class 3 phantom for multiple-frequency data at 1.5, 2.0, and 2.5 GHz. a) Background includes exact skin region and homogeneous interior of the volume-averaged tissue properties, b) Background is the homogeneous immersion medium. Axes in cm. Figure 5.14 Fidelity of the e s estimate versus SNR for estimation for the Class 3 phantom with and without the exact skin model included in the background media. Multiple-frequency data at 1.5, 2.0, and 2.5 GHz are used in the estimation. The analyses presented in this paper exploit the availability of the exact total field over the entire unknown imaging region in order to attain an evaluation of an upper bound on the ability of the microwave scattering data to image the breast. field term, In practice, the total in the scattering relation of Eqn. (5.2) is unknown and dependent on the unknown contrast. Solutions to the inverse scattering problem must therefore approximate the total field term. The Born approximation [16] is commonly used, for which the total field is replaced the by incident field in the background region. The Born approximation 123 is known to be inaccurate for objects with high contrast relative to the background, which is certainly the case for heterogeneous biological tissues. The error due to the use of the Born approximation in the microwave breast imaging application can be evaluated by the TSVD analysis. A multiple-frequency scattering matrix is formulated with the exact total fields replaced by the fields measured in the background. The background includes the skin region and the homogeneous, volume-averaged properties in the breast interior. The effect of the Born field approximation on the achievable image fidelity is demonstratedin Fig. 5.15 for the static permittivity and static conductivity Debye parameters for data from 820 channels at 1.5, 2.0, and 2.5 GHz. The comparison underscores the limited utilization of information in the data when the Born approximation is employed. Iterative methods of nonlinear optimization such as the Born iterative and distorted-Born iterative methods are often used to overcome this limitation. In the context of such methods, Fig. 5.15 gives an expectation of performance for the first iteration. The performance of the conductivity estimation is significantly more degraded than the permittivity estimation and suggests the use of an iterative technique that favors permittivity solutions in early iterations. The lower optimal truncation index of the Born approximation formulation suggests that low-pass regularization should be increased during early iterations. e , exact fields s s a , exact fields e s , Born approx. a s , Born approx. 0.4 0.2 1000 2000 3000 Truncation Index 4000 Figure 5.15 Fidelity of the £s and a s estimates versus truncation index with and without the Born field approximation for multiple-frequency data at 1.5, 2.0, and 2.5 GHz. 124 5.6 Discussion A number of insights are provided by the results of Section 5.5. In the single-frequency formulations, we see the tradeoff between noise performance and image resolution. The measurement sensitivity required to image experimental or clinical data can be estimated by analysis of the data over SNR such as in Figs. 5.4b and 5.8. The single-frequency performance relative to the singular value truncation index in Fig. 5.4a shows that as frequency is lowered, fewer and fewer of the singular solution vectors corresponding to the higher indices are informed by the data. Since the spatial frequencies of the singular vectors generally increases with singular value index, this trend suggests that low frequency formulations should be increasingly regularized by low-pass filtering to avoid noise components in the part of the singular spectrum for which the data has no information. Furthermore, this trend may suggest that a frequency-hopping method of imaging may be effective, in which the solution for each frequency of operation is filtered accordingly. The results of Figs. 5.6 and 5.8a reveal that the conductivity properties of the unknown objects are much more difficult to estimate from the scattering data than the permittivity properties. Given the rough correlation between the permittivity and conductivity properties of biological tissues, a higher fidelity estimate of the conductivity might be obtained from the permittity estimate rather than from direct estimation from the data. The goal of the analyses presented in this paper is to provide a sense of the upper bound of the imaging performance that can be expected from a microwave breast imaging system. However, it is important when speaking of upper bounds on imaging performance to clarify that the TSVD analysis optimizes the fidelity metic over the i 2 norm. Cost functions and regularizations based on the i 2 norm are known to have a low-pass filtering effect on the spatial frequencies of the solution. The TSVD is itself a method of inversion and regularization, with an effect similar to the Tikhonov method of regularization [96]. Alternative formulations and techniques of inversion intended to improve resolution and edge preservation have been investigated elsewhere and have potential to out-perform the results illustrated in Section 125 5.5. However, all approaches are equally limited with respect to the information available from the data. The TSVD analysis is adept at the comparative evaluation of that information. It is therefore a powerful tool for the design of the measurement system and problem formulation, irrespective of any particular choice of imaging technique used to recover the unknown object. There are a number of idealizations used in our analysis in service of the goal of evaluating upper-bound imaging performance. Practical performance expectations are tempered by some of these idealizations. The exact total fields in the imaging region are available in a computational study, but in practice they must be estimated as discussed in the analysis of the Born approximation in Section 5.5.4. The Hertzian current sources simplify polarization analysis, but exhibit none of the mutual coupling or modeling errors arising with the use of actual antennas. The inclusion of the exact skin model in the background allows us to isolate the ability of the microwaves to resolve realistic tissue structures in the breast interior, but the skin tissue presents a scattering response that is quite large compared to interior scattering and therefore small errors in the modeling or imaging of the skin layer can result in significant image distortion, as seen in Section 5.5.4. Although idealizations such as these may be seen as overly contrived from a practical perspective, the effect of each idealization and source of error can be analyzed by the method and their performance cost can be objectively quantified. In addition, the idealized images provide a better benchmark for practical imaging results than comparison to the exact profile of the numerical phantoms. Many important issues in microwave imaging design are not included in this investigation, but most can be analyzed in a similar manner. For example, the shape and population of the antenna array can be optimized using the TSVD analysis. The properties of the immersion medium has significant implications on the antenna design and the coupling of illumination into the breast interior. The TSVD analysis can be used to evaluate the imaging performance for various selections of immersion media. Although inverse scattering techniques for breast imaging typically operate at frequencies in the UHF band (300 MHz to 3 GHz), tissue-penetrating radar techniques (e.g. [15,114]) are often designed in the UWB band (3.1 126 to 10.6 GHz). Since both methods are based on the same scattering phenomena it is a natural extension of the TSVD analysis to consider the usable information content of data in the UWB band, particularly with respect to measurement sensitivity. These and many other considerations in system design and problem formulation are readily analyzed by the techniques employed herein. 5.7 Conclusion A method of analysis using the TSVD has been presented in which the scattering data may be systematically evaluated independently from the development of imaging methods to be applied to the data. The method was used to compare the information available from the data for frequencies in the UHF band, for single- and multiple-frequency formulations, for different selections of measurement channel locations, and for errors due to measurement noise, field approximations, and modeling inaccuracies. The results illustrated the trade-off between field penetration and image resolution over frequency and indicated the measurement sensitivity necessary to utilize the data at each frequency. Limitations in data utilization were shown when using the scalar field approximation or Born approximation in the formulation of the breast imaging problem. In addition to the considerations investigated in this paper, the TSVD analysis of the scattering data and inverse formulation can be used to systematically optimize and evaluate a microwave imaging system, from the design of the test fixture to the analysis of noise and modeling errors. Any microwave breast imaging method stands to benefit from the betterment of the information in the data available from a system developed in this manner. 127 Chapter 6 Future Work 6.1 Summary of research At the outset of this research, the UWCEM laboratory published the first fully 3-D microwave breast imaging results for realistic breast models [21]. Others has previously imaged 3-D clinical data, but only in 2-D slices. Prior experimental and numerical tissue models were typically simple shapes with inaccurate dielectric properties. The UWCEM contribution was the first imaging work to use anatomically realistic numerical breast phantoms, and the first to use realistic tissue properties as reported by UWCEM [1,18]. A single spherical tumor model was added to each phantom. The numerical models allowed perfect registration of the resulting images with the known spatial profiles and known dielectric properties of tissue structures in the breast phantoms. However, the low resolution of the resulting images made their interpretation ambiguous. It was not clear which parts of the reconstruction were related to scattering from the healthy tissues, and which parts were related to the malignant tumor. Furthermore, a custom spatial basis was used to reduce the degrees of freedom of the solution. The custom basis intentionally limited the available image resolution according to the wavelength of the illumination in the breast tissues, so it was not clear whether the image resolution was limited artificially by the basis or fundamentally by the interrogating waves. After refactoring our forward solution and inverse solution code, our 3-D numerical breast phantom solutions were no longer restricted by computational hardware limitations. The memory and storage efficiency improvements allowed for an imaging solution on a relatively 128 fine voxel grid with increments on the order of a tenth of the shortest wavelengths. On this voxel basis, the resolution of the image solution is likely to be either fundamentally limited by the interrogating wavelengths or limited by the method of optimization and choice of regularization of the inverse solution. The numerical breast phantoms were imaged in Chapter 3 without any malignant tissues in order to remove any ambiguity in the interpretation of tissue types in the reconstructions. However, the resolution of the resulting images could not be conclusively attributed to a single causal factor. The improvement in the resolution could be due to either the higher resolution of the solution basis or improvements in the techniques of optimization. The limitation of the resolution also could not be conclusively associated with either the interrogating wavelengths or the low-pass behavior of the regularization technique used in the inverse solution. Although the resolution was still too low to recover the smallest millimeter-scale details of the tissue structures, the imaging performance was markedly improved and reconstructed areas of fibroglandular density clearly matched the true phantom profiles. These results showed that, at minimum, the microwave images were a reliable representation of density of the breast, an important factor in cancer risk assessment. The revelation of the dielectric properties studies [1,18] of reduced contrast expectations for malignant tissues further complicated the state of the research in microwave breast imaging. Many prior investigations throughout the community had assumed not only inaccurate dielectric properties for the constituent breast tissues, but had grossly overestimated the contrast between malignant and normal glandular tissues and underestimated the contrast between malignant glandular tissues and adipose tissues. The correction of these tissue properties in numerical models made the imaging and detection of tumors in the breast more difficult by reducing their scattered signal amplitude while increasing the cluttering signals scattered by the healthy fibroglandular tissue structures. In Chapter 4 we showed for a spherical scatterer of malignant tissue surrounded by healthy tissue that the scattered signal level is a function of both tumor size and contrast. The signal scattered from the sphere rolls off sharply when the diameter falls below one half-wavelength of the illumination in the 129 background medium, suggesting a fundamental limitation of the imaging resolution available from scattering observations. We then investigated the use of exogenous contrast agents to perturb the scattering response of the spherical tumor models placed in locations of dense tissue in numerical breast phantoms of each density classification. The tumors were detected in each phantom by differentially imaging the phantoms before and after applying the assumed effect of a contrast agent to the tumor properties. Resolution of the differential tumor image was improved in the dense phantom cases where the tumor was placed in a dense background, resulting in a larger magnitude of the reconstructed contrast. However, the noise performance was better for the adipose-dominated phantoms in which the higher contrast between tumor and fat results in higher scattered signal strength. Despite the qualified success of this preliminary contrast-enhanced detection study, the imaging results further underscored the deficit of resolution available from the microwave illumination, as the differential tumor reconstruction revealed a fairly large point spread function in a realistic breast interior. Furthermore, the same measurement design and inverse solution were used as in the study of normal tissue phantoms of Chapter 3, and thus there remained the same ambiguity in interpreting the limitations of resolution and image fidelity. In response to the deficit of resolution from the microwave inverse scattering solution, researchers in the community increasingly turned to the pursuit of inverse techniques capable of improving resolution and preserving sharp edges at the transitions between different tissue types. What often remain fixed in these investigations was the design of the set of scattering data being imaged. However, if there is an optimal performance to be striven for from a given imaging method, the natural question then becomes whether there is an optimal data set that can be evaluated as independently as possible from any particular imaging technique. In Chapter 5 we explored the fundamental information content of scattering data using a linear algebraic analysis of the scattering matrix and scattered field data vector. The analysis was applied comparatively to various selections of scattered data and formulations of the scattering matrix. The effective rank and estimated resolution of single-frequency and multiple-frequency formulations were shown for frequencies across the UHF band. The effect 130 of various field approximations and data errors were evaluated. In particular, the imaging performance was charted over signal-to-noise ratios, showing that useful imaging performance depends on high measurement sensitivity. In addition, the estimations of upper-bound image performance provide a much more appropriate benchmark for imaging system development than comparison to either exact profiles or prior performance. The method of analysis is also quite useful in the systematic evaluation of system design parameters, data and modeling errors, and formulations of the inverse problem. Perhaps most importantly, the estimation of the capability of a system provides an immediate evaluation of the clinical utility of its best-case imaging performance. Our preliminary analyses are not yet conclusive as to the efficacy of UHF-band microwave imaging systems, but they suggest that much can be done to select more information-rich and efficient sets of scattering data. 6.2 Directions for future research Going forward, there remains much to be learned about performance and limits of microwave imaging of the breast, and short-range microwave imaging in general. Although the developments to date have been largely inconclusive as to the viability of the modality, they justify continuation of general directions in breast imaging research. Alternative methods of optimization and regularization must continue to be pursued by the community with the goal of resolution enhancement and edge preservation using a priori information about the spatial structure of the breast. Better use should be made of anatomically realistic numerical and experimental breast phantoms, and careful image registration techniques should be developed for pre-clinical studies. Using an MRI as a reference image and a numerical phantom derived from the MRI, clinical images can be spatially registered against the MRI, and registered against the numerical image to corroborate ex-vivo tissue properties with in-vivo estimates from a quantitative image. The work presented in this dissertation leaves room for extensions of the studies and development of the methods, and has motivated new research directions related to the application. The following sections highlight some of the possible topics for future research. 131 6.2.1 System Design via Scattering Information Analyses The evaluation of scattering data in Chapter 5 only scratched the surface of the useful analyses to be made by the technique. The study should be extended to the higher frequencies of the UWB band by acquiring fields from phantoms on a 0.5 mm grid using FDTD. The design of a measurement system should be thoroughly optimized for considerations including channel selection, immersion medium, antenna design, array shape and location, modeling errors, and measurement errors. The formulation of the inverse problem should be thoroughly optimized for considerations including unknown parameter selection, frequency-dependence constraints, data weighting, scattering matrix preconditioning, and pre-regularization of the solution basis. Furthermore, some issues pertaining to specific inverse techniques can be analyzed. In Gauss-Newton techniques, the initial guess for the background media in the forward solution can be evaluated for first-iteration optimization performance. The utility of updating the Green's functions by recomputing them as the object is reconstructed can be shown by iterating the analysis once to use the computed Green's functions of the object estimated by the first analysis. 6.2.2 Custom Spatial Basis for Pre-Regularization Pre-regularization of the inverse problem was previously investigated for the application by using a custom spatial basis [21]. However, the design of the basis used a guess about the resolution and arbitrarily-shaped generating functions. A better custom basis can be constructed using the propagation and scattering phenomena in a model of the breast. If a residual scattering term, ||Ax — b||, is used in a cost function, the product A x is limited by the span of A, so x need not have more degrees of freedom than the rank of A. The space spanned by A could therefore be used as a custom basis for x. For example, such a basis is formed by a truncated set of right-singular vectors with cardinality equal to the effective rank of A. The practical difficulty in this approach is the unavailability of the exact object fields needed for the exact formulation of A. The fields may approximated to generate an approximation to the span of the exact scattering matrix for a breast, but the breast interior 132 must be approximated with a medium having permittivity near that of the highest percentile fibroglandular tissue, and conductivity near that of the lowest percentile adipose tissue. The resulting field approximations give an upper bound of spatial frequency and lower bound on conductive loss in the tissues. Given the eigen-decomposition AAH = W A W H , a rank R basis consisting of a set = AffWA-1/2. of right-singular vectors is found as $ = The basis $ replaces the original voxel basis in the inversion problem by expressing the unknown profile update in terms of unknown basis coefficients, (p, and transforming the linear system accordingly. The estimated coefficients are then transformed back to the basis of the forward solution in order to update the background profile. (e* — eb) j ^g* = <f>cp gb^ [J$] ip (e* - eb) jpmeas jjjcalc = Emeas - Ecalc = A preliminary investigation of this technique showed that the Gauss-Newton imaging results for a numerical breast phantom were substantially similar when using the voxel basis or the custom eigenvector basis. Cross-sections of the two reconstructions are compared in Fig. 6.1. The similarity in imaging performance of the two bases demonstrates that the custom basis offers a substantial savings in the computational complexity of the inverse problem at little cost to the quality of the solution. The pre-regularization achieved by the custom basis may even improve imaging performance since the inverse solution can be more aggressive in its approach to regularization. It worth noting the strong similarities between the formulation of the problem on an eigenvector spatial basis and the formulation of the problem in the Fourier spatial frequency domain. Microwave imaging has previously been treated in the Fourier space as a diffraction tomography problem using linear arrays of sources and receivers [33,34], It may be of future 133 interest to extend the diffraction tomography treatment to a cylindrical array, and to solve the inverse scattering problem in the Fourier spatial frequency domain. 0 2 4 6 8 10 12 (a) 0 2 4 6 8 10 12 (b) Figure 6.1 Reconstructions of static permittivity, e s , the Class 2 (ID 010204) numerical breast phantom for (a) the 2 mm voxel basis of over 65,000 basis functions, and (b) the custom eigenvector basis of fewer than 2,300 basis functions. 6.2.3 Generalization of Radar and Inverse Scattering The microwave imaging problem for the breast has been separately posed as a shortrange radar problem and as an inverse scattering problem. Work in both areas is ongoing in the community, with radar techniques typically operating in the UWB band and inverse techniques typically operating in the UHF band. Given that both problems are based on the same scattering phenomena, it is unsurprising that they share many attributes. For example, radar techniques generally require an approximate propagation model for focusing scattered fields to their source locations in the breast. This is analogous to the Green's function used in inverse scattering. The weights selected for the beamformer comprise a type of inverse method in which the system is diagonalized to localize source locations. In fact, the two methods are so similar that it would likely be of use to the community to generalize the microwave imaging problem with a common theoretical framework for both radar and inverse scattering. This is of particular importance in the comparison of their 134 performance given that radar methods are not iterative and therefore cannot overcome the nonlinearity of the wave scattering phenomena. The frequency domain radar form of the integral electric field scattering relationship in Eqn. 6.1 pits the bistatic propagation matrix, P , against the Green's function - total field product, G b E \ of Eqn. 6.2. E m e a s ^ _ E i n c ^ = k2 J (6.1) P(r, r')S(r')dr' V E™a°(r) - Eb(r) = k20 J G\r, r')E\r') (ec{r') - ebc{r')) dr' (6.2) v The source function, S, is the qualitative unknown scattering strength in the radar solution, while the dielectric contrast function, £ c (r') — £b{r'), is the unknown in the quantitative inverse scattering solution. One might further note the similarity of these equations to the direct and source equations in contrast-source inversion (CSI) methods (e.g. [45,80]). The breast hyperthermia problem [154] is also closely related to these equations in that it is the one-way analog to the round-trip imaging problem. Given the potential for a single test fixture to both administer hyperthermia treatments and image the breast interior, there is motivation to formalize their mathematical relationship. 6.2.4 Data Acquisition via Beamforming The most immediate application of the generalization between radar and inverse scattering is to incorporate the beamforming techniques of radar methods in the measurement strategy for inverse scattering methods. Although beamforming can be done synthetically by applying weights to a set of single-source channel measurements, it may be advantageous to the signal-to-noise ratio of the data to acquire data by scanning the beamformer over the imaging volume. Consider the difference between two single-frequency, phase-only beamforming scenarios for focused illumination at point rQ in the imaging domain, with observation at point Tj. In the first scenario, antennas numbered i — 1 , . . . , N of an iV-element array are sourced simultaneously and are phase-adjusted with source terms such that 135 the total fields due to each antenna are in-phase at r0. In the second scenario, no phase term is added to the sources and each antenna is sourced individually for a total of N independent measurements at r,-. Beamforming is then performed numerically by applying the phase factors to each measurement in post-processing. The goal of each beamformer is to selectively measure the scattering from r0 by coherently adding the contributions due to all N source antennas. Define to be the portion of the scattered field arriving at a receiver at Tj from the object contrast at r0 due to a source antenna at V{. Gb(rj, r0)Et(r0, r*) (k\r0) - k2b(r0)) (6.3) Let n be a random sample of complex white Gaussian noise in each measurement, and let be an error term representing the scattering from r ^ rQ due to imperfect focusing. If the phased-array beamforming process is employed during the measurement, then the data at receiver j, d°, resulting from focusing the array to r0 contains a single noise sample. N (6.4) If numerical beamforming is employed during post-processing, then the data resulting from focusing of N measurements to r0 contains an independent noise sample for each measurement. N i=i Clearly the noise performance of (6.4) is superior to (6.5), ceteris paribus. The use of beamforming during the measurement could improve the SNR by as much as 101og10(7V), provided that the system dynamic ranges at the source power used in each method are equivalent. 6.2.5 Generalization of Sourcing Functions The commonality between radar and inverse scattering approaches to microwave imaging motivate the use of excitation strategies that aren't limited to the one-at-a-time antenna sourcing of inverse scattering measurements. In the general scattering relation, any number 136 of simultaneous sources can be used, just as in beamforming for radar or hyperthermia techniques. Each channel measurement then becomes defined as the scattered field observed for a particular source configuration. The source configurations should be generalized as source functionals, defined as either discrete sets of points or continuous functions over a surface. In reality, an antenna is not a point source but a continuous source function defined over the active area of the antenna. Unless the antenna is electrically small, the point- source assumption for a practical antenna may actually introduce errors since the sourced (or received) signal is actually related to a surface integral of the currents (or fields) over the active area. This potential error is avoided by inverse methods that compute the Green's functions from forward solutions that include antenna models. However, it worth considering whether the averaging of the field over the active area effected by the integration kernel actually degrades the information available from the measurement. It may be advantageous to use electrically small antennas that better approximate point sampling of the scattered fields. However, the reduced gain of electrically small antennas reduces the measurement sensitivity and thus degrades the potential imaging performance of the system. A novel approach would be to excite and observe source functionals, or modes, over the entire surface of the source space. Upon examination of the scattering equation, e.g. (6.2), it becomes apparent that the span of the scattering matrix will be related to the mode-diversity of fields excited over the imaging region. The mode-diversity for a scattering system could be analyzed by the eigen-factorization or Fourier representation of the steadystate field distributions set up in the imaging volume for each source configuration. (Of course, this is exactly analogous to beamforming, in which a linear combination of these modes is used to focus the fields to a spot.) Ultimately, a continuous source functional approach does not translate well to practical microwave hardware. However, the short range, inward-looking antenna arrays used in microwave breast imaging should be considered in the general perspective of source functionals over the array surface. Furthermore, commonly used concepts such as far-field radiation and beam patterns do not have useful meaning in this 137 type of array configuration. Instead of evaluating isolated antenna performance in these terms, it is far more appropriate to evaluate the modes and amplitudes of the interrogating fields excited in the imaging volume in the presence of the complete array and physical test fixture. 6.3 Conclusion At the conclusion of this dissertation, though there remain some unanswered questions for basic research, microwave breast tomography has entered a translational research phase. Contributions must continue be made to system design, algorithmic technique, and clinical evaluation, and the theory and analysis of performance expectations or limitations must be be advanced. It is clear from the overview in Chapter 2 of the state of the art in microwave breast imaging that there remain significant steps in bringing the modality to clinical relevance. The dielectric contrast between malignant and normal glandular tissues is smaller than previously assumed, but it will require a well-registered comparison of images from more mature microwave imaging technology to answer such concerns more definitively. The resolution demonstrated by early clinical results and images of highly realistic numerical phantoms is insufficient for the reconstruction of the fine detail of tissue structures - detail that may be vital to the subjective interpretation of images for diagnostic evaluation. Although resolution can be further improved, the full dynamics and trade-offs involved in the factors influencing achievable resolution are not yet fully understood in the context of microwave breast imaging. While there is work to be done in improving image resolution, early results have supported the value of microwave imaging in classifying breast density, an important indicator of cancer risk. In particular, the 3-D nature of the images offers a significant improvement over the 2-D X-ray projections upon which current density assessments are based. Furthermore, microwave illumination is potentially better suited to imaging radiographically dense breasts than X-rays, given the known weaknesses of the latter. Finally, the dependence of the microwave dielectric properties of tissues on changing physiological factors such as thermal 138 conditions, vascularization, water content, or the accumulation of exogenous contrast agents points to the value of microwave imaging in tracking or detecting such changes. 139 LIST OF REFERENCES [1] M. Lazebnik, D. Popovic, L. McCartney, C. B. Watkins, M. J. Lindstrom, J. Harter, S. Sewall, T. Ogilvie, A. Magliocco, T. M. Breslin, W. Temple, D. Mew, J. H. Booske, M. Okoniewski, and S. C. Hagness, "A large-scale study of the ultrawideband microwave dielectric properties of normal, benign, and malignant breast tissues obtained from cancer surgeries," Phys. Med. Biol., vol. 52, pp. 6093-6115, 2007. [2] Cancer Facts & Figures 2009. Atlanta, GA: American Cancer Society, 2009. [3] "Surveillance, epidemiology, and end results http://www.seer.cancer.gov, National Cancer Institute, Research Program, Cancer Statistics Branch, 2009. (SEER) DCCPS, program," Surveillance [4] J. S. Michaelson, M. Silverstein, J. Wyatt, G. Weber, , R. Moore, E. Halpern, D. B. Kopans, and K. Hughes, "Predicting the survival of patients with breast carcinoma using tumor size," Cancer, vol. 95, no. 4, pp. 713-723, 2002. [5] "The world health report 2008," http://www.who.int, World Health Organization, Geneva, Switzerland, 2008. [6] Breast Cancer Facts & Figures 2009. Atlanta, GA: American Cancer Society, 2009. [7] "Breast cancer screening (PDQ)," http://www.cancer.gov/cancertopics/pdq/, tional Cancer Institute. Na- [8] P. M. Meaney, M. W. Fanning, D. Li, S. P. Poplack, and K. D. Paulsen, "A clinical prototype for active microwave imaging of the breast," IEEE Trans. Microw. Theory Tech., vol. 48, no. 11, pp. 1841-1853, Nov. 2000. [9] Cancer Prevention & Early Detection Facts & Figures. Atlanta, GA: American Cancer Society, 2009. [10] W. Teh and A. Wilson, "The role of ultrasound in breast cancer screening. A consensus statement by the European Group for Breast Cancer Screening," European Journal of Cancer. 140 [11] E. Fear and M. Stuchly, "Microwave detection of breast cancer," IEEE Trans. Microw. Theory Tech., vol. 48, no. 11, pp. 1854-1863, Nov 2000. [12] E. Fear, P. Meaney, and M. Stuchly, "Microwaves for breast cancer detection?" Potentials, vol. 22, no. 1, pp. 12-18, Feb/Mar 2003. IEEE [13] E. C. Fear, X. Li, S. C. Hagness, and M. A. Stuchly, "Confocal microwave imaging for breast cancer detection: Localization of tumors in three dimensions," IEEE Trans. Biomed. Eng., vol. 49, no. 8, pp. 812-822, Aug. 2002. [14] X. Li, E. J. Bond, B. D. V. Veen, and S. C. Hagness, "An overview of ultrawideband microwave imaging via space-time beamforming for early-stage breast cancer detection," IEEE Antennas and Propagation Magazine, vol. 47, no. 1, pp. 19-34, Feb. 2005. [15] J. M. Sill and E. C. Fear, "Tissue sensing adaptive radar for breast cancer detection - experimental investigation of simple tumor models," IEEE Trans. Microw. Theory Tech., vol. 53, no. 11, pp. 3312- 3319, Nov. 2005. [16] W. C. Chew, Waves and Fields in Inhomogeneous Media. Press, 1995. Piscataway, NJ: IEEE [17] D. Colton and R. Kress, Inverse Acoustic and Electromagnetic 2nd ed. New York: Springer-Verlag Berlin Heidelberg, 1998. Scattering Theory, [18] M. Lazebnik, L. McCartney, D. Popovic, C. B. Watkins, M. J. Lindstrom, J. Harter, S. Sewall, A. Magliocco, J. H. Booske, M. Okoniewski, and S. C. Hagness, "A largescale study of the ultrawideband microwave dielectric properties of normal breast tissue obtained from reduction surgeries," Phys. Med. Biol., vol. 52, pp. 2637-2656, 2007. [19] E. Zastrow, S. K. Davis, M. Lazebnik, F. Kelcz, B. D. Van Veen, and S. C. Hagness, "Development of anatomically realistic numerical breast phantoms with accurate dielectric properties for modeling microwave interactions with the human breast," IEEE Trans. Biomed. Eng., vol. 55, no. 12, pp. 2792-2800, Dec. 2008. [20] C. J. D'Orsi, L. W. Bassett, and W. A. Berg, Breast Imaging Reporting and Data System: ACR BI-RADS-Mammography, 4th ed. Reston, VA: American College of Radiology, 2003. [21] D. W. Winters, J. D. Shea, P. Kosmas, B. D. V. Veen, and S. C. Hagness, "Threedimensional microwave breast imaging: Dispersive dielectric properties estimation using patient-specific basis functions," IEEE Trans. Med. Imag., vol. 28, no. 7, pp. 969981, July 2009. [22] P. M. Meaney, M. W. Fanning, T. Raynolds, C. J. Fox, Q. Q. Fang, C. A. Kogel, S. P. Poplack, and K. D. Paulsen, "Initial clinical experience with microwave breast imaging in women with normal mammography," Academic Radiology, vol. 14, no. 2, pp. 207-218, Feb. 2007. 141 [23] A. W. Guy, "History of biological effects and medical applications of microwave energy," IEEE Trans. Microw. Theory Tech., vol. 32, no. 9, pp. 1182-1200, Sep. 1984. [24] 0 . P. Gandhi, "Biological effects and medical applications of RF electromagnetic fields," IEEE Trans. Microw. Theory Tech., vol. 30, no. 11, pp. 1831-1847, Nov 1982. [25] H. I. Vargas, W. C. Dooley, R. A. Gardner, K. D. Gonzalez, R. Venegas, S. H. HeywangKobrunner, and A. J. Fenn, "Focused microwave phased array thermotherapy for ablation of early-stage breast cancer: Results of thermal dose escalation," Annals of Surgical Oncology, vol. 11, no. 2, pp. 139-146, Feb. 2004. [26] C. Simon, D. Dupuy, and W. Mayo-Smith, "Microwave ablation: Principles and applications," Radiographics, vol. 25, no. Sp. Iss. SI, pp. S69-S84, Oct. 2005. [27] M. Converse, E. Bond, B. Veen, and C. Hagness, "A computational study of ultrawideband versus narrowband microwave hyperthermia for breast cancer treatment," IEEE Trans. Microw. Theory Tech., vol. 54, no. 5, pp. 2169-2180, May 2006. [28] A. Rosen, M. Stuchly, and A. Vander Vorst, "Applications of RF/microwaves in medicine," IEEE Trans. Microw. Theory Tech., vol. 50, no. 3, pp. 963-974, Mar 2002. [29] L. E. Larsen and J. H. Jacobi, "Microwave interrogation of dielectric targets. Part I: By scattering parameters," Medical Physics, vol. 5, no. 6, pp. 500-508, 1978. [30] J. H. Jacobi and L. E. Larsen, "Microwave interrogation of dielectric targets. Part II: By microwave time delay spectroscopy," Medical Physics, vol. 5, no. 6, pp. 509-513, 1978. [31] L. E. Larsen and J. H. Jacobi, "Microwave scattering parameter imagery of an isolated canine kidney," Medical Physics, vol. 6, no. 5, pp. 394-403, 1978. [32] J. H. Jacobi and L. E. Larsen, "Microwave time delay spectroscopic imagery of isolated canine kidney," Medical Physics, vol. 7, no. 1, pp. 1-7, 1980. [33] M. Slaney, A. Kak, and L. Larsen, "Limitations of imaging with first-order diffraction tomography," IEEE Trans. Microw. Theory Tech., vol. 32, no. 8, pp. 860-874, Aug 1984. [34] L. E. Larsen and J. H. Jacobi, "Microwave imaging with first order diffraction tomography," Medical Applications of Microwave Imaging, edited by L. E. Larsen and J. H. Jacobi, pp. 184-212, 1986. [35] J. H. Jacobi and L. E. Larsen, "Linear FM pulse compression radar techniques applied to biological imaging," Medical Applications of Microwave Imaging, edited by L. E. Larsen and J. H. Jacobi, pp. 138-147, 1986. 142 [36] J. Bolomey, A. Joisel, and G. Gaboriaud, "Applications of microwave active imaging in medicine: past and future," in Application of Microwaves in Medicine, IEE Colloquium on, Feb 1995, pp. 3/1-3/5. [37] T. Guo, W. Guo, and L. Larsen, "A local field study of a water-immersed microwave antenna array for medical imagery and therapy," IEEE Trans. Microw. Theory Tech., vol. 32, no. 8, pp. 844-854, Aug 1984. [38] , "Recent developments in microwave medical imagery - phase and amplitude conjugations and the inverse scattering theorem," Medical Applications of Microwave Imaging, edited by L. E. Larsen and J. H. Jacobi, pp. 167-183, 1986. [39] J. Bolomey, A. Izadnegahdar, L. Jofre, C. Pichot, G. Peronnet, and M. Solaimani, "Microwave diffraction tomography for biomedical applications," IEEE Trans. Microw. Theory Tech., vol. 30, no. 11, pp. 1998-2000, 1982. [40] C. Pichot, L. Jofre, G. Peronnet, and J. C. Bolomey, "Active microwave imaging of inhomogeneous bodies," IEEE Trans. Antennas Propag., vol. 33, no. 4, pp. 416-425, 1985. [41] J. Bolomey, C. Pichot, and G. Garboriaud, "Planar microwave imaging camera for biomedical applications - critical and prospective analysis of reconstruction algorithms," Radio Science, vol. 26, no. 2, pp. 541-549, MAR-APR 1991. [42] D. Ghodgaonkar, O. Gandhi, and M. Hagmann, "Estimation of complex permittivities of three-dimensional inhomogeneous biological bodies," IEEE Trans. Microw. Theory Tech., vol. 31, no. 6, pp. 442-446, Jun. 1983. [43] L. Jofre, M. S. Hawley, A. Broquetas, E. de los Reyes, M. Ferrando, and A. Elias-Fuste, "Medical imaging with a microwave tomographic scanner," IEEE Trans. Biomed. Eng., vol. 37, no. 3, pp. 303-312, March 1990. [44] R. F. Remis and P. M. van den Berg, "On the equivalence of the Newton-Kantorovich and Distorted Born methods," Inverse Problems, vol. 16, no. 1, pp. 1-4, 2000. [45] P. M. van den Berg and R. E. Kleinman, "A contrast source inversion method," Inverse Problems, vol. 13, no. 6, pp. 1607-1620, 1997. [46] M. M. Ney, A. M. Smith, and S. S. Stuchly, "A solution of electromagnetic imaging using pseudoinverse transformation," IEEE Trans. Med. Imag., vol. 3, no. 4, pp. 155162, Dec. 1984. [47] T. C. Guo and W. W. Guo, "Physics of image formation by microwave scattering," in Soc. of Photo-Opt. Instr. Eng. (SPIE) Conf Series, ser. Presented at the Soc. of Photo-Opt. Instr. Eng. (SPIE) Conf., S. J. Dwyer III k R. H. Schneider, Ed., vol. 767, Jan. 1987, p. 30. 143 [48] S. Caorsi, G. Gragnani, and M. Pastorino, "Equivalent current density reconstruction for microwave imaging purposes," IEEE Trans. Microw. Theory Tech., vol. 37, no. 5, pp. 910-916, May 1989. [49] W. Chew and Y. Wang, "Reconstruction of 2-dimensional permittivity distribution using the distorted born iterative method," IEEE Trans. Med. Imag., vol. 9, no. 2, pp. 218-225, JUN 1990. [50] P. Meaney, K. Paulsen, and T. Ryan, "Two-dimensional hybrid element image reconstruction for TM illumination," IEEE Trans. Antennas Propag., vol. 43, no. 3, pp. 239-247, Mar 1995. [51] A. Devaney and M. Dennison, "Inverse scattering in inhomogeneous background media," Inverse Problems, vol. 19, no. 4, pp. 855-870, AUG 2003. [52] N. Joachimowicz, C. Pichot, and J. P. Hugonin, "Inverse scattering: an iterative numerical method for electromagnetic imaging," IEEE Trans. Antennas Propag., vol. 39, no. 12, pp. 1742-1753, Dec. 1991. [53] P. Meaney, K. Paulsen, A. Hartov, and R. Crane, "Microwave imaging for tissue assessment: initial evaluation in multitarget tissue-equivalent phantoms," IEEE Trans. Biomed. Eng., vol. 43, no. 9, pp. 878-890, Sep. 1996. [54] S. Y. Semenov, R. H. Svenson, A. E. Bulyshev, A. E. Souvorov, A. E. Souvorov, V. Y. Borisov, Y. Sizov, A. N. Starostin, K. R. Dezern, G. P. Tatsis, and V. Y. Baranov, "Microwave tomography: Two-dimensional system for biological imaging," IEEE Trans. Biomed. Eng., vol. 43, no. 9, pp. 869-877, Sep. 1996. [55] A. Souvorov, A. Bulyshev, S. Semenov, R. Svenson, A. Nazarov, Y. Sizov, and G. Tatsis, "Microwave tomography: a two-dimensional newton iterative scheme," IEEE Trans. Microw. Theory Tech., vol. 46, no. 11, pp. 1654-1659, Nov 1998. [56] S. Y. Semenov, R. H. Svenson, A. E. Bulyshev, A. E. Souvorov, A. G. Nazarov, Y. E. Sizov, A. V. Pavlovsky, V. Y. Borisov, B. A. Voinov, G. I. Simonova, A. N. Starostin, V. G. Posukh, G. P. Tatsis, and V. Y. Baranov, "Three-dimensional microwave tomography: experimental prototype of the system and vector Born reconstruction method," IEEE Trans. Biomed. Eng., vol. 46, no. 8, pp. 937-946, Aug. 1999. [57] A. E. Bulyshev, A. E. Souvorov, S. Y. Semenov, R. H. Svenson, A. G. Nazarov, Y. E. Sizov, and G. P. Tatsis, "Three-dimensional microwave tomography, theory and computer experiments in scalar approximation," Inverse Problems, vol. 16, p. 863875, 2000. 144 [58] S. Semenov, A. Bulyshev, A. Abubakar, V. Posukh, Y. Sizov, A. Souvorov, P. van den Berg, and T. Williams, "Microwave-tomographic imaging of the high dielectric-contrast objects using different image-reconstruction approaches," IEEE Trans. Microw. Theory Tech., vol. 53, no. 7, pp. 2284-2294, July 2005. [59] A. Abubakar, P. M. van den Berg, and J. J. Mallorqui, "Imaging of biomedical data using a multiplicative regularized contrast source inversion method," IEEE Trans. Microw. Theory Tech., vol. 50, no. 7, pp. 1761-1771, july 2002. [60] A. Abubakar, W. Hu, P. M. van den Berg, and T. M. Habashy, "A finite-difference contrast source inversion method," Inverse Problems, vol. 24, no. 6, p. 065004 (17pp), 2008. [61] E. C. Fear, S. C. Hagness, P. M. Meaney, M. Okoniewski, and M. A. Stuchly, "Enhancing breast tumor detection with near-field imaging," IEEE Microwave Magazine, pp. 48-56, Mar. 2002. [62] A. H. Barrett and P. C. Myers, "Subcutaneous temperatures: A method of noninvasive sensing," Science, vol. 190, no. 4215, pp. 669-671, 1975. [63] M. Lazebnik, M. C. Converse, J. H. Booske, and S. C. Hagness, "Ultrawideband temperature-dependent dielectric properties of animal liver tissue in the microwave frequency range," Phys. Med. Biol, vol 51, pp. 1941-1955, 2006. [64] J. LIN, "Frequency optimization for microwave imaging of biological tissues," Proc. IEEE, vol. 73, no. 2, pp. 374-375, 1985. [65] R. Kruger, K. Miller, H. Reynolds, W. Kiser, D. Reinecke, and G. Kruger, "Breast cancer in vivo: Contrast enhancement with thermoacoustic CT at 434 MHz - Feasibility study," Radiology, vol. 216, no. 1, pp. 279-283, Jul. 2000. [66] L. Nie, D. Xing, Q. Zhou, D. Yang, and H. Guo, "Microwave-induced thermoacoustic scanning CT for high-contrast and noninvasive breast cancer imaging," Medical Physics, vol. 35, no. 9, pp. 4026-4032, Sep. 2008. [67] A. Mashal, J. H. Booske, and S. C. Hagness, "Towards contrast-enhanced microwaveinduced thermoacoustic imaging of breast cancer: An experimental study of the effects of microbubbles on simple thermoacoustic targets," Phys. Med. Biol, vol. 54, pp. 641650, 2009. [68] P. Kosmas and C. Rappaport, "A matched filter FDTD-based time reversal algorithm for microwave breast cancer detection," IEEE Trans. Antennas Propag., vol. 54, no. 4, pp. 1257-1264, Apr. 2006. 145 [69] S. Hagness, A. Taflove, and J. Bridges, "Three-dimensional FDTD analysis of a pulsed microwave confocal system for breast cancer detection: design of an antenna-array element," IEEE Trans. Antennas Propag., vol. 47, no. 5, pp. 783-791, May 1999. [70] S. Davis, E. Bond, X. Li, S. Hagness, and B. Veen, "Microwave imaging via space-time beamforming for early detection of breast cancer: Beamformer design in the frequency domain," Journal of Electromagnetic Waves and Applications, vol. 17, no. 2, p. 357, 2003. [71] E. J. Bond, S. C. Hagness, and B. D. Van Veen, "Microwave imaging via space-time beamforming for early detection of breast cancer," IEEE Trans. Antennas Propag., vol. 51, no. 8, pp. 1690-1705, 2003. [72] M. Klemm, I. J. Craddock, J. A. Leendertz, A. Preece, and R. Benjamin, "RadarBased Breast Cancer Detection Using a Hemispherical Antenna Array-Experimental Results," IEEE Trans. Antennas Propag., vol. 57, no. 6, pp. 1692-1704, Jun. 2009. [73] Y. Xie, B. Guo, L. Xu, J. Li, and P. Stoica, "Multistatic adaptive microwave imaging for early breast cancer detection," IEEE Trans. Biomed. Eng., vol. 53, no. 8, pp. 1647-1657, AUG 2006. [74] A. E. Bulyshev, S. Y. Semenov, A. E. Souvorov, R. H. Svenson, A. G. Nazarov, Y. E. Sizov, and G. P. Tatsis, "Computational modeling of three-dimensional microwave tomography of breast cancer," IEEE Trans. Biomed. Eng., vol. 48, no. 9, pp. 10531056, Sep. 2001. [75] Q. Fang, P. M. Meaney, S. D. Geimer, A. V. Streltsov, and K. D. Paulsen, "Microwave image reconstruction from 3-D field coupled to 2-D parameter estimation," IEEE Trans. Med. Imag., vol. 23, no. 4, pp. 475-484, Apr. 2004. [76] Q. Fang, P. M. Meaney, and K. D. Paulsen, "Microwave image reconstruction of tissue property dispersion characteristics utilizing multiple-frequency information," IEEE Trans. Microw. Theory Tech., vol. 52, no. 8, pp. 1866-1875, 2004. [77] Z. Q. Zhang and Q. H. Liu, "Three-dimensional nonlinear image reconstruction for microwave biomedical imaging," IEEE Trans. Biomed. Eng., vol. 51, no. 3, pp. 544548, 2004. [78] T. Rubaek, P. M. Meaney, P. Meincke, and K. D. Paulsen, "Nonlinear microwave imaging for breast-cancer screening using Gauss-Newton's method and the CGLS inversion algorithm," IEEE Trans. Antennas Propag., vol. 55, no. 8, pp. 2320-2331, Aug. 2007. [79] J. E. Johnson, T. Takenaka, and T. Tanaka, "Two-dimensional time-domain inverse scattering for quantitative analysis of breast composition," IEEE Trans. Biomed. Eng., vol. 55, no. 8, pp. 1941-1945, Aug. 2008. 146 [80] C. Gilmore, A. Abubakar, W. Hu, T. Habashy, and P. van den Berg, "Microwave biomedical data inversion using the finite-difference contrast source inversion method," IEEE Trans. Antennas Propag., vol. 57, no. 5, pp. 1528-1538, May 2009. [81] T. Rubaek, O. Kim, and P. Meincke, "Computational validation of a 3-D microwave imaging system for breast-cancer screening," IEEE Trans. Antennas Propag., vol. 57, no. 7, pp. 2105-2115, July 2009. [82] N. Irishina, D. Alvarez, O. Dorn, and M. Moscoso, "Structural level set inversion for microwave breast screening," Inverse Problems, vol. 26, no. 3, Mar 2010. [83] J. D. Shea, P. Kosmas, B. D. Van Veen, and S. C. Hagness, "Contrast-enhanced microwave imaging of breast tumors: A computational study using 3D realistic numerical phantoms," Inverse Problems, vol. 26, no. 7, Sp. Iss. SI, July 2010. [84] J. D. Shea, P. Kosmas, S. C. Hagness, and B. D. Van Veen, "Three-dimensional microwave imaging of realistic numerical breast phantoms via a multiple-frequency inverse scattering technique," Medical Physics, vol. 37, no. 8, pp. 4210-4226, Aug 2010. [85] P. Meaney, K. Paulsen, and J. Chang, "Near-field microwave imaging of biologicallybased materials using a monopole transceiver system," IEEE Trans. Microw. Theory Tech., vol. 46, no. 1, pp. 31-45, Jan 1998. [86] P. Meaney, K. Paulsen, B. Pogue, and M. Miga, "Microwave image reconstruction utilizing log-magnitude and unwrapped phase to improve high-contrast object recovery," IEEE Trans. Med. Imag., vol. 20, no. 2, pp. 104-116, Feb 2001. [87] P. M. Meaney, E. Demidenko, N. K. Yagnamurthy, D. Li, M. W. Fanning, and K. D. Paulsen, "A two-stage microwave image reconstruction procedure for improved internal feature extraction," Medical Physics, vol. 28, no. 11, pp. 2358-2369, 2001. [88] P. M. Meaney, N. K. Yagnamurthy, and K. D. Paulsen, "Pre-scaling of reconstruction parameter components to reduce imbalance in image recovery process," Phys. Med. Biol., vol. 47, pp. 1101-1119, 2002. [89] S. P. Poplack, T. D. Tosteson, W. A. Wells, B. W. Pogue, P. M. Meaney, A. Hartov, C. A. Kogel, S. K. Soho, J. J. Gibson, and K. D. Paulsen, "Electromagnetic breast imaging: Results of a pilot study in women with abnormal mammograms," Radiology, vol. 243, no. 2, pp. 350-359, May 2007. [90] Q. Fang, P. Meaney, and K. Paulsen, "Viable three-dimensional medical microwave tomography: Theory and numerical experiments," Antennas and Propagation, IEEE Transactions on, vol. 58, no. 2, pp. 449-458, Feb. 2010. 147 [91] C. Gilmore, P. Mojabi, A. Zakaria, M. Ostadrahimi, C. Kaye, S. Noghanian, L. Shafai, S. Pistorius, and J. LoVetri, "A wideband microwave tomography system with a novel frequency selection procedure," Biomedical Engineering, IEEE Transactions on, vol. 57, no. 4, pp. 894-904, april 2010. [92] C. Gilmore, P. Mojabi, A. Zakaria, S. Pistorius, and J. LoVetri, "On super-resolution with an experimental microwave tomography system," Antennas and Wireless Propagation Letters, IEEE, vol. 9, pp. 393 -396, 2010. [93] M. Ali and M. Moghaddam, "3D nonlinear super-resolution microwave inversion technique using time-domain data," Antennas and Propagation, IEEE Transactions on, vol. 58, no. 7, pp. 2327-2336, july 2010. [94] C. A. Balanis, Advanced Engineering Electromagnetics. and Sons, 1989. Hoboken, NJ: John Wiley [95] J. Nocedal and S. J. Wright, Nonlinear Optimization, 2nd ed. New York: Springer, 2006. [96] P. C. Hansen, Rank-Deficient and Discrete Ill-Posed Problems: Numerical Aspects of Linear Inversion. Philadelphia, PA: SIAM, 1998. [97] Y. Saad, Iterative Methods for Sparse Linear Systems, 2nd ed. SIAM, 2003. Philadelphia, PA: [98] A. N. Tikhonov and V. T. Arsenin, Solutions of Ill-Posed Problems. D.C.: Winston, 1977. Washington, [99] R. Acar and C. Vogel, "Analysis of bounded variation penalty methods for ill-posed problems," Invers Problems, vol. 10, no. 6, pp. 1217-1229, Dec 1994. [100] P. Lobel, C. Picbota, L. Blanc Feraud, and M. Barlaud, "Conjugate gradient algorithm with edge-preserving regularization for image reconstruction from experimental data," in Antennas and Propagation Society International Symposium, 1996. AP-S. Digest, vol. 1, Jul 1996, pp. 644-647 vol.1. [101] P. van den Berg, A. Abubakar, and J. Fokkema, "Multiplicative regularization for contrast profile inversion," Radio Science, vol. 38, no. 2, APR 22 2003. [102] E. Resmerita, "Regularization of ill-posed problems in Banach spaces: convergence rates," Inverse Problems, vol. 21, no. 4, pp. 1303-1314, Aug 2005. [103] C. Gilmore, P. Mojabi, and J. LoVetri, "Comparison of an enhanced distorted born iterative method and the multiplicative-regularized contrast source inversion method," IEEE Trans. Antennas Propag., vol. 57, no. 8, pp. 2341-2351, Aug. 2009. 148 [104] J. A. Harvey and V. E. Bovbjerg, "Quantitative assessment of mammographic breast density: Relationship with breast cancer risk," Radiology, vol. 230, pp. 29-41, 2004. [105] K. Kerlikowske, L. Ichikawa, D. L. Miglioretti, D. S. M. Buist, R M. Vacek, R. SmithBindman, B. Yankaskas, P. A. Carney, and R. Ballard-Barbash, "Longitudinal measurement of clinical mammographic breast density to improve estimation of breast cancer risk," Journal of the National Cancer Institute, vol. 99, no. 5, pp. 386-395, Mar 7 2007. [106] P. C. Stomper, D. J. DSouza, P. A. DiNitto, and M. A. Arredondo, "Analysis of parenchymal density on mammograms in 1353 women 25-79 years old," American Journal of Roentgenology, vol. 167, no. 5, pp. 1261-1265, Nov. 1996. [107] N. F. Boyd, H. Guo, L. J. Martin, L. Sun, J. Stone, E. Fishell, R. A. Jong, G. Hislop, A. Chiarelli, S. Minkin, and M. J. Yaffe, "Mammographic density and the risk and detection of breast cancer," New England Journal of Medicine, vol. 356, no. 3, pp. 227-236, Jan. 2007. [108] G. A. Colditz, W. A. Willett, D. J. Hunter, M. J. Stampfer, J. E. Manson, C. H. Hennekens, B. A. Rosner, and F. E. Speizer, "Family history, age, and risk of breastcancer - prospective data from the nurses health study (vol 270, pg 338, 1993)," Journal of the American Medical Association, vol. 270, no. 13, p. 1548, Oct. 1993. [109] K. E. Martin, M. A. Helvie, C. Zhou, M. A. Roubidoux, J. E. Bailey, C. Paramagul, C. E. Blane, K. A. Klein, S. S. Sonnad, and H. Chan, "Mammographic density measured with quantitative computer-aided method: Comparison with radiologists' estimates and BI-RADS categories," Radiology, vol. 240, no. 3, pp. 656-665, Sep. 2006. [110] S. van Engeland, P. R. Snoeren, H. Huisman, C. Boetes, and N. Karssemeijer, "Volumetric breast density estimation from full-field digital mammograms," IEEE Trans. Med. Imag., vol. 25, no. 3, pp. 273-282, Mar. 2006. [111] C. K. Glide-Hurst, N. Duric, and P. Littrup, "A new method for quantitative analysis of mammographic density," Medical Physics, vol. 34, no. 11, pp. 4491-4498, Nov. 2007. [112] K. Nie, J. Chen, S. Chan, M. I. Chau, H. J. Yu, S. Bahri, T. Tseng, O. Nalcioglu, and M. Su, "Development of a quantitative method for analysis of breast density based on three-dimensional breast MRI," Medical Physics, vol. 35, no. 12, pp. 5253-5262, Dec. 2008. [113] W. A. Berg, J. D. Blume, A. M. Adams, R. A. Jong, R. G. Barr, D. E. Lehrer, E. D. Pisano, W. P. Evans, M. C. Mahoney, L. Hovanessian Larsen, G. J. Gabrielli, and E. B. Mendelson, "Reasons Women at Elevated Risk of Breast Cancer Refuse Breast MR Imaging Screening: ACRIN 66661," Radiology, vol. 254, no. 1, pp. 79-87, 2010. 149 [114] X. Li, S. K. Davis, S. C. Hagness, D. van der Weide, and B. D. V. Veen, "Microwave imaging via space-time beamforming: Experimental investigation of tumor detection in multi-layer breast phantoms," IEEE Trans. Microw. Theory Tech., vol. 52, no. 8, pp. 1856-1865, Aug. 2004. [115] M. Klemm, I. Craddock, J. Leendertz, A. Preece, and R. Benjamin, "Experimental and clinical results of breast cancer detection using UWB microwave radar," IEEE International Symposium on Antennas and Propagation, July 2008. [116] M. Lazebnik, E. L. Madsen, G. R. Prank, and S. C. Hagness, "Tissue-mimicking phantom materials for narrowband and ultrawideband microwave applications," Phys. Med. Biol., vol. 50, pp. 4245-4258, 2005. [117] D. Li, P. M. Meaney, and K. D. Paulsen, "Conformal microwave imaging for breast cancer detection," IEEE Trans. Microw. Theory Tech., vol. 51, no. 4, pp. 1179-1186, Apr 2003. [118] P. Mojabi and J. LoVetri, "Microwave biomedical imaging using the multiplicative regularized Gauss-Newton inversion," Antennas and Wireless Propagation Letters, IEEE, vol. 8, pp. 645-648, 2009. [119] http://uwcem.ece.wisc.edu, University of Wisconsin Computational Electromagnetics Laboratory, UWCEM Numerical Breast Phantom Repository. [120] S. Gabriel, R. W. Lau, and C. Gabriel, "The dielectric properties of biological tissues: III. Parametric models for the dielectric spectrum of tissues," Phys. Med. Biol, vol. 41, pp. 2271-2293, 1996. [121] A. Taflove and S. C. Hagness, Computational Electrodynamics: The Finite-Difference Time-Domain Method, 3rd ed. Norwood, MA: Artech House, 2005. [122] Z. Q. Zhang, Q. H. Liu, C. Xiao, E. Ward, G. Ybarra, and W. T. Joines, "Microwave breast imaging: 3-D forward scattering simulation," IEEE Trans. Biomed. Eng., vol. 50, no. 10, pp. 1180-1189, Oct. 2003. [123] A. Bulyshev, A. Souvorov, S. Semenov, V. Posukh, and Y. Sizov, "Three-dimensional vector microwave tomography: theory and computational experiments," Inverse Problems, vol. 20, no. 4, pp. 1239-1259, Aug 2004. [124] T. J. Cui and W. C. Chew, "Inverse scattering of two-dimensional dielectric objects buried in a lossy earth using the distorted Born iterative method," IEEE Trans. Geosci. Remote Sens., vol. 39, no. 2, pp. 339-346, 2001. [125] T. J. Cui, Y. Qin, G. Wang, and W. C. Chew, "Low-frequency detection of twodimensional buried objects using high-order extended born approximations," Inverse Problems, vol. 20, p. S41S62, 2004. 150 [126] M. Lazebnik, M. Okoniewski, J. H. Booske, and S. C. Hagness, "Highly accurate debye models for normal and malignant breast tissue dielectric properties at microwave frequencies," Microwave and Wireless Components Letters, IEEE, vol. 17, no. 12, pp. 822-824, Dec 2007. [127] A. Abubakar, T. M. Habashy, V. L. Druskin, L. Knizhnerman, and D. Alumbaugh, "2.5D forward and inverse modeling for interpreting low-frequency electromagnetic measurements," Geophysics, vol. 73, no. 4, pp. F165-F177, Jul-Aug 2008. [128] D. Calvetti, G. Landi, L. Reichel, and F. Sgallari, "Non-negativity and iterative methods for ill-posed problems," Inverse Problems, vol. 20, pp. 1747-1758, 2004. [129] D. W. Winters, E. J. Bond, B. D. V. Veen, and S. C. Hagness, "Estimation of the frequency-dependent average dielectric properties of breast tissue using a time-domain inverse scattering technique," IEEE Trans. Antennas Propag., vol. 54, no. 11 (2), pp. 3517-3528, Nov. 2006. [130] T. Rubaek and V. Zhurbenko, "Prototype of microwave imaging system for breastcancer screening," in 13th International Symposium on Antenna Technology and Applied Electromagnetics and the Canadian Radio Sciences Meeting, Banff, AB, Feb. 2009. [131] D. W. Winters, J. D. Shea, E. L. Madsen, G. R. Frank, B. D. Van Veen, and S. C. Hagness, "Estimating the breast surface using UWB microwave monostatic backscatter measurements," IEEE Trans. Biomed. Eng., vol. 55, no. 1, pp. 247-256, Jan. 2008. [132] T. Williams, J. Sill, and E. Fear, "Breast surface estimation for radar-based breast imaging systems," IEEE Trans. Biomed. Eng., vol. 55, no. 6, pp. 1678-1686, June 2008. [133] D. Popovic, L. McCartney, C. Beasley, M. Lazebnik, M. Okoniewski, S. C. Hagness, and J. H. Booske, "Precision open-ended coaxial probes for in vivo and ex vivo dielectric spectroscopy of biological tissues at microwave frequencies," IEEE Trans. Microw. Theory Tech., vol. 53, no. 5, pp. 1713-1722, May 2005. [134] J. Shea, S. Hagness, and B. Van Veen, "Hardware acceleration of FDTD computations for 3-D microwave breast tomography," in Proceedings of the IEEE International Symposium on Antennas and Propagation, Charleston, SC, june 2009. [135] J. D. Shea, P. Kosmas, S. C. Hagness, and B. D. Van Veen, "Contrast-enhanced microwave breast imaging," in 13th International Symposium on Antenna Technology and Applied Electromagnetics and the Canadian Radio Sciences Meeting, Banff, AB, Feb. 2009. [136] F. Chen and W. C. Chew, "Experimental verification of super resolution in nonlinear inverse scattering," Applied Physics Letters, vol. 72, no. 23, pp. 3080-3082, 1998. 151 [137] T. Cui, W. Chew, X. Yin, and W. Hong, "Study of resolution and super resolution in electromagnetic imaging for half-space problems," IEEE Trans. Antennas Propag., vol. 52, no. 6, pp. 1398-1411, jun 2004. [138] C. K. Kuhl, "MRI of breast tumors," Eur. Radiol., vol. 10, no. 1, pp. 46-58, Jan. 2000. [139] H. Maeda, "The enhanced permeability and retention (EPR) effect in tumor vasculature: The key role of tumor-selective macromolecular drug targeting," in Advances in Enzyme Regulation, Vol J^l, ser. Advances in Enzyme Regulation, G. Weber, Ed., 2001, vol. 41, pp. 189-207. [140] Z. Liu, W. Cai, L. He, N. Nakayama, K. Chen, X. Sun, X. Chen, and H. Dai, "In vivo biodistribution and highly efficient tumour targeting of carbon nanotubes in mice," Nature Nanotechnology, vol. 2, no. 1, pp. 47-52, Jan. 2007. [141] A. Mashal, B. Sitharaman, J. H. Booske, and S. C. Hagness, "Dielectric characterization of carbon nanotube contrast agents for microwave breast cancer detection," in AP-S International Symposium, Charlotte, NC, Jun. 2009. [142] A. Mashal, B. Sitharaman, X. Li, P. Avti, A. V. Sahakian, J. H. Booske, and S. C. Hagness, "Toward carbon-nanotube-based theranostic agents for microwave detection and treatment of breast cancer: Enhanced dielectric and heating response of tissuemimicking materials," IEEE Trans. Biomed. Eng., vol. (in press), 2010. [143] A. Ishimaru, Electromagnetic Wave Propagation, Radiation, and Scattering. wood Cliffs, NJ: Prentice Hall, 1991. Engle- [144] J. S. Michaelson, S. Satija, D. Kopans, R. Moore, M. Silverstein, A. Comegno, K. Hughes, A. Taghian, and B. S. S. Powell, "Gauging the impact of breast carcinoma screening in terms of tumor size and death rate," Cancer, vol. 98, no. 10, pp. 2114-2124, 2003. [145] M. Lazebnik, J. H. Booske, and S. C. Hagness, "Dielectric-properties contrast enhancement for microwave breast cancer detection: Numerical investigations of microbubble contrast agents," in URSI General Assembly, Chicago, IL, Aug. 2008. [146] Q. Fang, P. M. Meaney, and K. D. Paulsen, "Singular value analysis of the Jacobian matrix in microwave image reconstruction," IEEE Trans. Antennas Propag., vol. 54, no. 8, pp. 2371-2380, Aug. 2006. [147] P. M. van den Berg and R. E. Kleinman, "A total variation enhanced modified gradient algorithm for profile reconstruction," Inverse Problems, vol. 11, no. 3, pp. LS-L10, 1995. 152 [148] P. Mojabi and J. LoVetri, "Overview and classification of some regularization techniques for the gauss-newton inversion method applied to inverse scattering problems," Antennas and Propagation, IEEE Transactions on, vol. 57, no. 9, pp. 2658-2665, Sept. 2009. [149] , "Enhancement of the Krylov subspace regularization for microwave biomedical imaging," Medical Imaging, IEEE Transactions on, vol. 28, no. 12, pp. 2015-2019, Dec. 2009. [150] O. Bucci and G. Franceschetti, "On the spatial bandwidth of scattered fields," IEEE Transactions on Antennas and Propagation, vol. 35, no. 12, pp. 1445-1455, Dec 1987. [151] , "On the degrees of freedom of scattered fields," IEEE Transactions on Antennas and Propagation, vol. 37, no. 7, pp. 918-929, July 1989. [152] O. Bucci and T. Isernia, "Electromagnetic inverse scattering: Retrievable information and measurement strategies," Radio Science, vol. 32, no. 6, pp. 2123-2137, Nov-Dec 1997. [153] M. Hori, "Inverse analysis method using spectral decomposition of Green's function," Geophysical Journal International, vol. 147, no. 1, pp. 77-87, Oct 2001. [154] E. Zastrow, S. C. Hagness, and B. D. V. Veen, "3D computational study of noninvasive patient-specific microwave hyperthermia treatment of breast cancer," Physics in Medicine and Biology, vol. 55, no. 13, p. 3611, 2010.

1/--страниц