Забыли?

?

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

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

###### Документ
Категория
Без категории
Просмотров
0
Размер файла
3 856 Кб
Теги
sdewsdweddes
1/--страниц
Пожаловаться на содержимое документа