ALGORITHMS AND CALIBRATION TECHNIQUES FOR HIGH RESOLUTION MICROWAVE BREAST IMAGING by Timothy J. Colgan A dissertation submitted in partial fulfillment of the requirements for the degree of Doctor of Philosophy (Electrical and Computer Engineering) at the UNIVERSITY OF WISCONSIN–MADISON 2015 Date of final oral examination: 4/10/15 The dissertation is approved by the following members of the Final Oral Committee: Barry D. Van Veen, Professor, Electrical and Computer Engineering Susan C. Hagness, Professor, Electrical and Computer Engineering Nader Behdad, Associate Professor, Electrical and Computer Engineering Stephen J. Wright, Professor, Computer Sciences Walter F. Block, Professor, Medical Physics ProQuest Number: 10128906 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. ProQuest 10128906 Published by ProQuest LLC (2016). Copyright of the Dissertation is held by the Author. All rights reserved. This work is protected against unauthorized copying under Title 17, United States Code Microform Edition © ProQuest LLC. ProQuest LLC. 789 East Eisenhower Parkway P.O. Box 1346 Ann Arbor, MI 48106 - 1346 i acknowledgments I want to thank my advisors, Drs. Barry Van Veen and Susan Hangess, for their support, patience, and guidance throughout my graduate studies. I would also like to thank my research committee members, Drs. Nader Behdad, Steven Wright, and Wally Block, for their support with my research. Thanks to all the current and past graduate students that have helped make my time here both productive and enjoyable, especially Luz María Neira, Owen Mays, James Sawicki, Zach Hertzberg, Prathyusha Sharma, Dr. Jacob Shea, Dr. Suzette Aguilar, Dr. Steve Kennedy, and Dr. Matthew Burfeindt. I also want to thank the faculty at the University of Wyoming, especially Drs. Jerry Hamann, Robert Kubichek, John Pierre, and Cameron H. G. Wright; without the exceptional education and guidance you provided me I would not have made it this far. To all my friends and family, your endless support and overwhelming positive influence in my life has given me a lifelong passion for learning and the determination to chase my dreams. ii contents Contents ii List of Tables iv List of Figures v Abstract x 1 Introduction 1 2 Background 4 2.1 Introduction 4 2.2 Microwave Breast Imaging 4 2.3 Level Set Imaging 8 2.4 Model Error in Microwave Imaging 12 3 A Level Set Method for 3D Microwave Breast Imaging 14 3.1 Introduction 14 3.2 Dielectric Properties of Breast Tissue 15 3.3 Level Set Method for Microwave Imaging 18 3.4 Algorithm Implementation Details 23 3.5 Results 26 3.6 Discussion 31 3.7 Conclusion 33 4 Model Error in Model-Based Inverse Scattering Algorithms for Microwave Breast Imaging 34 4.1 Introduction 34 4.2 Model Error Metric 35 4.3 Experimental and Numerical Testbeds 40 iii 4.4 4.5 4.6 Results 44 Discussion 52 Conclusion 53 5 Calibration and Imaging with Model Error 55 5.1 Introduction 55 5.2 Simple Phantom Calibration 55 5.3 Multiple Object Calibration 59 5.4 Model Error Reduction using Resonance Matching 63 5.5 Oil and Substrate Imaging 63 5.6 TMM Phantom Imaging 65 5.7 Improvements to reduce the effect of model error 69 6 Future Work and Conclusion 74 6.1 Future Work 74 6.2 Conclusion 76 Bibliography 78 iv list of tables 3.1 Average absolute error in Debye parameters with the two reconstruction techniques . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Volumetric percent density estimates for the nine phantoms and the two reconstruction techniques . . . . . . . . . . . . . . . . . . . . . . . . . . 29 4.1 4.2 Dielectric Properties at 2 GHz . . . . . . . . . . . . . . . . . . . . . . . . SER improvement with calibration . . . . . . . . . . . . . . . . . . . . . 43 49 5.1 Two pole Debye parameters for the immersion medium and imaging object materials. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . SER (4.5) values for calibration of the rectangular prism phantom. . . . Dielectric properties at 2 GHz for calibration materials . . . . . . . . . . SER (4.5) values for a multiple object calibration of the 3D printed breast shell filled with TMMs. . . . . . . . . . . . . . . . . . . . . . . . . . . . . Array design for reducing model error due to 10% error in immersion medium. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Reconstruction quality vs imaging domain size with 5% error in immersion medium. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Reconstruction quality vs imaging domain size with antenna position errors. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.2 5.2 5.3 5.4 5.5 5.6 5.7 27 56 57 61 62 70 71 72 v list of figures 2.1 2.2 2.3 2.4 2.5 3.1 3.2 Measured dielectric properties for healthy breast tissue [1] and malignant breast tissue [2]. Solid red lines denote the 50th percentile of malignant tissue, the dash-dot lines represent the 50th percentile of fibroglandular tissue, and the dashed line is the 50th percentile of adipose tissue. Permittivity is plotted in (a) and conductivity is in (b). . . . . . . 6 Sequential coronal cross-sections, spaced every 8 mm, of the static permittivity for a phantom from the UWCEM repository [3] are shown in the first row, and the second row shows sequential cross-sections from a DBIM reconstruction of the static permittivity [4]. Both rows are plotted on the same color scale shown at the bottom of the figure. . . . . . . . . 7 2D example of a single level set φ segmenting a sample domain into distinct regions. Sections with (φ > 0) will have dielectric properties that correspond to Region A, and sections with (φ 6 0) will have dielectric properties that correspond to Region B. . . . . . . . . . . . . . . . . . . 9 Result from [5] where the far left is the phantom being simulated with three distinct tumors, the middle is a pixel by pixel reconstruction that serves as the initial guess, and on the right is the final reconstructed image. 11 Photograph of the UWCEM experimental imaging setup. . . . . . . . . 13 Linear Debye-parameter relationship of (3.2), shown by the black curve, compared to the Debye parameters that correspond to the 25th, 50th, and 75th percentile values [4] of measured dielectric properties [1]. (a) Infinite permittivity and static conductivity versus static permittivity of adipose tissue; (b) Infinite permittivity and static conductivity versus static permittivity of fibroglandular tissue. . . . . . . . . . . . . . . . . . 2D example of a single level set φ segmenting a sample domain into distinct regions. Sections with (φ > 0) have dielectric properties that correspond to fibroglandular tissue, and regions with (φ 6 0) have dielectric properties that correspond to adipose tissue. . . . . . . . . . . 17 19 vi 3.3 (a) Coronal cross section of static permittivity from a DBIM reconstruction. (b) Thresholded version of the DBIM reconstruction where fibroglandular tissue is represented by s = 49.1 (i.e. red) and adipose tissue is represented by s = 4.8 (i.e. dark blue). (c) Initial guess of the level set obtained using the DBIM reconstruction. Voxels in the level set are assigned a value based on the distance from the boundary between tissue types, with negative distance used for the region corresponding to adipose tissue. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.4 Coronal cross sections of static permittivity from a Class 3 phantom (top) and the corresponding cross sections from the 3D DBIM reconstruction (middle) and the 3D level set method reconstruction (bottom). The coronal cross sections are taken every 8 mm between the top and bottom rings of antennas with the cross sections closest to the phantom base shown at the far left. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.5 Coronal cross sections of static permittivity from a Class 4 phantom (top) and the corresponding cross sections from the 3D DBIM reconstruction (middle) and the 3D level set method reconstruction (bottom). The coronal cross sections are taken every 8 mm between the top and bottom rings of antennas with the cross sections closest to the phantom base shown at the far left. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.6 Volume of fibroglandular density in 2 mm thick coronal cross sections as a function of cross-section position (distance from the base of the phantom). The markers represent the mean value in the reconstructed cross section. Error bars represent the standard deviation of the estimates over 10 noise realizations. (a) Class 1 phantom, (b) Class 2 phantom, (c) Class 3 phantom, and (d) Class 4 phantom. . . . . . . . . . . . . . . . . 4.1 (a) Photograph of the experimental patch array. (b) Simulated model of the patch array. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 25 27 28 31 41 vii 4.2 4.3 4.4 4.5 4.6 Side view (a) and top view (b) of the 3D printed breast surface used in calibration and imaging results. (c) Picture of 3D printed breast surface filled with TMMs; the two cylinders mimic fibroglandular tissue and the remaining material mimics adipose tissue. . . . . . . . . . . . . . . Sagittal, axial, and coronal cross sections of the relative permittivity (left) and effective conductivity (right) of a simulated version of the 3D printed breast surface filled with TMMs. . . . . . . . . . . . . . . . . . . (a) DBIM imaging performance, assessed with the cosine metric, versus SER for different types of model error with a Class 3 phantom. (b) Extended domain DBIM imaging performance versus SER for different types of model error. The baseline performance of both algorithms is indicated by a dashed line and markers denote 10% error in the dielectric properties of the material. . . . . . . . . . . . . . . . . . . . . . . . . . . (a) DBIM imaging performance, assessed with the cosine metric, versus SER for different types of model error with a Class 2 phantom. (b) Extended domain DBIM imaging performance versus SER for different types of model error. The baseline performance of both algorithms is indicated by a dashed line and markers denote 10% error in the dielectric properties of the material. . . . . . . . . . . . . . . . . . . . . . . . . . . (a) Sagittal, axial, and coronal cross sections of the relative permittivity (left) and effective conductivity (right) for a reconstruction of the TMM filled breast shell (Fig. 4.3) using simulated data without model error. (b) Sagittal, axial, and coronal cross sections of a reconstruction of the TMM filled breast shell (Fig. 4.2c) from experimental measurements with 5 dB SER. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Side view (a) and top view (b) of the experimental 3D printed plastic rectangular prism. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5.2 Sagittal, axial, and coronal cross sections of the relative permittivity (left) and effective conductivity (right) of a simulated model of the 3D printed plastic rectangular prism filled with a 90% glycerol 10% water mixture. 42 43 46 47 51 5.1 56 57 viii 5.3 5.4 5.5 5.6 5.7 5.8 5.9 5.10 5.11 Sagittal, axial, and coronal cross sections of the relative permittivity (left) and effective conductivity (right) of a simulated reconstruction of the rectangular prism filled with 90% glycerol 10% water. . . . . . . . . . . Sagittal, axial, and coronal cross sections of the relative permittivity (left) and effective conductivity (right) of an experimental reconstruction of the rectangular prism filled with 90% glycerol 10% water that was calibrated with a 80% glycerol 20% water filled rectangular prism. . . . Sagittal, axial, and coronal cross sections of the relative permittivity (left) and effective conductivity (right) of the initial guess for a simulated reconstruction of the TMM breast shell that is compared to an experimental reconstruction of the same phantom. . . . . . . . . . . . . Sagittal, axial, and coronal cross sections of the relative permittivity (left) and effective conductivity (right) of the first iteration of a simulated reconstruction of the TMM breast shell that is compared to an experimental reconstruction of the same phantom. . . . . . . . . . . . . Sagittal, axial, and coronal cross sections of the relative permittivity (left) and effective conductivity (right) of the first iteration of a simulated reconstruction of the TMM breast shell that is compared to an experimental reconstruction of the same phantom. . . . . . . . . . . . . Sagittal, axial, and coronal cross sections of the relative permittivity (left) and effective conductivity (right) of the first iteration of a simulated reconstruction of the TMM breast shell that is compared to an experimental reconstruction of the same phantom. . . . . . . . . . . . . Sagittal, axial, and coronal cross sections of the relative permittivity (left) and effective conductivity (right) of the first iteration of a simulated reconstruction of the TMM breast shell that is compared to an experimental reconstruction of the same phantom. . . . . . . . . . . . . Mean-squared error (MSE) vs iteration for reconstructing the TMM filled plastic shell in simulation and experimentally. . . . . . . . . . . . . . . Dipole array surrounding a breast phantom. . . . . . . . . . . . . . . . 58 59 65 66 67 67 68 68 69 ix 5.12 (a) Breast interior imaging domain of the dipole array. (b) Imaging domain surrounding the ring of antennas. . . . . . . . . . . . . . . . . 5.13 (a) DBIM imaging performance, assessed with the cosine metric, versus SER for different types of model error with a Class 3 phantom. (b) Extended domain DBIM imaging performance versus SER for different types of model error including antenna position errors. The baseline performance of both algorithms is indicated by a dashed line and markers denote 10% error in the dielectric properties of the material. . . . . . . 71 73 x abstract Microwave breast imaging is being investigated as a low cost, non-ionizing, and widely available imaging modality for breast cancer risk assessment, treatment monitoring, and tumor identification. The dielectric properties of the breast, which are linked to the tissue composition, are iteratively reconstructed using model-based inverse scattering algorithms. However, conventional microwave imaging methods have significant limitations in terms of resolution and experimental performance. Common microwave imaging methods have moderate-to-low resolution that distorts boundaries between different tissue types. Also, any errors in the simulated model of the array setup can corrupt the iterative process of these algorithms and degrade experimental imaging performance. Numerical 2D microwave imaging studies demonstrate that the use of a level set method preserves the natural boundary between tissue types resulting in a higher resolution reconstruction of breast structure. Previously proposed level set imaging algorithms have a high computational cost and are impractical in 3D. I developed a computationally tractable 3D microwave imaging algorithm based on level sets and assessed its efficacy for microwave imaging by investigating its performance with evaluating breast density. The density estimates from my level set algorithm are more accurate than those of conventional microwave imaging, and the accuracy is greater than that reported for mammographic density estimation. These results demonstrate the feasibility of high resolution 3D microwave breast imaging using a level set method. I also developed a quantitative metric for determining the amount of model error when modeling an experimental array and predicting the impact model error will have on microwave breast imaging. I demonstrate the use of this metric to evaluate experimental imaging performance, experimental calibration, and the performance of different imaging algorithms in the presence of model error. These results demonstrate the efficacy of this metric for improving the performance of microwave imaging in a clinical environment. Finally, I present some of my experimental calibration and imaging results xi for microwave imaging. These include an analysis of experimental calibration techniques, imaging results with different experimental phantoms, and imaging algorithms designed to perform better experimentally. 1 1 introduction In 2015 the American Cancer Society estimates there will be 231,840 new cases of women diagnosed with breast cancer and 40,290 deaths from breast cancer [6]. High volumetric breast density, defined as a large percentage of fibroglandular tissue, is one of the strongest predictors of a patient’s breast cancer risk [7]. Breast density estimation can play an important role in assessing cancer risk and monitoring preventative interventions designed to lower breast density. Another important application of breast imaging is to determine if a tumor will become invasive. This can prevent over-treatment and over-diagnosis of an indolent condition since not all forms of breast cancer will spread from the breast and be fatal [8]. However, current imaging modalities have several shortcomings in terms of 3D breast density evaluation and distinguishing between indolent and invasive tumors. Overdiagnosis and over-treatment has been estimated to cause anywhere from one to three deaths for every breast cancer death that is avoided using current screening techniques [9]. Improved breast imaging can further reduce the mortality rate, lead to better preventative treatments, and reduce the rate of over-diagnosis with breast cancer. The most common imaging modality for breast cancer screening is X-ray mammography. Mammography involves painful compression of the breast and exposure to harmful ionizing radiation. Repeated exposure to ionizing radiation is especially problematic and undesirable when monitoring the impact of preventative interventions, such as drug trials or lifestyle changes to lower breast density, for high risk patients. It is also difficult to accurately infer volumetric density from the 2D projection of the breast volume provided by a mammogram. Estimating 3D tissue distributions from 2D mammography results in a relative error of around 10-15% [10] or an absolute error of 5-10 percentage points [11]. It is not possible to determine which tumors will metastasize using X-ray mammography [8]. Using mammography to monitor the growth of a tumor to see if it will become invasive is undesirable as well since it requires repeated exposure to ionizing radiation. Magnetic resonance imaging (MRI) is another breast imaging modality. MRI 2 generates 3D images of the breast, which enables excellent tissue contrast in fibroglandular tissue and an excellent ability to distinguish adipose from fibroglandular tissue. MRI has high sensitivity ( 90% in a meta-analysis) for detecting breast cancer and moderate specificity (72%) [12]. MRI image contrast can be tuned to be weighted to proton density, but considerable development would be needed to make it a quantitative method of absolute proton density. Further development would also be necessary to make MRI imaging protocols correlative to variations in mammography density. Even if this development was complete, MRI exams are currently far too long and expensive in efforts to either detect or monitor dense breast tissue in a clinical application. Microwave tomography is a promising alternative breast imaging modality. Microwave tomography involves non-invasive measurements in the microwave frequency range and an inverse scattering algorithm to reconstruct the breast’s 3D distribution of dielectric properties. In the microwave frequency range there is a large contrast between the dielectric properties of fibroglandular and adipose tissue allowing for differentiation between healthy tissue types [1] and consequent estimation of volumetric density. There is a moderate contrast between healthy fibroglandular tissue and malignant tissue [2], and contrast agents are being developed to increase this contrast [13, 14]. Microwave imaging does not involve painful compression and ionizing radiation, unlike X-ray mammography. Also sequential microwave scans can safely measure the spatio-temporal changes of a tumor since microwave imaging does not involve ionizing radiation. The equipment used in microwave imaging is portable and low-cost, thus microwave tomography has the potential to be more widely available and cost effective than MRI for both estimating density and detecting invasive tumors. The performance of microwave tomography for breast density estimation depends heavily on the inverse scattering algorithm used to reconstruct the dielectric properties of the breast. Many common microwave tomography algorithms are based on iteratively solving a series of least squares problems with Tikhonov-Philips regularization [15], which results in smoothed edges in the reconstructed dielectric properties. This smoothing or blurring produces reconstructions that distort the 3 size and shape of the fibroglandular and adipose tissue, which results in inaccurate estimates of breast density and makes tumor identification difficult as well. A simulated forward solution that generates the measured data for a numerical model of the array and the current estimate of the dielectric properties in the breast is required at each iteration. Mismatch between simulated results and measured results due to errors between the assumed model and actual measurement system is a key factor influencing the performance of experimental microwave tomography. Model mismatch, or model error, can degrade imaging performance if it is larger than the scattered signals from the breast interior. The objectives of my research were to develop edge preserving techniques for microwave breast imaging and investigate methods to reduce model error in the forward solver of experimental microwave breast imaging. More accurate reconstructions of the edges between adipose, fibroglandular, and any malignant tissue will improve microwave tomography for breast density evaluation, tumor detection, and tumor monitoring. Reduction of model error or the sensitivity of imaging algorithms to model error will improve the performance of experimental microwave imaging. The following specific aims were pursued to achieve these goals: • Develop a microwave imaging algorithm based on the level set method and investigate its performance for breast density estimation using numerical testbeds. • Investigate calibration strategies and imaging algorithm modifications to mitigate the effects of mismatch between simulation and measurement in model-based microwave imaging algorithms. 4 2 2.1 background Introduction This chapter describes the state of current research involving microwave breast imaging. An overview of the current state of microwave imaging algorithms is given, followed by a section on incorporating level sets into inverse problems and a section on current microwave imaging techniques that use level sets. A brief history of how spatial priors have been incorporated into medical imaging modalities is given, as well as an existing technique that incorporates a spatial prior into microwave imaging. This chapter ends with a discussion on how model error impacts experimental microwave imaging. 2.2 Microwave Breast Imaging Microwave breast imaging typically involves surrounding the breast with an array of antennas that are designed to operate in the ultra-high frequency (UHF) range (300 MHz - 3 GHz), then one of the antennas is sourced and the electric fields are measured throughout the array. The scattered field that emanates from the breast is a function of the unknown dielectric properties as expressed in the electric field integral equation (EFIE) [16]. The scalar field formulation of the EFIE [16] relates the scattered electric field Es (r) at a given measurement point r due to an unknown object of permittivity (·) as follows: Es (r) = Et (r) − Ei (r) Z 2 b 0 t 0 0 b 0 = ω µ0 G (r|r )E (r ) (r ) − (r ) dr0 . (2.1) (2.2) V The frequency of interest is ω, and µ0 is the permeability of free space. The background medium’s permittivity is b (r0 ) and the Green’s function in that medium is Gb (r|r0 ). An analytic solution for the Green’s function exists for simple imaging do- 5 mains but must be computed numerically for more complicated imaging scenarios. The total electric field inside the unknown object, Et (r0 ), is unknown and must be approximated. The electric field at the measurement point r due to an excitation without the object present is the incident electric field Ei (r). The scattered field cannot be directly measured, only the total field, but (2.1) relates the total field at a measurement point to the known incident field and the scattered field due to the unknown object. Dielectric Properties of Breast Tissue In 2007, a large-scale study done by Lazebnik et al. measured the dielectric properties of healthy breast tissue from reduction surgeries [1]. Lazebnik et al. reported a large contrast ratio between the dielectric properties of high water content tissue (fibroglandular) and low water content tissue (adipose). Lazebnik et al. also conducted a large scale study using cancer surgeries to measure the dielectric properties of malignant tissue and reported a moderate contrast between fibroglandular tissue and malignant tissue [2]. Smaller studies in the past have reported a higher contrast between healthy and malignant tissue [17, 18], but similar ex-vivo measurements to Lazebnik et al. were observed in a recent study by Halter et al. [19]. The dispersive dielectric properties from Lazebnik et al. are shown in Fig. 2.1. Microwave Breast Imaging Algorithms Most microwave-based algorithms for breast cancer screening and detection fall into one of two groups: 1) radar based tumor detection and 2) microwave tomography. Radar based algorithms rely on detecting regions of strong backscatter due to a large contrast in dielectric properties [20–25]. The efficacy of radar based tumor detection is limited due to the moderate contrast between healthy and malignant tissue and the heterogeneity of the healthy tissues inside the breast. Contrast enhancement of the malignant tissue might improve radar imaging and make it more practical for tumor detection, but radar based methods are not designed to estimate breast density. Microwave tomography involves voxel by voxel reconstruction of the 6 60 30 median malignant median fibroglandular median adipose 40 30 20 10 0 0 25 Effective Conductivity Relative Permittivity 50 20 15 10 5 5 10 Frequency (GHz) (a) 15 20 0 0 5 10 Frequency (GHz) 15 20 (b) Figure 2.1: Measured dielectric properties for healthy breast tissue [1] and malignant breast tissue [2]. Solid red lines denote the 50th percentile of malignant tissue, the dash-dot lines represent the 50th percentile of fibroglandular tissue, and the dashed line is the 50th percentile of adipose tissue. Permittivity is plotted in (a) and conductivity is in (b). dielectric properties of the breast. The distribution of tissue can be inferred from the reconstructed dielectric properties. The breast density and the presence of malignant tissue are estimated from the resulting tissue distribution. This makes microwave tomography more versatile than radar based methods since it can be used for risk assessment and tumor detection. In recent years algorithms for microwave tomography have been proposed to solve the nonlinear ill-posed inverse scattering equations [4, 15, 26–33]. Most of these algorithms use an iterative approach that simulates the fields with the current estimate of the dielectric properties and then updates the estimate based on the difference between the simulated and measured fields. Some algorithms [33] use the Born iterative method (BIM) [16] to update the reconstructed properties, while others [4] use the distorted Born iterative method (DBIM) [16]. The DBIM has a faster convergence rate with high contrast scatterers, like the breast, than the BIM [34]. Another imaging technique [27, 30] is to use the contrast source inversion (CSI) method [35], which avoids simulating the fields every iteration by 7 solving for updates to the fields and unknown permittivity. CSI has been shown to exhibit similar performance to DBIM [36] but with less computational complexity. The algorithms in [15, 29, 31, 32] take a Gauss-Newton based approach [37] to find the estimated properties that minimize the squared error between measured and simulated fields. The Gauss-Newton approach is equivalent to the DBIM [38]. Overall, most of the microwave imaging algorithms under investigation are either based on the Gauss-Newton approach (or the equivalent DBIM), or exhibit performance similar to the Gauss-Newton approach. The techniques used to stabilize the DBIM (or the equivalent Gauss-Newton) results in Tikhonov-Philips regularization [15], which produces images with blurred edges or low spatial resolution. This blurring of edges reduces the effectiveness of microwave tomography for density estimation, tumor detection, monitoring changes in breast tissue, and monitoring tumor growth. An example of an image obtained using the DBIM as described in Shea et al. [4] is shown in Fig. 2.2. 5 10 15 20 25 30 35 40 45 50 55 Figure 2.2: Sequential coronal cross-sections, spaced every 8 mm, of the static permittivity for a phantom from the UWCEM repository [3] are shown in the first row, and the second row shows sequential cross-sections from a DBIM reconstruction of the static permittivity [4]. Both rows are plotted on the same color scale shown at the bottom of the figure. Several different ways to preserve edges in the dielectric properties reconstruction will be studied in this thesis. The first method is developing a level set based algorithm that is designed to preserve edges in microwave imaging. The efficacy of 8 incorporating a priori information about the tissue structure into a level set algorithm will also be investigated. Another proposed technique to preserve edges in microwave imaging is the use of total variation regularization, instead of TikhonovPhilips regularization, in microwave imaging. Experimental microwave imaging performance is negatively impacted by large model error. If the mismatch between measured and simulated fields is larger than the scattered fields from the breast interior, an imaging algorithm will attempt to account for the model error by changing properties in the imaging region. This will produce inaccurate reconstructions of the dielectric properties and inhibit the performance of every experimental imaging algorithm. Physical Dimensions and Dielectric Properties of Skin One of the first steps of any microwave imaging algorithm is to determine the location of the skin. The surface of the skin can be conformed to a known surface, physically measured, or determined with an imaging algorithm. Machining or 3D printing tolerances can be less than a millimeter and measurement techniques, like using a 3D digitizer, can be accurate within one to two millimeters. There are several microwave imaging algorithms that determine the outer surface of the skin within several millimeters [39, 40]. The skin thickness of the breast has been measured to range between 0.5-2.7 mm [41], and the dielectric properties of skin are well characterized from previous studies [42, 43]. Consequently the skin’s location, thickness, and dielectric properties are often assumed to be known in many imaging algorithms. 2.3 Level Set Imaging Level set techniques were originally used for computer vision and computer and segmentation applications, but in recent years they have been used for imaging modalities like X-ray computerized tomography (CT) and positron emission tomography (PET). A level set approach to solving inverse problems was first laid 9 out by Santosa [44]. Several studies have used a level set approach for microwave imaging of homogeneous objects in 2D [45–48], and there was an experimental study done with microwave tomography and a level set algorithm in 2D [49]. One research group has developed an inverse scattering algorithm based on level sets for imaging inhomogeneous objects [50–54], and applied it to 2D tumor detection for breast cancer screening with promising results [5]. In the context of microwave tomography a level set is a real valued function that is defined everywhere in the imaging domain. With a single level set any voxel that has a positive value will have dielectric properties corresponding to Region A, and any location with a negative or zero value will have dielectric properties corresponding to Region B. An example of a single level set segmenting a domain into distinct regions is shown in Fig. 2.3. φ>0 φ>0 φ60 Figure 2.3: 2D example of a single level set φ segmenting a sample domain into distinct regions. Sections with (φ > 0) will have dielectric properties that correspond to Region A, and sections with (φ 6 0) will have dielectric properties that correspond to Region B. 10 This setup produces the following model for permittivity in the imaging domain: (r) = 1 (r)H φ(r) + 2 (r) 1 − H φ(r) . (2.3) Here φ(r) is the level set, H(·) is the Heaviside or unit step function, 1 (r) is the permittivity inside Region A, and 2 (r) is the permittivity inside Region B. The permittivities 1 (r) and 2 (r) are spatially varying inside the regions, and can incorporate a priori information about the statistical distribution of dielectric properties for the corresponding tissue type. Irishina et al. [5] developed a level set model of permittivity using three level sets to represent four regions corresponding to different tissue types (skin, adipose, fibroglandular, malignant) as follows: (x) = 1 1 − H(φ1 ) + H(φ1 ) 2 1 − H(φ2 ) + H(φ2 ) 3 1 − H(φ3 ) + 4 H(φ3 ) . (2.4) The permittivity throughout the unknown object is (x) and it is a function of the permittivity inside each of the four regions 1 , 2 , 3 , 4 and the value of the level sets φ1 , φ2 , φ3 . When φ1 6 0 the permittivity is 1 and is a region that corresponds to skin. Sections where φ1 > 0 and φ2 6 0 correspond to a tumor with permittivity 2 . When φ1 > 0 and φ2 > 0 the level set φ3 differentiates between adipose and fibroglandular tissue. Sections with φ3 6 0 correspond to fibroglandular tissue with permittivity 3 , and φ3 > 0 corresponds to fatty tissue with permittivity 4 . Irishina et al. used the least squares cost function shown in (2.5), to determine how close the current reconstructed properties (φ1 , φ2 , φ3 , 1 , 2 , 3 , 4 ) are to a set of measurements. 1 J (φ1 , φ2 , φ3 , 1 , 2 , 3 , 4 ) = kR (φ1 , φ2 , φ3 , 1 , 2 , 3 , 4 ) k2 2 (2.5) The difference between the reconstructed and measured electric fields is the residual R (φ1 , φ2 , φ3 , 1 , 2 , 3 , 4 ) . Irishina et al. calculate the derivative of the cost function with respect to the evolution time t to determine descent directions for 11 each of the level set parameters as follows: * + 3 4 X X ∂ ∂ dJ dφ d 0 v v = < R []∗ R[], + . dt ∂φv dt ∂v dt v=1 v=1 (2.6) P Here < is the real part of the inner product in the parameter space denoted as h, iP . 0 The expression R []∗ R[] represents the Fréchet derivative of the residual and is calculated using an adjoint formulation [52]. The adjoint formulation requires two forward simulations every iteration. One of the forward simulations is computes the electric fields given a set of parameters and the other computes an adjoint simulation of the fields with different source parameters [52]. The adjoint formulation can result in a considerable computation burden, especially in 3D. Figure 2.4: Result from [5] where the far left is the phantom being simulated with three distinct tumors, the middle is a pixel by pixel reconstruction that serves as the initial guess, and on the right is the final reconstructed image. One of the results from Irishina et al. [5] is shown in Fig. 2.4. This result has more distinct boundaries between tissue types than the results shown in Fig. 2.2. All three tumors are easily distinguishable and are reconstructed in the proper location. Each iteration of the Irishina et al. algorithm is more computationally intensive than a Gauss-Newton iteration due to the adjoint formulation. Also, the Irishina et al. 2D reconstructions were obtained after 100 iterations, whereas the 3D reconstructions in Fig. 2.2 were obtained after only five iterations. The results 12 presented in Irishina et al. avoid the regularization that create blurred edges as with Gauss-Newton based methods, but at a high computational cost that does not extend well to 3D microwave tomography. 2.4 Model Error in Microwave Imaging Errors between the simulation and measurement of a known object can impact the reconstruction accuracy of most microwave imaging algorithms. Any errors in the electric field due to modeling errors of the antennas, the array, or materials in the domain will create artifacts when imaging an unknown object. Model error has not been sufficiently addressed by the microwave imaging community, but Burfeindt [55] proposed a metric to quantitatively assess the amount of model error is as follows: Smeas (i) . SMER = 20 log mediani (2.7) Smeas (i) − Ssim (i) The signal to model error ratio (SMER) is the median, over all channels, of the measured S-parameters for a channel, Smeas (i), divided by the difference between measured and simulated S-parameters, Ssim (i), for that channel. The experimental SMER was measured to be as low as 5-10 dB for an FDTD simulation model and the experimental array shown in Fig. 2.5. This indicates a significant amount of error between simulation and measurement that is deteriorating the experimental imaging performance. A single cause of this model error has not been found so I will focus on developing a quantitative metric for determining the amount of model error when modeling a system and for predicting the impact model error will have on microwave imaging. 13 Figure 2.5: Photograph of the UWCEM experimental imaging setup. 14 3 3.1 a level set method for 3d microwave breast imaging Introduction In the microwave frequency range there is a moderate to large contrast between the dielectric properties of different breast tissue types [1, 2]. This contrast reflects both the type and physiological state of the tissue. Hence, the dielectric properties distribution within the breast conveys clinically relevant anatomical and functional information. Reconstructions of the dielectric properties are potentially useful for several clinical applications including density estimation, tumor detection, and treatment monitoring. However, the inverse scattering problem is ill-posed, so microwave imaging performance depends heavily on the algorithm used to reconstruct the dielectric properties of the breast. Many common microwave imaging algorithms are based on iteratively solving a least squares problems with TikhonovPhilips regularization (e.g. [15]) which results in over-smoothing or blurring of the boundaries between tissues of contrasting dielectric properties. Level set algorithms [5,51,56] for microwave imaging have yielded high-resolution reconstructions of dielectric properties. An algorithm for imaging several 3D, homogeneous, non-overlapping targets is described in [56], but the only level set algorithms suitable for imaging the spatially complex and heterogeneous tissue structure of the breast are applied to 2D numerical breast phantoms [5, 51]. These 2D level set algorithms have a high computational cost and are impractical for 3D breast imaging. They involve a large number of iterations (over 100) [5], with each iteration requiring multiple runs of a computational electromagnetics simulation that can be time-consuming for large domain sizes [51]. In this paper, we present a computationally feasible level set method for 3D microwave breast imaging. Our algorithm requires only one electromagnetics simulation per iteration and fewer overall iterations than previously proposed level set approaches. The effectiveness of our algorithm is demonstrated by comparing our 3D level set reconstructions of anatomically realistic MRI-derived numerical phantoms [57] to those from the distorted Born iterative method [4]. 15 We quantitatively evaluate the performance of our level set method in the context of breast density estimation – a promising application for microwave imaging. Current imaging modalities for 3D breast density evaluation have shortcomings in terms of accuracy and cost. High volumetric breast density, defined as a large percentage of fibroglandular tissue, is one of the strongest predictors of breast cancer risk [7]. Hence, breast density estimation can play an important role in assessing breast cancer risk, as described by Harvey and Bovbjerg [58], as well as monitoring preventative interventions designed to lower breast density and the consequent cancer risk. We use density estimation as a biologically relevant performance metric for the reconstructions obtained with our level set method. The high resolution images produced by our level set method yields more accurate breast density estimates than those obtained with Tikhonov-regularized GaussNewton methods. The accuracy of our density estimates also exceeds that reported previously for mammographic density estimation [10, 11]. The next section presents the physical relevance of the level set technique for representing the dielectric properties of breast tissue. We also introduce a single parameter dispersion model that significantly reduces computational requirements. Our level set model for dielectric properties and the optimization algorithm we use for imaging are presented in Section 3.3. Section 3.4 describes the implementation details of our algorithm and Section 3.5 presents the results from our imaging study using numerical phantoms. A discussion of our results is found in Section 3.6 and the conclusions are presented in Section 3.7. Throughout this work boldface symbols represent vectors or matrices and superscript T represents matrix or vector ∂ transpose. The symbol ∇ represents the Fréchet derivative, whereas ∂x represents the partial derivative with respect to variable x. A preliminary version of this work was reported in [59]. 3.2 Dielectric Properties of Breast Tissue Previous studies have shown that there exists a large contrast in the dielectric properties of different types of healthy breast tissues at microwave frequencies [1]. 16 The dielectric properties of fibroglandular tissue, which consists of epithelial and connective tissue, are significantly higher than those of adipose (fatty) tissue. The contrast ratio of fibroglandular tissue to adipose tissue is large in both permittivity and effective conductivity [1]. This contrast in properties suggests using level sets to represent the distinctly different adipose and fibroglandular tissue types. The microwave inverse scattering problem is ill-posed, as the number of unknowns is much larger than the number of electric field measurements. The illposedness of the microwave breast imaging problem is reduced by using a dispersive model for the tissue dielectric properties and making measurements at multiple frequencies. Throughout the frequency range of interest (0.5-3 GHz) the dispersive nature of breast tissue is well-approximated by a single-pole Debye model for the complex permittivity [60]: (ω) = ∞ + ∆ σs + . 1 + jωτ jω0 (3.1) The only parameters in the Debye model that are assumed to be spatially varying are ∆, ∞ , and σs . The parameter ∆ is the difference between static permittivity (s ) and infinite frequency permittivity (∞ ). The static conductivity of the material is σs . A fixed, spatially invariant relaxation time constant (τ) of 15 ps is used for all tissue types. This Debye model results in three unknown parameters per voxel in the imaging domain. We further reduce the number of unknowns by exploiting the following linear relationship between the three unknown Debye parameters, ∞ = 0.3265s + 1.6326 (3.2a) σs = 0.0151s − 0.0365. (3.2b) The linear model above well approximates the underlying dispersion curves for both adipose and fibroglandular tissue as demonstrated in Fig. 3.1. Figure 3.1 compares the linear relationship in (3.2) to the Debye parameters that correspond to the 25th, 50th, and 75th percentile values [4] of measured dielectric properties 17 of breast tissue [1]. The resulting linear model reduces the number of unknowns in the imaging problem to one parameter, s , per voxel. Reduction to a single unknown parameter at each voxel significantly reduces the complexity of the multifrequency microwave breast imaging problem. Once s is reconstructed throughout the imaging domain, the frequency dependent dielectric properties of each voxel may be calculated according to the relationship in (3.2) and the Debye model of (3.1). 5 30 25 ε∞ ε ∞ 4 25th 50th 75th 3 2 4 5 6 εs 7 20 15 10 35 8 0.1 σs σs 40 45 50 55 45 50 55 εs 1.25 1 0.05 25th 50th 75th 0 4 25th 50th 75th 5 6 εs (a) 7 25th 50th 75th 0.75 0.5 8 0.25 35 40 εs (b) Figure 3.1: Linear Debye-parameter relationship of (3.2), shown by the black curve, compared to the Debye parameters that correspond to the 25th, 50th, and 75th percentile values [4] of measured dielectric properties [1]. (a) Infinite permittivity and static conductivity versus static permittivity of adipose tissue; (b) Infinite permittivity and static conductivity versus static permittivity of fibroglandular tissue. 18 3.3 Level Set Method for Microwave Imaging Mathematical Model for Permittivity The breast is composed of distinct healthy tissue types, namely adipose (fatty) tissue and fibroglandular tissue. On a microscopic (cellular) scale, there are clearcut natural boundaries between these different tissues. On a macroscopic scale, these distinct boundaries persist, albeit in a spatially complex manner due to the inherent heterogeneity of the the mammary network. Nevertheless, level sets can accommodate such spatial complexity. Thus level sets are indeed well defined for this problem. In this paper we use a single level set to represent the distribution of healthy tissue throughout the breast. A level set is a real valued function that is defined everywhere in the imaging domain. An example of a single level set segmenting a domain into distinct regions is shown in Fig. 3.2, where sections with (φ > 0) have dielectric properties that correspond to fibroglandular tissue, and regions with (φ 6 0) have dielectric properties that correspond to adipose tissue. This representation gives rise to the following model for complex permittivity in the imaging domain, as follows: a (r) = (r)H φ(r) + (r) 1 − H φ(r) . f (3.3) Here φ(r) is our level set, H(·) is the Heaviside or unit step function, f (r) is the complex permittivity inside the region corresponding to fibroglandular tissue, and a (r) is the complex permittivity inside the region corresponding to adipose tissue. The static permittivity in the fibroglandular and adipose regions, fs (r) and as (r) respectively, are used to calculate the complex permittivity in each region by substituting (3.2) into (3.1). The level set function segments the imaging domain into arbitrarily shaped regions of fibroglandular and adipose tissue. We note that this inherent segmentation makes breast density estimation a trivial calculation based on the level set function itself. The complex permittivities f (r) and a (r) 19 φ>0 φ>0 φ60 Figure 3.2: 2D example of a single level set φ segmenting a sample domain into distinct regions. Sections with (φ > 0) have dielectric properties that correspond to fibroglandular tissue, and regions with (φ 6 0) have dielectric properties that correspond to adipose tissue. may be spatially varying inside the regions; they may also incorporate any a priori information about the statistical distribution of the dielectric properties of tissue. Level Set Optimization We use an iterative approach to reconstruct the level set and corresponding permittivities of an unknown object from a set of electric fields measured with the unknown object present. The following least squares cost function is used to assess how similar the measured electric field Etm (rm , ωf ) is to the total electric field 20 Etr rm , ωf , (r) computed for the reconstructed properties, (r): 2 F M 1 X X t t r , ω , (r) E (r , ω ) − E C() = m f m f r 2 f=1 m=1 m FM 2 1 X = Ri . 2 i=1 (3.4) The location of the receiver for one of the M transmit-receive pairs is denoted by rm , while ωf is one of the F frequencies being used. The residual Ri is defined as the difference between the measured and reconstructed electric fields for a transmit-receive pair at a given frequency. We iteratively optimize the least squares cost function with respect to (r) as follows. An initial guess for (r) is made, then a forward solution is computed to obtain the corresponding electric fields. A search direction is identified, based on the difference between the reconstructed and measured electric fields, and is used to update each dielectric properties voxel. The search direction used in other level set based microwave imaging algorithms [5] is based on an adjoint solution [51]. This adjoint solution incurs a considerable computational burden, on the order of another forward simulation. We avoid this computational cost by calculating the Jacobian matrix instead. The Jacobian matrix of the least squares cost function in (3.4) is calculated using only the Green’s function and total field simulated inside a known permittivity distribution, as derived in [32], as follows: Gr (r|r0 ) ◦ Etr (r0 ) = J(ffl) ∇RT1 ∇RT2 = .. . ∇RTMF . (3.5) 21 Here, Gr (r|r0 ) is the MF × K Green’s function matrix, where K is the number of imaging voxels. Each row in this matrix represents the Green’s function, as it varies across the K imaging voxels, for a particular channel and freq ωf for a medium with permittivity . In the MF × K total field matrix Etr (r0 ), each row represents the total electric field for a particular channel calculated at each of the K imaging voxels for a specific frequency ωf . The Jacobian matrix for the least squares cost function evaluated for a specific permittivity is the element-by-element multiplication (Hadamard product) of the Green’s function matrix and the total field matrix. Each row of the Jacobian matrix is interpreted as the Fréchet derivative of the residuals in (3.4) with respect to the complex permittivity . The finite-difference time-domain method (FDTD) is used to compute the total electric fields necessary to evaluate the Jacobian at each iteration. The Green’s function matrix is calculated from FDTD simulations of the total field [61] assuming each source excitation is a z-directed current source, Iz , of known length, L. Each row in the Green’s function matrix is calculated as follows: Gr (rm |r, ω) = j Et |z (r|rm , ω). ωµ0 LIz r (3.6) This Green’s function calculation is described in detail in [4]. The Fréchet derivative of our cost function with respect to permittivity is calculated using the Jacobian matrix from (3.5) and the cost function in (3.4) as follows [37]: ∇C() = J()T R(). (3.7) Calculating the Fréchet derivative in this manner greatly reduces the computational complexity compared to the adjoint method [52]. Using our method requires only one FDTD simulation, whereas the adjoint method requires two electromagnetics simulations to obtain this derivative at every iteration. Therefore, our method of calculating the derivative using the Jacobian reduces our computational burden by half at every iteration. 22 The partial derivatives of (3.3) with respect to φ, as , and fs are computed as follows: ∂ f a = (r) − (r) δ φ(r) (3.8a) ∂φ f ∂ ∂ = H φ(r) (3.8b) f ∂s ∂fs a ∂ ∂ = 1 − H φ(r) (3.8c) ∂as ∂as where ∂f /∂fs and ∂a /∂as are derivatives of the complex permittivity in (3.1) with respect to the static permittivity with the linear model of (3.2) substituted for the ∞ and σs parameters. The partial derivatives of the cost function with respect to the level set parameters are obtained using the chain rule to combine (3.8) with the Fréchet derivative in (3.7), as shown below: ∂ ∂C = ∇C() ∂φ ∂φ ∂C ∂ = ∇C() ∂fs ∂fs ∂ ∂C = ∇C() a . a ∂s ∂s (3.9a) (3.9b) (3.9c) These partial derivatives are used to obtain a gradient-based descent optimization approach, as follows: ∂C ∂φ ∂C f,n+1 = f,n s s − αf ∂fs ∂C a,n+1 = a,n − αa a s s ∂s φn+1 = φn − αφ (3.10a) (3.10b) (3.10c) where n is the iteration number. The step sizes αφ , αf , and αa are chosen inde- 23 pendently to control the descent of each parameter while taking into account any a priori information about the tissue types and their structure. 3.4 Algorithm Implementation Details Extension Velocity The partial derivative with respect to the level set is only defined at the boundary between tissue types due to the delta function in (3.8a). Choosing a suitable region around the boundary to extend the derivative is often referred to as choosing an extension velocity [52]. There are many choices for extension velocities but here we simply extend the derivative to the entire imaging domain. Consequently, the level set function is updated everywhere in the imaging domain at each iteration, according to (3.10a). The permittivites for the two different regions do not need an extension velocity and are updated everywhere in the imaging domain according to (3.10b), (3.10c). Step Size Selection We determined appropriate step sizes αφ , αf , and αa in (3.10) based on a priori information about the statistical distribution of dielectric properties for the two tissue types. The step sizes are also chosen to take into account the general behavior of the optimization algorithm. The step size of the level set function αφ is chosen to allow a fixed percentage of the estimated number of fibroglandular voxels to change from adipose to fibroglandular (negative value to positive value) or vice versa. This step size prevents the algorithm from taking erroneously large steps that add too many fibroglandular or adipose voxels in one step. The step size for the region corresponding to fibroglandular tissue αf is larger than the step size for the region corresponding to adipose tissue αa to account for the larger variation in complex permittivity of fibroglandular tissue than adipose tissue as reported by Lazebnik et al. [1]. These step sizes can be increased or 24 decreased to encourage more or less heterogeneity within tissue types. At each iteration the static permittivity fs of the fibroglandular region is restricted to the range 20 to 65, and the static permittivity as of the adipose region is constrained to be in the range 2 to 10 based on the expected values for the measured tissue properties [1]. Multiple Frequency Approach We use several frequencies independently in our algorithm in a frequency hopping fashion similar to the algebraic reconstruction technique (also known as the Kaczmarz method) [62]. One frequency is used to update every parameter (φ, fs , as ) for multiple iterations and then we switch to the next frequency and so forth until the algorithm converges. This strategy of iterating through the frequencies is an attempt to avoid any local minima that might exist at one frequency but not at other frequencies. This approach prevents the steepest descent algorithm from getting stuck in the nearest local minimum associated with any individual frequency. We terminate the iterative process when the change in the cost function is less than 1% of the initial cost, after ensuring that each frequency is given the same number of iterations. Initial Parameter Generation Many different biologically relevant initial guesses may be used with the level set method. The initial guess we chose in this work is based on a distorted Born iterative method (DBIM) reconstruction obtained using the approach outlined in Shea et al. [4]. The boundary of the level set is chosen by thresholding the DBIM reconstruction with the threshold set midway between the 50th percentile for static permittivity of fibroglandular tissue and the 50th percentile for static permittivity of adipose tissue [1]. Inside each region the level set function takes on a value corresponding to the signed Euclidean distance from the boundary. Level set voxels corresponding to fibroglandular tissue assume a value of the Euclidean distance from the boundary, whereas those corresponding to adipose tissue are the 25 negative of the distance to the boundary as shown in Fig. 3.3. Using the Euclidean 10 20 30 40 50 (a) 60 10 20 30 (b) 40 50 60 −1 −0.5 0 0.5 1 (c) Figure 3.3: (a) Coronal cross section of static permittivity from a DBIM reconstruction. (b) Thresholded version of the DBIM reconstruction where fibroglandular tissue is represented by s = 49.1 (i.e. red) and adipose tissue is represented by s = 4.8 (i.e. dark blue). (c) Initial guess of the level set obtained using the DBIM reconstruction. Voxels in the level set are assigned a value based on the distance from the boundary between tissue types, with negative distance used for the region corresponding to adipose tissue. distance effectively reduces the likelihood of the algorithm changing voxels in the interior of the two initial guess regions. This strategy was chosen because DBIM reconstructions approximate the rough shape of the different regions well, but tend to blur boundaries, reducing confidence in the initial guess near those boundaries. The initial values for the permittivity in the region corresponding to fibroglandular tissue fs (·) and adipose tissue as (·) are homogeneous and set to the 50th percentiles as reported by Lazebnik et al. [1]. Numerical Phantoms We tested our level set method using nine MRI-derived numerical breast phantoms from the UWCEM online repository [3]. These 3D phantoms span the four breast density classes defined by the American College of Radiology Breast Imaging Reporting and Data System (BI-RADS) Atlas. The BI-RADS denotes Class 1 as mostly adipose, Class 2 as scattered fibroglandular, Class 3 as heterogeneously 26 dense, and Class 4 as very dense. The phantoms are surrounded by an oil immersion medium and a cylindrical array of 40 dipole antennas, arranged in five rings of eight antennas, as described in Shea et al. [4]. The location of the skin, its thickness, and its dielectric properties were assumed known for all reconstructions. The antennas were excited with an ideal wideband current source, and we record data at four frequencies: 1.0, 1.5, 2.0, and 2.5 GHz. The fields received at the antennas were computed using a 2 mm FDTD grid. Gaussian noise was added to obtain a 30 dB SNR, where SNR is defined as the ratio of total non-monostatic received power to total noise power. 3.5 Results We begin by comparing reconstructions obtained with our level set method to those obtained with the DBIM algorithm. Figures 3.4 and 3.5 present reconstructions of Class 3 and Class 4 phantoms, respectively. The first row shows the phantom properties, the second row shows the reconstruction using the DBIM algorithm described in Shea et al. [4], and the third row shows the reconstruction after 80 iterations of our level set method. The cross sections shown represent regions of the breast contained within the span of the rings of the antenna array. The eight coronal cross sections are spaced every 8 mm for both the Class 3 and Class 4 phantom, where the base of the breast phantom is nearest to the first column of images (far left). All the cross sections are of static permittivity. Images of the other Debye parameters (static conductivity and infinite permittivity) are similar due to the linear relationship of (3.2) and are omitted for brevity. The proposed level set method more accurately reconstructs breast tissue structure than conventional microwave imaging methods. For example the DBIM reconstruction misses fibroglandular tissue in the fifth and sixth columns of Fig. 3.4 that is accurately reconstructed by the level set method. The level set method also identifies the small regions of fibroglandular tissue extending from the bottom of the large fibroglandular tissue region in the second column of Fig. 3.5; this feature is not present in the DBIM reconstruction. The average absolute errors of the re- 27 10 20 30 40 50 60 Figure 3.4: Coronal cross sections of static permittivity from a Class 3 phantom (top) and the corresponding cross sections from the 3D DBIM reconstruction (middle) and the 3D level set method reconstruction (bottom). The coronal cross sections are taken every 8 mm between the top and bottom rings of antennas with the cross sections closest to the phantom base shown at the far left. constructed Debye parameters using DBIM and the level set method are shown in Table 3.1 for the adipose and fibroglandular regions of the Class 3 and 4 phantoms (Figs. 3.4 and 3.5). While the average absolute errors in adipose tissue properties are comparable between the two methods, the absolute error in fibroglandular properties are consistently smaller for the level set method. Table 3.1: Average absolute error in Debye parameters with the two reconstruction techniques DBIM Class 3 Class 4 Level Set Method s ∞ σs s ∞ σs Adipose 6.3 2.3 0.09 6.3 2.1 0.09 Fibroglandular 16.2 5.4 0.29 15.0 4.7 0.24 Adipose 7.2 2.7 0.11 7.2 2.4 0.11 Fibroglandular 16.4 5.6 0.27 16.3 5.1 0.26 28 Next we compare the density of the phantom to the densities estimated from the DBIM and level set method reconstructions. All three densities are determined by comparing the static permittivity at each voxel to a threshold. Static permittivity is chosen because it exhibits the largest difference between fibroglandular and adipose tissue types. The threshold is chosen to account for the variance in measured static permittivity by fitting Weibull distributions [63] to the adipose and fibroglandular properties, respectively, reported in Lazebnik et al. [1]. The threshold is chosen as the value giving equal area under the two distributions. Adipose tissue was modeled as a Weibull distribution with λ = 7 and K = 1.63, and fibroglandular tissue was modeled as a Weibull distribution with a λ = 51.82 and a K = 4.29. We used the coefficient of determination (R2 value) to assess the fit of these models to the 25th, 50th, and 75th percentiles of the Lazebnik et al. study [1]. The coefficient of determination was 0.91 for the fit to adipose tissue and 0.95 for the fit to fibroglandular tissue. These models lead to a threshold of 18.8 for static 10 20 30 40 50 60 Figure 3.5: Coronal cross sections of static permittivity from a Class 4 phantom (top) and the corresponding cross sections from the 3D DBIM reconstruction (middle) and the 3D level set method reconstruction (bottom). The coronal cross sections are taken every 8 mm between the top and bottom rings of antennas with the cross sections closest to the phantom base shown at the far left. 29 permittivity. Any static permittivity values greater than 18.8 were classified as fibroglandular tissue and any values less than or equal to 18.8 were classified as adipose tissue. The breast density estimates (volumetric percentages) for the nine phantoms of the UWCEM repository [3] are shown in Table 3.2. Applying the static permittivity threshold of 18.8 to the level set method reconstructions gives the exact same density estimate as the level set segmentation itself. This occurs because the level set method reconstruction forces all voxels in the fibroglandular region to have a static permittivity greater than 20 and all voxels in the adipose region less than 10, as noted in Section 3.4. Table 3.2: Volumetric percent density estimates for the nine phantoms and the two reconstruction techniques Breast Density (%) Class 1 Class 2 Class 3 Class 4 Phantom Level Set Method DBIM Phantom 1 2.6 1.6 1.7 Phantom 2 2.9 1.4 0.6 Phantom 1 5.0 2.8 2.1 Phantom 2 6.4 3.8 4.7 Phantom 3 4.5 2.9 2.6 Phantom 1 29.5 28.7 36.5 Phantom 2 24.2 22.9 33.2 Phantom 3 23.0 19.9 25.9 Phantom 1 33.4 33.4 42.7 The quantitative density estimates shown in Table 1 reveal that the level set method outperforms DBIM for the higher density classes, and performs similarly to DBIM for the lower density phantoms. The level set density estimates for the Class 3 and Class 4 phantoms, the two density classes with the highest cancer risk, have an average error of 1.3 percentage points or a 5% relative error. This performance is much better than the average error of 7.1 percentage points or 25% relative error 30 of the DBIM reconstructions. The density estimation performance of our level set method is also better than the 5-10 percentage points [11] or 10-15% relative error [10] reported for mammographic density estimation. We calculated density estimates in each cross section to further investigate imaging performance and evaluate whether the reconstructed fibroglandular tissue was being spatially localized correctly. Fig. 3.6 compares the fibroglandular tissue volume in sequential 2 mm cross sections of the phantom to that of DBIM and level set method reconstructions. The error bars in Fig. 3.6 represent the standard deviation of the density estimates over 10 different noise realizations with the same SNR. The sensitivity of our level set method to errors in the assumed skin properties was also tested. This represents a realistic error that might occur in a clinical application. We introduced 10% error in the dielectric properties of the reconstructed skin properties and compared the density estimates to reconstructions without errors. A 10% error in skin properties resulted in less than a 0.5% change in the overall density estimates reported in Table 3.2. Furthermore, the volume of fibroglandular tissue in individual cross sections changed less than 0.2 cm3 from the mean values shown in Fig. 3.6. This demonstrates our level set method is robust to noise and assumed skin properties, both of which are potential sources of error in clinical application. The results shown in Fig. 3.6 indicate that for the highest density classes (Fig. 3.6c and Fig. 3.6d) the level set method accurately reconstructs the tissue throughout the breast. DBIM overestimates the volume of fibroglandular tissue near the base and underestimates it in the apical coronal cross sections. Overall the fibroglandular tissue distribution of the level set method tracks the profile of the phantom confirming the spatial accuracy of the reconstruction. 31 Phantom DBIM Reconstruction LS Reconstruction 6 4 2 0 3 8 fibroglandular tissue volume (mL) fibroglandular tissue volume (mL) 8 6 4 2 0 4 5 6 7 8 9 10 distance from base of breast (cm) Phantom DBIM Reconstruction LS Reconstruction 2 3 4 5 6 7 distance from base of breast (cm) (a) (b) Phantom DBIM Reconstruction LS Reconstruction 6 4 2 0 2 3 4 5 6 distance from base of breast (cm) (c) 8 fibroglandular tissue volume (mL) fibroglandular tissue volume (mL) 8 8 Phantom DBIM Reconstruction LS Reconstruction 6 4 2 0 2 3 4 5 6 distance from base of breast (cm) (d) Figure 3.6: Volume of fibroglandular density in 2 mm thick coronal cross sections as a function of cross-section position (distance from the base of the phantom). The markers represent the mean value in the reconstructed cross section. Error bars represent the standard deviation of the estimates over 10 noise realizations. (a) Class 1 phantom, (b) Class 2 phantom, (c) Class 3 phantom, and (d) Class 4 phantom. 3.6 Discussion Current level set imaging algorithms for 2D microwave breast imaging have a high computational cost that make them impractical for 3D imaging. Level set algorithms like the one developed by Irishina et al. [5] use an adjoint method [52] that requires two electromagnetics simulations to calculate the Fréchet derivative at each iteration. 32 The algorithm of [5] also follows a staged approach that does not update every parameter at each iteration, which results in a large number of iterations (over 100) in 2D imaging. Our level set method utilizes the Jacobian matrix to calculate the Fréchet derivative which only requires one electromagnetics simulation per iteration, effectively cutting our computational burden in half. Our algorithm also updates every level set parameter each iteration, utilizes frequency hopping across different iterations, and incorporates a priori information to determine step sizes and bounds for the properties of the different tissue regions. This optimization strategy gives our algorithm a faster convergence rate than other level set algorithms and results in fewer iterations (≈80) in 3D imaging. We selected DBIM as our benchmark instead of other microwave breast imaging algorithms for several reasons. DBIM, a well known technique that is equivalent to the Gauss-Newton algorithm [38], has Tikhonov regularization when implemented using conjugate gradient least squares [64] as described in Shea et al. [4]. It is able to localize coarse features of the breast, which motivates the use of DBIM for generating an initial guess of the level set. Robustness of the level set method to the initial guess was not investigated in this study. Certainly other initial guesses (for example, Kurrant et al. [65]) could be employed for the level set method. In this study, we used breast density evaluation to quantitatively assess the efficacy of our level set method for 3D microwave imaging. Conventional microwave imaging algorithms often use a minimum norm regularizer to select a unique solution to the ill-posed inverse problem. Minimum norm regularizers are know to blur the tissue structure of the breast and produce misleading density estimates as illustrated by the DBIM results in Table 3.2. This blurring inhibits the detection of small tissue features, especially in regions with less fibroglandular tissue, such as the apical coronal cross sections of the breast. This blurring effect also erroneously extends the size of large tissue features like those near the base of the breast. Both of these effects are apparent in Figs. 3.6c and 3.6d, and their impact on the overall density estimate is clearly evident in Table 3.2. Our level set method does not suffer from these blurring effects since it incorporates a priori information about the tissue properties and preserves distinct boundaries between tissues. Requiring 33 the solution to contain two tissue types constrains the set of possible solutions and forces the algorithm to produce physically meaningful results. This results in a more accurate boundary between tissue types. Consequently, our level set method improves upon the reconstructed tissue structure of DBIM and results in better localization and identification of the fibroglandular tissue as shown in Fig. 3.6 and Table 3.2. The level set method presented here was formulated with a single level set to accommodate two types of tissue: adipose and fibroglandular tissue. This formulation is relevant to imaging healthy breast tissue, as is the case for breast density estimation. An additional level set can be added to the model in (3.3) to represent cancerous tissue in a manner similar to that presented by Irishina et al. [5]. This addition would allow the application of the algorithm to tumor detection. 3.7 Conclusion We have developed a level set method that requires fewer iterations and fewer electromagnetics simulations per iteration than other level set based implementations and leads to a feasible level of computational complexity for full 3D imaging. Reconstructions of anatomically realistic numerical breast phantoms demonstrate the ability of the proposed algorithm to reconstruct the tissue structure better than conventional microwave imaging approaches. The entire-breast density estimates as well as the sequence of coronal cross section density estimates demonstrate the high accuracy that our level set method achieves in reconstructing the actual tissue structure of the phantom. 34 4 model error in model-based inverse scattering algorithms for microwave breast imaging 4.1 Introduction Model-based optimization algorithms used to solve the inverse scattering problem iteratively minimize the difference between measured scattered electric fields from the actual (unknown) object and simulated scattered electric fields from an estimate of the unknown object’s properties. Any errors in the simulated model of the array environment can corrupt this iterative process by introducing differences in the electric fields due to errors in the model of the known environment instead of errors in the dielectric properties estimate of the unknown object. Consequently, model error can reduce the image quality of microwave breast imaging. Potential sources of model error include the difference between the simulated and experimental performance of various antennas used in microwave imaging systems [66, 67]. Prior investigations have yielded imaging algorithms that are designed to be robust to noise [68] and other sources of model error [69]. There has also been work to analyze how system design parameters can affect model error [70] and how immersion medium properties can impact microwave imaging performance [71]. However, to date there has been no common mathematical framework for analyzing different types of model error or comparing the performance of different arrays, calibration techniques, or imaging algorithms in the presence of model error. Furthermore, to the best of our knowledge the previous studies have been conducted in the absence of a suitable metric for quantifying model error in the context of imaging. In this paper we propose a metric for quantifying the amount of model error in the forward solver of an inverse scattering algorithm and the impact it has on microwave breast imaging. We develop a mathematical framework for the effect of model error on microwave imaging and demonstrate the use of the proposed metric for analyzing an experimental microwave imaging system. 35 The rest of this paper is organized as follows: Section 4.2 presents a brief overview of model-based inverse scattering algorithms used in microwave imaging and describes a mathematical framework and quantitative metric for evaluating the significance of model error, Section 4.3 describes the experimental and simulated testbeds used in this study, and Section 4.4 presents the use of the proposed metric to assess system calibration, experimental imaging performance, and the performance of different imaging algorithms in the presence of several types of model error. 4.2 Model Error Metric The theoretical underpinning for most model-based inverse scattering algorithms is the electric field integral equation [16], which relates the scattered field to the scattering object’s properties as follows: Z s 2 E (r) = ω µ h i Gb (r|r0 ) Et (r0 ) (r0 ) − b (r0 ) dr0 . (4.1) V Here (·) is the spatial distribution of dielectric properties for the unknown object and b (·) is the dielectric properties of the background environment. The total electric field is Et (·) and the background Green’s function is Gb (·). The frequency of the scattered field is ω, and µ is the permeability of free space. The measurement location is denoted as r and the integral is over the volume of the unknown object denoted V. The total electric field in (4.1) depends on the dielectric properties of the object, which are unknown. Model-based inverse scattering algorithms typically approximate the total electric field with a simulated background electric field via the Born approximation [16]. There is no closed form expression for the Green’s function when the background is heterogeneous, so it must be computed via numerical simulation. The simulated scattered electric fields outside the unknown object, the simulated background electric fields inside V, and the background Green’s function electric 36 fields are calculated using a numerical solver that models the physical setup, the antenna or source excitation, and the estimated properties for the unknown object. Differences between the measured and simulated scattered fields could be due to errors in the estimated properties of the object and the objective of model-based inverse scattering algorithms is to minimize those differences by improving the estimate of the object’s dielectric properties. However, they could also be due to errors in the modeling of the system setup. Discretizing the integral in (4.1) for K voxels in the imaging domain results in the following expression: Es = Gb ◦ Et − b . (4.2) Here, Gb and Et are 1 × K vectors, and b are K × 1 vectors, and ◦ is the element by element multiplication (Hadamard product). The terms in front of the integral as well as any scale factors produced by the discretization of the integral are absorbed into the GreenâŁ™s function vector Gb for notational simplicity. Suppose we use the EFIE to simulate the scattered electric field associated with the true dielectric properties of the object in the presence of model error. We may write Es = Gb + α Gb ◦ Et + β Et − b − δb (4.3) where Gb , Et , and b are the true Green’s function, true total electric fields, and true baseline dielectric properties. The nonlinear function α Gb represents model error in the simulated Green’s function. This error could result from errors in the dielectric properties of the antenna substrate or immersion medium, the locations of antennas, or source properties. Note that the error is a function of the true Green’s function as it depends on the true underlying properties of the system. For example, an error in the assumed dielectric properties of the immersion medium will have different effects with different antenna arrays. The nonlinear function β (Et ) represents errors between the true total electric field and the simulated one, 37 while the δb term represents errors in the assumed dielectric properties for the baseline conditions. The expression in (4.3) may be rewritten to combine all the terms associated with model error into an error term Ese , as shown by Es = Gb ◦ Et ( − b ) + Ese = Eso + Ese . (4.4a) (4.4b) Here Eso is the true scattered field, i.e. the scattered field in the absence of model error due solely to the difference between the object properties and the baseline properties. Hence, the simulated scattered field in the presence of model error is expressed as the sum of the scattered field in the absence of model error, which we term “signal”, and an error term. The absolute power level of Eso and Ese are not of great importance, since they will change with array design, different objects, and even input power to the system; instead we are interested in their relative size. We define the signal-to-error ratio SER = 10 · log |Eso |2 |Ese |2 ! (4.5) as a metric for characterizing the impact of model error on the image reconstruction problem. This metric represents the size of the signal used to reconstruct the feature of interest relative to model error in describing the system. This metric is similar in concept to signal to noise ratio. The SER metric is based on the power of the signal associated with the object of interest. The term in the denominator reflects the size of the model error. Some model-based inverse scattering algorithms will perform better at a given SER than others algorithms, similar to the case where some algorithms perform better at a given SNR than others. 38 Signal power calculation The numerator of the SER may be approximated with S-parameters from two scenarios, one containing a baseline scenario and another containing the baseline plus the object of interest. That is, we approximate |Eso |2 as |Eso |2 M X F X So (i, f) − Sb (i)2 . ∝ (4.6) i=1 f=1 The index i is used to sum over the M channels of the array (transmit-receive pairs), and the index f is used to sum over the F frequencies of interest used in the reconstruction. The subscript on the S-parameters denotes the object of interest or baseline scenario from the desired testbed (either simulated or measured). Accordingly, So (·) refers to the S-parameters for the object, and Sb (·) refers to the measured S-parameters for the baseline scenario in the same testbed. One example for microwave breast imaging would be when reconstructing the tissue structure of the breast, the full phantom is the entire breast phantom and the baseline phantom is the same phantom profile but the interior is replaced by a homogenous material. Another example is when attempting tumor detection, the full phantom is the entire breast including the tumor and the baseline phantom is only the healthy tissues of the breast with no tumor present. This method of estimating the signal power ensures that the calculated difference is due to the feature we are interested in reconstructing and not any other signal present in the system, e.g. reflections from other objects, and paths between antennas that do not interact with the feature. This quantity is proportional to the SER numerator since S-parameters are calculated from normalized electric fields, but the ratio of terms in (4.5) will cause these normalization factors to cancel. This calculation can be impractical to conduct experimentally with spatially complex phantoms and will result in small changes to the setup or phantom when adding or removing the feature of interest. Using simulations would allow us to avoid changes in the setup or phantom that might occur when attempting experimentally, as well as determine this quantity with spatially complex phantoms. 39 However, this simulated difference will not be identical to the true experimental difference since the true Green’s function, total field, and baseline dielectric properties are unknown. Subsequently, the agreement between the measured and simulated signal power will depend on the match between the true and simulated properties. A comparison between the simulated and experimental signal power will be presented in this study for an experimental phantom. Error calculation We determine the denominator of the SER in (4.5) as the difference between experimental and numerical S-parameters for the phantom with the object or feature of interest. Hence, the SER denominator is proportional to |Ese |2 M X F X So,m (i, f) − So,s (i, f)2 . ∝ (4.7) i=1 f=1 Here i is an index of the M channels, and f is summed over the F frequencies of interest. The first subscript on the S-parameters denotes the object of interest or baseline scenario, and the second subscript denotes that the S-parameters are simulated or measured. The power of the scattered electric field due to model error is estimated as the squared difference between an experimental measurement of the object So,m (·) and a simulated version of the object So,s (·). This difference includes model error, background noise present in the experimental measurement, co-registration errors between simulation and measurement, and errors in the phantom model. Meticulous co-registration and measuring the properties of the phantom can reduce the impact of mismatch between the experimental and simulated phantoms, but errors in the phantom model will be due to the same inaccuracies that cause model error of the system setup itself. Also, other measurement techniques can be used to assess the amount of background noise present in the system and determine whether model error or background noise contributes more to this quantity. 40 For applications to simulated scenarios this quantity could also be calculated as |Ese |2 ∝ M X F X So,s (i, f) − So,ŝ (i, f)2 . (4.8) i=1 f=1 The power of the scattered electric field due to model error is estimated as the squared difference between a simulation of the object So,s (·) and a simulated version of the object So,ŝ (·) where the system properties were changed to introduce model error. For example, to test the impact of errors in the immersion medium of a simulated array the first simulation (So,m (·)) contains the object and assume background properties, while the second simulation So,ŝ (·)) is identical except with a perturbation in the dielectric properties of the immersion medium. This allows the metric to be calculated before an experimental array is built and assesses the sensitivity of a proposed array or imaging algorithm to common sources of model error. 4.3 Experimental and Numerical Testbeds Miniaturized Patch Antenna Array We developed an array, based on the design from Burfeindt et al. [72], that contains 32 miniaturized patch antennas on a Rogers RO4360G2 substrate to investigate the effects of model error experimentally. The patch antennas have resonances at 1.36 GHz and 2.31 GHz and the antenna layout is shown in the digital model in Fig. 4.1b. As shown in Fig. 4.1a, the patch array is connected to a ScenarioTek switch matrix and data for each channel is collected using an Agilent PNA. The digital model shown in Fig. 4.1b and the finite-difference time-domain (FDTD) method were used to simulate the electric fields with various objects in the array. A two-pole Debye model was used to represent dispersive media. A graded mesh, varying from 0.5 mm to 1 mm, was used to model the fine features of the substrate and the coarser features of the interior. The domain was terminated with 41 (a) (b) Figure 4.1: (a) Photograph of the experimental patch array. (b) Simulated model of the patch array. a convolutional perfectly matched layer absorbing boundary. The antenna feeds were modeled as lumped-element resistive voltages sources and the electric fields were recorded at the two resonant frequencies. Tissue Mimicking Phantom A phantom with dielectric properties similar to healthy breast tissue was created to study the effect of model error on calibration and imaging. The outer surface of a numerical phantom from the UWCEM breast phantom repository [3] was converted to a facet model and printed using a 3D printer, similar to the procedure described in [73]. As shown in Figures 4.2a and 4.2b, six mm thick plastic was used to replace the skin of the numerical phantom, while the interior of the numerical phantom was 42 replaced with a void. A plastic slab was added to the top of the phantom to allow for positioning within the array. The void in the 3D-printed surface was then filled with tissue mimicking materials (TMMs) that were designed using the procedure in [74]. The bulk of the phantom was filled with an adipose tissue mimicking material and then two test tubes of different sizes, that contained a fibroglandular mimicking material, were added as shown in Fig. 4.2c. A simulated version of the phantom is shown in Fig. 4.3, where the dielectric properties of the TMMs were measured using an Agilent 85070E slim form probe. (a) (b) (c) Figure 4.2: Side view (a) and top view (b) of the 3D printed breast surface used in calibration and imaging results. (c) Picture of 3D printed breast surface filled with TMMs; the two cylinders mimic fibroglandular tissue and the remaining material mimics adipose tissue. The 3D printed surface in Fig. 4.2a was also filled with different materials to test calibration of the TMM phantom. Dielectric properties of the materials of the phantom and properties of potential calibration objects are shown in Table 4.1. Once again the dielectric properties of the liquids were measured using an Agilent 85070E slim form probe, and the plastic was characterized using the technique specified in Burfeindt et al. [73]. 43 0 2 cm 5 10 15 0 0.5 1 Figure 4.3: Sagittal, axial, and coronal cross sections of the relative permittivity (left) and effective conductivity (right) of a simulated version of the 3D printed breast surface filled with TMMs. Table 4.1: Dielectric Properties at 2 GHz Material Relative permittivity Conductivity (S/m) Oil immersion 2.8 0.03 3D printed plastic 2.3 0.004 Adipose TMM 10.2 0.3 Fibroglandular TMM 34.3 1.0 Glycerol 6.1 0.4 50/50 Glycerol and Water 43 2.1 Water 76 0.9 44 Dipole Array A simulated dipole array was used to investigate the effects of model error on microwave imaging. This array is identical to the setup described in Shea et al. [4]. It consists of 40 electrically short dipoles; 5 rings of 8 antennas that surround a numerical phantom from the UWCEM repository [3]. The FDTD method is used to simulate array measurements of the numerical phantoms and data is collected at 1.0, 1.5, 2.0, and 2.5 GHz. The simulations were conducted on a 2.0 mm grid and the domain was terminated with uniaxial perfectly matched layer boundary conditions. The dispersive tissues of the numerical breast phantoms were represented using a single-pole Debye model. Imaging Algorithms Two different imaging algorithms are used to reconstruct the phantoms from experimental or numerical measured fields. The interior domain imaging algorithm is identical to the DBIM algorithm described in [4]. The exterior domain algorithm is similar except the imaging domain is not only the inside of the breast, but also includes a region outside the breast. The extended domain DBIM method can change properties both interior and exterior to the breast. 4.4 Results Three different scenarios are employed to illustrate the utility of the SER metric. In the first subsection we explore the relationship between SER, various types of model error, and image quality using simulated data acquired with the dipole array. In the second subsection we use SER to evaluate the effectiveness of different calibration objects for reducing experimental model error in the patch array. Finally, the third subsection evaluates the effect of model error with an experimental image of the phantom in Fig. 4.1 using the patch array. 45 SER and Imaging Performance in the Presence of Model Error We begin by evaluating the relationship between SER and imaging performance with different levels and types of model error. These results are based on simulated data from the dipole array described in Section 4.3 using class 2 and class 3 anatomically realistic numerical phantoms. Use of simulations enables exploration of a wide variety of model errors in a controlled fashion. Figures 4.4 and 4.5 illustrate the performance of both the interior and exterior domain algorithms with different types of model error and different levels of SER for a Class II and Class III phantom from the UWCEM repository [3]. The SERs in Figures 4.4 and 4.5 were calculated in simulation with the complete numerical breast phantom as the full phantom and a homogenous phantom as the baseline phantom using (4.6). The denominator of (4.5) estimates the model error power using (4.8) for the specific level and type of model error being studied. The performance of these algorithms was tested with varying levels of three different types of model error: noise, immersion medium dielectric properties error, and skin dielectric properties error. The noise model error results are based on adding white Gaussian noise to the simulated data and averaging imaging performance over five realizations of the noise at each level. The skin dielectric properties are varied from 55% below to the true value and the and immersion medium dielectric properties are varied from 15% below to the true value to obtain the range of SER values shown. The markers on the skin and immersion properties curves depict a 10% error between properties used for data acquisition and image reconstruction. The dashed line depicting baseline performance represents the value of the cosine metric when there is no error between the numerical phantom and model used for image reconstruction. 46 0.84 Cosine Metric 0.82 0.8 0.78 0.76 0.74 Noise Immersion Medium Error Skin Error Baseline Performance 0.72 0.7 0.68 0 10 20 30 SER (a) 0.84 Cosine Metric 0.82 0.8 0.78 0.76 0.74 Noise Immersion Medium Error Skin Error Baseline Performance 0.72 0.7 0.68 0 10 20 30 SER (b) Figure 4.4: (a) DBIM imaging performance, assessed with the cosine metric, versus SER for different types of model error with a Class 3 phantom. (b) Extended domain DBIM imaging performance versus SER for different types of model error. The baseline performance of both algorithms is indicated by a dashed line and markers denote 10% error in the dielectric properties of the material. 47 0.84 Cosine Metric 0.82 0.8 0.78 0.76 0.74 Noise Immersion Medium Error Skin Error Baseline Performance 0.72 0.7 0.68 0 10 20 30 SER (a) 0.84 Cosine Metric 0.82 0.8 0.78 0.76 0.74 Noise Immersion Medium Error Skin Error Baseline Performance 0.72 0.7 0.68 0 10 20 30 SER (b) Figure 4.5: (a) DBIM imaging performance, assessed with the cosine metric, versus SER for different types of model error with a Class 2 phantom. (b) Extended domain DBIM imaging performance versus SER for different types of model error. The baseline performance of both algorithms is indicated by a dashed line and markers denote 10% error in the dielectric properties of the material. 48 Evaluating Experimental Calibration Objects using SER The experimental array shown in Fig. 4.1 and the SER metric were used to quantitatively assess the ability of calibration to remove model error from the phantom shown in Fig. 4.2. Numerical simulations and data collected with the experimental phantom were used to calculate the SER of this array and phantom. The location of the plastic surface and array dimensions were measured, and the dielectric properties of the plastic, oil immersion, and substrate were measured or supplied by the manufacturer. In this imaging scenario the feature of interest is the two fibroglandular mimicking cylinders embedded in the adipose mimicking material. In terms of (4.6), the object measurement contains the two cylinders whereas the baseline measurement is homogenous adipose TMM inside the plastic breast surface. The model error power was calculated using the difference between a measurement of the phantom with both cylinders and a simulation of the phantom with both cylinders according to (4.7). The only calibration performed was correcting any phase errors and attenuation due to the cables, switch matrix, and PNA. The SER for this scenario is -22 dB. This means that the signal used to reconstruct the two fibroglandular mimicking cylinders is 22 db below the agreement between our simulations and measurements of this phantom. The level of background noise present in the system was determined using back to back measurements of the data and was estimated to be 60 db below the signals used in reconstruction. Therefore noise was not a significant source of error in this array. The numerator of the SER could also be calculated in simulation instead of experimentally. Using simulated data collected with the phantom shown in Fig. 4.3 and the same phantom with the fibroglandular TMM material removed, the SER was -23 dB instead of -22 dB obtained with an experimental numerator. Therefore, there is excellent agreement between these two methods of estimating signal power with this phantom and with our match between simulation and experimental. Using simulations to estimate signal power allows the metric to be more broadly used, but with more mismatch between simulation and experimental these quantities would not agree as well. 49 We investigated the efficacy of using a calibration object to increase the SER above -22 dB, or equivalently to decrease this error between simulation and measurement. This calibration method involves calculating channel by channel complex correction factors using a measurement and simulation of a known object. The correction factors are the ratio of the complex S-parameter for a simulation to the measured complex S-parameter, where both S-parameters are calculated with the calibration object in the array. This correction factor will then be multiplied by the measured S-parameters for the phantom in Fig. 4.2 to correct any magnitude and phase differences between measurement and simulation. The improvement in SER using different calibration objects is shown in Table 4.2. A simulation and measurement of each calibration object in Table 4.2 was used to generate calibration factors. These calibration factors were applied to the TMM phantom with two cylinders, and the effectiveness of the calibration was determined using the SER metric. Table 4.2: SER improvement with calibration Calibration object SER (dB) Oil filled array -9 50/50 glycerol-water mixture filled breast surface -5 Water filled breast surface -5 Glycerol filled breast surface -4 Adipose TMM filled breast surface 5 These calibration results combined with the SER metric demonstrate that using a calibration object that is drastically different than the imaging object (like the oil filled array) reduces the model error much less than an object that is similar to the imaging object. Likewise, the calibration objects that have the exact same shape but are filled with different material, like the breast surface filled with glycerol, water, or mixtures of both, also reduce model error much less than using the same material that is the bulk of the imaging object, i.e. the adipose TMM filled breast surface. This is an unrealistic calibration scenario with human subjects, but it demonstrates 50 how similar a calibration object needs to be to the imaging object in order to improve model error significantly. Imaging Results with an Experimental Array The Distorted Born Iterative Method (DBIM), as described in [4], was used to image the phantom and determine the impact of our experimental SER on imaging performance. The adipose TMM filled breast surface was used as a calibration object, as described in the previous section, since it had the largest improvement in SER. Consequently, the calibrated SER was 5 dB for this imaging scenario and was calculated using (4.5), where the S-parameters for the measured phantom were multiplied by the complex calibration factors. The experimental imaging result is compared to a simulated imaging scenario with no model error in Fig. 4.6. In Fig. 4.6 the simulated imaging scenario is done without noise or model error and demonstrates the high quality image that can be obtained with the DBIM for the phantom shown in Fig. 4.3. The experimental imaging result shown is the best image we could obtain with the given level of model error. The experimental imaging result has much lower resolution, larger errors in the dielectric properties, and localizes the two cylinders worse than the simulated imaging scenario without model error. Modifying the algorithms optimization parameters or using any additional iterations did not improve the experimental imaging results at all. 51 0 2 cm 5 10 15 0 0.5 1 0.5 1 (a) 0 2 cm 5 10 15 0 (b) Figure 4.6: (a) Sagittal, axial, and coronal cross sections of the relative permittivity (left) and effective conductivity (right) for a reconstruction of the TMM filled breast shell (Fig. 4.3) using simulated data without model error. (b) Sagittal, axial, and coronal cross sections of a reconstruction of the TMM filled breast shell (Fig. 4.2c) from experimental measurements with 5 dB SER. 52 4.5 Discussion Several trends are noteworthy in Figures 4.4 and 4.5. Noise effects the interior and exterior domain DBIM algorithms comparably. In contrast, the exterior domain algorithm is much more robust to immersion medium errors. This is consistent with the fact that the exterior domain algorithm adjusts the dielectric properties in the region surrounding the breast. The exterior domain algorithm is also slightly more robust to errors in skin dielectric properties. Differences in the cosine metric between the class II and III phantoms are not meaningful as the class II phantom has a much lower degree of heterogeneity. The higher performance of the exterior-domain algorithm is likely directly due to the fact that it adjusts the dielectric properties in the region surrounding the antennas. Line of sight propagation paths that do not pass through the breast interior do not provide significant information about the interior of the object, but do provide information about the immersion medium and exterior environment. The interior-domain algorithm can only account for errors in these propagation paths by changing the properties of the breast interior. This leads to significant reductions in image quality for small errors in the assumed immersion medium properties. In contrast, the exterior-domain algorithm can easily compensate for errors in the immersion medium without disturbing the properties of the breast interior. The SER associated with a 10% error in skin properties is significantly greater than the SER for a 10% error in the immersion medium, especially in the class III phantom. The impact of skin dielectric properties errors for both algorithms is similar. These observations support the conclusion that errors in skin properties have less an impact on image quality. This is likely due to the much smaller volume occupied by skin and consequent reduced influence on the fields. Figures 4.4 and 4.5 clearly show that image reconstruction quality is proportional to SER. A threshold phenomenon is evident. Once the SER drops below a certain threshold, the performance begins to rapidly deteriorate. The threshold SER varies for different algorithms. It also varies for different types of error with 53 the interior-domain imaging algorithm. These observations are consistent with prior expectations. High quality imaging generally results in these examples for SERs greater than 20 dB. The strong dependence on calibration object shown in Table 4.2 can be explained by the non-linearity of electromagnetic scattering and model error as described by (4.3). The model error is a nonlinear function of the object and errors, so a single correction factor calculated with a calibration object will not significantly reduce the error of another object unless the two are very similar. As shown here, the SER metric can be used to evaluate the efficacy of using different calibration objects to remove model error and potentially other calibration techniques as well. The experimental imaging result shown in Fig. 4.6 has much lower resolution, larger errors in the dielectric properties, and localizes the two cylinders worse than the simulated imaging scenario without model error shown in Fig. 4.3. This degradation in imaging performance is predicted by low SER because the signal used to reconstruct the interior of the phantom is not much larger than the model error. Eventually the difference in scattered electric fields due to model error dominates, and the algorithm is updating properties based on model error instead of differences due to the interior properties. 4.6 Conclusion In this paper we presented a mathematical framework for evaluating the impact of model error on microwave breast imaging. The proposed metric relates the power of the received signal due to the desired feature, to the power of the received signal due to model error. This metric was used to assess the performance of experimental calibration techniques for removing model error. We presented the SER for an experimental prototype microwave breast imaging array and compared the experimental imaging performance to an analogous simulated array without any model error. Model error degrades the experimental imaging performance compared to the model error free simulated results, even when using an excellent calibration object. In this paper the SER metric was also used to assess the imaging 54 performance in a simulated array with different types and varying levels of model error. An imaging algorithm that can modify the background dielectric properties in an extended imaging domain is more robust to model errors than an algorithm that is restricted to imaging only inside the breast. Consequently, analysis employing the proposed SER metric can be used to influence array design, calibration techniques, and imaging algorithms in order to improve the performance of experimental microwave breast imaging. 55 5 5.1 calibration and imaging with model error Introduction I investigated several calibration techniques and imaging algorithms to improve microwave imaging performance in the presence of model error. I determined the efficacy of using calibration methods to reduce model error and investigated imaging algorithms that are robust to model error using the metric described in the previous chapter. 5.2 Simple Phantom Calibration I tried using a single calibration object to correct the channel by channel differences between simulation and measurement with the array described in Fig. 4.1. The choice of calibration object must be similar to the unknown object due to the nonlinearity, with respect to the dielectric properties of the object, of the scattered fields. Measured and simulated fields of a known object are used to generate complex calibration coefficients for each channel. Then these coefficients are applied to the fields of a test object, and the SER with different calibration objects is compared. Dielectric properties for the materials in this calibration study are shown in Table 5.1. The dispersive nature of these materials is described using a two-pole Debye model as follows: (ω) 1 1 σs = ∞ + ∆ + + . 0 1 + jωτ1 1 + jωτ2 jω0 (5.1) Here ∞ , ∆, and σs describe the same Debye parameters as (3.1). The parameters τ1 and τ2 correspond to the two relaxation time constants of the Debye Model. I determined how similar a calibration object and an imaging object should be in order to increase SER, and therefore reduce model error. I investigated channel by channel calibration of the 3D printed rectangular prism, shown in Figure 5.1. The 56 Table 5.1: Two pole Debye parameters for the immersion medium and imaging object materials. Material ∞ ∆ σs τ1 (ps) τ2 (ps) Air 1 0 0 — — Oil immersion 2.35 0.33 0 20 202 3D printed Plastic 2.25 0 0.0036 — — 90% Glycerol 10% water mixture 6.6 17.6 0 187 744 80% Glycerol 20% water mixture 8.0 18.9 0.12 132 541 70% Glycerol 30% water mixture 9.24 23.7 0 58 154 (a) (b) Figure 5.1: Side view (a) and top view (b) of the experimental 3D printed plastic rectangular prism. measured dielectric properties of the 3D printed rectangular prism are shown in Figure 5.2. Several calibration scenarios with this simple phantom were tested. Data was collected with the prism filled with one material and then channel by channel calibration coefficients were calculated using a simulation of the same object. The resulting SER values after using different calibrations scenarios is reported in Table 5.2. 57 15 15 10 10 5 5 0 150 0 150 10 10 10 5 5 0 0 10 0 0 10 15 15 10 10 5 10 5 0 2 0 cm 0 5 0 0 10 10 15 0 10 0.5 1 Figure 5.2: Sagittal, axial, and coronal cross sections of the relative permittivity (left) and effective conductivity (right) of a simulated model of the 3D printed plastic rectangular prism filled with a 90% glycerol 10% water mixture. Table 5.2: SER (4.5) values for calibration of the rectangular prism phantom. Imaging Object Air filled Calibration 90% Glycerol 10% Water Calibration 80% Glycerol 20% Water Calibration 70% Glycerol 30% Water Calibration 90% Glycerol 10% Water 2.4 — 11.6 10.9 80% Glycerol 20% Water 1.5 11.3 — 16.1 70% Glycerol 30% Water 0.7 10.4 15.7 — 58 Calibrating one glycerol/water mixture with another did raise SER significantly. However, the calibrations using air did not perform well since the dielectric properties of air are very different than the glycerol/water mixtures. The relationship between these SER values and image quality was determined using simulated and experimental image reconstructions. A simulated reconstruction of the 90% glycerol 10% water filled rectangular prism (Fig. 5.2) is shown in Fig. 5.3. The experimental reconstruction of the same phantom is shown in Fig. 5.4. Experimentally this phantom was calibrated channel by channel with an 80% glycerol 20% water filled rectangular prism, which resulted in an SER of 11.3 dB. 15 15 10 10 5 5 0 150 0 150 10 10 10 5 5 0 0 10 0 0 10 15 15 10 10 5 10 5 0 2 cm 0 0 0 0 10 5 10 15 0 10 0.5 1 Figure 5.3: Sagittal, axial, and coronal cross sections of the relative permittivity (left) and effective conductivity (right) of a simulated reconstruction of the rectangular prism filled with 90% glycerol 10% water. The experimental imaging performance is worse than simulation due to the presence of model error. The outer edges and the overall dielectric properties are similar between the simulated and experimental reconstructed images. However, the experimental reconstruction cannot localize the bottom of this simple phantom 59 15 15 10 10 5 5 0 150 0 150 10 10 10 5 5 0 0 10 0 0 10 15 15 10 10 5 10 5 0 2 cm 0 0 0 0 10 5 10 15 0 10 0.5 1 Figure 5.4: Sagittal, axial, and coronal cross sections of the relative permittivity (left) and effective conductivity (right) of an experimental reconstruction of the rectangular prism filled with 90% glycerol 10% water that was calibrated with a 80% glycerol 20% water filled rectangular prism. and it has a large error in the dielectric properties in the interior of the phantom due to the level of model error. 5.3 Multiple Object Calibration Channel by channel calibration with a single reference object, as discussed in sections 4.4 and 5.2, does not prevent imaging artifacts due to model error. Consequently, I investigated mitigating the effects of model error using calibration with multiple objects. Experiments were conducted using the 3D printed breast shell, Figures 4.2a and 4.2b, that was filled with varying mixtures of glycerol and water. Those measurements were then used to calibrate the 3D printed breast shell filled with TMMs, shown in Fig. 4.2c. Determining the calibration coefficients involved 60 solving a system of equations for each transmitter as follows: (5.2a) ··· Object N ··· A = GÂ α1 βk ... = βk αL Object 1 Object 1 Object N (5.2b) Here A is an L × N matrix, where L is the number of antennas and N is the number of calibration objects. The ith column of A is the data measured at the L antennas for the ith calibration scenario of a specific transmitter. Similarly Â is an L × N matrix, where the ith column is the simulated data at the L antennas for the ith calibration scenario of the same transmitter. The calibration matrix G is L × L and will correct magnitude and phase errors in bistatic channels (α) as well as mutual coupling errors (β) for each transmitter. In order for a unique solution of (5.2a) to exist, the number of calibration scenarios N must be greater than or equal to the number of sensors L. The optimal calibration matrix G of (5.2a) is the solution with the minimum Frobenius norm as described by Pierre and Kaveh [75] as follows: min kA − GÂkF (5.3) G G = AÂH ÂÂH −1 . (5.4) The Hermitian transpose is represented by H and the matrix inverse is denoted by −1 . For the inverse to exist the matrix ÂÂH must be full rank. If the calibration objects do not lead to a full rank matrix several techniques can approximate the inverse, e.g. using a pseudoinverse. Using the pseudoinverse to approximate the solution to (5.2a) has the potential to produce a high quality calibration with fewer calibration objects than the number of sensors as well. The number of sensors L in the clinical prototype under development is on the order of 32, so collecting that 61 many unique calibration scenarios is impractical if not impossible. Calculating the SER with the calibrated array will demonstrate if this technique will achieve the desired calibration goals. Data was collected using ten different calibration scenarios and then the calibration matrices were estimated using the pseudoinverse. The dielectric properties of the ten materials used to fill the 3D printed breast surface are shown in Table 5.3. Then the SER was computed using the pseudoinverse with different ranks (keeping fewer eigenvalues in the inverse), and the SER after calibration is shown in Table 5.4. Table 5.3: Dielectric properties at 2 GHz for calibration materials Material Relative permittivity Effective conductivity Oil 2.4 0.02 Pure Glycerol 6.1 0.38 85% Glycerol 15% Water 8.2 0.67 75% Glycerol 25% Water 19.3 1.9 60% Glycerol 40% Water 31.3 1.9 50% Glycerol 50% Water 42.8 2.0 40% Glycerol 60% Water 56.3 2.0 30% Glycerol 70% Water 61.0 1.6 15% Glycerol 85% Water 72.1 1.2 Pure Water 75.9 0.83 The baseline SER without any calibration for this object was -23 dB, so the multiple object calibration did raise the SER. Keeping only five eigenvalues resulted in the best SER, but this was still much less than using only the adipose TMM calibration object as presented in Section 4.4. Even a wide range of properties in the multiple object calibration does not improve the SER as much as using an object with similar properties. Therefore, it appears the only satisfactory calibration is when the calibration object is closely matched to the imaging object. Several 62 Table 5.4: SER (4.5) values for a multiple object calibration of the 3D printed breast shell filled with TMMs. Pseudoinverse rank SER (dB) 10 -16.4 9 -15.8 8 -15.3 7 -14.3 6 -14.0 5 -13.3 4 -13.6 3 -14.5 2 -15.2 1 -15.7 calibration scenarios that are close to the true object might improve the results but was not tested in this study. 63 5.4 Model Error Reduction using Resonance Matching I attempted to reduce model error by changing the simulated substrate and immersion medium dielectric properties of the miniaturized patch array (Fig. 4.1) to make the resonant frequencies match between simulation and experiment. The basis for this idea was the assumption that the sources of model error are also the reason the resonances do not match between simulation and experiment. Simulations of a single array panel in air were compared to a measurement of an empty array. The simulated substrate properties (dielectric properties, grid sizes, and thickness) were adjusted to make the simulated resonant frequencies closer to the experimental resonant frequencies. Then simulations with an oil filled array were compared to a measurement of an oil filled array and the simulated dielectric properties of oil were adjusted to make the resonant frequencies match again. Unfortunately, a perfect match was not obtained since changing any parameter made both resonances move. One resonant frequency would be getting closer to matching, while the other resonanct frequency was moving further apart. Also, the idea that the performance of a single panel in air was identical to an air filled array was not proven to be accurate. The single panel simulation was only used to save time and speed up the process. After this entire process the SER did not change dramatically. Using the old substrate and oil properties, the SER was -7 dB with a glycerol filled shell calibrating the TMM phantom of Fig. 4.2. The same scenario had an SER of -8 dB after the properties were adjusted to improve the resonance match. 5.5 Oil and Substrate Imaging I also attempted to use an imaging algorithm to determine the dielectric properties of the substrate and oil in the array shown in Fig. 4.1. This was done using an enhanced version of DBIM that was developed for microwave imaging using a 64 spatial prior by Luz María Neira. This involved gathering experimental data with an array filled with only oil, then the enhanced DBIM algorithm was used to determine the Debye parameters of the substrate and the oil. This required recording electric fields in both the substrate and the oil. The enhanced DBIM algorithm solves the following system of equations: b = H · x. (5.5) Here b is the difference between the measured electric fields and simulated electric fields for the estimated properties of substrate and oil. Similar to regular DBIM, the matrix H is the Hadamard product of the Green’s function and the total electric field with the current estimated properties. In this case the full vectorized fields were used (Ex , Ey , and Ez ). The vector x is the update for the Debye parameters of the substrate and oil as shown below. oil ∞ oil ∆ σoil st x = sub (5.6) ∞ ∆sub σsub st The H matrix is then setup according to the EFIE (2.2). This is broken into two sums, one over voxels containing oil and the other over voxels containing substrate. This results in an overdetermined system of equations where we have 496 equations to solve for six unknowns. To solve for x the pseudoinverse of H is multiplied by b. However, I could not get this algorithm to actually minimize the difference between simulated and experimental fields. Only an incredibly small step size that barely changed the properties would decrease the error and any significant updates to the properties would always increase cost. This most likely occurred because small changes to the substrate or oil properties in front of the antennas would drastically change the performance of the antennas. Consequently, this is a highly nonlinear 65 cost function and will be difficult to descend to a good solution. 5.6 TMM Phantom Imaging This section contains detailed results of the experimental reconstruction of the 3D printed breast shell filled with TMMs in Fig. 4.2, that was calibrated using the 3D printed breast shell filled completely with adipose TMM. In the following plots a simulated reconstruction will be on the left side of the figure, and the experimental reconstruction will be on the right. The initial guess is identical for both algorithms and is shown in Fig. 5.5. Initial Guess 15 Simulated 15 15 10 10 10 10 5 5 5 5 0 0 15 0 0 15 0 0 15 10 10 Experimental 15 0 0 15 10 10 10 10 10 5 5 5 5 0 0 0 0 10 0 0 10 0 0 10 15 15 15 15 10 10 10 10 5 5 5 0 0 0 cm 0 0 2 0 cm 0 5 10 15 0 10 0.5 1 10 5 0 10 10 5 2 0 0 10 10 15 0 10 0.5 1 Figure 5.5: Sagittal, axial, and coronal cross sections of the relative permittivity (left) and effective conductivity (right) of the initial guess for a simulated reconstruction of the TMM breast shell that is compared to an experimental reconstruction of the same phantom. The first iteration of both reconstructions is shown in Fig. 5.6. Both algorithms perform similarly at this point but further iterations of the experimental reconstruction are corrupted by model error whereas the simulated reconstruction keeps improving every iteration. 66 Iteration 1 15 Simulated 15 15 10 10 10 10 5 5 5 5 0 0 15 0 0 15 10 Experimental 15 0 0 15 10 0 0 15 10 10 10 10 10 5 5 5 5 0 0 0 0 10 0 0 10 0 0 10 15 15 15 15 10 10 10 10 5 5 0 0 0 cm 0 5 0 2 0 cm 0 5 10 15 0 10 0.5 1 10 5 0 10 10 5 2 0 0 10 10 15 0 10 0.5 1 Figure 5.6: Sagittal, axial, and coronal cross sections of the relative permittivity (left) and effective conductivity (right) of the first iteration of a simulated reconstruction of the TMM breast shell that is compared to an experimental reconstruction of the same phantom. Figures 5.7 through 5.9 are the subsequent iterations of both algorithms. The simulated reconstruction improves the overall dielectric properties and the shape of the two cylinders, while the experimental reconstruction is adding only artifacts to correct for model errors in the system. A plot of mean squared error (MSE) versus iteration of both algorithms is shown in Fig. 5.10. In simulation each iteration is decreasing MSE while it is improving image quality. Conversely, the experimental reconstruction increases cost on some iterations and barely decreased the MSE back to the cost of the initial guess. 67 Iteration 2 15 Simulated 15 15 10 10 10 10 5 5 5 5 0 0 15 0 0 15 10 Experimental 15 0 0 15 10 0 0 15 10 10 10 10 10 5 5 5 5 0 0 0 0 10 0 0 10 0 0 10 15 15 15 15 10 10 10 10 5 5 0 0 0 cm 0 5 0 2 0 cm 0 5 10 15 0 10 0.5 1 10 5 0 10 10 2 5 0 0 10 10 15 0 10 0.5 1 Figure 5.7: Sagittal, axial, and coronal cross sections of the relative permittivity (left) and effective conductivity (right) of the first iteration of a simulated reconstruction of the TMM breast shell that is compared to an experimental reconstruction of the same phantom. Iteration 3 15 Simulated 15 15 10 10 10 10 5 5 5 5 0 0 15 0 0 15 0 0 15 10 10 Experimental 15 0 0 15 10 10 10 10 10 5 5 5 5 0 0 0 0 0 0 10 10 0 0 10 15 15 15 15 10 10 10 10 5 5 5 0 0 0 cm 0 0 2 0 cm 0 5 10 15 0 10 0.5 1 10 5 0 10 10 5 2 0 0 10 10 15 0 10 0.5 1 Figure 5.8: Sagittal, axial, and coronal cross sections of the relative permittivity (left) and effective conductivity (right) of the first iteration of a simulated reconstruction of the TMM breast shell that is compared to an experimental reconstruction of the same phantom. 68 Iteration 4 15 Simulated 15 15 10 10 10 10 5 5 5 5 0 0 15 0 0 15 10 Experimental 15 0 0 15 10 0 0 15 10 10 10 10 10 5 5 5 5 0 0 0 0 10 0 0 10 0 0 10 15 15 15 15 10 10 10 10 5 5 0 0 0 cm 0 5 0 2 0 cm 0 5 10 15 0 10 0.5 1 10 5 0 10 10 5 2 0 0 10 10 15 0 10 0.5 1 Figure 5.9: Sagittal, axial, and coronal cross sections of the relative permittivity (left) and effective conductivity (right) of the first iteration of a simulated reconstruction of the TMM breast shell that is compared to an experimental reconstruction of the same phantom. 1.4 Normalized MSE 1.2 1 Simulated Experimental 0.8 0.6 0.4 0.2 0 0 1 2 Iteration 3 4 Figure 5.10: Mean-squared error (MSE) vs iteration for reconstructing the TMM filled plastic shell in simulation and experimentally. 69 5.7 Improvements to reduce the effect of model error This section discusses several improvements I have found that will make microwave imaging systems more robust to model error. This includes changes to array design that will limit the sensitivity to common sources of model error. As well as improvements to imaging algorithms that improve performance in the presence of model error. Array design One feature of microwave array design that I suspected of increasing model error was whether or not the array was enclosed with metal walls. Metal walls enclosing an array will increase the time it takes for reflections to dissipate, which is hypothesized to lead to increased error due to modeling the immersion medium. To test this we compared the SER of the array in Fig. 4.1 to variations of the array shown in Fig. 5.11, which is described in Section 4.3. Figure 5.11: Dipole array surrounding a breast phantom. In simulation I introduced error into the immersion medium and calculated the SER for the miniaturized patch array, the dipole array, and the dipole array with 70 metal walls surrounding it. Model error was calculated as the difference between two simulations where one had a 10% error introduced into the immersion medium, and the signal of interest was calculated as the difference between a numerical breast phantom and a homogenous phantom. The SER with these different configurations is shown in Table 5.5. Table 5.5: Array design for reducing model error due to 10% error in immersion medium. Array type SER (dB) Patch array -13 Enclosed dipole array 0.5 Dipole array 12.5 With 10% error in the immersion medium the two enclosed arrays had much worse SER than the open dipole array. The extra reflections for the metal walls appear to increase the amount of model error. Also the more complicated patch array appears to be more sensitive to this type of error and it will potentially have additional sources of model error due to substrate errors that the dipole arrays will not have. Imaging algorithms In this section I analyze the impact of the imaging domain size on the performance of imaging algorithms in the presence of model error. This will be analyzed with the dipole array without metal walls from the last section. Typically the imaging domain is restricted to the interior of the breast, but I analyzed the impact on performance of extending the imaging domain outside the breast. The location and properties of the skin are assumed known, but the imaging domain is extended outside the breast as shown in Fig. 5.12. DBIM was then performed with different imaging domain sizes and 5% error was introduced into the immersion medium. DBIM was allowed to update the parameters anywhere in the different imaging domains, but different constraints 71 (a) (b) Figure 5.12: (a) Breast interior imaging domain of the dipole array. (b) Imaging domain surrounding the ring of antennas. were placed on voxels outside the breast than inside the breast. Voxels outside the breast are not constrained to dielectric properties corresponding to tissue; instead they have relaxed constraints that keep them in a reasonable range of properties for the immersion medium. Table 5.6 contains the reconstruction quality for different imaging domains with 5% error introduced into the immersion medium. Table 5.7 contains the reconstruction quality for different imaging domains with random 2 mm position errors introduced into antenna positions. Table 5.6: Reconstruction quality vs imaging domain size with 5% error in immersion medium. Imaging Domain Reconstruction quality (cosine metric) Breast interior only 0.76 Inside antenna ring 0.79 Surrounding ring of antennas 0.81 Tables 5.6 and 5.7 suggest that increasing the imaging domain size improves the reconstruction quality in the presence of model error. Compared to imaging only the interior of the breast, the imaging performance is slightly improved when 72 Table 5.7: Reconstruction quality vs imaging domain size with antenna position errors. Imaging Domain Reconstruction quality (cosine metric) Breast interior only 0.67 Inside antenna rings 0.70 Surrounding rings of antennas 0.78 the imaging domain is outside the breast but inside of the rings of antennas. However, the imaging performance is drastically improved when the imaging domain surrounds the rings of the antenna. This is due to the imaging domain containing every channel path in the array. When a channel does not pass through the breast it provides information about model error, and the channels that do pass through the breast provide information about the breast interior. When these channels are not included in the imaging domain any errors in these channel paths get reflected back into the imaging domain. The robustness of this algorithm to different types of imaging algorithms is further demonstrated in Fig. 5.13 that also contains antenna position errors unlike the results from the previous chapter. 73 0.82 Cosine Metric 0.8 0.78 0.76 0.74 Noise Immersion Medium Error Antenna Position Error Skin Error Baseline Performance 0.72 0.7 0.68 0 10 20 30 SER (dB) (a) 0.82 Cosine Metric 0.8 0.78 0.76 0.74 Noise Immersion Medium Error Antenna Position Error Skin Error Baseline Performance 0.72 0.7 0.68 0 10 20 30 SER (dB) (b) Figure 5.13: (a) DBIM imaging performance, assessed with the cosine metric, versus SER for different types of model error with a Class 3 phantom. (b) Extended domain DBIM imaging performance versus SER for different types of model error including antenna position errors. The baseline performance of both algorithms is indicated by a dashed line and markers denote 10% error in the dielectric properties of the material. 74 6 6.1 future work and conclusion Future Work The previous chapters demonstrate a high resolution microwave breast imaging algorithm, a quantitative way to assess the impact of model error on microwave imaging, and several techniques to mitigate the effects of model error on experimental microwave imaging performance. However, measurement techniques and manufacturer specifications are not precise enough to sufficiently characterize the dielectric materials in an array so the biggest challenge facing microwave breast imaging is still model error. Future research in microwave imaging should focus on experimental validation of the extended domain DBIM imaging results, reducing the impact of model error on imaging performance by investigating array designs and imaging algorithms that are not sensitive to potential model error sources, and investigating a hybrid extended domain level set imaging algorithm. Experimental exploration of a dipole array This work has presented an array that is less sensitive to common types of model error as well as an imaging algorithm that will perform better in the presence of model error. Section 5.7 demonstrated that an open dipole array is not as sensitive to errors in the immersion medium as an enclosed array. And Sections 4.4 and 5.7 demonstrated that an imaging domain that surrounds the antennas will be robust to model error. These simulation should be extended to include a realistic dipole array. This includes adding a support structure for the dipole antennas as well as feed lines for these dipoles. A 3D printed box could be used to support the necessary feed lines that run from outside the array and then to each antenna in the cylindrical dipole array. Also, simulations should be done to verify whether or not dipoles can operate in air, which will eliminate model error due to immersion medium. These simulations should also investigate whether or not errors in the positioning of the 75 dipole feeds can affect the performance of the extended domain DBIM imaging algorithm. If these simulations still suggest that the extended domain DBIM imaging algorithm is still robust to any types of model error, these results should be validated experimentally. The same investigations of calibration and imaging performance with different phantoms could be completed as was done with the miniaturized patch array. If this array and algorithm performs as well as expected it could be directly applied to clinical studies and demonstrate superior performance than any other microwave imaging studies. The SER metric to improve array design and imaging algorithms As demonstrated in Section 5.7 the SER metric can be used to improve array design and imaging algorithms. Array features, like metal enclosing walls, that increase the sensitivity to different types of model error can be identified and limit the amount of model error present in experimental imaging. Potential features might include antenna feed models, effective width of metal in FDTD, graded meshing to include small features, or errors in the physical dimensions of the array. Other imaging algorithm improvements that reduce the effects of model error on microwave imaging can also be investigated with the SER metric. Algorithms like the contrast source inversion method [69] could be evaluated with the SER metric similar to the way the extended domain DBIM technique was evaluated in Section 5.7. Other changes to imaging algorithms, such as regularization methods or other optimization techniques that utilize the correlated nature of model error, can be evaluated with the SER metric to assess the potential improvement of experimental imaging. Extended domain level set imaging The structured nature and high levels of model error will affect the performance of any imaging algorithm that is not specifically designed to perform well in its presence, including the level set method presented in Chapter 3. Consequently, the SER 76 metric could be used to determine how much the high resolution reconstructions and accurate density estimates will be degraded by the presence of model error. Since the level set method is similar to DBIM, in terms of imaging domain setup and the use of the Jacobian for updates, it will most likely be as sensitive to model error as DBIM. The high resolution imaging algorithm of Chapter 3 and the algorithm robust to model error presented in Chapter 4 could be combined to produce high resolution microwave images. The imaging domain of the level set method could be extended to the exterior of the breast to ensure that it is robust to model error as was shown with DBIM. The exterior properties could be designated as another level set that does not change shape but the dielectric properties can vary to account for model error. This method will then provide the same higher resolution reconstruction of the breast, while allowing the exterior properties to change and reduce the effects of model error experimentally. 6.2 Conclusion Microwave imaging is potentially a low cost, non-ionizing, and widely available imaging modality for breast cancer risk assessment, treatment monitoring, and tumor identification. However, conventional microwave imaging methods have significant limitations in terms of resolution and experimental performance. Common microwave imaging methods have moderate-to-low resolution that distorts boundaries between different tissue types. Also, any errors in the simulated model of the array setup can corrupt the iterative process of these algorithms and degrade experimental imaging performance. This document described several techniques to address these concerns and improve the performance of microwave breast imaging and prepare it for the clinical application stage. Chapter 3 presented a computationally tractable 3D microwave breast imaging algorithm based on level sets and assessed its efficacy for microwave imaging by investigating its performance with evaluating breast density. The density estimates from my level set algorithm are more accurate than those of conventional microwave 77 imaging, and the accuracy is greater than that reported for mammographic density estimation. These results demonstrate the feasibility of high resolution 3D microwave breast imaging using a level set method. Next, chapter 4 described a quantitative metric (SER) for determining the amount of model error when modeling an experimental array and predicting the impact model error will have on microwave breast imaging. I demonstrated the use of SER to evaluate experimental imaging performance, experimental calibration, and the performance of different imaging algorithms in the presence of model error. Model error was demonstrated to significantly degrade experimental imaging performance, and calibration was shown to be too dependent on the calibration object to remove a satisfactory amount of model error. However, the SER metric was used to demonstrate an imaging algorithm that is robust to model error. These results demonstrate the efficacy of using SER for improving the performance of microwave imaging in a clinical environment. Chapter 5 further explored the use of the SER metric to improve experimental microwave imaging performance. The reduced imaging performance when simulating an experimental patch array was presented with several different imaging objects. The ineffectiveness of many calibration techniques to reduce the amount of model error was also demonstrated. Enclosed array designs were shown to be more sensitive to immersion medium errors than open array designs. And the extended domain DBIM imaging algorithm was shown to be robust to many types of model error as well as large amounts of model error. Future work involves investigating techniques to mitigate the effects of model error on experimental imaging performance. Experimental validation of the extended domain DBIM imaging algorithm will determine if this technique will perform well in a clinical setting. The SER metric could also be used to design arrays and imaging algorithms that are not sensitive to common sources of model error. The level set method and extended domain DBIM algorithms can be also combined to provide high resolution reconstructions that are robust to model error and improve upon any existing experimental or clinical microwave imaging results. 78 bibliography [1] 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 large-scale study of the ultrawideband microwave dielectric properties of normal breast tissue obtained from reduction surgeries,” Phys. Med. Biol., vol. 52, pp. 2637–2656, April 2007. [2] M. Lazebnik, D. Popovic, L. McCartney, C. Watkins, M. J. Lindstrom, J. Harter, S. Sewall, T. Ogilvie, A. Magliocco, T. 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. [3] UWCEM numerical breast phantom repository. [Online]. Available: http: //uwcem.ece.wisc.edu [4] J. D. Shea, P. Kosmas, S. C. Hagness, and B. D. Van Veen, “Three-dimensional microwave imaging of realistic numerical breast phantoms via a multiplefrequency inverse scattering technique,” Med. Phys., vol. 37, no. 8, pp. 4210– 4226, August 2010. [5] N. Irishina, D. Alvarez, O. Dorn, and M. Moscoso, “Structural level set inversion for microwave breast screening,” Inverse Probl., vol. 26, 035015, 2010. [6] “Cancer Facts & Figures 2013,” American Cancer Society, Inc. [7] K. Kerlikowske, L. Ichikawa, D. L. Miglioretti, D. S. M. Buist, P. M. Vacek, R. Smith-Bindman, B. Yankaska, P. A. Carney, and R. Ballard-Barbash, “Longitudinal measurement of clinical mammographic breast density to improve estimation of breast cancer risk,” J. Natl. Cancer Inst., vol. 99, no. 5, pp. 386–395, 2007. 79 [8] National Institutes of Health State-of-the-Science Conference Statement: Diagnosis and Management of Ductal Carcinoma In Situ, vol. 102, no. 3, Semptember 22-24 2009. [9] M. Baum, “Harms from breast cancer screening outweigh benefits if death caused by treatment is included,” BMJ, vol. 346, pp. 385–386, 2013. [10] S. van Engeland, P. R. Snoeren, H. Huisman, C. Boetes, and N. Karssemeijer, “Volumetric breast density estimation from full-field digital mammograms,” IEEE Transactions on Medical Imaging, vol. 25, no. 3, pp. 273–282, March 2006. [11] O. Pawluczyk, B. J. Augustine, M. J. Yaffe, D. Rico, J. Yang, and G. E. Mawdsley, “A volumetric method for estimation of breast density on digitized screen-film mammograms,” Med. Phys., vol. 30, no. 3, pp. 352–364, March 2003. [12] N. H. Peters, I. H. Borel Rinkes, N. P. Zuithoff, W. P. Mali, K. G. Moons, and P. H. Peeters, “Meta-analysis of MR imaging in the diagnosis of breast lesions,” Radiology, vol. 246, no. 1, pp. 116–24, 2008, epub 2007/11/21. doi: 10.1148/radiol.2461061298. PubMed PMID: 18024435. [13] A. Mashal, B. Sitharaman, X. Li, P. K. 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 tissue-mimicking materials,” IEEE Transactions on Biomedical Engineering, vol. 57, no. 8, pp. 1831–1834, August 2010. [14] G. Bellizzi, O. M. Bucci, and I. Catapano, “Microwave cancer imaging exploiting magnetic nanoparticles as contrast agent,” IEEE Transactions on Biomedical Engineering, vol. 58, no. 9, pp. 2528–2536, September 2011. [15] T. Rubaek, P. M. Meaney, P. Meincke, and K. Paulsen, “Nonlinear microwave imaging for breast-cancer screening using Gauss-Newton’s method and the CGLS inversion algorithm,” IEEE Transactions on Antennas and Propagation, vol. 55, no. 8, pp. 2320–2331, August 2007. 80 [16] W. C. Chew, Waves and Fields in Inhomogeneous Media. Press, 1995. Piscataway, NJ: IEEE [17] W. T. Joines, Y. Z. Dhenxing, and R. L. Jirtle, “The measured electrical properties of normal and malignant human tissues from 50 to 900 MHz,” Med. Phys., vol. 21, pp. 547–550, April 1994. [18] S. S. Chaudhary, R. K. Mishra, A. Swarup, and J. M. Thomas, “Dielectric properties of normal and malignant human breast tissues at radiowave and microwave frequencies,” Indian J. Biochem. Biophys., vol. 21, pp. 7679–7679, 1984. [19] R. J. Halter, T. Zhou, P. M. Meaney, A. Hartov, R. J. Barth Jr., K. M. Rosenkranz, W. A. Wells, C. A. Kogel, A. Borsic, E. J. Rizzo, and K. D. Paulsen, “The correlation of in vivo and ex vivo tissue dielectric properties to validate electromagnetic breast imaging: initial clinical experience,” Physiol. Meas., vol. 30, pp. 121–136, 2009. [20] E. Fear and M. Stuchly, “Microwave detection of breast cancer,” IEEE Transactions on Microwave Theory and Techniques, vol. 48, no. 11, pp. 1854–1863, November 2000. [21] E. Fear, X. Li, S. Hagness, and M. Stuchly, “Confocal microwave imaging for breast cancer detection: Localization of tumors in three dimensions,” IEEE Transactions on Biomedical Engineering, vol. 49, no. 8, pp. 812–822, August 2002. [22] E. Bond, S. Hagness, and B. Van Veen, “Microwave imaging via space-time beamforming for early detection of breast cancer,” IEEE Transactions on Antennas and Propagation, vol. 51, no. 8, pp. 1690–1705, August 2003. [23] X. Li, S. Davis, S. Hagness, D. van der Weide, and B. Van Veen, “Microwave imaging via space-time beamforming: Experimental investigation of tumor detection in multilayer breast phantoms,” IEEE Transactions on Microwave Theory and Techniques, vol. 52, no. 8, pp. 1856–1865, August 2004. 81 [24] S. Davis, H. Tandradinata, S. Hagness, and B. Van Veen, “Ultraphantoms microwave breast cancer detection: A detection-theoretic approach using the generalized likelihood ratio test,” IEEE Transactions on Biomedical Engineering, vol. 52, no. 7, pp. 1237–1250, July 2005. [25] P. Kosmas and C. Rappaport, “Time reversal with the FDTD method for microwave breast cancer detecction,” IEEE Transactions on Microwave Theory and Techniques, vol. 53, no. 7, pp. 2317–2323, July 2005. [26] 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 Transactions on Microwave Theory and Techniques, vol. 48, no. 11, pp. 1841–1853, November 2000. [27] A. Abubakar, P. van den Berg, and J. Mallorqui, “Imaging of biomedical data using a multiplicative regularized contrast source inversion method,” IEEE Journal of Oceanic Engineering, vol. 50, no. 7, pp. 1761–1771, July 2002. [28] A. Bulyshev, A. Souvorov, S. Semenov, V. Posukh, and Y. Sizov, “Threedimensional vector microwave tomography: Theory and computational experiments,” Inverse Probl., vol. 20, pp. 1239–1259, 2004. [29] Q. Fang, P. S. Geimer, A. Streltsov, and K. Paulsen, “Microwave Image reconstruction from 3-D fields coupled to 2-D parameter estimation,” IEEE Transactions on Medical Imaging, vol. 23, no. 4, pp. 475–484, April 2004. [30] Z. Zhang and Q. Liu, “Three-dimensional nonlinear image reconstruction for microwave biomedical imaging,” IEEE Transactions on Biomedical Engineering, vol. 51, no. 3, pp. 544–548, March 2004. [31] P. Meaney, Q. Fang, T. Rubaek, E. Demidenko, and K. Paulsen, “Log transformation benefits parameter estimation in microwave tomographic imaging,” Med. Phys., vol. 34, pp. 2014–2023, 2007. 82 [32] A. Franchois and C. Pichot, “Microwave imaging - Complex permittivity reconstruction with a Levenberg-Marquardt method,” IEEE Transactions on Antennas and Propagation, vol. 2, no. 2, pp. 203–215, February 1997. [33] M. Haynes, S. Clarkson, and M. Moghaddam, “Electromagnetic inverse scattering algorithm and experiment using absolute source characterization,” IEEE Transactions on Antennas and Propagation, vol. 60, no. 4, pp. 1854–1867, April 2012. [34] W. C. Chew and Y. M. Wang, “Reconstruction of two-dimensional permittivity distribution using the distorted born iterative method,” IEEE Transactions on Medical Imaging, vol. 9, pp. 218–225, 1990. [35] P. M. van den Berg and R. E. Kleinman, “Contrast source inversion method,” Inverse Probl., vol. 13, pp. 1607–1620, 1997. [36] C. Gilmore, P. Mojabi, and J. LoVetri, “Comparison of an enhanced distorted born iterative method and the multiplicative-regularized contrast source inversion method,” IEEE Transactions on Antennas and Propagation, vol. 57, no. 8, pp. 2341–2351, August 2009. [37] J. Nocedal and S. Wright, Nonlinear Optimization, 2nd ed. Springer, New York, 2006. [38] R. F. Remis and P. M. van den Derg, “On the equivalence of the NewtonKantorovich and distorted Born methods,” Inverse Probl., vol. 16, no. 1, pp. L1–L4, February 2000. [39] D. Winters, J. Shea, E. Madsen, G. Frank, B. D. Van Veen, and S. C. Hagness, “Detection the breast surface using UWB microwave monostatic backscatter measurements,” IEEE Transactions on Biomedical Engineering, vol. 55, no. 1, pp. 247–256, 2008. 83 [40] T. Williams, J. Sill, and E. Fear, “Breast surface estimation for radar-based breast imaging systems,” IEEE Transactions on Biomedical Engineering, vol. 55, no. 6, pp. 1678–1686, 2008. [41] T. L. Pope Jr., M. E. Read, T. Medsker, A. J. Buschi, and A. N. Brenbridge, “Breast skin thickness: normal range and causes of thickening shown on film-screen mammography,” J. Can. Assoc. Radiol., vol. 35, no. 4, pp. 365–368, December 1984. [42] S. Gabriel, R. W. Lau, and C. Gabriel, “The dielectric properties of biological tissue: II. Measurements on the frequency range 10 Hz to 20 GHz.” Phys. Med. Biol., vol. 41, no. 11, pp. 2251–2269, November 1996. [43] ——, “The dielectric properties of biological tissue: III. Parameteric models for the dielectric spectrum of tissues,” Phys. Med. Biol., vol. 41, no. 11, pp. 2271–2293, Novemeber 1996. [44] F. Santosa, “A level-set approach for inverse problems involving obstacles,” ESAIM: Control, Optimisation and Calculus of Variations, vol. 1, pp. 17–33, January 1996. [45] R. Ferrayé, J.-Y. Dauvignac, and C. Pichot, “An inverse scattering method based on contour deformations by means of a level set method using frequency hopping technique,” IEEE Transactions on Antennas and Propagation, vol. 51, no. 5, pp. 1100–1113, May 2003. [46] M. R. Hajihashemi and M. El-Shenawee, “Shape reconstruction using the level set method for microwave applications,” IEEE Antennas and Wireless Propagation Letters, vol. 7, pp. 92–96, 2008. [47] M. Benedetti, D. Lesselier, M. Lambert, and A. Massa, “A multi-resolution technique based on shape optimization for the reconstruction of homogeneous dielectric objects,” Inverse Probl., vol. 25, p. 015009, 2009. 84 [48] M. El-Shenawee, O. Dorn, and M. Moscoso, “An adjoint-field technique for shape reconstruction of 3-D penetrable object immersed in lossy medium,” IEEE Transactions on Antennas and Propagation, vol. 57, pp. 520–534, 2009. [49] D. A. Woten, M. R. Hajihashemi, A. M. Hassan, and M. El-Shenawee, “Experimental microwave validation of level set reconstruction algorithm,” IEEE Transactions on Antennas and Propagation, vol. 58, no. 1, pp. 230–233, January 2010. [50] O. Dorn, H. Bertete-Aguirre, J. G. Berryman, and G. C. Papanicolaou, “A nonlinear inversion method for 3D electromagnetic imaging using adjoint fields,” Inverse Probl., vol. 15, pp. 1523–1558, 1999. [51] O. Dorn, E. L. Miller, and C. M. Rappaport, “A shape reconstruction method for electromagnetic tomography using adjoint fields and level sets,” Inverse Probl., vol. 16, pp. 1119–1156, 2000. [52] O. Dorn and D. Lesselier, “Level set methods for inverse scattering,” Inverse Probl., vol. 22, pp. R67–R131, 2006. [53] N. Irishina, O. Dorn, and M. Moscoso, “A level set evolution strategy in microwave imaging for early breast cancer detection,” Computers and Mathematics with Applications, vol. 56, pp. 607–618, 2008. [54] N. Irishina, D. Alvarez, O. Dorn, and M. Moscoso, “Detecting and imaging dielectric objects from real data: A shape-based approach,” Mathematical and Computer Modelling, vol. 50, pp. 743–749, 2009. [55] M. J. Burfeindt, “Methods for mitigating the effect of noise, interference, and model error on microwave breast imaging,” Ph.D. dissertation, University of Wisconsin-Madison, 2013. [56] M. R. Hajihashemi and M. El-Shenawee, “Level Set Algorithm for Shape Reconstruction of Non-Overlapping Three-Dimensional Penetrable Targets,” IEEE 85 Transactions on Geoscience and Remote Sensing, vol. 50, no. 1, pp. 75–86, January 2012. [57] 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 Transactions on Biomedical Engineering, vol. 55, no. 12, pp. 2792–2800, December 2008. [58] J. A. Harvey and V. E. Bovbjerg, “Quantitative assesment of mammographic breast density: Relationship with breast cancer risk,” Radiology, vol. 230, pp. 29–41, January 2004. [59] T. J. Colgan, S. C. Hagness, and B. D. V. Veen, “3-D Microwave Imaging using a Level Set Method for Breast Density Evaluation,” presented at USNC-URSI National Radio Science Meeting, Boulder, CO, January 2013. [60] 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,” IEEE Microwave and Wireless Components Letters, vol. 17, no. 12, pp. 822–824, December 2007. [61] T. J. Cui, Y. Qin, G.-L. Wang, and W. C. Chew, “Low-frequency detection of twodimensional buried objects using high-order extended Born approximations,” Inverse Probl., vol. 20, pp. S41–S62, 2004. [62] K. T. Ladas and A. J. Devaney, “Generalized ART algorithm for diffraction tomography,” Inverse Probl., vol. 7, no. 1, pp. 109–125, 1991. [63] A. Papoulis and S. U. Pillai, Probablity, Random Variable, and Stochastic Processes, 4th ed. Boston: McGraw-Hill, 2002. [64] P. C. Hansen, Rank-deficient and Discrete Ill-posed Problems: Numerical Aspects of Linear Inversion. 3600 Market Street, Floor 6, Philadelphia, PA 19104: Philadelphia, Pa: Society for Industrial and Applied Mathematics, 1998. 86 [65] D. Kurrant and E. Fear, “Regional estimation of the dielectric properties of inhomogeneous objects using near-field reflection data,” Inverse Probl., vol. 28, 075001 (27pp), 2012. [66] S. M. Aguilar, M. A. Al-Joumayly, M. J. Burfeindt, N. Behdad, and S. C. Hagness, “Multiband Miniaturized Patch Antennas for a Compact, Shielded Microwave Breast Imaging Array,” IEEE Transactions on Antennas and Propagation, vol. 62, no. 3, pp. 1221–1231, March 2014. [67] N. R. Epstein, P. M. Meaney, and K. D. Paulsen, “Microwave Tomographic Imaging Utilzing Low-Profile, Rotating, Right Angle-Bent Monopole Antennas,” International Journal of Antennas and Propagation, vol. 2014, 2014. [68] M. Burfeindt, J. Shea, B. V. Veen, and S. Hagness, “Beamforming-enhanced inverse scattering for microwave breast imaging,” IEEE Transactions on Antennas and Propagation, vol. 62, no. 10, pp. 5126–5132, October 2014. [69] P. M. van den Berg and A. Abubakar, “Contrast source inversion method: state of art,” Progress In Electromagnetics Research, vol. 34, pp. 189–218, 2001. [70] J. LoVetri, P. Mojabi, A. Zakaria, M. Ostadrahimi, and I. Jeffrey, “System and formulation options for biomedical microwave imaging,” in General Assembly and Scientific Symposium (URSI GASS), 2014 XXXIth URSI, 2014. [71] C. Gilmore, A. Zakaria, J. LoVetri, and S. Pistorius, “A study of matching fluid loss in a biomedical microwave tomography system,” Medical Physics, vol. 40, no. 2, p. 023101, 2013. [72] M. J. Burfeindt, N. Behdad, B. D. V. Veen, and S. C. Hagness, “Quantitative Microwave Imaging of Realistic Numerical Breast Phantoms Using and Enclosed Array of Multiband, Miniaturized Patch Antennas,” IEEE Antennas and Wireless Propagation Letters, vol. 11, pp. 1626–1629, 2012. [73] M. J. Burfeindt, T. J. Colgan, R. O. Mays, J. D. Shea, N. Behdad, B. D. V. Veen, and S. C. Hagness, “MRI-Derived 3-D-Printed Breast Phantom for Microwave 87 Breast Imaging Validation,” IEEE Antennas and Wireless Propagation Letters, vol. 11, pp. 1610–1613, 2012. [74] M. Lazebnik, E. L. Madsen, G. R. Frank, and S. C. Hagness, “Tissue mimicking phantom materials for narrowband and ultrawideband microwave applications,” Physics in Medicine and Biology, vol. 50, pp. 4245–4258, 2005. [75] J. W. Pierre and M. Kaveh, “Experimental evaluation of high-resolution direction-finding algorithms using a calibrated sensor array testbed,” Digital Signal Processing, vol. 5, pp. 243–254, 1995.

1/--страниц