close

Вход

Забыли?

вход по аккаунту

?

Contributions to cost reduction and sensitivity improvement of microwave breast cancer detection

код для вставкиСкачать
CONTRIBUTIONS TO COST REDUCTION AND SENSITIVITY IMPROVEMENT OF
MICROWAVE BREAST CANCER DETECTION
by
Min Zhao
A dissertation submitted in partial fulfillment of
the requirements for the degree of
Doctor of Philosophy
(Electrical Engineering)
at the
UNIVERSITY OF WISCONSIN-MADISON
2009
UMI Number: 3400124
All rights reserved
INFORMATION TO ALL USERS
The quality of this reproduction is dependent upon the quality of the copy submitted.
In the unlikely event that the author did not send a complete manuscript
and there are missing pages, these will be noted. Also, if material had to be removed,
a note will indicate the deletion.
UMT
Dissertation Publishing
UMI 3400124
Copyright 2010 by ProQuest LLC.
All rights reserved. This edition of the work is protected against
unauthorized copying under Title 17, United States Code.
ProQuest LLC
789 East Eisenhower Parkway
P.O. Box 1346
Ann Arbor, Ml 48106-1346
o
ID
W
O
3"
O
u
Q.
C
o
3
3
O
O
o
u
?
c
3
ai
**•
(A
IQ
*
3
a
IB
o
o
3
5'
a
(A
(0
*
*
*
*
*
*
**
*
*
*
*
ft
»
»
t
**
*
»
*
**
s J
3
CO
IQ
3
01
<
—
o
•a
>
Committee's Page. This page is not to be hand-written except for the signatures
(Q
3
5"
• • *
01
3
5*
CQ
C
(A
>
c
K)
O
O
•O
c
</>
3
w. >
*<
c
o
O
(D
3
o-
<D
o
fl>
a
a.
CD
Q.
at
2L
O
SL
m
x
(0
3
O
o
I
<v
-<
a
u
B»
2
o
Q
O
D
M
a-
3
(D
3
CD
CD
>< en
3
CD
•g3" I~
OT
O
- * CD
If
-»
8" f
D "
o • 9.
O
-*
a.
3 5
0)
3-
CD
CD
O
s
3"
CD
3"
o
o
o
CO
0)
c
0)
a.
Q
5 o
C
O
~
ccr
C 3
3
3:
<' QCD.
in
<
—^
o
o
CO
O
C
D
Q
D
CO
CD
CO
D
o"
CD | i
a
—-s
D
o a
0
o
Q
8 CV
CD
^
CO
o o
o O
o
^
!§
c
3 cr
CD
CD
O
3
a
CD
<D
3
3
o
0)
it*
CD
Q.
Committee's Page. This page is not to be hand-written except for the signatures
1
To my parents Qin and Zhao.
11
ACKNOWLEDGMENTS
It was quite a journey pursuing the Ph.D., and I would not have been able to complete this
journey without the support of many people over the past several years.
I must first express my gratitude towards my advisors, Professor Susan Hagness and Professor
Daniel van der Weide. They not only helped me keep growing technically and professionally, but
their leadership, vision, and strategic thinking also set an example I hope to match some day.
I wish to thank those faculties and technical staff at UW for their consulting and assistance in
technical areas: Alan Bettermann, Robert Blick, Timothy Hall, Steve Limbach, Robert Marsland,
Charles Paulson, Barry Van Veen, and Tomy Varghese.
I would like to thank the graduate students I have worked with on different projects in UWCEM
and vdW groups: Min-Ki Choi, Anuj Dron, Hong-Joon Kim, Al Mashal, and Jake Shea. I look
forward to future collaboration with any of them.
I also thank some of my fellow Ph.D. students: Suzette Aguilar, Matthew Burfemdt, Zhen Ji,
Steve Kennedy, Mariya Lazebnik, Keely Willis, and Earl Zastrow. They each helped make my
time in the Ph.D. program more fun and interesting.
Last but not least, millions of thanks to my parents Qin and Zhao who have been supporting
every decision I made over the years.
iii
TABLE OF CONTENTS
Page
LIST OF TABLES
v
LIST OF FIGURES
vi
ABSTRACT
1
2
1
1.1
1.2
1.3
1.4
1.5
1
2
3
4
5
Background of Current Breast Cancer Detection Techniques
Overview of Active Microwave Imaging Techniques for Breast Cancer Detection .
Motivation to Develop a Low-Cost Microwave Refiectometry
Motivation to Investigate an Acoustic and Microwave Hybrid Modality .
Contributions
Ultrawideband Compact Mixer-based Reflctometer (CMR)
Background Knowledge
Theory and Design
Calibration
Fabrication and Measurement
Results
.
Motivation and Background
Free-Space CMR-Antenna Measurement System
Free-Space Calibration Method
Results
Concept of Acoustic and Microwave Hybrid Modality
4.1
4.2
7
.
Calibrated Free-Space Measurements with A CMR-Antenna Subsystem
3.1
3.2
3.3
3.4
4
xi
Introduction
2.1
2.2
2.3
2.4
2.5
3
.:....
Theory
. .
Methods of Inducing Tissue Deformation
7
8
10
11
12
15
15
16
16
20
26
26
27
IV
Page
5
Global-Acoustic and Microwave Hybrid Modality
5.1
5.2
5.3
5.4
6
7
8
9
2-D Numerical Simulation
Discussion of Limitations
3-D Numerical Simulation with a Bulk Translation Removal Method
Discussion of Limitations
29
. 29
37
37
46
Focused-Acoustic and Microwave Hybrid Modality
6.1
Local Excitation with a Focused Acoustic Transducer
6.2
Overview of the 3-D Numerical Techniques
49
. 49
49
3-D Simulation of Focused-Acoustic Induced Tissue Deformation
52
7.1
Acoustic Beamformer Design and Modeling
52
7.2
Induced Deformation Results
54
3-D Simulation of Microwave Interaction with a Deforming Tissue Target
8.1 Derivation of 3-D SBCs
8.2 3-D FDTD-SBC Implementation
8.3 Validation
3-D Computational Investigation of Focused-Acoustic and Microwave Hybrid Modality
......
9.1 Results - Microwave Scattering from a Canonically Deforming Target
9.2 Results - Microwave Scattering from a Focused-Acoustic Induced Deforming Target
9.3 Results - Impact of Dielectric and Elastic Properties Variations on the Microwave
Doppler Scattering
-....'
60
60
66
68
74
74
75
78
10 Conclusion and Future Work
83
LIST OF REFERENCES . . . . . . .
86
APPENDICES
Appendix A:
Appendix B:
Appendix C:
Theoretical Analysis of the Doppler Reconstruction Method
List of Authored Publications
Developed Codes for FDTD-SBC Implementation
93
96
98
V
LIST OF TABLES
Table
5.1
Page
Introduced errors in calculating the bulk translation
/
47
VI
LIST OF FIGURES
Figure
Page
2.1
Signal flow analysis of a reflectometer
7
2.2
Mixer-based 6dB bridge-T attenuator circuit
8
2.3
Mixer-based reflectometer system with two external signal sources
9
2.4
A prototype of the cascaded 6dB bridge-T attenuators
2.5
(a) Magnitude and (b) Phase of the reflection coefficients of the DUTs measured by
the reflectometer prototype as well as by a commercial PNA
13
2.6
Simulated and measured directivity and coupling factor for the reflectometer prototype
14
3.1
CMR-antenna system: The compact CMR and a UWB horn antenna are connected
through a semi-rigid cable
16
Prototype: The compact CMR and a UWB horn antenna are connected through a
semi-rigid cable
17
3.2
. .
10
3.3
The controlled experiment setup used in the free-space calibration and measurement
procedure
. 19
3.4
Raw time-domain responses from the standards and MUT
21
3.5
Time-gated and calibrated time-domain responses from the standards and MUT . . . .
22
3.6
Measured and calibrated magnitude of complex reflection coefficients of a silicone
rubber sheet (er = 3, thickness = 1.2cm). The measurement results obtained from
CMR-antenna system are compared with those obtained using commercial PNA-antenna
setup as well as with the theoretical values
22
Vll
Figure
3.7
Page
Measured and calibrated phase of complex reflection coefficients of a silicone rubber
sheet (er = 3, thickness = 1.2cm). The measurement results obtained from CMRantenna system are compared with those obtained using commercial PNA-antenna
setup as well as with the theoretical values.
23
Measured reflection coefficients of (a) a Silicone rubber sheet (er = 3, thickness =
1.2cm) and (b) a Skycast material (er = 12, thickness = 2.6cm) using the free-space
CMR-antenna prototype. Results (circles) compare favorably with the theory values
(dotted lines)
24
Measured reflection coefficients of (c) an FR4 board (er = 4.34, thickness = 0.2cm)
and (d) a Plexiglas material (er = 2.7, thickness = 2.4cm) using the free-space CMRantenna prototype. Results (circles) compare favorably with the theoretical values
(dotted lines)
25
Potential ways of inducing tissue deformation in breast: (a) global excitation with a
compression plate, and (b) local excitation with a focused acoustic transducer
28
Comparison of the amplitudes of the normalized Doppler component computed via
FDTD-SBCs (dashed line) and the normalized difference between stationary scattered
field solutions computed via FDTD-REP (solid line). The 2D scattering object is a
deforming 4-mm-diameter circular cylinder undergoing a peak radial displacement of
0.45 mm
33
Phantom configuration for the hybrid simulations, with axis of cylindrical inclusion
into the page .
33
Deformed inclusion boundaries, under a 10% peak applied strain (Bulk translation
is removed). Peak deformation of boundary point (0, 2) is 0.24 mm for the normal
fibroglandular inclusion and 0.07 mm for the malignant tumor inclusion
35
Computed Doppler components in the backward (0 = 180°) and forward (</> = 0°)
scattered fields under different peak applied strains, for malignant and normal fibroglandular scatterers
36
3-D numerical phantom testbed for the (a) mechanical simulation and (b) electromagnetic simulation of the global-acoustic and microwave hybrid modality
39
5.6
An illustration of the 'null' detection procedure
41
5.7
An illustration of the bulk translation calculation procedure
41
3.8
3.9
4.1
5.1
5.2
5.3
5.4
5.5
Vlll
Appendix
Figure
5.8
5.9
Calibrated elastic boundary deformation of the malignant tumor and the normal fibroglandular tissue clutter. Bulk translation for both cases are removed with a synthetic tracking method through the use of an antenna array. The normal fibroglandular
tissue clutter experiences about three times the boundary deformation of the malignant
tumor.
7.1
7.2
43
The Doppler scattering response from the deforming malignant tumor and deforming
normal fibroglandular tissues. Despite the dielectric properties variation in the normal
fibroglandular tissue, the Doppler scattered field from the deformed malignant tumor
is observed to be 5-10 dB weaker than that from the deformed normal fibroglandular
tissue clutter
. 44
5.10 The Doppler scattering response from a deforming malignant tumor and a deforming
normal fibroglandular tissue clutter, when different levels of errors in approximating
the bulk translation are introduced in the testbed
6.1
Page
45
Computational steps of simulating microwave scattering response from a deforming
target under focused-acoustic excitation
51
Acoustic intensity (W/cm2) in the (a) lateral and (b) axial focal plane generated by a
focused square transducer in a normal fibroglandular tissue background: The achieved
focal region is roughly 1 cm x 1 cm x 3.5 cm
55
3-D computational testbed for simulating the microwave Doppler scattering response
from a tissue target excited by a focused acoustic transducer. The target is either
a malignant tumor or a normal fibroglandular tissue clutter, surrounded by normal
breast tissues
56
7.3
Focused-acoustic induced time-varying radial deformation of a tumor boundary, ipiaterai
is the angle with respect to the y-axis in the y-z plane and ipaxiai is the angle with
respect to the x-axis in the x-z plane. The group of solid-line curves represents deformation for different snapshots on the acoustic time scale. Its largest boundary
deformation is 0.07 mm
57
7.4
Focused-acoustic induced time-varying radial deformations of a normal fibroglandular
tissue clutter boundary, ipiaterai is the angle with respect to the y-axis in the y-z plane
and ipaxiai is the angle with respect to the x-axis in the x-z plane. The group of solidline curves represents deformation for different snapshots on the acoustic time scale.
Its largest boundary deformation is 0.23 mm
58
ix
Appendix
Figure
8.1
8.2
8.3
8.4
8.5
8.6
8.7
9.1
-..
Page
Sheet impedances Zs and admittances Ys for a 3-D dielectric object are derived by
enforcing equivalence in the scattered field solutions for a) a perturbed (r = a + 5)
dielectric sphere, and b) an unperturbed (r = a) dielectric sphere with SBCs applied
at the boundaries.
61
3-D numerical testbed for validating the FDTD-SBC implementation: simulate the
microwave scattering response from a deforming malignant tumor surrounded by a
normal fibroglandular tissue background
69
The computed total scattered fields from a slightly bloated spherical tumor with a
uniform radial perturbation across its boundary, ip is the angle with respect to the
y-axis in the y-z plane, ip = —90° is the monostatic angle
69
The computed total scattered fields from a perturbed spherical tumor with a nonuniform radial perturbation across its boundary, ip is the angle with respect to the
y-axis in the y-z plane, ip — —90° is the monostatic angle.
70
The computed total scattered fields from a perturbed ellipsoidal tumor with a nonuniform radial perturbation across its boundary, -p is the angle with respect to the
y-axis in the y-z plane. tp_ = —90° is the monostatic angle
70
The computed Doppler scattered fields from a sinusoidally vibrating spherical tumor
with a uniform peak radial deformation across its boundary, ip is the angle with respect
to the y-axis in the y-z plane, tp = —90° is the monostatic angle
72
The computed Doppler scattered fields from a sinusoidally vibrating spherical tumor
with a non-uniform peak radial deformation across its boundary, ip is the angle with
respect to the y-axis in the y-z plane, ip = —90° is the monostatic angle
73
Simulated magnitude of Doppler scattering responses of a canonically deforming target, tpiaterai is the angle with respect to the y-axis in the y-z plane and ipiaterai = —90°
is the monostatic angle. There is a 6-8 dB contrast in the Doppler signal strength between the uniformly-deforming malignant tumor (5=0.1 mm) and normal fibroglandular tissue clutter (5=0.3 mm) cases. This contrast is larger and inverted, in comparison
with the contrast in fundamental scattering from the two types of targets as in pure
microwave detection techniques
76
gure
Pag
2
Simulated magnitude of Doppler scattering responses of a focused-acoustic induced
deforming target, ipiaterai is the angle with respect to the y-axis in the y-z plane and
^lateral = —90° is the monostatic angle. The Is* Doppler scattering contrast between the deforming malignant tumor (5max=Q.01 mm) and normal fibroglandular
tissue clutter (5max=Q.2?> mm) cases is larger and inverted, in comparison with the
contrast in the fundamental scattering from the two types of targets as in pure microwave detection techniques.
77
3
Simulated magnitude of Doppler scattering responses of a focused-acoustic induced
deforming target which varies in its dielectric properties. The assigned dielectric properties are: the 50th percentile data of normal fibroglandular tissue (e r -44, CT=4.16), the
75th percentile data of normal fibroglandular tissue (er=48, (7=4.71), and the 50th percentile data of malignant tumor (er=50, <r=4.91). ^lateral is the angle with respect to
the y-axis in the y-z plane and ipiaterai = —90° is the monostatic angle
. . . 81
4
Simulated magnitude of Doppler scattering responses of a focused-acoustic induced
deforming target which varies in its elastic properties. The assigned elastic properties
are: E=12 kPa, 9 kPa, and 7.5 kPa that represents different types of malignant tumors,
and E=1.82 kPa that represents a normal fibroglandular tissue clutter, v = 0.495 is
assigned to all these target cases, ipiaterai is the angle with respect to the y-axis in the
y-z plane and ^; atera ; =—90° is the monostatic angle. . .
82
XI
ABSTRACT
This dissertation focuses on technology development that contributes to cost reduction and
sensitivity improvement in microwave breast cancer detection systems.
A microwave refiectometer measures scattered signals from breast tissue, which allows subsequent signal processing to obtain breast images for screening and diagnosis purposes. In current
systems, a commercial vector network analyzer (VNA) functions as a refiectometer, which connects to either a single antenna or an antenna array via switches. An ideal configuration would be a
fully-populated reflectometer-antenna array with each antenna connected to its own refiectometer.
The VNA is too bulky and expensive for this approach. To address this problem, a low-cost mixerbased reflectometer-antenna subsystem is developed. A free-space calibration method is developed
to eliminate the need for connectorized mechanical standards, compensate for antenna reverberation, and allow flexibility in the choice of reference plane position in measurements. Measurement
results using the prototype compare favorably with those using a VNA. This free-space calibrated
reflectometer-antenna subsystem offers potential for cost reduction in a full-array based microwave
breast imaging system.
Heterogeneity in the dielectric properties of normal breast tissue poses a challenge for current microwave techniques to distinguish malignant from normal fibroglandular tissue. Hence, an
acoustic and microwave hybrid modality that exploits both dielectric and elastic properties contrasts in breast tissue is proposed. Two methods of inducing tissue deformation in breast — global
excitation with a compression plate and local excitation with a focused acoustic transducer — are
Xll
investigated. 3-D numerical techniques are developed to efficiently model this multi-physics problem. An acoustic beamformer is designed and the finite-element-method based acoustic-structural
simulation is conducted to compute the induced tissue deformation. The finite-difference timedomain based numerical method with the use of sheet boundary conditions is adopted and extended
to 3-D to simulate the microwave scattering response of the acoustically excited tissue. The numerical study demonstrates a larger and inverted microwave Doppler scattering contrast between
the malignant and normal fibroglandular tissues in comparison with their fundamental scattering
contrast of a pure microwave detection scheme. This reveals the potential of the hybrid modality
for improved sensitivity in detecting a malignant tumor in a heterogeneously dense breast environment.
1
Chapter 1
Introduction
1.1 Background of Current Breast Cancer Detection Techniques
Breast cancer is the second most common cause of death from cancer and the most common
non-skin-related malignancy among women in the United States, according to the US Center of
Disease Control and Prevention. In 2005, over 180,000 new female and over 1,700 new male cases
of invasive breast cancer were diagnosed; in the same year, more than 40,000 women and more
than 370 men died from the disease. Early detection of malignancy, based on regular screening
and diagnosis programs, and thus timely medical treatment is very important in saving lives of the
breast cancer patients.
X-ray breast mammography is currently the primary screening tool for detecting non-palpable
early-stage breast cancer, and is recommended for women aged 35 and older in the United States.
Technical advances such as digital mammography, as well as advances in radiological expertise,
have greatly improved the overall quality and interpretation of mammographic images. However, the sensitivity-related limitations of mammography result in a relatively high false negative
rate (4-34%) [1], particularly in patients with radiographically dense breast tissue such as young
women [2]. Furthermore, the relatively high false positive rate (70%) [3] of mammography leads
to many unnecessary biopsies. Other drawbacks of X-ray mammography include uncomfortable
or painful breast compression and exposure to low levels of ionizing radiation.
Ultrasound scanning, also called sonography, is one of the main alternative imaging techniques
used in clinics [4, 5]. The ultrasound breast exam does not use ionizing radiation, is non-invasive,
and captures images in real-time. The main use of the ultrasound procedure is determining if an
2
abnormality is solid (which may be a non-cancerous lump of tissue or a cancerous tumor) or fluidfilled (such as a benign cyst) or both cystic and solid. However, ultrasound scanning alone cannot
determine the nature of a solid abnormality.
Contrast-enhanced MRI is another main alternative imaging technique for breast cancer detection and is playing an important role in the diagnostic evaluation of mammographically detected
lesions [4, 5]. Contrast-enhanced MRI is minimally-invasive, has good resolution, and provides
valuable information for diagnostic purposes. Yet, this modality is quite costly for routine screening purposes, takes longer to perform than other modalities, and is not recommended for patients
with metallic medical implants.
The tremendous toll that breast cancer takes, combined with the limitations of the current methods, has motivated intense research efforts to advance alternative techniques as complimentary
diagnostic tools.
1.2
Overview of Active Microwave Imaging Techniques for Breast Cancer Detection
Active microwave techniques that exploit the dielectric-properties contrast in breast tissue have
gained increased research interests in the past decade for breast cancer detection. The fundamental
physics that serves as the basis for the active microwave techniques is the dielectric-properties contrast in breast tissue in the microwave frequency range suggested by several dielectric spectroscopy
studies over the past several decades [6, 7, 8, 9, 10]. In microwave breast imaging techniques,
tissue-dependent microwave scattered signals, received by antennas and measured by microwave
refiectometers, are processed to obtain breast images for screening and diagnosis purposes. Microwave tomography [11, 12, 13, 14, 15, 16] and microwave radar [17, 18, 19, 20, 21, 22] are the
two main classes of current active microwave imaging techniques.
Microwave tomography recovers the dielectric-properties profile of an object from measurement of the transmission of microwave energy through the object. This requires the solution of
an ill-conditioned nonlinear inverse-scattering problem via elaborate image reconstruction algorithms. Some promising initial results for breast imaging have been obtained under specialized
3
conditions in 2-D [13], and recent advances in 3-D microwave inverse scattering have been reported in literature [14, 15, 16].
Microwave radar imaging technique uses scattered signals to infer the locations of significant
scatterers in breast. Significant scattering arises from significant contrasts in the dielectric properties, for example, between normal adipose dominated tissues and malignant tumors. This modality
departures from the complicated image reconstruction techniques required by conventional microwave tomography, as it only seeks to identify the presence and location of strong scatterers in
the breast, rather than to completely reconstruct the dielectric-properties profile. Several promising
experimental studies have been reported [21, 22].
Both the active microwave imaging techniques exploit the dielectric properties contrast in
breast tissue, and the current prototypes in both techniques employ a commercial VNA to function
as the refiectometer for measuring microwave scattered signals [13, 21, 22].
1.3
Motivation to Develop a Low-Cost Microwave Reflectometry
A microwave refiectometer, a key hardware in microwave breast cancer detection systems,
can measure the tissue-dependent microwave scattered signals, which allows subsequent signal
processing to obtain breast images for screening, diagnosis and monitoring purposes. In current
systems, one single commercial VNA is employed to function as a high performance refiectometer, which connects either to a mechanically-scanned single antenna or to an antenna array via
switches. However, both approaches have their limitations. The effect of cable movements during
the antenna scanning process is found to be a significant source of measurement non-repeatability.
For the latter approach, additional calibration errors would be introduced, and the switching implementation would limit the employment of multi-channel transmit-receive strategies that offer
potential improvement in resolution and robustness to artifacts in the imaging systems. An ideal
configuration instead would be a fully-populated refiectometer-antenna array with each antenna
connected to its own refiectometer. VNA is too bulky and expensive for this approach, which
4
motivates the development of a low-cost and compact reflectometer that can be paired with an ultrawideband (UWB) antenna as a rigid subsystem for receiving and measuring microwave scattered
signals.
Different methods were proposed to implement a microwave reflectometer. A low-cost reflectometer was constructed using planer directional couplers which are composed of cascaded
Wilkinson power dividers [23]. However, this is a narrowband design and therefore is not suitable
for an ultrawideband application. A six-port architecture based on power detection was proposed
for microwave network analysis, but it suffers from the problem of low dynamic range which is
a side effect of employing power detection instead of coherent detection, limited bandwidth and
difficulty of calibration [24, 25, 26, 27]. Wheatstone bridge was also used for reflection coefficient
measurement and has its application in mm-wave instrumentation area [28, 29].
1.4
Motivation to Investigate an Acoustic and Microwave Hybrid Modality
As mentioned in Section 1.2, the fundamental physics that serves as the basis for active mi-
crowave techniques is the dielectric-properties contrast in breast tissue in the microwave frequency
range. Among those reported dielectric spectroscopy studies [6, 7, 8, 9, 10], early studies [6, 7, 8]
characterized the dielectric properties for a general group of normal tissue, whereas more recent
studies [9, 10] broke down the normal tissue characterization into groups with different adipose
content. It is known from the studies that the dielectric properties contrast between malignant
and normal adipose dominated tissues is as large as on the order of 10:1; However, the contrast
between malignant and normal fibroglandular tissues is small, about 10%. A human breast is a
dielectrically-heterogeneous tissue environment, and very often tumors grow in or near normal
fibroglandular tissues. Therefore, researchers are seeking novel ways to enhance the contrast between malignant and normal fibroglandular tissues in order to achieve better detection sensitivity.
In the past decades, researchers have also conducted experimental studies to measure the elastic modulus of human breast tissue. Reported measurements from different research groups were
often conducted following different procedures which do not necessarily produce exactly the same
results on the same type of tissue. In addition, different tissue handling techniques have been
5
used in generating the data, which is also a source of uncertainty in evaluating the results found
in literature. But overall, despite the above variabilities, previous studies [30, 31, 32, 33, 34] indicate a significant elasticity contrast between malignant tumor and normal breast tissues (both
adipose-dominated and fibroglandular types). As a brief summary, the Young's modulus of different carcinomas are found to be much larger than that of either normal fibroglandular tissue or
normal adipose tissue [30, 31]. Shear modulus and shear viscosity in breast cancer are found to be
higher than those in surrounding normal breast tissues [32, 33]. Invasive and infiltrating ductal carcinoma is found to be significantly stirrer than normal glandular tissues [34]. Such elastic contrast
has motivated research effeorts in various elastography applications for breast imaging, including
breast sonoelastography [35, 36], breast magnetic resonance elastography (MRE) [37], and breast
vibro-acoustography [38]. The author is thus motivated to investigate an acoustic and microwave
hybrid modality that exploits both dielectric and elastic properties contrasts in breast tissue.
Microwave-induced thermoacoustic tomography was proposed for breast cancer detection as
in [39,40,41,42]. However, the physics basis of the microwave-induced thermoacoustic tomography is very different from the proposed acoustic and microwave hybrid approach presented in this
dissertation. In microwave-induced thermoacoustic tomography, microwave energy is selectively
absorbed by the high conductivity tissues in breast and induces acoustic waves by means of thermoelastic expansion. Ultrasound transducers are used to detect the thermoacoustic signals which
are then processed for image reconstruction.
1.5
Contributions
In this dissertation, a low-cost mixer-based reflectometer-antenna subsystem is developed. The
novel mixer-based reflectometer design allows significant cost and size reduction. A free-space calibration method is developed to improve the measurement repeatability by eliminating the need for
contacting mechanical standards, compensating for antenna reverberation, and allowing flexibility
in the choice of reference plane positions in measurements. Measurement results using the prototype compare favorably with those using a VNA. This free-space calibrated reflectometer-antenna
6
subsystem offers potential for cost reduction in a full-array based microwave breast imaging system.
In this dissertation, an acoustic and microwave hybrid modality that exploits both dielectric
and elastic properties contrasts in breast tissue is also proposed. Two methods of inducing tissue
deformation in breast — global excitation with a compression plate and local excitation with a
focused acoustic transducer — are investigated. 3-D numerical techniques are developed to efficiently model this multi-physics problem. An acoustic beamformer is designed and finite-elementmethod (FEM) based acoustic-structural simulations are conducted to compute the induced tissue
deformation. The finite-difference time-domain (FDTD) based numerical method with the use of
sheet boundary conditions (SBCs) is adopted and extended to 3-D to simulate the microwave scattering response of the acoustically excited tissue. From our numerical studies, we found that the
local excitation approach has its limitations in order to be used in this hybrid modality for breast
cancer detection such as the presence of the bulk translation and the difficulty of removing the bulk
translation with a practical antenna array. The focused-acoustic and microwave hybrid modality,
however, overcomes this problem. The numerical simulation demonstrates a larger and inverted
microwave Doppler scattering contrast between the malignant and normal fibroglandular tissues
in comparison with their fundamental scattering contrast of a pure microwave detection scheme.
This reveals the potential of the hybrid modality for improved sensitivity in detecting a malignant
tumor in a heterogeneously dense breast environment.
7
Chapter 2
Ultrawideband Compact Mixer-based Reflctometer (CMR)
This chapter presents the design, fabrication, calibration and measurement results of an ultrawideband compact mixer-based reflectometer (CMR).
2.1 Background Knowledge
\
k
L
D-directivity
C-coupling
3
2
UUT
2
Figure 2.1 Signal flow analysis of a reflectometer
A reflectometer is a circuit that uses a wave separator to isolate and sample the incident and reflected waves from a mismatched load. It measures the reflection coefficient of a one-port network
and, in a more general configuration, the S parameter of a two-port network. Figure 2.1 shows the
signal flow analysis of a reflectometer [43]. In operation, the separator provides a sample Vi of the
incident wave and a sample Vr of the reflected wave. A realistic reflectometer has finite directivity
D, which means both the incident wave and the reflected wave contribute to both Vi and Vr, and
therefore an error in the reflection coefficient measurement is introduced. Assume a unit incident
wave, the expressions of Vi and Vr can be obtained from the signal flow analysis as
Vi = C
C, 9
TeP ,
D
C_
+ CTe?*
D
Vr
(2.1)
where T is the complex reflection coefficient of the device under test (DUT). 9 and 4> are the
phase delay difference through the reflectometer circuit, C and D are the coupling factor and the
directivity of the reflectometer respectively. The calibrated ratio Vr/Vi can represent the reflection
coefficient of the DUT as
I v. I
I r l +4
(2.2).
V
The measurement uncertainty due to finite directivity is approximately ± (1+ | F |) /D, which
implies that higher directivity results in better measurement accuracy.
2.2 Theory and Design
R1
Zin
ZS
*«^/\^u~~»
R25
Zout:
-WW
< ^ ^ >
Figure 2.2 Mixer-based 6dB bridge-T attenuator circuit
A resistive attenuator circuit is employed as a single wave separator in this reflectometer design to achieve compactness. Based on the simulation using Agilent Advanced Design System
(ADS) software, a 6 dB bridge-T attenuator is chosen because it maintains good port characteristic
impedance over the 1-12 GHz frequency range used in this ultrawideband microwave imaging application. Furthermore, in a novel arrangement, the normal 50 Q, resistors to ground are replaced by
frequency mixers with a similar input impedance instead as shown in Figure 2.2. This reflectometer
enables direct detection and measurement of the complex reflection coefficients while absorbing
9
the parasitic resistance of the active elements into part of the attenuator circuit to dramatically reduce the effects of parasitic. Because the broadband mixers chosen for construction do not have a
characteristic impedance of Z0 at the RF input port, a 3 dB attenuator is placed between the resistive part and the RF input port of the mixer to increase isolation for better broadband performance.
VvV
•AM-
OUT
(V,RF
3dBAttn
3dB Attn
Computer
{
N ' IF1
"lo
\'S} IF2
Reference
1
Loefeirv.
Amplifier
Figure 2.3 Mixer-based reflectometer system with two external signal sources
Two cascaded bridge-T attenuator circuits function to separate and detect samples of the incident and reflected waves. In this coherent detection system as shown in Figure 2.3, two synchronized external sources provide RF and LO signals and a lock-in-amplifier measures the complex
IF signals detected by the mixers. The expressions for the two IF signals can be obtained quantitatively through microwave theory analysis as
[-V +
-Vrei29}Aei2e\
IF1
=
IF2
= [^Vejd +
-VTe^9]Bej29"
4
4
(2.3)
where V represents the input RF signal, T is the complex reflection coefficient of the DUT, 9 is
the electrical length of the transmission line connected between two attenuator circuits, Ae^26' and
Bej2e" represent the contribution of LO signal to IF output signals. The signals IF1 and IF2 both
depend on the incident wave as well as the reflected wave; however, IF1 depends relatively more on
the incident wave whereas IF2 depends more on the reflected wave. Thus, the difference between
the two IF signals after calibration can represent the complex reflection coefficient T:
IF1-K-IF2
oc T
(2.4)
10
where K is the complex calibration factor. This complex reflection coefficient is referred as the
raw reflection coefficient Traw in the following discussions.
2.3 Calibration
A
Figure 2.4 A prototype of the cascaded 6dB bridge-T attenuators
Three types of measurement errors are introduced in the microwave refiectometer system: systematic errors, random errors, and drift errors.
Systematic errors are caused by imperfection in the refiectometer setup, for example, the imperfect broadband performance of the bridge-T attenuators, the mismatches at the RF input port,
DUT port and IF output ports, the reflection between the resistive circuit and the input port of
frequency mixer, and the imperfection due to multiple use of adapters and connectors, etc. If these
systematic errors do not vary over time, assuming that a one-time measurement is short enough in
time, then the systematic errors can be removed mathematically through a calibration procedure.
Random errors vary as a function of time and cannot be removed mathematically since they are
unpredictable. In the refiectometer system, the main sources of random errors are the instrument
noise (phase noise of the source signal and sample noise in the lock-in-amplifier) and the connector
n
non-repeatability. Random errors yet can be reduced by increasing the input source power, narrowing the system IF bandwidth and averaging over multiple sweeps. Trial experiment with the
reflectometer prototype has been conducted to find the optimal driven LO and RF power as well as
the appropriate IF bandwidth and averaging without sacrificing too much measurement time.
Drift errors occur when the system performance changes after a calibration has been performed,
which is mainly caused by temperature variation and mechanical instability. The thermal characteristics of multiple cables used in the reflectometer system make the system sensitive to drift errors,
especially when the system becomes heated after being in operation for hours. Efforts are made to
ensure the mechanical stability of the whole experimental setup to further reduce the drift errors.
A one-port error network models the systematic errors which are categorized into three main
terms: ED, ERT, and Es. ED is the directivity error due to signal leakage, ERT is the frequency
response error caused by reflection and transmission tracking within the system, and Es is the
error due to source or impedance mismatches. These three error terms can be derived from the
general equation below by establishing three equations with three known DUTs:
TA = ED + ERT •
F
X
™
III
(2.5)
1 raw
To establish the equations, three known calibration standards need to be measured. In the
reflectometer system, an open, a short and a load (50 fi), which is known as the SOL calibration
standards, are used. Solving these equations then yields the systematic error terms and makes
it possible to derive the actual reflection coefficient r ^ of a DUT from measured raw reflection
coefficient Traw.
2.4 Fabrication and Measurement
A mixer-based reflectometer measurement system is shown in Figure 2.3. The cascaded 6 dB
bridge-T attenuators as shown in Figure 2.4 are fabricated on a Rogers RO3010 printed circuit
board (PCB) (er = 10.2) with surface mount chip resistors. This microstrip line design has an
approximately 2 cm x 2 cm planar size and thus can accommodate the size of SMA connectors.
A Marki Ml-0212 frequency mixer which covers the RF frequency range 2-12 GHz is selected for
12
frequency conversion. An SRS844 RF lock-in-amplifier receives and detects the complex value of
IF1 and IF2 signals. For this low power detection application, a +9 dbm RF input signal is chosen
to feed the reflectometer which, after attenuation through the cascaded 6 dB bridge-T circuits,
produces an approximately -3 dbm power signal ready to be transmitted to the DUT. LO driven
power for the mixer is chosen to be +16 dbm for the optimal linearity performance. Mechanical
open, short and load standards from Agilent 85033D calibration kit are used to perform the oneport SOL calibration. Two different DUTs are used for measurement: one is an open-ended 3 dB
attenuator, and the other is a DUT composed of several connectors and fixed-value termination.
Since the exact reflection coefficient of the second DUT is not known intuitively, it is referred as
an 'unknown DUT' in the following discussions.
2.5 Results
The measured complex reflection coefficient results for the two DUTs are summarized in Figure 2.5. The measurement results using the reflectometer prototype show a fairly good agreement with the results measured by a commercial Agilent E8364A Performance Network Analyzer
(PNA). In addition, directivity and coupling factor, as the important performance indices of a reflectometer, has also been simulated in ADS and measured in the lab. As shown in Figure 2.6, the
coupling factor of the reflectometer is 7 dB ± 1.5 dB, and its directivity is around 6-12 dB in the
wide frequency range of 1-12 GHz.
In the next chapter, this developed compact mixer-based reflectometer technique is integrated
into an ultrawideband CMR-antenna subsystem for free space microwave measurements.
13
-5 vw3?£f&*££?&>'
§-10
-15
-20
o 3dB atteunator DUT (Our Prototype)
« Unknown DUT
(Our Prototype)
3dB attenuator DUT (PNA)
Unknown DUT
(PNA)
4
6
8
Frequency (GHz)
10
12
10
12
(a)
4
6
8
Frequency (GHz)
(b)
Figure 2.5 (a) Magnitude and (b) Phase of the reflection coefficients of the DUTs measured by
the reflectometer prototype as well as by a commercial PNA
14
16
\A
Directivity Measurement
—«— Directivity Simulation
Coupling factor Measurement
Coupling factor Simulation
121S10;
6
8
frequency (GHz)
Figure 2.6 Simulated and measured directivity and coupling factor for the refiectometer prototype
15
Chapter 3
Calibrated Free-Space Measurements with A CMR-Antenna Subsystem
This chapter presents a free-space calibration method and the free-space calibrated microwave
measurements using an ultrawideband CMR-antenna subsystem.
3.1 Motivation and Background
Current microwave imaging systems under development typically employ antennas in combination with a single reflectometer, mostly a commercial VNA, to receive and measure scattered
signals [13, 21, 22]. Traditional calibration procedures used in these systems require repeated
mechanical contact and do not compensate for antenna dispersion, mechanical artifacts or reverberation, which motivates work to investigate a free-space calibration method that is suitable for
the CMR-antenna subsystem in order to avoid the above problems.
A variety of free space measurement systems and calibration techniques have been reported
for use in the microwave regime [44, 45, 46, 47, 48]. The measurement systems in [44, 45] are
based on a two-port thru-reflect-line (TRL) calibration. A modified reduced-cost measurement
system which follows this approach is reported in [46]. Reference [47] constructs a two-parameter
model to approximate antenna-media coupling and uses shifted metal plates as calibration standards. Reference [48] discusses free space calibration procedures using two-term and three-term
error models, as well as free space calibration standards and errors associated with the calibration
methods.
16
3.2 Free-Space CMR-Antenna Measurement System
ref
Lock-in
PC
t IF1
^ S
f-C^Zyd
JVW
Figure 3.1 CMR-antenna system: The compact CMR and a UWB horn antenna are connected
through a semi-rigid cable.
A calibrated free-space measurement procedure is developed for a monostatic system comprised of a microwave CMR and a UWB horn antenna. The calibration procedure involves a precalibration using the empty room standard case, a time-gating process, and a three-term one-port
error model. Multiple versions of the low-cost CMR circuit block, as shown in Figure 3.1, can be
arranged in parallel for array-based systems. This contactless measurement system eliminates the
need for connectorized mechanical standards at the reference planes, thus providing better measurement reliability and repeatability. The calibration method compensates for the UWB antenna
dispersion, the reverberation from the antenna, and reflections between the antenna and the material under test (MUT). It also allows for flexibility in the choice of reference plane positions.
Figure 3.2 shows a combined CMR and UWB horn antenna prototype.
3.3 Free-Space Calibration Method
The measurement and calibration procedure used to compensate for the error sources and measurement artifacts inherent in the whole system is described below.
First, extract the broadband pre-calibrated measured reflection coefficient Traw from the raw
outputs IF1 and IF2 of the CMR. X is a pre-calibration coefficient determined by forcing the Fraw
of the empty room case to be zero.
rP„,„ = X-IF1-
IF2
(3.1)
17
^
I
I
Figure 3.2 Prototype: The compact CMR and a UWB horn antenna are connected through a
semi-rigid cable.
18
Second, apply a time-gating step to obtain a smoothed reflection coefficient FQ- Only the timedomain reflection contributions from the standards or MUT is of interest, thus a window function
H(t) is selected which excludes any other reverberation or background noise and clutter that fall
outside the window. This window is applied to the inverse discrete Fourier transformed timedomain data ^ r _ 1 {r r . o l u }. As a result, a time-gated response TG is obtained by Fourier transforming
the time-domain gated information back into the frequency domain as
YG =
T{T-l{Yraw}-H{t)}
(3.2)
Finally, apply a one-port error model to extract the free space system error terms ei, e2, e3, and
to obtain the calibrated reflection coefficient Tc-
rc = -
Tc-ei
(3.3)
e2 + e 3 • ( r G - ei)
In the free-space measurement system, e\ represents the return loss present in the CMR circuit
and the mismatch between the CMR and antenna; e2 accounts for any UWB antenna dispersion,
multipath effects, semi-rigid cable loss and free space path loss; e 3 takes into account the multiple
reflections returning from the MUT. To find the error terms, three known calibration standards are
employed: 'Empty room' (load), a flat copper sheet placed at the reference plane (reference short),
and the same copper sheet shifted a distance Ad away from the reference plane (offset short), as
illustrated in Figure 3.3. The reference plane is placed in the far field region for this application.
Finally, the three error terms are obtained and used to determine the unknown material reflection
coefficient calibration, TER, TRS> and Tos are the time gated reflection coefficients of the three
standards.
1
0
0
ei
0
— {FRS — FER)
-1
e2
FRS—
e3
TQS
0 _e-3P2Ad{Tos
__
TER)
_e-jf32Ad
FER
_
FER
FER
(3.4)
19
y y y y y y y v j v yV y v y v v
a)
CMR
K]
S
<
<
<
i
i
<
<
^AiVsAiiiiai^
b)
£sSn*!\£iA
1 Copper sheet<]
Cft/IR
yMi^iUi^Sw^^^^
7^?V?;7"^
r-.^v'v^TVVvvf.v-^
M
c)
CIViR
/vVSfTv^
C)
CIMR
Reference
plane
Figure 3.3 The controlled experiment setup used in the free-space calibration and measurement
procedure
20
3.4 Results
Compact (2 cm x 2 cm) microstrip-line based bridge-T attenuator circuits as described in
Chapter 2 were fabricated, and a UWB horn antenna was connected to the output port of circuits
through a 30.5 cm semi-rigid SMA cable. The antenna used was the Tera-X Model 530, a 4.2 cm
x 2.6 cm modified horn antenna that is based in part on work in [49]. The RF input source was
swept from 1 to 12 GHz with a step size of 50 MHz, yielding 221 data points. The RF power
was optimized at +8 dbm for better wideband performance. As shown in Figure 3.3, the free-space
calibration and measurement area was shielded with absorbing foam (ECCOSORB-AN) to provide
an anechoic environment (61 cm x 76.2 cm x 76.2 cm). Optical rails were also used to achieve
precisely controlled and repeatable standard placement. The desirable reference plane position is
in the far-field region for measurement. However, the reference plane position needs to be closer
enough to the transmit antenna in order to achieve a reasonable signal-to-noise (SNR) ratio. With
this trade-off considered, we chose the reference plane to be 10 cm away from the antenna aperture.
A 30.5 cm x 30.5 cm copper sheet was used as the reference short and the offset short standards;
whereas a 30.5 cm x 30.5 cm flat silicone rubber sheet (er = 3, thickness = 1.2 cm, Durometer =
70 A) served as the MUT.
The calibrated time-domain responses in Figure 3.4 and Figure 3.5 show a vivid picture of the
effectiveness of the free-space calibration. The inverse discrete Fourier transformed time-domain
return pulses from the standards as well as the MUT, both prior to calibration as shown in Figure 3.4, i.e., .F -1 {rV ot0 }, and after time-gating and free-space calibration as shown in Figure 3.5,
i.e., Jr~1{Tc}
are shown. A bandpass Gaussian window is used in the transformation. Without
calibration, the original reference plane for the time-domain observation is at the output port of
the CMR. The pulses complete a round trip, traveling through the semi-rigid cable and over the
free-space path, then being reflected and returning. As in Figure 3.4, they arrive approximately 4.3
ns after passing the CMR reference plane. In Figure 3.5, the calibrated return pulses are observed
21
immediately after a time delay defined for the modulated Gaussian source pulse because the reference plane has been relocated by the free-space calibration procedure. It is apparent that the pulse
broadening due to dispersion in the circuit and antenna is effectively compensated.
1]
5
Time (ns)
Figure 3.4 Raw time-domain responses from the standards and MUT
Figure 3.6 and Figure 3.7 show the calibrated complex reflection coefficients of a silicone rubber sheet measured using the CMR-antenna system. This result compares favorably with those obtained using an Agilent E8364A PNA-antenna setup in the frequency range of 3 to 12 GHz. For the
experiments, both instrument configurations employed the same free-space calibration procedure
and controlled measurement setup. The result also compares well with the theoretical reflection
coefficient. As expected, the response below 3 GHz does not merit consideration due to antenna
gain limitations at lower frequencies.
Additional reflection coefficient measurement of materials that have different dielectric properties was conducted with this free-space CMR-antenna prototype. A group of dielectric materials,
including the silicone rubber sheet, a Skycast material (er = 12, thickness = 2.6cm), an FR4
board (er = 4.34, thickness = 0.2cm), and a Plexiglas material (er = 2.7, thickness = 2.4cm), was
selected for measurements. The results, as summarized in Figure 3.8 and Figure 3.9, show good
agreement with the theorectical values.
22
t
3
3CD
<n
c
o
Q.
W
CD
DC
0.5
—
i'n
Ref short
Offset short
MUT
Ac :
fit :
* ^ t l "yr^^lmll j i T T l >||ii<llil>
0
„
(
HIT
CO
£
o
Q
-0.5
0)
1 l': _:
I i". :
11: i
l |
E
Time-Gated & Calibrated -
'-'•
• i ::
•l ~
i
-1
-
0.5
1
Time (ns)
1.5
Figure 3.5 Time-gated and calibrated time-domain responses from the standards and MUT
DQ
-5
o
f
1-15
o
O
c -20
o
1-25
1-30
-35.
%
5
o CMR/antenna
PNA/antenna Theory
6
7
8
9
Frequency (GHz)
10
11
12
Figure 3.6 Measured and calibrated magnitude of complex reflection coefficients of a silicone
rubber sheet (er = 3, thickness = 1.2cm). The measurement results obtained from CMR-antenna
system are compared with those obtained using commercial PNA-antenna setup as well as with
the theoretical values.
23
o CMR/antennaPNA/antenna
Theory
200
0
T3
ts
c
100
0
o
o
O
0-
••G - 1 0 0
0
0
!^
-200
5
6
7
8
9
Frequency (GHz)
10
11
12
Figure 3.7 Measured and calibrated phase of complex reflection coefficients of a silicone rubber
sheet (er = 3, thickness = 1.2cm). The measurement results obtained from CMR-antenna system
are compared with those obtained using commercial PNA-antenna setup as well as with the
theoretical values.
24
Fmynmcf
|GH2|
(a)
I1
0
-10
i
-20
J;
f
" |i
*
£ -30
'•}
'*
ii
\
;;
!
;
'<
';f ~
:'
i
J
i
-
-40
•SO '
;:'
f
i
—J
i
1
— - i ~
-
7
B
Ffetjfl?Has)
(b)
Figure 3.8 Measured reflection coefficients of (a) a Silicone rubber sheet
(er = 3, thickness = 1.2cm) and (b) a Skycast material (er = 12, thickness = 2.6cm) using the
free-space CMR-antenna prototype. Results (circles) compare favorably with the theory values
(dotted lines).
25
Freq(GHz)
(C)
6
7
8
Frequency (GMt)
(d)
Figure 3.9 Measured reflection coefficients of (c) an FR4 board (er = 4.34, thickness = 0.2cm)
and (d) a Plexiglas material (er = 2.7, thickness = 2.4cm) using the free-space CMR-antenna
prototype. Results (circles) compare favorably with the theoretical values (dotted lines).
26
Chapter 4
Concept of Acoustic and Microwave Hybrid Modality
This chapter presents the concept of the acoustic and microwave hybrid modality for breast
cancer detection.
4.1 Theory
In the acoustic and microwave hybrid modality, an acoustic wave is applied to induce tissuedependent deformation in breast interior, and low-power microwave signals also illuminate the
breast through an antenna array which in turn measures the scattered microwave signals. The
induced tissue deformation modulates the scattered microwave signals, contributing to the microwave Doppler scattering response detected by the receiver. This modulation effect of the tissue
deformation on microwave scattering is illustrated with a simple 1-D example. Assume that a halfspace dielectric interface, displacing at an acoustic frequency um, is illuminated by a continuous
plane wave at a microwave frequency LO0, the scattered microwave field from the interface is then
written as
Es(z,t)
= ^{Te-j20ld{t)e+j(}lZej^t}
(4.1)
where F is the reflection coefficient of the dielectric interface, /3i = y- is the microwave wavenumber in the medium, and the interface displacement d(t) is
d(t) = Am sm(ujmt)
(4.2)
The acoustic-modulated scattered field is then expressed as
+00
Es{t)\z=o = Tcos(u0t - 2/31Amsm(umt))
= F ^
.
Jn{
.
m
Ai
) cos ({u0 + num)t)
(4.3)
27
where Jn{x) is referred to as the Bessel function of the first kind of order n. Assume that the
displacement is much smaller than the microwave wavelength, i.e. Am <C Ai which is true for
the proposed application in this dissertation, and simplify the modulated scattered field using a
small-argument Bessel approximation:
t
1
JU(K)
= { ±/c/2
for
n =0
for
\n\ = 1
0
for
\n\ > 1
The acoustic-modulated scattered field is approximated as
Ea(t)\z=0
2irA
2irA
= Tcos(wot) + T — ^ cos((o;o + um)t) - T — ^ cos((u;o - uJm)t)
Ai
(4.4)
Ai
The first term in this expression is the fundamental scattered component which depends on the
dielectric properties of the scatterer. The second and third terms represent the Doppler scattered
components which depend on both dielectric and elastic properties of the scatterer. This theoretical
analysis demonstrates the advantage of introducing acoustic modulation into microwave scattering,
i.e., utilizing both dielectric and elastic properties contrasts in breast tissue instead of just one.
4.2 Methods of Inducing Tissue Deformation
There are two main methods of coupling acoustic power into human tissue, as used in various
biomedical applications: (a) a global excitation approach [50, 51, 52, 53, 54, 55], and (b) a local
excitation approach [56, 57, 58, 59, 60, 61, 62, 63, 64, 65, 66, 67, 68, 69, 70, 71, 72]. For breast
cancer detection, potential ways of inducing tissue deformation in breast are illustrated in Figure 4.1. External acoustic or mechanical excitation with a compression plate is a contact method.
Local excitation with a focused acoustic transducer, however, can be non-contacting, distribute
acoustic power deep in breast, and be configured for real-time beam steering as needed. Regarding
the two options, the local excitation approach with the use of a focused acoustic transducer usually
generates smaller displacement in human tissue (less than a fraction of mm), whereas the global
excitation approach with the use of a compression plate can generate relatively larger deformation. The induced deformation in breast tissue directly affects the signal strength of the microwave
28
Doppler scattering. Numerical studies of the hybrid modality using both acoustic excitation approaches are presented in the following chapters.
antenna
compression
plate _ f
antenna
breast j > -
focused
transducer
(a)
(b)
Figure 4.1 Potential ways of inducing tissue deformation in breast: (a) global excitation with a
compression plate, and (b) local excitation with a focused acoustic transducer
29
Chapter 5
Global-Acoustic and Microwave Hybrid Modality
This chapter presents the numerical study of a global-acoustic and microwave hybrid modality
which employs an external mechanical or acoustic excitation approach for inducing tissue deformation in breast. In addition, a synthetic means of tracking and removing bulk translation through
the use of an antenna array and predictive signal processing is investigated. Finally, the limitations
of using this global excitation approach in the hybrid modality are discussed.
5.1 2-D Numerical Simulation
Modeling the acoustic (or mechanical) wave and microwave interaction in the hybrid modality requires a multiphysics numerical technique. Two distinct time scales are involved in this
problem: the shorter time scale associated with high-frequency microwave illumination and the
longer time scale associated with low-frequency acoustic (or mechanical) excitation. Simulating
such wave interaction with a time-domain numerical technique such as the finite-difference timedomain (FDTD) method [73] would require a prohibitive number of time steps in computation.
Thus, an efficient approach, the quasi-stationary (QS) approximation [74], is employed to solve
this problem. In this multiphysics problem, the tissue target, under mechanical or acoustic excitation (~kHz), deforms sufficiently slow in comparison with the high-frequency (~GHz) microwave
illumination. Thus, the microwave scattered fields should consist of, at any time t, the field pattern corresponding to the instantaneous deformed tissue target at any time Tj on the acoustic time
scale, calculated as if the target is stationary. This field pattern is termed quasi-stationary. This QS
method remains valid under the condition that ujm <C UJQ, where um is the angular frequency of the
30
low-frequency wave that induces tissue deformation and UJQ is the angular frequency of the highfrequency wave that interacts with the deforming tissue target which results in modulated scattered
fields.
Furthermore, the induced tissue deformation in breast is much smaller in comparison with the
microwave wavelength. Resolving such small deformation in the numerical grid also results in
a heavy computational load for FDTD. Thus, sheet boundary conditions (SBCs) [75, 76, 77] are
adopted to efficiently model the microwave interaction with the deforming tissue target excited by
acoustic (or mechanical) wave. With the SBC method, the excited object is modeled as a stationary
object with time-varying sheet boundary conditions which account for the object perturbation on
the relatively slow acoustic (or mechanical) time scale. SBCs were first postulated and derived for
a 2-D dielectric object by Lawrence et al. [75], wherein SBCs were also implemented with the
Method of Moments (MoM) in 2-D. More recently, Buerkle et al. [77] implemented SBCs with
FDTD in 2-D.
Overall, the numerical modeling of the global-acoustic and microwave hybrid modality involves two types of simulations. First, a mechanical simulation, based on a finite element method
(FEM) solution of a plain-strain partial differential equation (PDE), is conducted to compute tissue
deformation under mechanical excitation which is then imported into the subsequent electromagnetic simulation. Next, electromagnetic simulations, based on a FDTD solution of Maxwell's
equations with the use of SBCs, model the microwave interactions with stationary snapshots of the
the tissue target. In FDTD-SBC implementation, the total scattered field is decomposed into the
unperturbed field scattered by the unperturbed object, E 0 , and a perturbation component which
can be reconstructed as the Doppler scattering component, Ed, associated with the vibration of the
object. The steps of the electromagnetic simulations are as follows: (1) An FDTD simulation is
conducted to compute the scattered fields from the unperturbed object illuminated by the incident
electromagnetic wave. (2) The equivalent surface currents for the SBCs are calculated using the
unperturbed field solutions and peak boundary displacements computed from the mechanical simulation. Detailed expressions for the currents are given in [76] and are valid for small displacements
31
relative to the microwave wavelength — a regime that is indeed consistent with the proposed application with a wavelength of 1 cm in the tissue medium for the center microwave frequency 5 GHz.
(3) FDTD-SBC simulations are conducted to compute perturbation components using equivalent
surface currents on the stationary object as the only sources. (4) The perturbation components are
then demodulated to reconstruct the magnitude of the Doppler scattering component.
Implementation of the FDTD-SBC method for a 2D TM2 case is validated by comparing its results with those generated by a brute-force method in which the boundary of the scattering object
is repositioned. This benchmark approach, referred as FDTD-REP, requires two FDTD simulations in which the scattered fields are computed for the same microwave illumination but different
stationary scattering objects: an unperturbed object and the same object under peak deformation.
The FDTD-REP solution, AEs(t), is calculated as the difference between the two sets of scattered
field solutions. It can be shown analytically with the same 1-D problem described in Section 4.1
that
AEs{t)\z=0
= rcos(o;ot) - Tcos(io0t =
1fixAm)
r{cos(cu0t)(l — cos(—r—^)) — sin((j0i) sin(———)}
Ai
Ai
(5.1)
With the assumption that the displacement is much smaller than the microwave wavelength, i.e.
A-m <C Ai, the expression can be approximated as
4-7T A
AEa(t)\*=o
=
-r—-^sin(^t)
(5.2)
Ai
Recall from Section 4.1 that the theoretically derived expression for the Doppler component of the
acoustic-modulated microwave scattered field in this 1-D problem is
Ed{t)
= T
cos((cu0 + ujm)t) - r
Ai
=
—r——^ sin(u)0t) sm((dmt)
Ai
cos((o;o - wTO)t)
M
(5.3)
32
Then it is evident that the amplitude of AES is an excellent approximation to the amplitude of the
Doppler component, provided that the displacement is much smaller than the microwave wavelength. Therefore, the use of the amplitude of the FDTD-REP solution as a benchmark is appropriate for the small deformation case considered here.
The scattering object in the 2D TM2 validation study is a 4-mm-diameter cylinder (er = 50, as
= 4.91 S/m at 5 GHz) in a homogeneous background material (er = 44, as = 4.16 S/m at 5 GHz).
The radius of the inclusion is assumed to vary sinusoidally in time and uniformly across azimuthal
angle, with a peak radial deformation of 0.45 mm. The incident wave is launched by an electric
current line source 7.5 cm away from the inclusion center. Electric field observation points are 4
cm away from the inclusion center. Figure 5.1 shows the normalized amplitude of Ed and A E S
computed at 5 GHz using FDTD-SBCs and FDTD-REP, respectively. The excellent agreement
between these two set of results demonstrates the accuracy of the FDTD-SBC implementation.
In the 2-D numerical testbed of the global-acoustic and microwave hybrid modality, a scenario
that is challenging for a pure microwave detection scheme that is based solely on dielectric contrast
is considered, i.e., investigating microwave scattering from a cylindrical inclusion of malignant tissue (er = 50, as = 4.91 S/m at 5 GHz) surrounded by normal tissue that has significant fibroglandular content (er = 35, as - 3.09 S/m at 5 GHz). For reference, scattering from a fibroglandular mass
(er = 44, as = 4.16 S/m at 5 GHz) surrounded by the same normal tissue environment is also considered. In both cases, the scattering object (tumor or fibroglandular inclusion) is a 4-mm-diameter
circular cylinder. The dielectric properties values assigned to the different tissues are derived from
the data reported in [10]. While the tumor detection scenario described above involves a relatively
small dielectric contrast between the scattering object and the surrounding medium, there is a large
elastic properties contrast. The Young's modulus of the tumor is assigned as E = 12 kPa, while
the surrounding normal tissue is modeled with the Young's modulus reported for fibroglandular
tissue (E = 1.8 kPa) as reported in [30]. For the reference scenario, Young's modulus throughout
the entire phantom (both the fibroglandular inclusion as well as the surrounding normal tissue) is
assumed as E = 1.8 kPa. All tissue types are modeled using the same Poisson's ratio (y - 0.495).
33
-REP, [\Es|2/|Ei|2
SBC, |Ed|2/[Ei[2
-40
m
T3
-50
-60
180
4> (deg)
360
Figure 5.1 Comparison of the amplitudes of the normalized Doppler component computed via
FDTD-SBCs (dashed line) and the normalized difference between stationary scattered field
solutions computed via FDTD-REP (solid line). The 2D scattering object is a deforming
4-mm-diameter circular cylinder undergoing a peak radial displacement of 0.45 mm.
compression plate
normal fibroglandular tissue
1
6L
y"
L=10cm
— inclusion
W=12cm
Figure 5.2 Phantom configuration for the hybrid simulations, with axis of cylindrical inclusion
into the page
34
In the mechanical simulation, the inclusion is embedded 5 cm below the upper surface of the
12 cmxlO cm domain, as shown in Figure 5.2. A mechanical compression plate is placed on
the upper surface to apply an excitation with a peak strain of 5L. In simulation, applied peak
strains from 1 mm to 10 mm (i.e., SL/L = 1-10%) have been tested. The structural mechanics module of Comsol Multiphysics software is used to conduct this FEM-based simulation and
obtain the tissue-dependent deformation. The total deformation of the inclusion under mechanical compression consists of both elastic deformation and bulk translation. In order to isolate the
tissue-dependent elastic deformation of the inclusion only, the bulk translation of the inclusion
center 5L' is subtracted from the total simulated deformation. This calibrated deformation data
from the mechanical simulation is imported into the FDTD-SBC electromagnetic simulation. The
sampling density of the FDTD grid is N = Xi/Ax
= 40, where Ai is the wavelength at 5 GHz
in the surrounding normal tissue. The incident wave is launched by an electric current line source
7.5 cm away from the scatterer center. In this feasibility study, complicating effects due to scattering from the breast surface are avoided by filling the entire electromagnetic computational domain
(aside from the inclusion region) with the normal breast tissue background medium. Electric field
observation points are 4 cm away from the scatterer center. The mechanical simulation results,
reported in Figure 5.3, show the elastic deformation of inclusion boundaries for the malignant tumor or normal fibroglandular inclusion, when a 10% peak strain applies. Bulk translation for both
cases have been removed to separate out the tissue-dependent elastic deformation. The normal
fibroglandular tissue inclusion experiences more than three times the boundary deformation of the
malignant mass. The peak deformation of the top boundary of the inclusion (coordinates x' = 0,
y' = 2 mm) is 0.24 mm for the normal fibroglandular inclusion and 0.07 mm for the malignant
tumor.
The subsequent electromagnetic simulation results, summarized in Figure 5.4, show the Doppler
component in the backward ((f) = 180°) and forward (0 = 0°) scattered fields for the deforming
scatterers as a function of the peak applied strain. Doppler scattering from the deforming normal
fibroglandular inclusion is observed to be stronger than that from the deforming malignant tumor.
The harmonic backscattered signal from the normal fibroglandular inclusion is 8 dB higher than
\
35
2.
1.5
-
•
(mm)
1
0.5
.
\ ^»& * *' ^ * *^^ /
V*j
'•I
0 •if
i il
>%
\ n
-0.5
V*
\» i
original
• *.
1
- ' - ' deformed (normal)
- - - deformed (malignant)
11
Pi
\A
-1
"*4
-1.5
•
^^^^^^^J^^^^
-2
-2
-1
0
x ' (mm)
•
1
2
Figure 5.3 Deformed inclusion boundaries, under a 10% peak applied strain (Bulk translation is
removed). Peak deformation of boundary point (0, 2) is 0.24 mm for the normal fibroglandular
inclusion and 0.07 mm for the malignant tumor inclusion.
36
-45
-55
CD
5. -65
• -
-75Y
^
—©—normal, <>
| = 180°
X
-85 V
- • - normal, § = 0°
*
-95
0
1
malignant, i|) = 180°
- * - malignant, $ = 0°
2
3
4
5
6
7
8
peak applied strain (%)
10
11
Figure 5.4 Computed Doppler components in the backward {<j> = 180°) and forward (</) - 0°)
scattered fields under different peak applied strains, for malignant and normal fibroglandular
scatterers.
37
that from the malignant tumor. In comparison, the backscattered signal from the unperturbed normal fibroglandular mass (in which case only the dielectric properties play a role) is 3.5 dB lower
than that from the unperturbed malignant tumor. Hence, these results show the promise of using this hybrid modality to enhance the overall microwave scattering contrast for the purpose of
distinguishing cancerous tumors from normal fibroglandular tissue in the breast.
5.2 Discussion of Limitations
Despite showing some promising contrast enhancement results, the global excitation approach
has some limitations. One assumption made in the 2-D numerical study is that the induced elastic
deformation can be separated out from the total deformation by removing the bulk translation
of the inclusion target. However, it is difficult for a practical antenna array to track the bulk
translation of the target when external vibration applies to the breast. The bulk translation is not as
tissue-dependent as the elastic deformation, and contributes to a majority of the total deformation.
So with the bulk translation involved, there may not be much observed contrast in the microwave
Doppler scattering between the malignant tumor and normal fibroglandular tissue cases. Therefore,
in the following section, a synthetic means of tracking and removing the bulk translation through
the use of an antenna array and predictive signal processing is investigated and tested with 3-D
numerical simulations.
5.3 3-D Numerical Simulation with a Bulk Translation Removal Method
As in 2-D, the QS approximation method is used to model this 3-D multiphysics problem.
Two types of simulations are involved: (1) mechanical simulation and (2) electromagnetic simulation. However, for the preliminary 3-D work in this section, FDTD-REP method, as described in
Section 5.1, is used to compute the amplitude of the Doppler scattering response from an object
under mechanical excitation. In the 3-D numerical testbed as shown in Figure 5.5, a malignant
tumor or a reference normal fibroglandular tissue clutter, both surrounded by normal breast tissue,
is investigated. The same dielectric and elastic properties as those in the 2-D numerical testbed
38
in Section 5.1 are assigned to these phantom materials. The inclusion (tumor or fibroglandular
tissue) is a 1-cm-diameter sphere. In the mechanical simulation testbed as shown in Figure 5.5(a),
the inclusion is embedded 5 cm below the upper surface of the cylindrical phantom which has a
height of 10 cm and a base diameter of 10 cm. A stainless steel compression plate is placed on the
upper surface to apply an excitation with a 10% peak axial strain. The structural mechanics module of Comsol Multiphysics software is used to conduct this FEM-based simulation and obtain the
tissue-dependent deformation. In the electromagnetic simulation testbed as shown in Figure 5.5(b),
a multi-ring antenna array is modeled to transmit and receive microwave signals. The transmit antenna element is modeled as a A/2 dipole, and the receive antenna elements are modeled as electric
field observation points. All the antenna elements are located 5 cm away from the central axis of
the cylindrical phantom. The sampling density of the FDTD grid is N = Xi/Ax
= 100, where
Ai is the wavelength at 5 GHz in the surrounding tissue background. Similar to what is found in
the mechanical simulation results in previous 2-D work, the total boundary deformation of the inclusion consists of both elastic deformation and bulk translation. The center of both the tumor and
fibroglandular tissue targets experiences a 0.57 cm vertical bulk translation under the peak axial
strain. This bulk translation is much larger than the elastic deformation of the inclusion which is
on the order of a fraction of millimeter or less.
The bulk translation is difficult to track and be removed by a physical antenna which would
need to move synchronously with the breast phantom under the mechanical or acoustic excitation.
Hence, a synthetic means of tracking and removing the bulk translation through the use of an
antenna array and predictive signal processing is investigated as follow.
Figure 5.6 and 5.7 illustrate simplified 2-D views of the 3-D phantom testbed, which helps
explain the bulk translation removal procedure. First, as shown in Figure 5.6, stress Fi is applied
to the upper surface of the phantom via the compression plate. A transmit antenna located at Zi
launches microwave signals to the phantom, and the other antenna elements in the array receive
the scattered microwave signals. Mechanical simulations are conducted to compute the tissue
deformation under this applied stress. The spherical inclusion experiences both bulk translation
and elastic deformation. Yet the elastic deformation is much smaller than the bulk translation.
39
stainless steel
compression plate
*
-^
5cm
H=10em
#- ¥ -
D=10cm
(a)
<p=90°
/C;;;;;;;;;;;;;;!!:*
dipole 5.
antenna
...•••'
\
observation points
(b)
Figure 5.5 3-D numerical phantom testbed for the (a) mechanical simulation and (b)
electromagnetic simulation of the global-acoustic and microwave hybrid modality.
40
Electromagnetic simulations then compute the microwave scattered fields from the inclusion both
before and after the compression. We calculate the difference of the pre-compression and postcompression microwave scattered signals AE S for each receive antenna location. In this testbed,
the tissue background is a homogeneous medium and the tissue target is a spherical inclusion. With
these assumptions, at a certain receive antenna location Zn\, the EM wave propagation path from
the transmit antenna Zt to the tissue target and then back to the receive antenna Zn\ is equivalently
the same for the pre-compression and post-compression cases. Thus, the signals received by the
antenna at Zn\ locations are approximately the same for pre-compression and post-compression,
which results in a minimum value ('null') in the response AE S at the antenna location Zn\ of
the array. Finding this minimum AE S response is referred as 'null' detection in the rest of the
discussion. Next, we apply two different stresses, F\ and F2 as shown in Figure 5.7 to the tissue
phantom, and calculate the bulk translation of the tissue target under compression based on the
'null' detection method as follow. The 'null' location in the AE S response is detected for each
stress case. Then the distance L between the transmit antenna location Zi and the 'null' location
Zn is obtained for each stress case. Assuming a linear stress-strain relationship of the elastic
modulus of the tissue phantom, we then obtain the bulk translation Az of the tissue target under
each applied stress as
AZl
= (L2 - Lx) • -^—
(5.4)
r2 — -f*i
Az2^(L2-L1)--^-r
(5.5)
r2 — r\
After the bulk translation is approximated using the 'null detection' approach, the total deformation result of the inclusion from the mechanical simulation can be calibrated accordingly. The
elastic deformation of the inclusion is obtained by subtracting the calculated bulk translation from
the total deformation result. Figure 5.8 shows the elastic boundary deformation of the malignant
tumor and the normal fibroglandular tissue clutter in this simulation testbed. The peak boundary
deformation of the malignant tumor is 0.24 mm, in comparison with 0.73 mm of the normal fibroglandular tissue clutter case. The calibrated elastic deformation is then imported into the subsequent
electromagnetic simulations.
pre-compression
post-compression
stress F1
V
+n*
© i
Az1
=-n1
-n1
Figure 5.6 An illustration of the 'null' detection procedure
post-compression
stress F1
post-compression
stress F2
V
T~*
imAz1
L:
V
Az,
=-n1
^n2
gure 5.7 An illustration of the bulk translation calculation procedure
42
FDTD-REP simulations are conducted to compute the scattering solutions that approximate
the magnitude of the Doppler component. Figure 5.9 shows the Doppler scattered fields received
at different positions in the central antenna array which aligns with the inclusion center. <>
/ = 180°
represents the dipole source location. A range of dielectric properties (e r =35~50, cr=4.162 S/m),
representing the group of normal tissue most difficult to distinguish from tumor, is assigned to the
normal fibroglandular inclusion in this testbed. Results show that the Doppler scattered field from
the deforming malignant tumor is observed to be 5-10 dB weaker than that from the deforming
normal fibroglandular tissue clutter. This Doppler scattering contrast is not only larger but also
inverted in comparison with the fundamental scattering contrast of the pure microwave detection
scheme.
There is, however, a limit of the physical spacing between the antenna elements in a practical array. Errors in approximating the bulk translation using the 'null' detection method may be
introduced with a coarse antenna spacing. Therefore, a linear estimation method is employed to
interpolate and predict the microwave scattered signals received at locations between the physical
antenna elements in the array. Consider an antenna array which is composed of N physical antenna
elements evenly distributed in space. Let X{ (i=l,2,...N) represents the microwave scattered signal
received by the ith physical antenna, and let X = [x1, x 2 , ...x N ], where the vector element in the
matrix x1 = [xi,xi+1, ...XN,XI,X2,
...£J_I]T
(i=l,2,...N). Also denote v = [vi,v2, . . . « M ] T as the
microwave scattered signal at M locations between the physical antennas in the array. The goal
of the linear estimation is to predict vt(i = 1,2, ...M) based on Xi(i = 1, 2, ...N) received by the
physical antennas. The relationship of the detected and predicted signals can be represented as
v = X«
(5.6)
where a is the coordinate vector of v with respect to X. Define R = X f f X and consider the inner
product of v with X as
j9 = X v
(5.7)
a = R~1f3
(5.8)
The coordinate vector a can be obtained as
43
-,^-^-,—^p.
unperturbed
deformed normal
deformed malignant
^HW.,..
-5 - 4 - 3 - 2 - 1
•..,jM«rr
0 1
x (mm)
2
3
Figure 5.8 Calibrated elastic boundary deformation of the malignant tumor and the normal
fibroglandular tissue clutter. Bulk translation for both cases are removed with a synthetic tracking
method through the use of an antenna array. The normal fibroglandular tissue clutter experiences
about three times the boundary deformation of the malignant tumor.
44
-20
• malignant (e =50,0=4.910)
normal (e=50,o=4.910)
-30h
m
-a
-40 Y
normal (e : 35-50, a=4.162)
LU
**-«> - 5 0
LU
-60
-70
u v%w
^*.-,^%fe ^ H
90
100 110 120 130 140 150 160 170 180
*(deg)
Figure 5.9 The Doppler scattering response from the deforming malignant tumor and deforming
normal fibroglandular tissues. Despite the dielectric properties variation in the normal
fibroglandular tissue, the Doppler scattered field from the deformed malignant tumor is observed
to be 5-10 dB weaker than that from the deformed normal fibroglandular tissue clutter.
45
-20
-30
m
• error=0
•error=1 mm
error=2 mm
-40
normal fibroglandular clutter
uT "50
10
W
<
en
-60
malignant tumor
-70h.
-80
90
100
110
120
130 140
* (deg)
150
160
170
180
Figure 5.10 The Doppler scattering response from a deforming malignant tumor and a deforming
normal fibroglandular tissue clutter, when different levels of errors in approximating the bulk
translation are introduced in the testbed.
46
The microwave scattered fields at locations between the physical antennas can then be predicted
based on the coordinate vector a. and signals received by the physical antennas X. Both the received and predicted scattered signals are used in the 'null' detection method to approximate the
bulk translation of the inclusion. Different choices of physical antenna spacing are tested in simulation (Table 5.1). It is found that larger spacing between the physical antennas results in a larger
error in the approximated bulk translation. It also shows that the spacing between the physical
antennas needs to be smaller than 7 mm in order to obtain a relatively accurate bulk translation that
has a moderate level of error, less than 2 mm. Errors introduced in the bulk translation approximation directly affects the computed Doppler scattering response in this computational testbed, as the
total boundary deformation of the inclusion from the mechanical simulation is calibrated by subtracting the bulk translation before it is imported into the electromagnetic simulation. Figure 5.10
shows the change in the Doppler scattered fields from the malignant and normal fibroglandular
inclusions when different levels of errors in approximating the bulk translation (error: 1 mm or 2
mm) are introduced. When such an error exists, the receive antennas do not track the inclusion
center exactly and therefore the the signal strength of the received Doppler response decreases for
both the malignant and normal cases. However, the Doppler scattering contrast between the malignant tumor and normal fibroglandular tissue clutter remains distinguishable when only moderate
levels of such errors are present.
5.4 Discussion of Limitations
As pointed out earlier in Section 5.1, the global-acoustic and microwave hybrid modality has
its limitation due to the bulk translation issue, as it is difficult for a practical antenna array to track
the bulk translation of the inclusion when external vibration applies to the breast. With the bulk
translation involved, not much contrast in the microwave Doppler scattering between the malignant
tumor and normal fibroglandular tissue clutter cases may be observed.
In Section 5.3, a synthetic method of tracking and removing the bulk translation through the
use of an antenna array and predictive signal processing is investigated and numerically tested
with a 3-D phantom testbed. It is found that small spacing between the physical antennas in
47
Table 5.1 Introduced errors in calculating the bulk translation.
antenna spacing A (mm)
error in calculating the bulk translation (mm)
0.2
0
0.4
0
1
0
2
0
3
0.2
4
0.2
5
0.4
6
1.0
7
1.9
48
the array is required so that only a moderate level of error in approximating the bulk translation is
introduced. However, this desired small spacing may cause other problems such as mutual coupling
effects between the antennas in a practical setting. Furthermore, this bulk translation removal
method does not apply for arbitrarily-shaped target cases, heavily relies on the homogeneity of
the breast phantom, and may complicate the operating procedure of the imaging system based on
this hybrid modality. Considering these combined limitations, a focused-acoustic and microwave
hybrid modality is proposed which will be presented in the following chapters.
49
Chapter 6
Focused-Acoustic and Microwave Hybrid Modality
This chapter introduces the focused-acoustic and microwave hybrid modality and provides an
overview of the 3-D numerical techniques for simulating this multiphysics problem.
6.1 Local Excitation with a Focused Acoustic Transducer
In a focused-acoustic and microwave modality, a focused acoustic transducer deposits energy
to induce tissue deformation in a region of interest (ROI) in breast. Focused transducers have been
employed in many biomedical applications, from low-power uses in sonography and various elastography technologies [60, 61, 62, 63, 64, 65, 66, 67] to high-intensity uses in hyperthermia [68],
thermal ablation [69, 70], and drug delivery [71, 72]. The transducer can be either geometrically
focused or electronically focused. The latter type is considered in this investigation, as it provides
better depth of focus as well as flexibility in convenient beam-steering.
6.2 Overview of the 3-D Numerical Techniques
The numerical modeling of the focused-acoustic and microwave hybrid modality involves two
types of simulations as shown in Figure 6.1, which also employs a QS approximation method as
introduced in Chapter 5.
(1) 3-D acoustic-structural simulation
A coupled acoustic-structural simulation, based on the FEM solutions of the acoustic propagation PDE (Equation 6.1), the stress-strain governing equations (Equation 6.2 and 6.3), and
Newton's Law of motion (Equation 6.4), are conducted to compute the acoustically excited tissue
50
deformation. (In these equations, p is the acoustic pressure, p is the density of the medium, c is
the longitudinal acoustic wave speed in the medium, UJ is the angular acoustic frequency, q and
Q are the acoustic source terms, a and r are the normal and shear stresses respectively, e and 7
are the normal and shear strain respectively, E is the Young's modulus, v is the Poisson's ratio,
u = [u, v, w]T is the deformation, and F is the applied force (N/m) [78].) This simulation is conducted using the acoustic-structure interaction module of the Comsol Multiphysics software. In the
simulation, an acoustic wave is launched as the excitation source, and the excited tissue-dependent
deformation is computed and exported to the subsequent electromagnetic simulation.
v- (-- v p - q
^=Q
<*x
l - i /
Gy
V
&z
TXy
Tyz
T~xz
E
(1 + !/)(! --2v)
V
l-v
v
0
0
0
ex
v
0
0
0
e
0
0
0
tz
2
0
0
Ixy
0
l-2i/
2
0
lyz
0
l-2i/
V
V
\-v
0
0
u
0
0
0
(6.1)
pcz
p
0
0
0
0
v
du
dx
dv
dy
£z
dw
dz
f-x
e
du
Ixy
1 dv
2
v
(6.2)
Ixz
(6.3)
dy "^ dx
dv
lyz
1 dw
dz ^
Txz
2
dy
du
1 dw
dz
dx
-pu \x = F
(6.4)
(2) 3-D electromagnetic simulation
Electromagnetic simulations, based on a FDTD solution of Maxwell's equations (Equation 6.5
and 6.6) using SBCs, model the microwave interaction with the deforming tissue target and compute the Doppler scattered response. (In the equations, E and H are the electric and magnetic
51
fields, e, a, fj,, and a* are the dielectric properties of the medium, and J s and M s are the electric
and magnetic source terms.) In this dissertation, the SBCs for a 3-D dielectric object are derived
and SBC implementation with a 3-D FDTD algorithm is validated.
^!
at
=
<9E
¥
_ I v x E - i ( M s + <7*H)
(6.5)
1
„
1,,
= - V xH--(J l + aB)
(6.6)
(j,
^
Details of the two types of simulations are presented in the following two chapters respectively.
Elastic and
Acoustic
Properties
N
Acoustic-Structural
Mechanics Simulation
Dielectric
Properties
Figure 6.1 Computational steps of simulating microwave scattering response from a deforming
target under focused-acoustic excitation
52
Chapter 7
3-D Simulation of Focused-Acoustic Induced Tissue Deformation
This chapter presents the coupled acoustic-structural simulation for computing the focusedacoustic induced tissue deformation.
7.1 Acoustic Beamformer Design and Modeling
In the acoustic-structural simulation, the acoustic wave is launched by a focused acoustic transducer which is modeled based on an acoustic beamformer design presented as follow. An acoustic
beamformer is a spatial filter that can be used with an array of acoustic transducer elements to
focus energy at some desired location in a spatial field. Beamformer is comprised of one complex
weight in each transducer element channel. The amplitude and phase in each channel is chosen to
obtain constructive interference at the focus location and destructive interference elsewhere. The
goal of the beamformer design is to maximize the energy deposited in a desired ROI in breast
while minimizing energy deposited throughout the remainder of the breast region. Assume that
the focused transducer is a square array transducer composed of N x N small circular pistons,
the frequency response associated with the propagation from the nth piston element to the focus
location is denoted as r by an(um, r), where um is the angular acoustic frequency in operation.
Follow [79] and let the piston surface move uniformly with a speed of Uoelu'mt normal to the baffle,
the analytical propagation response at a far-field location r(r, 6) is obtained as
an(ujm,r) =j — U0kae
2Ji(kasm9)
ka sin 6
(7.1)
53
where a is the radius of the piston, p is the density of the background medium, c is the velocity in
the background medium, and k is the wavenumber. Also define
a(r) =
[ai(ujm,r),a2(uJm,r),...,aN2(iom,T)}T
and let the TV2 x 1 vector h = [hi, h2,..., hN2]T be the concatenation of the weight vectors from
each channel. When we apply the complex weights to the transducer elements, the observed signal
at a location r is denoted as Y(r) as
Y(r) = h H a ( r )
Then the energy received at the location r is
qe(r) = h i / a ( r ) a / / ( r ) h
Thus, the energy deposited in the ROI $ / is expressed as hHQfh, where
Qf=
/
a(r)a*(r)dr
The goal of the weights design is to maximize the total transmit energy deposited in the region $ f
and can be written as
h = arg max
h
hHh
(7.2)
h = maxeig{Qf}
(7.3)
One solution for Equation (7.2) is
where the operator maxeig means selecting the eigenvector that corresponds to the maximum
eigenvalue of the matrix Qf.
Assume that the square transducer, composed of 21 x 21 circular piston elements, has an aperture size of 2 cmx 2 cm and is immersed in a homogeneous background medium (p = 1000 kg/m3, c
= 1500 m/s) that mimics the normal fibroglandular breast tissue. The acoustic operating frequency
is 150 kHz. The desired lateral focal plane is 2 cm away from the transducer aperture. Follow the
design approach above, the transducer design is implemented using the acoustic-structure interaction module of Comsol Multiphysics Software. Figure 7.1 shows an achieved focused beam spot
54
in the (a) lateral and (b) axial focal planes. The focal size is roughly 1 cm x 1 cm in the lateral
plane and 1 cm x 3.5 cm in the axial plane, which is a reasonable spot dimension for targeting a
ROI to search for early stage tumors usually less than 2 cm.
7.2 Induced Deformation Results
Then, with the use of this transducer as the excitation source, 3-D acoustic structural simulations are conducted to compute the induced tissue deformation excited. In the 3-D simulation
testbed as shown in Figure 7.2, a 1-cm-diameter spherical tissue target surrounded by a homogeneous normal fibroglandular tissue background (E = 1.8 kPa, Poisson's ratio u = 0.495) is investigated. This target is either a malignant tumor (E = 12 kPa, v = 0.495) or a normal fibroglandular
tissue clutter (E = 1.82 kPa, v = 0.495) for comparison. Coupled acoustic-structural simulations are
conducted using the acoustic-structure interaction module of Comsol Multiphysics software, and
reveal a roughly 3:1 contrast in the maximal radial deformation 5 between normal fibroglandular
and malignant cases. Tumor deforms less as it is much stiffen Figure 7.3 and Figure 7.4 show the
time-varying boundary deformation S of the malignant tumor and the normal fibroglandular tissue
clutter respectively. Boundary deformation of the targets in both the lateral and axial focal planes
are shown, ipiaterai is the angle with respect to the y-axis in the y-z lateral plane and ipaxiai is the
angle with respect to the x-axis in the x-z axial plane. The group of solid-line curves represents
the deformation for different snapshots on the acoustic time scale. It is observed that the largest
boundary deformation occurs at ipaxiai = 0° for both targets due to the acoustic illumination pattern in ROI generated by the focused transducer. The largest boundary deformation is 0.07 mm for
the malignant tumor case and 0.23 mm for the normal fibroglandular tissue clutter case. Furthermore, no bulk translation is observed from the resulted animations of the boundary deformation
for both targets, which certainly shows an advantage of using the local instead of global excitation
approach.
It is critical to assess the safety of using this type of focused transducer to induce tissue deformation in human breasts. Applying high levels of acoustic power may cause cavitation-related
bioeffects in breast. Therefore, investigation on the indicator for the potential for cavitation-related
55
-1.5 -1 -0.5
0
0.5
1
1.5
y (cm)
Figure 7.1 Acoustic intensity (W/cm2) in the (a) lateral and (b) axial focal plane generated by a
focused square transducer in a normal fibroglandular tissue background: The achieved focal
region is roughly 1 cm x 1 cm x 3.5 cm.
56
1.5 cm
acoustic
force
Rx
\ antennas
<
i
target
2 cm
normal fibroglandular
tissue background
Tx
/
i
/ antenna
2 cm
Figure 7.2 3-D computational testbed for simulating the microwave Doppler scattering response
from a tissue target excited by a focused acoustic transducer. The target is either a malignant
tumor or a normal fibroglandular tissue clutter, surrounded by normal breast tissues.
57
E
E
I
to
45
90
135 180
Materal
E
E
90 135 180
V.axial
Figure 7.3 Focused-acoustic induced time-varying radial deformation of a tumor boundary.
^lateral is the angle with respect to the y-axis in the y-z plane and ipaxtai is the angle with respect
to the x-axis in the x-z plane. The group of solid-line curves represents deformation for different
snapshots on the acoustic time scale. Its largest boundary deformation is 0.07 mm.
58
0.3
0.2
£
E
0.1
0
2
M -0.1
to
-0.2
-0.3
0
45
90
135 180
^lateral
45
90
135 180
^axial
Figure 7.4 Focused-acoustic induced time-varying radial deformations of a normal fibroglandular
tissue clutter boundary, ipiaterai is the angle with respect to the y-axis in the y-z plane and ip axial is
the angle with respect to the x-axis in the x-z plane. The group of solid-line curves represents
deformation for different snapshots on the acoustic time scale. Its largest boundary deformation is
0.23 mm.
59
bioeffects, such as the mechanical index (MI) or thermal index (TI), in this focused-acoustic and
microwave hybrid modality needs to be addressed in future work.
60
Chapter 8
3-D Simulation of Microwave Interaction with a Deforming Tissue Target
This chapter presents the derivation of SBCs for a 3-D dielectric object and the validation of
SBC implementation with a 3-D FDTD algorithm.
8.1 Derivation of 3-D SBCs
SBCs were postulated for a penetrable object in [75] as
h x (Ei - E 2 ) = n x h x (ZfHi + Z | H a )
(8.1)
-n x (H x - H 2 ) = n x h x' (>?Ei + y 2 s E 2 )
(8.2)
where Zf and Yf are the sheet impedances and admittances on the exterior (i = 1) or interior
(i = 2) boundary of the object, n is the unit vector normal to the object boundary, E i and Hx are
the electromagnetic fields outside the object, and E 2 and H 2 are the fields inside the object. If Zf
= 0 and Yf = 0, the above expressions are simplified to the source-free (unperturbed) boundary
conditions as
h x (Eoi - E 0 2 ) = 0
(8.3)
-h x (Hoi - H 0 2 ) = 0
(8.4)
The sheet impedances and admittances for a 3-D dielectric object are derived by enforcing
equivalence in the scattered field solutions for a perturbed (r = a + 8) dielectric sphere, and an
61
Figure 8.1 Sheet impedances Zs and admittances Ys for a 3-D dielectric object are derived by
enforcing equivalence in the scattered field solutions for a) a perturbed (r = a + 8) dielectric
sphere, and b) an unperturbed (r = a) dielectric sphere with SBCs applied at the boundaries.
62
unperturbed (r — a) dielectric sphere with SBCs applied at the boundaries, as illustrated in Figure 8.1. First, consider an unperturbed dielectric sphere (radius a) illuminated by a time-harmonic
plane wave
E 1 = xe~jkircose
(8.5)
H 1 = y^e-jkircose
Z\
(8.6)
where Z\ and k\ are the intrinsic impedance and wavenumber, respectively, of the background
medium. As described in [80], the incident fields can be represented as a superposition of TM r
11
and TE r modes expressed in terms of potentials A$ and Frr*
+oo
A
r = -^TaMkir)P^cos9)
(8.7)
+ 0O
Fj
=
^J2anJn(k1r)P^(cos9)
uZi
,
(8.8)
n=l
where u is the angular electromagnetic frequency, an = j n " ( „"^ , and Jn(x) =
•y/^Jn+i(x).
The eigenfunction solution for the scattered fields outside the sphere can be written as
COS(
Ars = —^V6„^(M^(cos0),
r>a
(8.9)
F/ = ^Vcn^(fc1r)Pn1(cos0),
r>a
(8.10)
n—l
The eigenfunction scattering solution inside the sphere can be written in a similar manner as
, +oo
A
,
=
S^YtdnJn(k2r)Pfcoa6),
.
r<a
(8.11)
r<a
(8.12)
. +oo
s
Fr = ^ V e „ J „ ( f c 2 r ) ^ ( c o s ^ ,
LUZ/2 —,
n—l
where Z2 andfc2are the intrinsic impedance and wavenumber, respectively, of the sphere. Then,
the electric and magnetic field solutions can be determined from the above potentials using the
following relations:
1
rj 2
63
=
EQ
Erh =
1 ld2Ar
1 1 dFr
JLUfj,e r drdO
e r sin 9 d(f>
(8.14)
_L
1 d2Ar |
lldFr
jujfxe r sin 9 drd<f>"*" e r 89
(8.15)
Hr
i
jufie
d2
dr 2
1 1 dAr
(i r sin 9 d<fi
#9
k2)Fr
(8.16)
1 l<9 2 F r
jw[ie r drd9
(8.17)
1 02FP
— r sin 9 drd(f)
[ir d9 + juj/ie
11 dA.
Hj, =
(8.18)
Enforcing the continuity of the tangential components of the fields at r = a yields the following
mode coefficients:
-y/^2fHJn'(kia)Jn(k2a)
,
y/efiliH^
(kia)Jn(k2a)
-yJe2HiJn(kia)Jn(k2a)
y/eifj,2Jn(kia)Jn'{k2a)
- ^ex\i2H\
(8.19)
'{ha)Jn'{k2a)
y/e2iJi\H\ '(kia)Jn'(k2a)
^e1fi2Jn'(k1a)Jn(k2a)
2
• ar,
- ^/e^Hk £r( )'/ (kia)Jn{k2a)
y/e2/j,iHk
- ^jr2H\
r(2)
dr.. =
+
+
(8.20)
(8.21)
(ha)Jn(k2a)
(ha)
Jn'(k2a)
(8.22)
C-T).
yJa-ilXiH^, '(ha)Jn'(k2a)
- ^/e1fi2Hn
(k1a)Jn(k2a)
where ei,//i are the dielectric properties of the medium outside the sphere, and e2, fj,2 are the
dielectric properties of the sphere.
When the sphere is perturbed with a small radial amount 8, radius a in Equations (8.19) and
(8.20) is replaced with a + S. Assuming that 5 is much smaller than the microwave wavelength and
the sphere radius, we apply the first-order Taylor series expansion of each term about 5 = 0. The
expressions then represent the perturbed mode coefficients bn and cn:
-yftm.
Jn'(ha)
+ hSJn"(ha)
Jn(k2a) +
+v/iT^i Jn(ha) + h8Jn'(ha)
\ftm
r(2)"
Hk> (ha) + h5H(2)
(ha)
-y^iM2
m2)(ha)
k25Jn'(k2a)
- \Jn'(k2a) +
k2SJn"(k2a)
an
\Jn(k2a) +
2}
+ h8Hkr(2)'
(ha)
k25Jn'(k2a)
Jn'(k2a) +
k25Jn"(k2a)
(8.23)
64
V£2Mi Jn(ha)
+ h5Jn'(ha)
+y/eifj,2 Jn'(ha)
Jn'(k2a) +
+ hdJn"(ha)
k25Jn"(k2a)
Jn(k2a) + k26 Jn'(k2a)
(8.24)
Cn —
v^I
H{2)(kxa) + h5H{2)'(ha)~\
H{2)\ha)
\jn'(k2a)
+ h5H{2)"(ha)}
+
k26Jn"(k2a)
\jn(k2a)
+
k2SJn'(k2a)
The SBCs (8.1) and (8.2) are also applied to the eigenfunction scattering solutions for the
unperturbed dielectric sphere of radius a to derive alternative expressions for the coefficients in
terms of Z{, Z2, Y{, and Y2. The scattered fields can be represented as a superposition of TM r
and TE r modes. For the TM r mode, applying the sheet boundary conditions to the E$ and He field
components at r = a yields the coefficient bn:
-y/m
[jn'(ha)
- if±ZiJn(ha)]
euT2 \jn(ha)
- JB-YfJ^ha)]
r(2)
(ha) - ^Z{H(2\ha)\
y/eZjIiW*
\Hk2)(ha)
-^jT2
[jn(k2a) +
- ^-Y{H{2)
jw«i
\jn'(k2a)
\Jn(k2a) +
(ha)\
JjfcYjJn'(k2a)
+
JUIC2
if±Z°2Jn(k2a)
^-Y2SJn'(k2a)
\Jn'(k2a) +
^Z}jn(k2a)
(8.25)
Similarly, for the TE r mode, applying the sheet boundary conditions to the E^ and He field components at r = a yields the coefficient <v
~^m
[Uha)
+^m
^T,
[Jn'(ha) - i^YfUha)}
\H^(ha)
^JT2
- ^-ZtJn'(ha)]
- ^-Z{H{2) r(2)' (ha)
jw/n
[jn'(k2a) +
^Y2°jn(k2a)
[jn(M +
\Jn'(k2a) +
r(2)
lit™ (ha) - ^Y{H%\ha)
^Z&'.foaj
™Y2°Jn(k2a)
\Jn(k2a) +
JUfl.2
-^-Z^(k2a
(8.26)
We equate the coefficient (8.23) with (8.25) to obtain the sheet impedances and admittances for
the TM r mode, assuming that 5 is much smaller than a and simplifying kS^j (x)
ks m2 i
Jn(x)
^J^TT
= kS [(^) 2 - l ] , where x = ka. In a similar way, we equate the coefficient (8.24) with
(8.26) to obtain the sheet impedances and admittances for the TE r mode. Consider the superposition of the TM r and TE r modes. The sheet impedances and admittances are derived to be:
Z{ =
-jhZ.S
(—?
K
ha>
(8.27)
65
Z\ =
+jk2Z25
Yl _
K
k2aJ
(8.29)
z1
1
jk26
Yf
(8.28)
- <£>''
(8.30)
where Z\ and Z2 are the intrinsic impedances of the medium outside the sphere and the sphere
itself, respectively. We note that the first-order expressions of 8.27-8.30 exhibit a large degree of
similarity to the first-order sheet impedance and admittance expressions derived in 2-D in [75, 76].
As discussed in [75, 76], the sheet impedances and admittances derived from the first-order
SBCs (8.1) and (8.2) depend on the mode number n. Second-order sheet impedances and admittances, independent of the mode number n and thus applicable to an arbitrarily shaped object, are
derived using the second-order SBCs presented below.
h x (E x - E 2 ) = n x h x
-n x (Hi — H 2 ) = n x n x
"' H ' + ;§< s 'i Hi)+ '* H2+ ! (B2 i H2)
7iEi+ ( Ei)+72E2+ ( E2)
^ ^l
l ^4
(8.31)
(8.32)
Jj is the surface tangent along the object boundary. We postulate that 77, B, 7, and A in 3-D
maintain the same degree of similarity to the corresponding expressions in 2-D (reported in [75,
76]) as observed between the first-order sheet impedances and admittances in 3-D (equations 8.278.30) and 2-D (reported in [75, 76]). Thus, our postulated second-order expressions are as follows:
-jhZiS
(8.33)
V2 = +Jk2Z2S
(8.34)
B'
(8.35)
??i =
Bs2
7is
=
' fci
jZ28
=
' k2
' Zi
72 = +
jk25
(8.36)
(8.37)
(8.38)
66
Ai = -
(8.39)
As —
A2
—
(8.40)
k2Z2
These expressions account for the time-varying boundary deformation as 5 changes on the acoustic
time scale. The second-order SBCs, independent of mode number n, can be applied to an arbitrarily
shaped 3-D dielectric object with arbitrary perturbation.
8.2 3-D FDTD-SBC Implementation
As described in Section 5.1, the total scattered field E from the deforming object is decomposed
into the unperturbed field scattered by the unperturbed object Eo and a perturbation component E.
For a perturbed, stationary object, E is at the same microwave frequency as the incident wave.
For a vibrating object, E can be reconstructed as the microwave Doppler component Ed. The
implementation of the SBCs with a 3-D FDTD algorithm is as follow.
First, an FDTD simulation is conducted to compute the scattered fields from the unperturbed
object illuminated by the incident electromagnetic wave. Then, the equivalent surface currents on
the object boundary which result from the object perturbation are calculated as follows:
J B = n x ( H 1 - H 3 ) = -n x n x 7iE 0 i + ^ ( A i — E o i ) +72E02 + ^ ( ^ 2 ^ E 0 2 )
—Ms = hx (Ei — E 2 ) =
hxnx r?iH01
|(i^H
0 1
) + ,2H02 + | ^ 2 | H
0 2
)
(8.41)
(8.42)
Let the unperturbed field solutions E 0 i = E 0 2 = E 0 and H 0 i = H 0 2 = H 0 , and insert the derived
sheet impedances and admittances (8.33) to (8.40) into (8.41) and (8.42) to obtain
d
1
n x nx ju(e1 - e2)6E0 - -os( ( Hi
M s = hx
fix juj(jii - M 2 ) £ H 0 - • 5 -((
OS
Use the dielectric properties e,
1' -i
d
)-^-E0)
H2 u os
)—£^-H0)
€\
62 . W
OS
- e\ + -?1 and ^ = /^o for both the background medium (i = 1)
and the deforming object (i = 2). The surface currents are then expressed as
J s = ft x h x [joj(e[ - e'2)6E0 + (CTI - a2)6E0]
67
M , = —n x h x
1
( eijw
T
+ ff
1
4(4Ho)
ds
t'2ju+f'dsK
Apply the Fourier Transform rules [81]. The time-domain expressions of the surface currents are
n x n x « - 4)<^o
= —n x n x
. 1 -2ft
e-!
+ (<n - <J2)8SQ
1 -2?t.
e' 2
'
. 5 ,_ d _ , , ,
'ds
ds
(8.43)
(8.44)
These ^-dependent surface currents account for the object perturbation.
Next, FDTD-SBC simulations are conducted to compute the perturbation component E associated with each snapshot on the acoustic time scale, using the equivalent surface currents on the
unperturbed object boundary as the only sources. These currents are introduced by adding the appropriate source terms in the FDTD update algorithm. The electric surface current or electric field
in Equation (8.43) is composed of the x-, y-, and ^-directed components which are not collocated
in a standard Cartesian FDTD grid. And for a curved object such as a sphere, taking the cross product on the electric field and the unit vector normal to the object boundary means that the f-directed
current component is contributed by all the three x-, y-, and ^-directed electric fields. So is the
y- and £-directed current components. Therefore, in order to obtain the <?0 value for calculating
the electric current component at a surface location, all the neighboring electric fields residing on
the same side of the boundary are averaged. For example, to obtain the £Q value for calculating
a x-directed surface current that locates on the interior boundary, all the neighboring y- and zdirected electric fields which are also inside the object boundary are averaged. Then the electric
current is calculated using the SQ value as in Equation (8.43). For calculating the magnetic surface
current at a surface location, the derivative of the magnetic field with respect to the surface tangent
in Equation (8.44) is calculated with a spatial central-difference method, using the difference of
two neighboring magnetic field values devided by an approximated area. This approximated area
is calculated using the radius, the latitude angle, and the azimuth angle of the neighboring field
locations. The computed perturbation component changes as the surface current sources change
on the acoustic time scale.
68
Finally, the computed perturbation component solution that varies on the acoustic time scale is
demodulated to reconstruct the Doppler scattering component E j . This is a two-step demodulation
process. The theoretical analysis of the Doppler reconstruction method is provided in Appendix
A. The computed perturbation component from the FDTD-SBC simulation is a time-domain pulse
response which corresponds to a wideband signal in the frequency domain, thus taking the Fourier
transform of this perturbation component can extract its phasor component at the center microwave
frequency associated with each time snapshot. This phasor component of the perturbation component which varies on the acoustic time scale is then demodulated via Fourier transform to obtain
the Doppler scattering component.
8.3 Validation
For validation purpose, a 3-D numerical testbed as shown in Figure 8.2 is considered. The scatterer in the validation study is a deforming malignant tumor (er -50,a = 4.91 at 5 GHz) in a normal
fibroglandular tissue background (er = 35, a - 3.09 at 5 GHz). Multiple canonically deforming
tumors of different shapes are chosen for the validation studies. The tumor is either perturbed and
stationary or vibrating. The incident electromagnetic wave is launched by a A/2 dipole 2 cm away
from the tumor center. 360 electric field observation points are selected to be 1.5 cm away from the
tumor center in a 2-D plane. The sampling density of the FDTD grid is TV = Ai/Ax = 100, where
Ai is the electromagnetic wavelength at 5 GHz in the normal fibroglandular tissue background,
which results in a computational grid size of 141 x341 x441. For the perturbed, stationary tumor
cases, the benchmark results are obtained by directly meshing the perturbed tumor and conducting
a standard FDTD simulation. For the vibrating tumor cases, the vibrating tumor is modeled with
two stationary snapshots on the acoustic time scale, i.e. the unperturbed tumor and the perturbed
tumor under peak deformation. We then directly mesh the tumor for both snapshots and conduct
standard FDTD simulations. With the assumption that the deformation is much smaller than the
microwave wavelength for a simple canonically vibrating case, the difference of these two stationary solutions can be calculated and used as a good approximation of the magnitude of Doppler
component.
69
normal fibroglandular
tissue background _
=
i.
m
^ Rx
antennas
',1.5 cm
2 cm
Tx
"antenna
Figure 8.2 3-D numerical testbed for validating the FDTD-SBC implementation: simulate the
microwave scattering response from a deforming malignant tumor surrounded by a normal
fibroglandular tissue background
-10
-15
•SBC
• benchmark
-20
•E,
-25
-45
0
45
V (deg)
135
180
Figure 8.3 The computed total scattered fields from a slightly bloated spherical tumor with a
uniform radial perturbation across its boundary, tp is the angle with respect to the y-axis in the y-z
plane, tp = —90° is the monostatic angle.
70
DO
M—
T3
(D
(D
•+—»
+-»
ro
o
co
"to
o
-50"
•
>
-180 -135 - 9 0
>
-45
•
•0
45
¥(deg)
90
135
180
Figure 8.4 The computed total scattered fields from a perturbed spherical tumor with a
non-uniform radial perturbation across its boundary, tp is the angle with respect to the y-axis in
the y-z plane, ip = —90° is the monostatic angle.
-50'
•
'
'
•
•-180 -135 -90 -45 ¥ (deg)
0
45
90
135 180
Figure 8.5 The computed total scattered fields from a perturbed ellipsoidal tumor with a
non-uniform radial perturbation across its boundary, ip is the angle with respect to the y-axis in
the y-z plane, ij: = —90° is the monostatic angle.
71
Figure 8.3 shows the total scattered fields E 0 + E computed at 5 GHz from a slightly bloated
spherical tumor with a uniform radial perturbation across its boundary. The total scattered fields in
this and the following figures are all normalized by the scattered field received at the tumor center
location but in a homogeneous normal fibroglandular tissue background (i.e., without the tumor
existence). The size of the unperturbed tumor is 1 cm in diameter. The uniform radial perturbation
is 5 = 0.3 mm. ip m this and the following figures is the angle with respect to the y-axis in the
y-z plane, representing different field observation locations. The monostatic angle is ip = —90°.
The FDTD-SBC result shows good agreement with the benchmark result. For reference, the total
scattered field for the unperturbed sphere is also included to demonstrate that E 0 + E differs from
EoFigure 8.4 shows the total scattered fields computed at 5 GHz from the same spherical tumor
which is perturbed with a non-uniform radial perturbation across its boundary. The radial perturbation across the boundary is 5(p,4>,ip) = <5mcos(20)| cos^l, where 5m = 0.3 mm and (f> is the
azimuth angle in the x-y plane. Figure 8.5 shows the total scattered fields computed at 5 GHz from
an ellipsoidal tumor which is perturbed with a non-uniform radial perturbation across its boundary.
2
2
2
The boundary of this ellipsoidal tumor can be represented as ^ + j - + ^ < 1, where a0 = 5 mm,
xl1
/2
z!2
bo = 4 mm, CQ = 3 mm. The perturbed boundary can be represented as-^- + ^- + -^- < 1, where a
= 5.4 mm, b = 4.3 mm, c = 3.2 mm. There is good agreement between FDTD-SBC results and the
benchmarks for both cases, which shows the validity of FDTD-SBC implementation for an object
under arbitrary boundary perturbation.
Figure 8.6 shows the computed Doppler scattered fields from a sinusoidally-vibrating spherical
tumor with a uniform peak radial deformation across its boundary. The Doppler scattered fields
in this and the following figures are all normalized by the scattered field received at the tumor
center location but in a homogeneous normal fibroglandular tissue background (i.e., without the
tumor existence). The size of the tumor is 1 cm in diameter. The uniform peak radial deformation
is 5(t) = 5mcos(u>mt), where 5m = 0.3 mm and u>m is the angular acoustic frequency. For a
simple canonically vibrating case such as this, FDTD-SBC simulation for only one snapshot on
the acoustic time scale, i.e. when the object is under peak deformation, needs to be conducted to
72
compute the magnitude of Doppler component. Again, good agreement between the FDTD-SBC
result and the benchmark is observed.
Figure 8.7 shows the computed Doppler scattered fields from the same spherical tumor which
is sinusoidally vibrating with a non-uniform peak radial deformation across its boundary. The
non-uniform peak radial deformation is 5(p, cj), ip, t) = Sm cos(2<f>)\ costp\ cos(ujmt), where 5m 0.3 mm and oom is the angular acoustic frequency. The agreement between the FDTD-SBC and
benchmark results in this case shows the validity of FDTD-SBC implementation for a vibrating
object with an arbitrary deformation across its boundary.
-?80 -135 -90 -45
0
45
¥ (cleg)
90
135 180
Figure 8.6 The computed Doppler scattered fields from a sinusoidally vibrating spherical tumor
with a uniform peak radial deformation across its boundary, ip is the angle with respect to the
y-axis in the y-z plane, ip — —90° is the monostatic angle.
73
-30
CD
U—
T3
(D
0
a
ro
o
w
-40
SBC
— benchmark
•
-50
\
-60
/
V
-70
i_
jy
o. =80
o.
o
Q
-90
-180 -135 -90 -45
0
45
¥(deg)
90
135 180
Figure 8.7 The computed Doppler scattered fields from a sinusoidally vibrating spherical tumor
with a non-uniform peak radial deformation across its boundary, ip is the angle with respect to the
y-axis in the y-z plane, ip = —90° is the monostatic angle.
74
Chapter 9
3-D Computational Investigation of Focused-Acoustic and Microwave
Hybrid Modality
Based on the numerical techniques presented in the previous two chapters, this chapter presents
a 3-D computational investigation of the focused-acoustic and microwave hybrid modality.
9.1 Results - Microwave Scattering from a Canonically Deforming Target
As a preliminary step, a canonically deforming target case in which the target is assumed to
be deforming sinusoidally with a uniform peak deformation across its boundary is first investigated. Consistent with the acoustic-structural simulation testbed in Chapter 7, the tissue target
under investigation is either a malignant tumor (er = 50, as = 4.91 S/m at 5 GHz) or normal fibroglandular tissue clutter. (er = 44, as = 4.16 S/m at 5 GHz). In both cases, the target (tumor
or fibroglandular tissue clutter) is a 1-cm-diameter sphere surrounded by normal tissue that has
significant fibroglandular content (er = 35, as = 3.09 S/m at 5 GHz). This represents a scenario
that is challenging for the pure microwave detection scheme solely based on dielectric contrast
to distinguish malignant tumor from normal tissue clutter. As shown in Figure 7.2, the incident
wave is launched by a A/2 dipole 2 cm away from the target center. 360 electric field observation
points are selected to be 1.5 cm away from the target center in a 2-D y-z plane. The sampling
density of the FDTD grid is N = X\/Ax = 100, where Ai is the wavelength at 5 GHz in the normal fibroglandular tissue background, resulting in a computation grid size of 141 x341x441. A
roughly 1:3 peak deformation contrast between malignant tumor and normal fibroglandular tissue
clutter from the acoustic-structural simulations is observed as in the Chapter 7; thus, the tumor
75
boundary is assumed to be deforming sinusoidally with a uniform peak deformation of 0.1 mm,
whereas the more elastic fibroglandular tissue clutter boundary is assumed to be deforming sinusoidally with a uniform peak deformation of 0.3 mm. For this simple canonically deforming case,
FDTD-SBC simulation for only one snapshot on the acoustic time scale, i.e., when the target is
under peak boundary deformation, is conducted to compute the magnitude of the Doppler scattering component. As reference, standard FDTD simulations are also conducted to compute the
fundamental scattering fields from unperturbed targets (tumor or fibroglandular tissue clutter), as
in the pure microwave detection approach. The simulation results are summarized in Figure 9.1
which shows the Doppler scattered fields from the two deforming malignant and normal tissue
targets as well as the fundamental scattered fields from both unperturbed targets, lateral is the
angle with respect to the y-axis in the y-z plane. Because the 1-cm-diameter spherical target under
the microwave illumination in this testbed is a very weak back-scatterer, the strongest responses in
both the fundamental and the Doppler scattering are observed to be at the forward scattered field
location (ipiaterai = 90°). The magnitude of Doppler component from the deforming malignant
tumor is observed to be lower than that from the deforming normal fibroglandular tissue clutter
in the same normal tissue background. There is a 6-8 dB contrast in the Doppler signal strength
between the two types of targets, which is observed to be larger and inverted in comparison with
the contrast in their fundamental microwave scattering signals.
9.2 Results - Microwave Scattering from a Focused-Acoustic Induced Deforming Target
Now the focused-acoustic excited tissue target case is addressed, wherein the simulated deformation results induced by the focused transducer, as shown in Figure 7.3 and Figure 7.4, is
imported into the 3-D electromagnetic simulations. The orientation of the acoustic transducer and
the microwave antennas is illustrated in Figure 7.2. The same dielectric properties as in the canonically deforming case are assigned to the targets investigated. FDTD-SBC simulations for a series
of 12 snapshots on the acoustic time scale are conducted to compute the perturbation scattering
component associated with each snapshot. The sampling number within one acoustic period is
76
-10
-20
•o
» -30
a>
"§ -40
i_
(U
(U
8 -50
-60
-180 -135 -90
-45
0
45
de
lateral < 9)
90
135
180
Figure 9.1 Simulated magnitude of Doppler scattering responses of a canonically deforming
target, ipiaterai is the angle with respect to the y-axis in the y-z plane and ipiaterai = —90° is the
monostatic angle. There is a 6-8 dB contrast in the Doppler signal strength between the
uniformly-deforming malignant tumor (5=0.1 mm) and normal fibroglandular tissue clutter
(5=0.3 mm) cases. This contrast is larger and inverted, in comparison with the contrast in
fundamental scattering from the two types of targets as in pure microwave detection techniques.
77
m
T3
to
O
U)
0
-10
-20
-30
-40
-50
-60
-70
-80
-90
• normal fibroglandular clutter
•malignant tumor
-10'
-180 -135 -90
fundamental
st
1 Doppler
-45
0
lateral
45
(cleg)
90
135
180
Figure 9.2 Simulated magnitude of Doppler scattering responses of a focused-acoustic induced
deforming target, ipiaterai is the angle with respect to the y-axis in the y-z plane and
^lateral = —90° is the monostatic angle. The Is* Doppler scattering contrast between the
deforming malignant tumor (Smax=0.07 mm) and normal fibroglandular tissue clutter (5max=0.23
mm) cases is larger and inverted, in comparison with the contrast in the fundamental scattering
from the two types of targets as in pure microwave detection techniques.
78
determined based on a trade-off of the requirement to reconstruct higher-order Doppler field components of interest and the computational economy. In our simulation, we are mainly interested
in the 1st Doppler scattering component. Then we demodulate the perturbation phasor component
that varies along the acoustic time scale via Fourier transform to reconstruct the Doppler scattering
components. A theoretical analysis of the Doppler reconstruction method is provided in Appendix
A. Figure 9.2 shows that the l 5 i Doppler scattering contrast between malignant tumor and normal
fibroglandular tissue clutter is again larger and inverted, in comparison with the contrast in the
fundamental scattering from the two types of targets as in pure microwave detection techniques,
which offers potential for using this hybrid modality in enhancing the overall microwave scattering contrast for the purpose of distinguishing cancerous tumors from normal fibroglandular tissue
clutter in breast.
9.3 Results - Impact of Dielectric and Elastic Properties Variations on the Microwave Doppler Scattering
In the hybrid modality, the signal strength of the focused-acoustic modulated Doppler scattering depends on both dielectric and elastic properties in breast tissue. Therefore, with this 3-D
focused-acoustic and microwave computational testbed, a series of simulations are conducted to
investigate the impact of dielectric properties variation and elastic properties variation of the deforming target on the Doppler scattering responses.
First, the dielectric properties of the target in this computational testbed are chosen to vary.
The large-scale characterization data reported in [10] suggests that the dielectric properties of the
normal fibroglandular tissue vary much more than those of the malignant tumor. Therefore, the
following dielectric properties are assigned to the target under investigation: the 50th percentile
data of normal fibroglandular tissue (er=44, CT=4.16 S/m), the 75th percentile data of normal fibroglandular tissue (er=48, cr=4.71 S/m), and the 50th percentile data of malignant tumor (er=50,
cr=4.91 S/m). For all the three cases, the elastic properties of the target are assigned as (E=1.82
kPa, ^=0.495). All other parameters remain the same as in the previous focused acoustic and microwave simulation testbed. Results in Figure 9.3 show that the Doppler signal strength increases
79
with the increase of the dielectric properties of the target. Since the same elastic properties are
assigned, the targets experience the same deforming behavior under the focused-acoustic excitation in all three cases. Therefore, the observed Doppler scattering contrast solely depends on the
dielectric properties contrast of the targets. It is reasonable to conclude from this result that the
Doppler scattering strength will become weaker for the group of normal breast tissue that has lower
dielectric properties such as the adipose-dominated tissue. However, this type of tissue is much
easier to distinguish from malignancy using the pure microwave detection scheme, as there is large
dielectric contrast between the two tissue types. The main potential of the hybrid modality exists in
situations such as distinguishing the normal breast tissue that has significant fibroglandular content
from the malignant tumor.
After examing variations of dielectric properties, the elastic properties of the target in this computational testbed are chosen to vary. The elastic properties of the malignant tumor vary widely
depending on its cancerous type and stage, resulting in a range of elastic contrast between malignant tumor and normal fibroglandular tissue as reported in the literature [30, 31, 32, 33, 34]. Thus
the following elastic properties are assigned to the target under investigation: E=12 kPa, 9 kPa, and
7.5 kPa, that represent the elastic properties of different types of malignant tumors, and E=l .82 kPa
that represents the stiffness of a normal fibroglandular tissue clutter; v = 0.495 is assigned to all
these target cases. The dielectric properties of the target under investigation is assigned as: e r =50,
cr=4.91 S/m. All other parameters remain the same as in the previous focused-acoustic and microwave simulation testbed. In the acoustic-structural simulations, the largest radial deformation
of the target across its boundary induced by the focused acoustic transducer is found to be: 0.07
mm (E=12 kPa target), 0.11 mm (E=9 kPa), 0.17 mm (E=7.5 kPa) and 0.23 mm (E=1.82 kPa). The
deformation of each target case is then imported to the subsequent electromagnetic simulations
to compute the Doppler scattering response. The targets, assigned with different elastic properties, experience different deforming patterns under the focused-acoustic excitation, which results
in different scattering patterns across the observation angle as shown in Figure 9.4. However, the
averaged Doppler signal strength detected across the observation angles is found to be: -78.6 dB
(E=12 kPa target), -75.8 dB (E=9 kPa), -71.9 dB (E=7.5 kPa), and -66.4 dB (E=1.82 kPa), which
80
shows that the overall Doppler signal strength also increases with the decrease of the elastic modulus of the target. The result also implies that this hybrid modality performs better in distinguishing
stiffer breast carcinomas from normal fibroglandular tissue clutter.
Overall, the computational investigation of the hybrid modality demonstrates a larger and inverted microwave Doppler scattering contrast between the malignant and normal fibroglandular
tissues in comparison with their fundamental scattering contrast of a pure microwave detection
scheme. This reveals the potential of the hybrid modality for improved sensitivity in detecting a
malignant tumor in a heterogeneously dense breast environment.
81
m
T3
2
CD
Q.
O
Q
180
Figure 9.3 Simulated magnitude of Doppler scattering responses of a focused-acoustic induced
deforming target which varies in its dielectric properties. The assigned dielectric properties are:
the 50th percentile data of normal fibroglandular tissue (er=44, cr=4.16), the 75th percentile data
of normal fibroglandular tissue (er=48, CT=4.71), and the 50th percentile data of malignant tumor
(er=50, (7=4.91). ipiaterai is the angle with respect to the y-axis in the y-z plane and
^lateral = —90° is the monostatic angle.
82
-30
CQ
T>
-nu
<>
/
-50
•a
<u
ip
- - - e =50,a=4.91
r
e =50,c=4.91
r
e =50,a=4.91
r
e=50,o=4.91
r
•
•
T)
P
CI)
-KO
S/m, E=12 kPa
S/m, E=9 kPa
S/m, E=7.5 kPa
S/m, E = 1 . 8 2 k P a
(0
u
-70
' ^
i_
\
CD
n
a.
o
Q
-HO
•
r
S
^
-
, ' '
—
1 «.
/,''
-. jt~
.*'.*..'*'—
'+
+
'
v
"
-
^ S _ _ ^ ^ ^
^—"^
s
-
s
-90
-100
-180 -135
,•
'v.'- .-V
-90
-45
.-''><-"-*».' '''•.,
.•• y
s
•••
-....•' /•
.'
Vs
.-;'
'
>~
\ V
~--'
45
, (deg)
90
135
180
Figure 9.4 Simulated magnitude of Doppler scattering responses of a focused-acoustic induced
deforming target which varies in its elastic properties. The assigned elastic properties are: E=12
kPa, 9 kPa, and 7.5 kPa that represents different types of malignant tumors, and E=1.82 kPa that
represents a normal fibroglandular tissue clutter, v — 0.495 is assigned to all these target cases.
^lateral is the angle with respect to the y-axis in the y-z plane and ipiaterai = —90° is the
monostatic angle.
83
Chapter 10
Conclusion and Future Work
This dissertation addresses two important issues in current microwave imaging systems for
breast cancer detection: cost reduction and sensitivity improvement.
The dissertation first presents the development of a low-cost CMR-antenna subsystem. The
specific contributions include:
(1) Cascaded mixer-incorporated bridge-T attenuators are used as the wave separators in this
refiectometer design. This approach allows significant cost and size reduction of a microwave
refiectometer.
(2) A three-step free-space calibration method is developed for the subsystem. Unlike the conventional mechanical calibration procedure, this contactless method eliminates the need for connecting mechanical standards repeatedly at the reference planes and compensates for the antenna
reverberation, thus providing better measurement repeatability. The method also allows flexibility
in the choice of reference plane positions.
(3) The free-space calibrated measurement results using the CMR-antenna prototype compare
favorably with those using a VNA-antenna setup. Therefore, multiple sets of the reflectometerantenna subsystem can be used to construct a low-cost array-based hardware system for microwave
breast imaging.
Following this track, future research involves the development of a fully-populated receiver
with the use of this free-space calibrated CMR-antenna subsystem. Specific tasks of potential
interest include the following:
84
(1) Conduct simulation to investigate the mutual coupling effects between neighboring antenna
elements in the array. Better understanding of the mutual coupling effects aids in the design of a
fully-populated configuration.
(2) Develop free-space calibration standards for use in a specific array. Depending on the array
configuration (circular, planar, etc.), different free-space calibration standards will be needed.
(3) Construct a fully-populated CMR-antenna array prototype, and investigate its performance
in measuring scattered response in both time and frequency domains.
This dissertation also presents an acoustic and microwave hybrid modality that exploits both
dielectric and elastic properties contrasts in breast tissue. The specific contributions include:
(1) Two methods of inducing tissue deformation in breast, global excitation with a compression
plate and local excitation with a focused acoustic transducer, are investigated.
(2) 3-D numerical techniques are developed to efficiently model this multi-physics problem.
The FEM-based acoustic-structural simulations are conducted to compute breast tissue deformation excited either by mechanical vibration or by focused-acoustic illumination. Simulation results
show that the malignant tumor which has higher elastic modulus deforms less than the normal
fibroglandular tissue target. The FDTD based numerical method with the use of sheet boundary
conditions is adopted and extended to 3-D to simulate the microwave scattering response of the
acoustically excited tissue.
(3) Despite showing some promising initial results, the numerical study of the global-acoustic
and microwave modality reveals the limitations of employing an external excitation approach,
including the difficulty of tracking and removing the bulk translation of the excited target via a
physical or synthetic antenna array.
(4) The numerical study of the focused-acoustic and microwave modality, however, shows that
it avoids the bulk translation problem. The computational investigation demonstrates a larger and
inverted microwave Doppler scattering contrast between the malignant and normal fibroglandular
tissues in comparison with their fundamental scattering contrast of a pure microwave detection
scheme. While the pure microwave imaging scheme creates a breast image of adipose-dominated
and fibroglandular tissues, the hybrid modality can provide additional contrast in the microwave
85
Doppler scattering response to help distinguish the malignant from normal fibroglandular tissues
and provide important diagnostic information. Therefore, the main role of this hybrid modality is
in detecting malignant tumors in a heterogeneously dense breast.
Following this track, future research tasks are as follows:
(1) Simulate the focused-acoustic modulated microwave scattering response with a realistic
3-D breast phantom testbed, wherein a heterogeneous breast interior (preferably an anatomically
realistic breast interior), skin layer, and coupling media layers for the acoustic transducer are included in the simulation testbed.
(2) Perform further parametric studies of the focused-acoustic modulated microwave Doppler
scattering strength. Parameters of interest include the size, the shape, and the density of the lesion,
the properties of the heterogeneous tissue background, the operating frequency and power of the
acoustic source, and the operating microwave frequency.
(3) Conduct safety assessment of the focused-acoustic excitation approach employed in the
hybrid modality.
86
LIST OF REFERENCES
[1] P. T. Huynh, A. M. Jarolimek, and S. Daye, "The false-negative mammogram," Radiographics, vol. 18, no. 5, pp. 1137-1154, 1998.
[2] V. P. Jackson, R. E. Hendrick, S. A. Reig, and D. B. Kopans, "Imaging of the radiographically
dense breast," Radiology, vol. 188, pp. 297-301, 1993.
[3] J. G. Elmore, M. B. Barton, V. M. Moceri, S. Polk, P. J. Arena, and S. W. Fletcher, "Ten-year
risk of false positive screening mammograms and clinical breast examinations," New Engl. J.
Med., vol. 338, no. 16, pp. 1089-1096, 1998.
[4] M. Sabel and H. Aichinger, "Recent developments in breast imaging," Phys. Med. Biol,
vol. 41, pp. 315-368, 1996.
[5] S. H. Heywang-Kobrunner, "Nonmammographic breast imaging techniques," Current Opinion in Radiology, vol. 4, pp. 146-154, 1992.
[6] A. J. Surowiec, S. S. Stuchly, J. R. Barr, and A. Swamp, "Dielectric properties of breast
carcinoma and the surrounding tissues," IEEE Trans. Biomed. Eng., vol. 35, pp. 257-263,
April 1988.
[7] W. T. Joines, Y. Zhang, C. Li, and R. L. Jirtle, "The measured electrical properties of normal
and malignant human tissues from 50 to 900 MHz," Med. Phys., vol. 21, pp. 547-550, April
1994.
[8] S. S. Chaudhary, R. K. Mishra, A. Swarup, and J. M. Thomas, "Dielectric properties of
normal and malignant human breast tissues at radiowave and microwave frequencies," Indian
J. Biochem. Biophys., vol. 21, pp. 76-79, February 1984.
[9] A. M. Campbell and D. V. Land, "Dielectric properties of female human breast tissue measured in vitro at 3.2 GHz," Phys. Med. Biol, vol. 37, pp. 193-210, 1992.
[10] 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. Hagess, "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.
87
[11] S. Y. Semenov, A. E. Bulyshev, A. E. Souvorov, R. H. Svenson, Y. E. Sizov, V. Y.
Borisov, V. G. Posukh, I. M. Kozlov, A. G. Nazarov, and G. P. Tatsis, "Microwave tomography:Theoretical and experimental investigation of the iteration reconstruction algorithm,"
IEEE Trans. Microw. Theory Tech., vol. 46, February 1998.
[12] A. E. Souvorov, A. E. Bulyshev, S. Y Semenonv, R. H. Svenson, A. G. Nazarov, Y E. Sizov,
and G. P. Tatsis, "Microwave tomography: A two-dimensional Newton iterative scheme,"
IEEE Trans. Microw. Theory Tech., vol. 46, November 1998.
[13] 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,
pp. 1841-1853, November 2000.
[14] J. D. Zaeytijd, A. Franchois, C. Eyraud, and J. M. Geffrin, "Full-wave three-dimensional
microwave imaging with a regularized Gauss-Newton method - theory and experiment," IEEE
Trans. Antennas Propag., vol. 55, no. 11, pp. 3279-3292, 2007.
[15] P. Kosmas, J. D. Shea, B. D. Van Veen, and S. C. Hagness, "Three-dimensional microwave
imaging of realistic breast phantom via an inexact Gauss-Newton algorithm," in Proc. Antennas Propag. Society Intl. Symp., July 2008.
[16] D. W. Winters, J. D. Shea, P. Kosmas, B. D. Van Veen, and S. C. Hagness, "Three-dimensional
microwave breast imaging: Dispersive dielectric properties estimation using patient-specific
basis functions," IEEE Trans. Med. Imag., vol. 28, pp. 969-981, July 2009.
[17] X. Li and S. C. Hagness, "A confocal microwave imaging algorithm for breast cancer detection," IEEE Microw. Wireless comp. Lett., vol. 11, March 2001.
[18] 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, 2002.
[19] S. K. Davis, E. J. Bond, X. Li, S. C. Hagness, and B. D. Van Veen, "Microwave imaging
via space-time beamforming for early detection of breast cancer: Beamformer design in the
frequency domain," 7. of Electromagnetic Wave and Applications, vol. 17, no. 2, 2003.
[20] S. K. Davis, H. Tandradinata, S. C. Hagness, and B. D. Van Veen, "Ultrawideband microwave
breast cancer detection: A detection-theoretic approach using the generalized likelihood ratio
test," IEEE Trans. Biomed. Eng., vol. 52, no. 7, 2005.
[21] X. Li, S. K. Davis, S. C. Hagness, D. W. van der Weide, and B. D. Van Veen, "Microwave
imaging via space-time beamforming: experimental investigation of tumor detection in multilayer breast phantom," IEEE Trans. Microw. Theory Tech., vol. 51, pp. 1856-1865, 2004.
88
[22] 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,
November 2005.
[23] A. Wangsanata, "A planar vector network analyzer," Master's thesis, University of WisconsinMadison, 2002.
[24] G. F. Engen, "The six-port reflectometer: An alternative network analyzer," IEEE Trans.
Microw. Theory Tech., vol. 25, pp. 1075-1080, 1977.
[25] G. F. Engen, "A (historical) review of the six-port measurement technique," IEEE Trans.
Microw. Theory Tech., vol. 45, pp. 2414-2417, 1997.
[26] F. R. de Sousa, B. Huyart, and R. C. S. Freire, "Low cost network analyzer using a six-port
reflectometer," in Proc. IEEE MTT-S Intl. Microw. Optoelectronics Conf., 2001.
[27] F. Wiedmann, B. Huyart, E. Bergeault, and L. Jallet, "New structure for a six-port reflectometer in monolithic microwave integrated-circuit technology," IEEE Trans. Instrum. Meas.,
vol. 46, pp. 527-530, 1997.
[28] B. Huyart, F. Wiedmann, L. Jallet, E. Bergeault, R. Benelbar, and R. G. Bosisio, "Microwave
measurement using Wheatsone's bridges," in Proc. European Microw. Conf., 1994.
[29] R. A. Marsland, C. J. Madden, D. W. van der Weide, M. S. Shakouri, and D. M. Bloom,
"Monolithic integrated circuits for mm-wave instrumentation," in Proc. IEEE Annual GaAs
ICSymp., 1990.
[30] A. Samani, J. Bishop, C. Luginbuhl, and D. B. Plewes, "Measuring the elastic modulus of ex
vivo small tissue samples," Phys. Med. Biol., vol. 48, pp. 2183-2198, 2003.
[31] A. Samani, J. Zubovits, and D. Plewes, "Elastic moduli of normal and pathological human
breast tissues: An inversion-technique-based investigation of 169 samples," Phys. Med. Biol,
vol. 52, pp. 1565-1576, January 2007.
[32] R. Sinkus, M. Tanter, T. Xydeas, S. Catheline, J. Bercoff, and M. Fink, "Viscoelastic shear
properties of in vivo breast lesions measured by MR elastography," Magnetic Resonance
Imaging, vol. 23, pp. 159-165, 2005.
[33] A. R Sarvazyan, "Shear acoustic properties of soft biological tissues in medical diagnostics,"
/. Acoust. Soc. Am. Proc, vol. 93, pp. 2329-2330, April 1993.
[34] T. A. Krouskop, T. M. Wheeler, F. Kallel, B. S. Garra, and T. Hall, "Elastic moduli of breast
and prostate tissue under compression," Ultrasound Imaging, vol. 20, pp. 260-274, 1998.
89
[35] A. Thomas, S. Kummel, F. Fritzsche, M. Warm, B. Ebert, B. Hamm, and T. Fischer, "Realtime sonoelastography performed in addition to B-mode ultrasound and mammography: Improved differention of breast lesions?," Acad. Radiology, vol. 13, pp. 1496-1504, December
2006.
[36] E. S. Burnside, T. J. Hall, A. M. Sommer, G. K. Hesley, G. A. Sisney, W. E. Svensson, J. P.
Fine, J. Jiang, and N. J. Hangiandreou, "Differentiating benign from malignant solid breast
masses with US strain imaging," Radiology, vol. 245, no. 2, pp. 401^10, 2007.
[37] R. Sinkus, J. Lorenzen, D. Schrader, M. Lorenzen, M. Dargatz, and D. Holz, "High-resolution
tensor MR elastography for breast tumor detection," Phys. Med. Biol, vol. 45, pp. 1649-1664,
2000.
[38] A. Alizad, M. Fatemi, D. H. Whaley, and J. F. Greenleaf, "Application of vibro-acoustography
for detection of calcified arteries in breast tissue," J. Ultrasound Med., vol. 23, pp. 267-273,
2004.
[39] R. A. Kruger, K. K. Kopecky, A. M. Aisen, D. R. Reinecke, G. A. Kruger, and W. L. Kiser,
"Thermoacoustic CT with radio waves: A medical imaging paradigm," Radiology, vol. 211,
pp. 275-278, 1999.
[40] R. A. Kruger, K. Stantz, and W. L. Kiser, "Thermoacoustic CT of the breast," in Proc. SPIE,
vol. 4682, pp. 521-525, 2002.
[41] G. Ku, B. D. Fomage, X. Jin, M. Xu, K. K. Hunt, and L. V. Wang, "Thermoacoustic and photoacoustic tomography of thick biological tissues toward breast imaging," Technol. Cancer
Res. Treatment, vol. 4, pp. 559-565, 2005.
[42] M. Xu, K. Geng, J. Xing, L. V. Wang, B. D. Fomage, and K. K. Hunt, "Breast cancer imaging
by microwave-induced thermoacoustic tomography," in Proc. SPIE, vol. 5697, pp. 45-48,
2005.
[43] D. M. Pozar, Microwave Engineering. Readings, Mass.: Addison-Wesley, 1990.
[44] D. K. Ghodgaonkar, V. V. Varadan, and V. K. Varadan, "A free-space method for measurement
of dielectric constants and loss tangents at microwave frequencies," IEEE Trans. Instrum.
Meas., vol. 37, pp. 789-793, June 1989.
[45] D. K. Ghodgaonkar, V. V. Varadan, and V. K. Varadan, "Free-space measurement of complex permittivity and complex permeability of magnetic materials at microwave frequencies,"
IEEE Trans. Instrum. Meas., vol. 39, pp. 387-394, April 1990.
[46] N. Gagnon, J. Shaker, L. Roy, A. Petosa, and P. Berini, "Low-cost free-space measurement of
dielectric constant at Ka band," in Proc. IEE Microw. Antennas Propag., vol. 151, pp. 271—
276, June 2004.
90
[47] M. Nakhkash, Y. Huang, W. Al-Nuaimy, and M. Fang, "An improved calibration technique for free-space measurement of complex permittivity," IEEE Trans. Geo. Remote Sens.,
vol. 39, pp. 453-455, February 2001.
[48] F. C. Smith, B. Chambers, and J. C. Bennett, "Calibration techniques for free space reflection coefficient measurements," in Proc. IEE A-Science Meas. Tech., vol. 139, pp. 247-253,
September 1992.
[49] X. Li, S. C. Hagness, M. K. Choi, and D. W. van der Weide, "Numerical and experimental
investigation of an ultrawideband ridged pyramidal horn antenna with curved launching plane
for pulse radiation," IEEE Antennas Wireless Propag. Lett, vol. 2, pp. 259-262, 2003.
[50] K. J. Parker, S. R. Huang, R. A. Musulin, and R. M. Lerner, "Tissue response to mechanical
vibration for sonoelasticity imaging," Ultrasound Med. Biol, vol. 16, pp. 241-246, 1990.
[51] Y. Yamakoshi, J. S ato, and T. Sato, "Ultrasonic imaging of the internal vibration of soft tissue
under force vibration," Ultrasound Ferroelectr. Freq. Control, vol. 37, pp. 45-53, 1990.
[52] R. Muthupillai, D. J. Lomas, P. J. Rossman, J. F. Greenleaf, A. Manduca, and R. L.
Ehman, "Magnetic resonance elastography by direct visulization of propagating acoustic
strain waves," Science, vol. 269, pp. 1854-1857, 1995.
[53] S. Catheline, F. Wu, and M. Fink, "A solution to diffraction biases in sonoelasticity: The
acoustic impulse technique," /. Acoust. Soc. Am., vol. 105, pp. 2941-2950, 1999.
[54] L. Sandrin, M. Tanter, S. Catheline, and M. Fink, "Shear modulus imaging with 2-D transient
elastography," IEEE Trans. Ultrason. Ferroelectr. Freq. Control, vol. 49, pp. 426^435, 2002.
[55] M. Zhao, J. D. Shea, S. C. Hagness, D. W. van der Weide, B. D. Van Veen, and T. Varghese,
"Numerical study of microwave scattering in breast tissue via coupled dielectric and elastic
contrasts," IEEE Antennas Wireless Propag. Lett., vol. 7, pp. 247-250, 2008.
[56] K. R. Nightingale, P. J. Kornguth, and G. E. Trahey, "The use of acoustic streaming in breast
lesion diagnosis: A clinical study," Ultrasound Med. Biol., vol. 25, pp. 75-87, September
1998.
[57] S. Calle, J. P. Remenieras, O. B. Matar, M. E. Hachemi, and F. Patat, "Temporal analysis of
tissue displacement induced by a transient ultrasound radiation force," J. Acoust. Soc. Am.,
vol. 118, no. 5, 2005.
[58] T. Sugimoto, S. Ueha, and K. Itoh, "Tissu hardness measurement using the radiation force of
focused ultrasound," in Proc. IEEE Intl. Ultrasounics Symp., pp. 1377-1380, 1990.
[59] W. F. Walker, "Internal deformation of a uniform elastic solid by acoustic radiation force," J.
Acoust. Soc. Am., vol. 105, pp. 2508-2518, 1999.
91
[60] A. P. Sarvazyan, O. V. Rudenko, S. D. Swanson, J. B. Fowlkes, and S. Y. Emelianov, "Shear
wave elasticity imaging: A new ultrasonic technology of medical diagnostics," Ultrasound
Med. Biol, vol. 24, no. 9, pp. 1419-1435, 1998.
[61] E. A. Barannik, "Acoustic remote palpation: The influence of tissue viscosity on the excitation and relaxation of local impulse shear strain," in Proc. WCU, pp. 7-10, September 2003.
[62] K. Nightingale, M. S. Soo, R. Nightingale, and G. Trahey, "Acoustic radiation force impulse
imaging: in vivo demonstration of clinical feasibility," Ultrasound Med. Biol, vol. 28, no. 2,
pp. 227-235, 2002.
[63] J. Bercoff, M. Tanter, S. Chaffai, and M. Fink, "Ultrafast imaging of beamformed shear wave
induced by the acoustic radiation force," in Proc. IEEE Int. Ultrasonics Symp., 2002.
[64] E. E. Konofagou, M. Ottensmeyer, S. Agabian, S. L. Dawson, and K. Hynynen, "Estimating
localized oscillatory tissue motion for assessment of the underlying mechanical modulus,"
Ultrasonics, vol. 42, pp. 951-956, 2004.
[65] M. Fatemi and J. F. Greenleaf, "Ultrasound-stimulated vibro-acoustic spectrography," Science, vol. 280, pp. 82-85, 1998.
[66] S. Calle, J. P. Remenieras, O. B. Matar, M. Defontaine, and F. Patat, "Application of nonlinear
phenomena induced by focused ultrasound to bone imaging," Ultrasound Med. Biol, vol. 29,
pp. 465^72, 2003.
[67] R. Souchon, R. Salomir, O. Beuf, L. Milot, D. Grenier, D. Lyonnet, J. Y. Chapelon, and
O. Rouviere, "Transient MR elastography (t-MRE) using ultrasound radiation force: Theory,
safety, and initial experiments in vitro," Magn. Res. Med., vol. 60, pp. 871-881, October
2008.
[68] C. J. Diederich and K. Hynynen, "Ultrasound technology for hyperthermia," Ultrasound Med.
Biol, vol. 25, pp. 871-887, July 1999.
[69] E. Kennedy, G. R. T. Haar, and D. Cranston, "High intensity focused ultrasound: Surgery for
the future?," Br. J. Radiol., vol. 76, pp. 590-599, September 2003.
[70] G. ter Haar, "Ultrasound focal beam surgery," Ultrasound Med. Biol, vol. 21, no. 9, pp. 10891100,1995.
[71] M. Kinoshita, N. McDannold, F. A. Jolesz, and K. Hynynen, "Noninvasive localized delivery of herceptin to the mouse brain by MRI-guided focused ultrasound-induced blood-brain
barrier disruption," in Proc. Natl Acad. Sci., vol. 103, pp. 11719-11723, 2006.
[72] F. Y. Yang, W. M. Fu, R. S. Yang, H. C. Liou, K. H. Kang, and W. L. Lin, "Quantitative
evalution of focused ultrasound with a contrast agent on blood brain barrier disruption," Ultrasound Med. Biol, vol. 33, pp. 1421-1427, September 2007.
92
[73] A. Taflove and S. C. Hagness, Computational Electrodynamics: The Finite-Difference TimeDomain Method. 3rd ed. Boston: Artech House, 2005.
[74] J. Van Bladel, "Electromagnetic fields in the presence of rotating bodies," in Proc. IEEE,
vol. 64, pp. 301-318, March 1976.
[75] D. E. Lawrence and K. Sarabandi, "Electromagnetic scattering from vibrating penetrable objects using a general class of time-varying sheet boundary conditions," IEEE Trans. Antennas
Propag., vol. 54, pp. 2054-2061, July 2006.
[76] A. M. Buerkle and K. Sarabandi, "Analysis of acousto-electromagnetic wave interaction using sheet boundary conditions and the finite-difference time-domain method," IEEE Trans.
Antennas Propag., vol. 55, pp. 1991-1998, July 2007.
[77] A. M. Buerkle and K. Sarabandi, "Analysis of acousto-electromagnetic wave interaction
using the finite-difference time-domain method," IEEE Trans. Antennas Propag., vol. 56,
pp. 2191-2199, August 2008.
[78] Comsol Multiphysics, Comsol Multiphysics Acoustic Module - User Guide.
[79] L. E. Kinsler, A. R. Frey, A. B. Coppens, and J. V. Sanders, Fundamentals of Acoustics. 4th
ed. New York: Wiley, 2000.
[80] R. F. Harrington, Time-Harmonic Electromagnetic Fields. New York: McGraw-Hill, 1961.
[81] A. V. Oppenheim and A. S. Willsky, Signals and Systems. 2nd ed. Englewood Cliffs, NJ:
Prentice Hall, 1996.
93
Appendix A: Theoretical Analysis of the Doppler Reconstruction
Method
A theoretical analysis of the Doppler reconstruction method used in the FDTD-SBC implementation is provided here, with the same 1-D example analyzed as in Section 4.1 and 5.1.
Consider that a half space dielectric interface, displacing at an acoustic frequency uim, is illuminated by a continuous plane wave at a microwave frequency UIQ. The scattered microwave field
from the interface can be written as
Es(z,t)
=
ft{re-j20ld(t)e+j0lZejuot}
(A.l)
where F is the reflection coefficient of the dielectric interface, j3\ is the microwave wave number
in the medium, and the interface displacement d(t) is
d{t) = Amsm(umt)
(A.2)
The acoustic-modulated scattered field can be expressed as
Ea(t)\z=0
= r-cos(w0t - 2/Mm am(umt))
(A.3)
If the dielectric interface is unperturbed, then this expression becomes the unperturbed scattered
field
£b(*)U=o = rcos(u;oi)
(A.4)
In FDTD-SBC method, the perturbation scattered component is defined as
E(t) = Es(t) - E0(t)
(A.5)
As in the QS FDTD-SBC implementation described in Section 8.2, FDTD-SBC simulations
are conducted to compute the perturbation components associated with a series of snapshots at r%
(i=l,2,...N) on the acoustic time scale. N is the sampling number within one acoustic period. Then
a two-step demodulation process is employed to reconstruct the Doppler component. The computed perturbation component from the FDTD-SBC simulation is a time-domain pulse response
which corresponds to a wideband signal in the frequency domain, thus taking the Fourier transform
94
of this perturbation component can extract its phasor component at the center microwave frequency
associated with each time snapshot. Next, this phasor component of the perturbation component
which varies on the acoustic time scale is demodulated via Fourier transform to obtain the Doppler
scattering component.
Assume that N snapshots are taken in one acoustic period. The displacement of the dielectric
interface associated with each snapshot Tj is
2m
d(ri) = Amsm(—),i
= l,2,...N
(A.6)
The total acoustic-modulated scattered field associated with the snapshot Tj is then expressed as
2m
Es(t, Ti)\z=0 = rcos(w 0 t - 2f31Amsin(—)),
i = 1,2, ...TV
(A.7)
The perturbation scattered component can then be obtained as
~
E(t,Ti)\z=0
2m
= r c o s ( w o * - 2 / M m S i n ( — ) ) - r C os(w 0 t),* = 1,2,...N
(A.8)
With the assumption that the displacement is much smaller than the microwave wavelength, i.e.
Am -C Ai which is true in the proposed application, the expression is then approximated as
~
E{t,Ti)\z=0
2m
= Tsm(uj0t)2f31Amsm(—),i
= l,2,...N
(A.9)
The phasor of the perturbation component at the central microwave frequency U>Q associated with
each snapshot Tj is extracted via Fourier transform as
27TZ
7T
N
j
E(r i ,w 0 ). = r 2 / 3 i ^ n s i n ( _ ) - - , i = l,2,..JV
(A.10)
Consider constructing an impulse train to represent the perturbation phasor component that varies
on the acoustic time scale
N
„
L
x(t) = r2f31Am(-)J2M^ )
.
+oo
J2 5^-^T--kT^
(A n
- >
where Tm = — is the period of the acoustic excitation. This time-varying perturbation phasor
component on the acoustic time scale is then demodulated via Fourier transform as
X{jku) = T2f31Am(-)J2sin^)uJm
J
i=l
E
e
fc=-oo
^(^)
(AJ2)
95
For the 1 Doppler component, where k = 1 and u = um, the expression above becomes
X(jwm)
=
vr.
T—v . ,27ri, _;2*i
r2/?!Am(-)o;m^sin(—)e~^
i=l
JV/2
iV/2-1
•^
z=l
= r2/3^ m (-)o; m V
=
.
sin(-—)e+^]
= r2/Mm(-)u;m[J>n(—)e"^+ £
i=0
-2|sin(^)|2
-Tl31Am^)umN
J
(A.13)
Therefore, the magnitude of the demodulated term is
\X(jum)\
= r01Am™mN
(A.14)
Recall that for this canonically deforming case, FDTD-SBC simulation for only one snapshot
on the acoustic time scale, i.e. when the dielectric interface is under peak displacement (d(ri) =
Am), needs to be conducted to approximate the magnitude of Doppler component. In that case, the
computed time-domain perturbation component is
E(t,T1)\z=0
= Tsm(uj0t)2p1Am
(A.15)
Based on which, the magnitude of the 1 s t Doppler component is approximated as
|E d (w m )| = r2/3 1 ^ m 7r
(A.16)
Compare Equations A.14 with A.16. It is evident that the expression obtained as in Equation A.14 calibrated by a normalization factor ^ ^ can be adopted to approximate the magnitude
of I s4 Doppler component as in Equation A.16. This factor ^jj
normalizes the metric (Equa-
tion A.14) by accounting for the chosen total sampling number within one acoustic period in the
reconstruction method.
96
Appendix B: List of Authored Publications
Peer-Reviewed Publications
[1] M. Zhao, S. C. Hagness, B. D. Van Veen, and D. W. van der Weide, "Three-dimensional
computational investigation of a focused-acoustic and microwave hybrid modality for breast cancer
detection and diagnosis," IEEE Trans. Biomed. Eng., to be submitted.
[2] M. Zhao, and S. C. Hagness, "Three-dimensional numerical modeling of electromagnetic
scattering from an acoustically excited object using the finite-difference time-domain method,"
IEEE Trans. Antennas Propag., to be submitted.
[3] M. Zhao, J. D. Shea, S. C. Hagness, D. W. van der Weide, and B. D. Van Veen, "Numerical
study of microwave scattering in breast tissue via coupled dielectric and elastic contrasts," IEEE
Antennas Wireless Propag. Lett., vol. 7, pp. 247-250, 2008.
[4] M. Zhao, J. D. Shea, S. C. Hagness, and D. W. van der Weide, "Calibrated free-space
microwave measurements with an ultrawideband reflectometer-antenna system," IEEE Microw.
Wireless Comp. Lett., vol. 16, no. 12, pp. 675-677, December 2006.
[5] M. K. Choi, M. Zhao, S. C. Hagness, and D. W. van der Weide, "Compact mixer-based
1-12 GHz reflectometer," IEEE Microw. Wireless Comp. Lett, vol
15, no. 11, pp. 781-783,
November 2005.
Conference Presentations
[1] M. Zhao*, S. C. Hagness, B. D. Van Veen, and D. W. van der Weide, "Computational
study of a focused acoustic and microwave hybrid sensing modality that exploits coupled dielectric
and elastic properties contrasts," IEEE APS International Symposium and URSI North American
Radio Science Meeting, Charleston, June 1-5, 2009.
[2] M. Zhao*, J. D. Shea, S. C. Hagness, B. D. Van Veen, D. W van der Weide, and T. Varghese, "Three-dimensional computational study of microwave scattering in breast tissue via coupled
dielectric and elastic contrasts," IEEE AP-S International Symposium and URSI North American
Radio Science Meeting, San Diego, July 5 - 1 1 , 2008.
97
[3] S. C. Hagness*, D. W. van der Weide, B. D. Van Veen, M. Zhao, and J. D. Shea"Recent
advances in microwave detection of breast cancer: a compact low-cost refiectometer and enhanced
sensing modality that couples dielectric and elastic properties contrasts," Department of Defense
Breast Cancer Research Program Era of Hope 2008 Meeting, Baltimore, June 25 - 28, 2008.
[4] M. Zhao*, J. D. Shea, S. C. Hagness, D. W. van der Weide, B. D. Van Veen, and T. Varghese,
"Contrast enhancement for breast cancer detection using acoustic-microwave coupling: a preliminary investigation," URSI North American Radio Science Meeting, Ottawa, Canada, July 22 - 26,
2007 (invited talk).
98
Appendix C: Developed Codes for FDTD-SBC Implementation
C.l
Determine the surface current locations
load iv.mat
load jv.mat
load kv.mat
ib=iv+l;
jb=jv+l;
kb=kv+l;
io=iv/2+l/2;
jo=jv/2+l/2;
ko=kv/2+l/2;
r_phy=5e-3;
load dx.mat
% dertermine the boundary voxels:
numVo=0;
iVo= [] ;
jVo=[];
kVo= [] ;
for ii=iv/2-60:iv/2+60
for jj=jv/2-60:jv/2+60
for kk=kv/2-60:kv/2+60
dr=r_phy-sqrt((ii-io)~2+(jj-jo)~2+(kk-ko)~2)*dx;
if dr>0 && dr<dx
numVo=numVo+l;
iVo(immVo)=ii;
jVo(numVo)=jj;
kVo(numVo)=kk;
end
end
end
end
% variables used to determine the surface locations
Dexl=r_phy-sqrt((iVo-io).~2+(jVo-0.5-jo).~2+(kVo-0.5-ko). ~2).*dx;
Dex2=r_phy-sqrt((iVo-io)."2+(jVo+0.5-jo).~2+(kVo-0.5-ko).~2).*dx;
Dex3=r_phy-sqrt((iVo-io).~2+(jVo-0.5-jo).~2+(kVo+0.5-ko)."2).*dx;
99
Dex4=r_phy-sqrt(-(iVo-io) . ~2+(jVo+0.5-jo) .~2+(kVo+0.5-ko)..~2) .*dx;
Deyl=r_phy-sqrt
Dey2=r_phy-sqrt
Dey3=r_phy-sqrt
Dey4=r_phy-sqrt
iVo-0.5-io
iVo+0.5-io
iVo-0.5-io
iVo+0.5-io
~2+(jVo-jo).~2+(kVo-0.5-ko)
~2+(jVo-jo).~2+(kVo~0.5-ko)
~2+(jVo-jo).~2+(kVo+0.5-ko)
~2+(jVo-jo).~2+(kVo+0.5-ko)
"2).*dx;
~2).*dx;
"2).*dx;
"2).*dx;
Dezl=r_phy-sqrt
Dez2=r_phy-sqrt
Dez3=r_phy-sqrt
Dez4=r_phy-sqrt
iVo-0.5-io
iVo+0.5-io
iVo-0.5-io
iVo+0.5-io
*2+(jVo-0.5-jo).~2+(kVo-ko)
*2+(jVo-0.5-jo).~2+(kVo-ko)
•2+(jVo+0.5-jo)."2+(kVo-ko)
'2+(jVo+0.5-jo).~2+(kVo-ko)
"2).*dx;
"2).*dx;
"2).*dx;
~2).*dx;
Dhxl=r_phy-sqrt
Dhx2=r_phy-sqrt
iVo-0.5-io
iVo+0.5-io
'2+(jVo-jo) . ~2+(kVo-ko)~2) *dx;
'2+(jVo-jo).~2+(kVo-ko) '2) *dx;
Dhyl=r_phy-sqrt
Dhy2=r_phy-sqrt
iVo-io).~2+(jVo-0.5-jo).~2+(kVo-ko) ~2) *dx;
iVo-io).~2+(jVo+0.5-jo).~2+(kVo-ko) ~2) *dx;
Dhzl=r_phy-sqrt
Dhz2=r_phy-sqrt
iVo-io).~2+(jVo-jo) '2+(kVo-0.5-ko) *2> *dx;
iVo-io).~2+(jVo-jo) *2+(kVo+0.5-ko) '2) *dx;
% x-dir interior surface locations (e)
iex_intl=iVo(find(Dexl>0))
iex_int2=iVo(find(Dex2>0))
iex_int3=iVo(find(Dex3>0))
iex_int4=iVo(find(Dex4>0))
iex_int=[iex_intl iex_int2 iex_int3 iex_int4];
jex_intl=jVo(find(Dexl>0))-0.5;
jex_int2=jVo(find(Dex2>0))+0.5;
j ex_int3=jVo(find(Dex3>0))-0.5;
j ex_int4=jVo(find(Dex4>0))+0.5;
jex_int=[jex_intl jex_int2 jex_int3 jex_int4];
kex_intl=kVo(find(Dexl>0))-0.5;
kex_int2=kVo(find(Dex2>0))-0.5;
kex_int3=kVo(find(Dex3>0))+0.5;
kex_int4=kVo(find(Dex4>0))+0.5;
kex_int=[kex_intl kex_int2 kex_int3 kex_int4];
% y_dir interior surface locations (e)
iey_intl=iVo(find(Deyl>0))-0.5;
iey_int2=iVo(find(Dey2>0))+0.5;
iey_int3=iVo(find(Dey3>0))-0.5;
iey_int4=iVo(find(Dey4>0))+0.5;
iey_int=[iey_intl iey_int2 iey_int3 iey_int4];
jey_intl=jVo(find(Deyl>0));
jey_int2=jVo(find(Dey2>0>);
jey_int3=jVo(find(Dey3>0));
jey_int4=jVo(find(Dey4>0));
jey_int=[jey_intl jey_int2 jey_int3 jey_int4];
key_intl=kVo(find(Deyl>0))-0.5;
key_int2=kVo(find(Dey2>0))-0.5;
key_int3=kVo(find(Dey3>0))+0.5;
key_int4=kVo(find(Dey4>0))+0.5;
key_int=[key_intl key_int2 key_int3 key_int4];
% z-dir interior surface locations (e)
iez_intl=iVo(find(Dezl>0))-0.5;
iez_int2=iVo(find(Dez2>0))+0.5;
iez_int3=iVo(find(Dez3>0))-0.5;
iez_int4=iVo(find(Dez4>0))+0.5;
iez_int=[iez_intl iez_int2 iez_int3 iez_int4];
jez_intl=jVo(find(Dezl>0))-0.5;
jez_int2=jVo(find(Dez2>0))-0.5;
jez_int3=jVo(find(Dez3>0))+0.5;
jez_int4=jVo(find(Dez4>0))+0.5;
jez_int=[jez_intl jez_int2 jez_int3 jez_int4];
kez_intl=kVo(find(Dezl>0));
kez_int2=kVo(find(Dez2>0));
kez_int3=kVo(find(Dez3>0));
kez_int4=kVo(find(Dez4>0));
kez_int=[kez_intl kez_int2 kez_int3 kez_int4];
'/, x-dir interior surface locations (h)
ihx_intl=iVo(find(Dhxl>0))-0.5;
ihx_int2=iVo(find(Dhx2>0))+0.5;
ihx_int=[ihx_intl ihx_int2];
jhx_intl=jVo(find(Dhxl>0));
jhx_int2=]Vo(find(Dhx2>0));
jhx_int=[jhx_intl jhx_int2];
khx_intl=kVo(find(Dhxl>0));
khx_int2=kVo(find(Dhx2>0));
khx_int=[khx_intl k h x _ i n t 2 ] ;
% y - d i r i n t e r i o r surface l o c a t i o n s (h)
ihy_intl=iVo(find(Dhyl>0));
ihy_int2=iVo(find(Dhy2>0));
ihy_int=[ihy_intl ihy_int2];
jhy_intl=jVo(find(Dhyl>0))-0.5;
jhy_int2=jVo(find(Dhy2>0))+0.5;
jhy_int=[jhy_intl jhy_int2];
khy_intl=kVo(find(Dhyl>0));
khy_int2=kVo(find(Dhy2>0));
khy_int=[khy_intl k h y _ i n t 2 ] ;
% z - d i r i n t e r i o r surface l o c a t i o n s (h)
ihz_intl=iVo(find(Dnzl>0));
ihz_int2=iVo(find(Dhz2>0));
ihz_int=[ihz_intl ihz_int2];
jhz_intl=jVo(find(Dhzl>0));
jhz_int2=jVo(find(Dhz2>0));
jhz_int=[jhz_intl jhz_int2];
khz_intl=kVo(find(Dhzl>0))-0.5;
khz_int2=kVo(find(Dhz2>0))+0.5;
khz_int=[khz_intl k h z _ i n t 2 ] ;
% adjust t o coordinates in FDTD g r i d
% x - d i r i n t e r i o r surface l o c a t i o n s (e)
102
iex_int_grid=iex_int;
jex_int_grid=jex_int+0.5;
kex_int_grid=kex_int+0.5;
% y-dir interior surface locations (e)
iey_int_grid=iey_int+0.5;
jey_int_grid=jey_int;
key_int_grid=key_int+0.5;
% z-dir interior surface locations (e)
iez_int_grid=iez_int+0.5;
j ez_int_grid=j ez_int+0.5;
kez_int_grid=kez_int;
% x-dir interior surface locations (h)
ihx_int_grid=ihx_int+0.5;
jhx_int_grid=jhx_int;
khx_int_grid=khx_int;
% y-dir interior surface locations (h)
ihy_int_grid=ihy_int;
jhy_int_grid=jny_int+0.5;
khy_int_grid=khy_int;
% z-dir interior surface locations (h)
ihz_int_grid=ihz_int;
jhz_int_grid=jhz_int;
khz_int_grid=khz_int+0.5;
% coordinate matrix in FDTD grid
°/0 interial surface location (e)
cex_int_grid=[iex_int_grid;j ex_int_grid;kex_int_grid]
cey_int_grid=[iey_int_grid;jey_int_grid;key_int_grid]
cez_int_grid=[iez_int_grid;jez_int_grid;kez_int_grid]
% interior surface location (h)
chx_int_grid=[ihx_int_grid;jhx_int_grid;khx_int_grid]
chy_int_grid=[ihy_int_grid;jhy_int_grid;khy_int_grid]
chz_int_grid=[ihz_int_grid;jhz_int_grid;khz_int_grid]
C.2
Import the interpolated displacement data from the acoustic-structural
simulation
°/0 import the Comsol file
load fem_app.mat
load fa.mat % acoustic frequency
wm=2*pi*fa;
T=l/fa;
M=12; % total snapshots
taxis=[0:T/N:T] ;
i i = l ; '/, snapshot*
ti=taxis(ii);
% convert surface location coordinates to
% coordinates in acoustic-structural simulations:
°/0 x-dir interior surface location (e)
xex_int=(iex_int-io).*dx+0.02;
yex_int=(jex_int-jo).*dx;
zex_int=(kex_int-ko).*dx;
% y-dir interior surface location (e)
xey_int=(iey_int-io).*dx+0.02;
yey_int=(j ey_int-j o).*dx;
zey_int=(key_int-ko).*dx;
7, z-dir interior surface location (e)
xez_int=(iez_int-io).*dx+0.02;
yez_int=(jez_int-jo).*dx;
zez_int=(kez_int-ko). *dx;
% x-dir interior surface location (h)
xhx_int=(ihx_int-io).*dx+0.02;
yhx_int=(jhx_int-jo).*dx;
zhx_int=(khx_int-ko).*dx;
% y-dir interior surface location (h)
xhy_int=(ihy_int-io).*dx+0.02;
yhy_int=(jhy_int-jo). *dx;
zhy_int=(khy_int-ko). *dx;
% z-dir interior surface location (h)
xhz_int=(ihz_int-io).*dx+0.02;
yhz_int=(jhz_int-jo). *dx;
zhz_int=(khz_int-ko). *dx;
% coordinate matrix for importing data
% from acoustic-structural simulations:
°/0 for interior surface locations (e)
cex_int= [xex_int;yex_int;zex_int];
cey_int=[xey_int;yey_int;zey_int] ;
cez_int= [xez_int;yez_int;zez_int];
% for interior surface locations (h)
chx_int=[xhx_int;yhx_int;zhx_int];
chy_int=[xhy_int;yhy_int;zhy_int];
chz_int=[xhz_int;yhz_int;zhz_int];
% import interpolated displacement
% for x-dir interior surface locations (e)
xm_ex_int=postinterp(fem_app,,u_amp_acsld',cex_int);
ym_ex_int=postinterp(fem_app,,v_amp_acsld,,cex_int);
H-
P
H-
P
J
J
<
B
o
l
P
13
13
CD
Hi
CO
M
^_^
ct
P
H-
X
O
(3*
P. P.
CO
M
P-
s_^
^-^
ct
s^
ct
^_^
ct
P
H-
p"
X
1
«n* *»
o *»
o
CO
M
O
CO
M
n
|
pX PJ p"
X
H- i
1
P H- Hct P
P
o
P- p .
h-
en1
<
P
P
l
1
1
1
1
13
13 pj Q) QJ
B
t3- P" 0
B
1
l
1 3 1 3 13
1
1
P P
n n P P P
P
13
13
B
1
P
13
13
B
B
3
1
p
13
13
1
P
13
13
CD
CD
CD
CD
Hi
Hi
/^^
Hi Hi
/•—N
/-"N
/•-N
H
13
H.
13
/~N
ct ct ct
CD CD CD
H
l-i H
13 1 3 1 3
P
H-
X
1
CJ*
0
x
ct
CD
P
P
H-
X
P
1
1
ct
CD
P
ct
P
H-
P
X
1
B
N
P P p
ct ct ct ct
II
II
II
II
II
1 3 1 3 13 1 3 1 3
O
O
O O
o
CO CO CO CO CO
c t c t c t ct c t
H- M - H - H - H -
-
P
X
1
1
P
X
i
i
X
13
<<:
13
p
-
• ~ N
CO
p
o
H-
p
ct
o
H-
P
Ct
ct
P
• - N
H-
H-
<<
H-
H-
*
CD
!->•
H-
t
h-'
00
H-
*-••
ct
CD
ct
CD
P
P
ct
CD
P
H-
ct
CO
ct
CD
P
H-
ct
CO
ct
CD
P
H-
ct
CO
Ct
CD
f-i
13
P
H-
ct
CO
H-
^ " N
y—\
•X
H-
^ - N
H-
CD
o
ct
P
H-
\_^
ct
P
H-
ct
v
p
H'
P
^^
ct
p
H-
H-
ct
**s
1
CD
i
CD
p.
CO
I-1
ct
^_x
P
H-
i
CD
P.
CO
M
«n "o "n
p-
CO
M
s
^-\
CD
^-^
CO
p
o
H-
ct
M
+ +
-X
•—\
/-N
H-
CD
H-
X
CD
ct
P
H-
X
1
H-
ct
P
H-
H-
ct
P
H-
ct
P
H-
ct
P
H-
ct
P
H-
P
ct
H-
X
1
X
CD
CD
^ <<: <<:
P
ct
H-
CD
X
1
X
^-^
^_1^
X
1
>-L
H-
P
H-
P
P
O
CO
M
P
n
1
1P3 1
1
<!»-
s~^
Hi
CO
M
-
ct
s_x
P
H-
X
1
O
CD
-
^-^
ct
p
H-
ct
'—•
P
H-
CD CD
X
X
1
1
«n «o-
ct
^-^
P
H-
X
CD
i
«n
p.
M
n
CO
v*
E|
1
1
13
P* g
1 13
P 1
n P
P
v*
Hi
• ~ \
+ P- P. P-
13 13
-X
ct
P
H-
1
CO
HM
13 1 3
1
|
X
1
~
E|
•X 1 3
^-N £3-
H-
^-N
^ "^ ^
H-
13
•X
ct
P
H-
X
1
CD
13
1
+ <<S+
N
H-
H- H P P
/~N / " N
Hi Hi
^
CD
1
^-
1
H-
CD CD CD CD
B B B B
1
1
l
l
P P
P
P
1 3 13 1 3 1 3
13 13 13 13
-
N
CD
\ '
H-
P
ct
•
•X
CD
X
13
1
N
CD
o
H-
1
H-
^ - ^ s_^ v_^
*~s
^s
\~s *. '
o
o
H-
B
ct ct ct ct
1
CD CD CD CD
CD H. H. l-i l-i
X 13 1 3 13 1 3
+
X
s; g E|
B
B
B
•X
-X
*
ct ct ct
^"N
H-
/"N
-* *
•~N
H-
13
v_^
1> >
X
13
CD
to> X
•X
P
ct
H-
X
1
1
CD
CD
+
< 1X
Ht o> P
ct
+
El
1
CD
> B"B
to
P
CD
V w /
s_^ v—'
v_^ **s *^x
/
^_x
h-'
00
!-»•
00
1
H-
ct
N^-X
P
H-
1
n
CD
h-1 M
p . P-
-
n
CD
P
Hi
i-i
p
CO
l-i
o
^<!
o o
H-
CD
O
P-
M
<< << ^
•X- •X • X
1 3 13 13
ct
P
P
ct
H-
1
H-
<<!
1
1
1
CD CD CD
+ + X+
N
13 ^
13 13
H-
^ - N
•X-
•—s
/ ^ N
H-
>-i
ct
CD
P
H-
CD
X
P P
P
P
CD c t ri- ri- riCO X
X X
II
ll
ll
ll
1
1
1 3 1 3 1 3 13
•£> 1
HH - H- O
H.
O
O
O
ct P
P P CO CO CO CO
• - v ct
ct ct ct ct ct ct
N
•—N
ll
P
M
^"N
P1
I—
P
ri-
CD
CD
H-
^<!
P
N
N
X
II 1 3 1 3 1 3 B
l-i 1
1
1
CD
CD
CD CD CD
P
X X X
X
H1 1
1
1
I
o
u .
•X-
^ ~ N
)
H.
H
H.
H
+
1 3 1^
3 13 1 3 1 3
< *<:
*<:
<<:
i
i
i
/
—
\
/
—
\
H i H i H i H i /->.
Hi
H - H - H - H i /^\
t o> P P P CD CD CD CD CD CD
ct ct ct B B B B B B
1
1
1
1
£+
*D C*D 13P 13P 1P3 1P3 13P 13p
CD C
NJ> X
X
X 1 3 13 1 3 1 3 1 3 1 3
\_^ 13 13 13
r—\
H - H - H - si
< P El < P
1
1
1
1
1
*
13
PJ 1P3 13
P" 0 P
13 1 3
^ 0^ B^ 1P 1P 1P 13
B
1
1
1
o
o
P
n
* c*t CO CO CO n Po Po
c*t ct
B
1
CD
H-
ct
CO
H-
ct
CO
HH
P.
<<
1
<
II
l-i
P. El
II
M
1
CD
X
!-»•
00 00
H-
P
ct
P
H-
1
N
s_^ *^s
p
ri
H-
1
N
^-^
CD
H-
O
l-i
Hi
5<
00
N^-X
ct
P
H-
I
N
«o »n* «n•
CD
CD
CD
hp.
CO
p
o
H-
ct
P
o
n
h-1
B
1
ct
P
+ + +
N
ct
P
ct
P
to l
1
CD (D
p
^ - N
H-
H-
w
3
1
CD CD
^1
l-J-
CD
N
ct
M
n
n
CO CO
1
B
13
1
P
13
1
P
1
QJ
B
P
n
CD
P
Hi
p
CO
H
O
^
H-
CD
P
HCt
H-
H-
H-
CD
B
N
<i <: ^1 *<:
•<
i
1 •<!
1 ^
1
^
N
X
1 3 13 1 3
1
1
1
CD CD CD
P
P
P P P P
ri- ri- ri- ri- ri- rill
ll
ll
ll
ll
ll
1 3 1 3 13 13 1 3 1 3
H- O
O
O O O O
y - >
N ^-\ X
CD ^
CD CD
<<! ^1 ^<
/^-s
II
l-i
CD
P1
h-
p
00
(-*•
N
CD
N
1
ct
B
1
p
13
13
CD
/*-^
Hi
CD
l-i
13
P
ri
H-
ct
i
P
13
13
P- p .
CO
M
B
l
0
Hi
/ ^ N
13
H
ct
CD
P
H-
ct
CO
f-i
ct
CO
•P
ll
P
ri-
H
f-i
CD
P1
h-
CD
P
M
<II
II
E:
00
P
H-
c t Ct
v—/ \~s
P
H-
13
1
B
1
PJ
P
13
13
n P
en o
h-1 M
p . P.
CO
o
B
l
CD
Hi
P
H-
ct
CO
HH
P.
N
1
H-
^1
p.
1
CD
^. -^ ~\
I-
v_^
^-x
1
X
CD
N
P
ct
v—'
P
H-
P.
CO
M
o
<
SI
P
1
1
1
1 3 1 3 13
P* tr P*
i
P
P P
B
1
P
13
13
CD
Hi
P
H-
ct
CO
CO
H-
O
M
Hi
o^
^^
*
ct
P
H-
o o o
00
•X
ct
P
H-
1 3 T) 1 3
H- H- H-
*
ct
P
H-
+ + +
H-
g
B
•if
£l
B
•X
ct
ct
/
N
H-
• - N
H-
P
13
CD
B
Hi
X 1*o
3 13
X)
* *
H-
•X
CD
X
B
1
PJ
CD
Hi
P
H-
ct
CO
ct c t ct c t
CD CD CD CD
H
H
H
H
1 3 13 13 1 3
/ — N s~^ ^ " N s~\
P
H-
ct
CO
N
X
13 13 13
1
1
1
o o n
CD CD CD CD CD CD
N
N
N N N N
1
1
1
1
1
H-
ct
•X
B
g
• - N
•X
13
h- 1
O
*
CD
to> X
£l
N
ct
P
«<!
^ - N
^
X
N
X
N
II
1 3 13 B
B
l-i 1
1
1
l
1
1
CD CD CD CD CD CD CD
N
N
N
N
P
N
N
1
1
1
1
p
N
X
P P P P P P
CD CD CD c t ri- ri- ri- ri- riN
N
N
II
ll
ll
ll
ll
ll
1
1
1 1 3 13 13 1 3 1 3 13
H- HH- O
O
O O
O O
• ~ \
H
CD
P>
<II
+ + X+
B
B
t o> l
^
CD CD CD
+
< 1N 1N 1N
HH- HIO) P
p
P
+ c t ct c t
P
^ ~ N
ct
H
CO
•T!
ct
II
P
P>
CD
O
CD
P>
H
Hi
P
CO
H.
O
I-!
H-
ct
CD
P
H-
HH
X
1
p.
l-i
o
^
II
p. E|
Hi
1
CD
N
1
H-
5-s
oca
o
zp_hx_int=postinterp(fem_app,,w_ph_acsld',chx_int);
u=real(xhx_int+xm_hx_int.*exp(i.*(wm*ti+xp_hx_int.*pi./180)))-xhx_int;
v=real(yhx_int+ym_hx_int.*exp(i.*(wm*ti+yp_hx_int.*pi./180)))-yhx_int;
w=real(zhx_int+zm_hx_iiit.*exp(i.*(wm*ti+zp_hx_int.*pi./180)))-zhx_int;
d_hx_int=sqrt(u.~2+v."2+w.~2);
°/0 for y-dir interior surface locations (h)
xm_hy_int=postinterp(fem_app,'u_amp_acsld',chy_int);
ym_hy_int=postinterp(fem_app,'v_amp_acsld',chy_int);
zm_hy_int=postinterp(fem_app)'w_amp_acsld',chy_int);
xp_hy_int=postinterp(fem_app,'u_ph_acsld',chy_int)
yp_hy_int=postinterp(fem_app,'v_ph_acsld',chy_int)
zp_hy_int=postinterp(fem_app,'w_ph_acsld',chy_int)
u=real(xhy_int+xm_hy_int.*exp(i.*(wm*ti+xp_hy_int.*pi./180)))-xhy_int;
v=real(yhy_int+ym_hy_int.*exp(i.*(wm*ti+yp_hy_int.*pi./180)))-yhy_int;
w=real(zhy_int+zm_hy_int.*exp(i.*(wm*ti+zp_hy_int.*pi./180)))-zhy_int;
d_hy_int=sqrt(u.~2+v.~2+w.~2);
% for z-dir interior surface locations (h)
xm_hz_int=postinterp(fem_app,,u_a^lp_acsld,,chz_int);
ym_,hz_int=postinterp(fem_app, 'v.amp.acsld' ,chz_int) ;
zm_hz_int=postinterp(fem_app,,w_amp_acsld',chz_int);
xp_hz_int=postinterp(fem_app,'u.ph.acsld',chz_int)
yp_hz_int=postinterp(fem_app,'v.ph.acsld',chz_int)
zp_hz_int=postinterp(fem_app,'w.ph.acsld',chz_int)
u=real(xhz_int+xm_hz_int.*exp(i.*(wm*ti+xp_hz_int.*pi./180)))-xhz_int;
v=real(yhz_int+ym_hz_int.*exp(i.*(um*ti+yp_hz_int.*pi./180)))-yhz_int;
w=real(zhz_int+zm_hz_int.*exp(i.*(wm*ti+zp_hz_int.*pi./180)))-zhz_int;
d_hz_int=sqrt(u.~2+v.~2+w.~2);
C.3
Create 'pointers' for checking the neighboring E fields on the same
of the boundary
% for x-dir interior location (e)
n_ex_int=size(iex_int,2);
In_ex_int=zeros(n_ex_int,8);
for kk=l:n ex int
% neighboring Ey
if (sqrt((iex_int(kk)-0.5-io)~2+(jex_int(kk)-0.5-jo)~2
+(kex_int(kk)-ko)~2)*dx)<r_phy
In_ex_int(kk,l)=l;
end
if (sqrt((iex_int(kk)+0.5-io)~2+(jex_int(kk)-0.5-jo)~2
+(kex_int(kk)-ko)"2)*dx)<r_phy
In_ex_int(kk,2)=1;
end
if (sqrt((iex_int(kk)-0.5-io)~2+(jex_int(kk)+0.5-jo)~2
+(kex.int(kk)-ko)~2)*dx)<r_phy
In_ex_int(kk,3)=1;
end
if (sqrt((iex_int(kk)+0.5-io)~2+(jex_int(kk)+0.5-jo)~2
+(kex_int(kk)-ko)"2)*dx)<r_phy
In_ex_int(kk,4)=l;
end
% neighboring Ez
if (sqrt((iex_int(kk)-0.5-io)~2+(jex_int(kk)-jo)~2
+(kex_int(kk)+0.5-ko)"2)*dx)<r_phy
In_ex_int(kk,5)=1;
end
if (sqrt((iex_int(kk)+0.5-io)~2+(jex_int(kk)-jo)~2
+(kex_int(kk)+0.5-ko)"2)*dx)<r_phy
In_ex_int(kk,6)=1;
end
if (sqrt((iex_int(kk)-0.5-io)~2+(jex_int(kk)-jo)~2
+(kex_int(kk)-0.5-ko)~2)*dx)<r_phy
In_ex_int(kk,7)=1;
end
if (sqrt((iex_int(kk)+0.5-io)~2+(jex_int(kk)-jo)~2
+(kex_int(kk)-0.5-ko)"2)*dx)<r_phy
In_ex_int(kk,8)=1;
end
end
'/, for y-dir interior location (e)
...
...
...
...
...
...
...
...
n_ey_int=size(iey_int,2);
In_ey_int=zeros(n_ey_int,8);
for kk=l:n_ey_int
% neighboring Ex
if (sqrt((iey_int(kk)-0.5-io)~2+(jey_int(kk)-0.5-jo)~2
+(key_int(kk)-ko)~2)*dx)<r_phy
In_ey_int(kk,1)=1;
end
if (sqrt((iey_int(kk)-0.5-io)~2+(jey_int(kk)+0.5-jo)~2
+(key_int(kk)-ko)~2)*dx)<r_phy
In_ey_int(kk,2)=l;
end
if (sqrt((iey_int(kk)+0.5-io)~2+(jey_int(kk)-0.5-jo)"2
+(key_int(kk)-ko)~2)*dx)<r_phy
In_ey_int(kk,3)=l;
end
if (sqrt((iey_int(kk)+0.5-io)~2+(jey_int(kk)+0.5-jo)~2
+(key_int(kk)-ko)"2)*dx)<r_phy
In_ey_int(kk,4)=l;
end
% neighboring Ez
if (sqrt((iey_int(kk)-io)~2+(jey_int(kk)-0.5-jo)~2
+(key_int(kk)+0.5-ko)~2)*dx)<r_phy
In_ey_int(kk,5)=1;
end
if (sqrt((iey_int(kk)-io)"2+(jey_int(kk)+0.5-jo)~2
+(key_int(kk)+0.5-ko)"2)*dx)<r_phy
In_ey_int(kk,6)=1;
end
if (sqrt((iey_int(kk)-io)~2+(jey_int(kk)-0.5-jo)~2
+(key_int(kk)-0.5-ko)"2)*dx)<r_phy
In_ey_int(kk,7)=1;
end
if (sqrt((iey_int(kk)-io)"2+(jey_int(kk)+0.5-jo)"2
+(key_int(kk)-0.5-ko)~2)*dx)<r_phy
In_ey_int(kk,8)=l;
end
end
...
...
...
...
...
...
...
...
/o for z-dir interior surface location (e)
n_ez_int=size(iez_int,2);
In_ez_int=zeros(n_ez_int,8);
for kk=l:n_ez_int
% neighboring Ex
if (sqrt((iez_int(kk)-0.5-io)~2+(jez_int(kk)-jo)~2
+(kez_int(kk)-0.5-ko)~2)*dx)<r_phy
In_ez_int(kk,1)=1;
end
if (sqrt((iez_int(kk)-0.5-io)~2+(jez_int(kk)-jo)~2
+(kez_int(kk)+0.5-ko)~2)*dx)<r_phy
In_ez_int(kk,2)=l;
end
if (sqrt((iez_int(kk)+0.5-io)~2+(jez_int(kk)-jo)~2
+(kez.int(kk)-0.5-ko)~2)*dx)<r_phy
In_ez_int(kk,3)=1;
end
if (sqrt((iez_int(kk)+0.5-io)~2+(jez_int(kk)-jo)~2
+(kez_int(kk)+0.5-ko)"2)*dx)<r_phy
In_ez_int(kk,4)=l;
end
...
...
...
...
% neighboring Ey
if (sqrt((iez_int(kk)-io)~2+(jez_int(kk)-0.5-jo)~2
+(kez_int(kk)-0.5-ko)"2)*dx)<r_phy
In_ez_int(kk,5)=1;
end
if (sqrt((iez_int(kk)-io)"2+(jez_int(kk)-0.5-jo)"2
+(kez_int(kk)+0.5-ko)"2)*dx)<r_phy
In_ez_int(kk,6)=1;
end
if (sqrt((iez_int(kk)-io)~2+(jez_int(kk)+0.5-jo)~2
+(kez_int(kk)-0.5-ko)"2)*dx)<r_phy
In_ez_int(kk,7)=1;
end
if (sqrt((iez_int(kk)-io)"2+(jez_int(kk)+0.5-jo)"2
+(kez_int(kk)+0.5-ko)"2)*dx)<r_phy
...
...
...
...
In_ez_int(kk,8)=1;
end
end
C.4 Calculate the electric surface currents Js
7„ the electric surface current calculation is done based on the
% 3-D FDTD code (Fortran) initially developed by Henri Tandradinata
ccc
ccc
ccc
ccc
local coordinates of the surface currents: iObsP,jObsP.kObsP
displacement imported from the acoustic-structural simulation: delta
object radius: a
'pointers' for checking neighbouring E filed: sidDerm
do n2=l,obsInPx
i=i0bsPx(n2)
j=j0bsPx(n2)
k=k0bsPx(n2)
cnt=sidDermx(n2,l)+sidDermx(n2,2)+sidDermx(n2,3)
1
+sidDermx(n2,4)+sidDermx(n2,5)+sidDermx(n2,6)
2
+sidDermx(n2,7)+sidDermx(n2,8)
1
2
3
4
5
6
7
8
9
imJxl = deltax(n2)*(epsStatM(l)-epsStatM(2))
*(sidDermx(n2Jl)*(Ey(i,j-l,k)-Eyold(i)j-l,k) )
+sidDermx(n2,2)*(Ey(i+l)j-l,k)-Eyold(i+l)j-l,k))
+sidDermx(n2,3)*(Ey(i,j,k)-Eyold(i,j,k))
+sidDermx(n2,4)*(Ey(i+l,j,k)-Eyold(i+l,j,k))
+sidDermx(n2,5)*(Ez(i)j,k)-Ezold(i)j,k))
+sidDermx(n2,6)*(Ez(i+l,j,k)-Ezold(i+l,j,k))
+sidDermx(n2,7)*(Ez(i,j,k-l)-Ezold(i,j,k-l))
+sidDermx(n2,8)*(Ez(i+l,j,k-i)-Ezold(i+l,j,k-l)))
/dt/cnt
1
2
imJx2=deltax(n2)*(sigM(l)-sigM(2))
*(sidDermx(n2,l)*Ey(i,j-l,k)
+sidDermx(n2,2)*Ey(i+l)j-l)k)+sidDermx(n2,3)*Ey(i,j,k)
H-
H-
H-
H-
<<
to
C_|
B
!->•
«<!
+
H-
C-i
0
CD
M
ct
ct
s_^
h-*
h-1
w
u
H
B
V_^
(_1.
s~\
P
/ - N
*<:
B
CD
i-i
a
1
C _ l .
• - N
v_^
M-
1
"sv
w
(_1.
w
s-\
H-
c_l.
H-
/ - N
* •*
m
N m
N
>—•
(V
<_i.
H-
/'—v
*
m
X
s—'
*oo
•
-J «
cn t o
w
P p
to to
^~N
<< <<:
H
B
v-^
I-"-
l_j.
/—N
•X-
v—'
^^
•X- ^^ ^^
* m
* m
tn m
N N X X
r~\
H- /->
H- H- M1
*
l_i. <
!-»•
+ + I+- ' **
«IV «IV « h-+
w
1
«SV
+
+ \—>
'—/ CH-O co
H+
cn
\n a- O.
O Ha
P CD CD P.
M-
'
H-
*«IV
t_l.
u
^-^
\ s
to
•~\
s
W
H-
h-*
|
|
CO
v-^
h-1
• " N
(TO
3
H-
/•"N
*
m
X
V
< • <
P- P- P- P.
O
PaD O
C
CD CD a
CD O P
H M H H CD
B B B B M /•"N
• • <
B P
<<!
^ " N
•""N <
•<
" ^: ^
/—N
to
/—N
P P p P ^
t o tO NO t o P ^-^
*
v
•
*00
* cn <.£.
* «t o t o C~ON
«
CO
H-
ct
H-
M-
M-
M-
M-
M-
H-
CO
ct
P
\n
to
• ™ \
^<
P p
to to
/—N
^
^- /
p
^_^
s—^
!->•
s—'
1
•™N
v
IV '
s
w
u
^ - N
^ - N
/—N
«
^_x
V_X
v—^
u
|_l
^-x
i
S
^-N
|-l
v_x
V '
v_x
l-i
^
^
IV
»
C_j.
»
v_^
^-^
S
^-\
tO
O c+
M CD
P, ct
p. ^
l«!
O
M
i m cn
in x w
tv
!->•
»
l_l.
-
H- CD
1 ct
I V ^-^ CD
^-^
i *d
-
-
U .
+
l-i.
1
I-1
^-^ HI V t - J . IV H1
"
/-N
1 W 1
PI w m
N 1 X
O ITI o
(-• >>! M
p . O P.
M
!->•
!V
+
t—i. (_t.
^-N
\ S
(_i.
s~\
•—N
P. (-"• P. H-
s_^
1
< _ l .
^ - N
s-^
^-\
H- ( _ J . <_!.
1
«
** + ^ +
(_!. tv1 h-»^-^ (_l.
I+
-' I!V
IV
+
«!V ^-^
1
**IV
p.
M
o
s_^
p.
(_i.
^ - N
^~N
+ ^tv h-+
«IV M-1 -IV
1
I - ' ^-^
1
1
tn
1 N m
N
o o
m
N M M
I-L
(_i.
^~N
^ ^
I-'
-X^-^ ^—'
* * * * * * ^*> ^* ^^
CD
^~\ ^^
m ^n3
m
m
m N
m m
N N
N X m
x m en
><! ^-s
X
W
^-\
H- H- H- H- H- /"^
H- H - / ^ c t
*~s
<<
P
ct
CD
M
P
^~s • - ^
to p
P
<<:
B
^ " N ^
tfs. 00 t o 00 ~J 05 ^cn
^ - N^^- N^^-
•™N
•—\
to to to
s~\
<<! ^ ^<!
p P ^
p~> P
<<j
P. P. P- P. P- P- P. HO p.
aD CaD a
o CDa CaD CD
D
C
CD CD
f-i
CD
^
M f-i
H
^
^
B B B B B B B 3
II
\p. co
+ W+ en
+ w+ w+ C+O ra
+ ^ *N
to
II
p.
*
•—\
M-
II
! - l
B
C-l
B
H-
C-l
H-
«)(XISO)01^UIOK
B
+ OT
+ C+O Cfl
+
co
M-
o
H-
oo to
p.
1=1
p.
CD
en ^
S~\
H-
P.
+O
C
H-
p
to
w
* ^ N
/•"N
P.
+
CO
H-
s—>
**_•
p
to
^<
B
CD
a
P.
CO
H-
+
00
P
to
CD
p P o
t o st. o l-i
s*
B
00 cn <<!
y—N
M
B
^ <<:
M
B
a
CD a
CD
P. P. t o
H-
\_/
+ C+O
w
>—'
N—S
to to
S^*
T) T)
I-*
s_>-
CO
I-'
P
II t o
HII
H-
O
P-
v_^
<<:
t o 13
-o
1
^ aCO
<< *< ^<! h-1
P P /~\
P P
CO
CO
<_1.
II II
IV C _ i .
o O O
a- cr cr
tv
P
to
^ - N
oD
P ^
p—\ C
t o t o l-i
<-^1
• -rf^ B^
v_/
/—N
B
l-i
^ <<:
B
l-i
a
CD o
CD
P- P-
H-
+ C+O
CO
H-
<<:
B
CD
t-i
o
P-
CO
H-
ct
II
o
1
to I- p
P.
O
P
P.
CD
B
X
to
C-l
B
+
H-
M
!-»•
<-i
H-
B
H-
\~^
Ji
v_x
to
t-i
B
X
+
+
(V IV
*
N
*
N
P P
to to
^i cn
a
M- + w+
^^ cn Hp. P.
O
CD CD
•-i Hi
B B><!
IV IV IV
+
* m*
t*n tn
N N ^<!
»*_•
to
P.
O
CD
H-
^-\
P ^-^
P
B
X
M
CD
00 0 )
p
to
•~\
B
X
CD
l-i
H-
P- P.
O O
+ CO+ CO+
CO
H-
cn rf=> oo
do n2=l,obsInPz
i=i0bsPz(n2)
j=j0bsPz(n2)
k=k0bsPz(n2)
cnt=sidDermz(n2,1)+sidDermz(n2,2)+sidDermz(n2,3)
1
+sidDermz(n2,4)+sidDermz(n2,5)+sidDermz(n2,6)
2
+sidDermz(n2,7)+sidDermz(n2,8)
1
2
3
4
5
imJzl = deltaz(n2)*(epsStatM(l)-epsStatM(2))
*(sidDermz(n2,l)*(Ex(i-lJj,k)-Exold(i-l,j,k) )
+sidDermz(n2,2)*(Ex(i-l,j,k+l)-Exold(i-l,j,k+l))
+sidDermz(n2,3)*(Ex(i,j,k)-Exold(i,j,k))
+sidDermz(n2,4)*(Ex(i, j ,k+l)-Exold(i,j,k+l))
+sidDermz(n2,5)*(Ey(i,j-l,k)-Eyold(i,j-l,k))
6
7
8
9
+sidDermz(n2)6)*(Ey(i)j-l,k+l)-Eyold(i,j-l»^+l))
+sidDermz(n2,7)*(Ey(i)j»k)-Eyold(i)j,k))
+sidDermz(n2)8)*(Ey(i,j)k+l)-Eyold(i,j,k+l)))
/dt/cnt
1
2
3
4
imJz2=deltaz(n2)*(sigM(l)-sigM(2))
*(sidDermz(n2,l)*Ex(i-l,j,k)
+sidDermz(n2,2)*Ex(i-1,j,k+1)+sidDermz(n2,3)*Ex(i,j,k)
+sidDermz(n2,4)*Ex(i)j)k+l)+sidDermz(n2,5)*Ey(i,j-l,k)
+sidDermz(n2,6)*Ey(i,j-l,k+l)+sidDermz(n2)7)*Ey(i)j,k)
5
+sidDermz(n2,8)*Ey(i,j,k+1))/cnt
imJz=imJzl+imJz2
end do
C.5
Determine the index matrix of the neighboring H surface fields
% spherical coordinates of the surface field locations
% interior surface hx
phi_hx_int=acos((ihx_int-io)./sqrt((ihx_int-io)."2 ...
+(jhx_int-jo)."2)); % in rad.
phi_hx_int(find((ihx_int-jo)<0))=2*pi-phi_hx_int(find((ihx_int-jo)<0));
the_hx_int=acos((khx_int-ko)./sqrt((ihx_int-io). ~2 ...
+(jhx_int-jo).~2+(khx_int-ko)."2));
% interior surface hy
phi_hy_int=acos((ihy_int-io)./sqrt((ihy_int-io)."2 ...
+(jhy_int-jo).~2)); % in rad.
phi_hy_int(f ind((jhy_int-j o)<0))=2*pi-phi_hy_int(f ind((j hy_int-jo)<0));
the_hy_int=acos((khy_int-ko)./sqrt((ihy_int-io)."2 ...
+(jhy_int-jo).~2+(khy_int-ko)."2));
% interior surface hz
phi_hz_int=acos((ihz_int-io)./sqrt((ihz_int-io)."2 ...
+(jhz_int-jo).~2)); % in rad.
phi_hz_int(find((jhz_int-jo)<0))=2*pi-phi_hz_int(find((jhz_int-jo)<0));
the_hz_int=acos((khz_int-ko)./sqrt((ihz_int-io)."2 ...
+(jhz_int-jo).~2+(khz_int-ko).~2));
% determine index matrix of neighboring H surface field:
% for interior surface hx
n_hx_int=size(ihx_int,2);
Ihx_int_nb= [] ;
Ihx_int_nb_l=[] ;
Ihx_int_nb_2=[] ;
for kk=l:n_hx_int
de=(ihx_int-ihx_int(kk).*ones(l,n_hx_int))."2+(jhx_int-jhx_int(kk)
.*ones(l,n_hx_int)).~2+(khx_int-khx_int(kk).*ones(l,n_hx_int))."2;
% avoid choose neighboring locations -> dthe=0 or dphi=0
de(find(the_hx_int==the_hx_int(kk)))=le9;
de(find(phi_hx_int==phi_hx_int(kk)))=le9;
[Y,I]=sort(de,2);
Ihx_int_nb(kk,l)=I(2);
for tmp=3:n_hx_int
if the_hx_int(I(tmp))~=the_hx_int(I(2)) && ...
phi_hx_int(I(tmp))~=phi_hx_int(I(2))
if (the_hx_int(I(tmp))-the_hx_int(kk))*(the_hx_int(I(2)) ...
-the_hx_int(kk))<0 && (phi_hx_int(I(tmp))-phi_hx_int(kk)) ..
*(phi_hx_int(1(2))-phi_hx_int(kk))<0
Ihx_int_nb(kk)2)=I(tmp);
break
end
end
end
13=1(tmp);
de=(ihx_int-ihx_int(I(2)).*ones(1,n_hx_int)).~2+(jhx_int-jhx_int(I(2))
.*ones(l,n_hx_int)).~2+(khx_int-khx_int(1(2)).*ones(l,n_hx_int))."2;
de(find(the_hx_int==the_hx_int(I(2))))=le9;
de(find(phi_hx_int==phi_hx_int(I(2))))=le9;
[Y_l,I_l]=sort(de,2);
Ihx_int_nb_l(kk,l)=I_l(2);
for tmp=3:n_hx_int
if the_hx_int(I_l(tmp))~=the_hx_int(I_l(2)) && ...
phi_hx_int(I_l(tmp))~=phi_hx_int(I_l(2))
if (the_hx_int(I_l(tmp))-the_hx_int(I(2)))*(the_hx_int(I_l(2)) ...
-the_hx_int(I(2)))<0 && (phi_hx_int(I_l(tmp))-phi_hx_int(1(2))) ...
*(phi_hx_int(I_l(2))-phi_hx_int(I(2)))<0
Ihx_int_nb_l(kk,2)=I_l(tmp);
break
end
end
end
de=(ihx_int-ihx_int(I3). *ones(l,n_hx_int)).~2+(jhx_int-jhx_int(I3)
.*ones(l,n_hx_int)).~2+(khx_int-khx_int(13).*ones(l,n_hx_int)). ~2;
de(find(the_hx_int==the_hx_int(I3)))=le9;
de(find(phi_hx_int==phi_hx_int(I3)))=le9;
[de_2,I_2]=sort(de,2);
Ihx_int_nb_2(kk,l)=I_2(2);
for tmp=3:n_hx_int
if the_hx_int(I_2(tmp))~=the_hx_int(I_2(2)) && . . .
phi_hx_int(I_2(tmp))~=phi_hx_int(I_2(2))
if (the_hx_int(I_2(tmp))-the_hx_int(I3))*(the_hx_int(I_2(2)) ...
-the_hx_int(I3))<0 kk (phi_hx_int(I_2(tmp))-phi_hx_int(I3)) ...
*(phi_hx_int(I_2(2))-phi_hx_int(I3))<0
Ihx_int_nb_2(kk,2)=I_2(tmp);
break
end
end
end
end
% for interior surface hy
n_hy_int=size(ihy_int,2);
Ihy_int_nb= [];
Ihy_int_nb_l=[];
Ihy_int_nb_2=[];
for kk=l:n_hy_int
de=(ihy_int-ihy_int(kk).*ones(l,n_hy_int)).~2+(jhy_int-jhy_int(kk)
.*ones(l,n_hy_int)).~2+(khy_int-khy_int(kk).*ones(l,n_hy_int)).~2;
% avoid choose neighboring locations -> dthe=0 or dphi=0
de(find(the_hy_int==the_hy_int(kk)))=le9;
de(find(phi_hy_int==phi_hy_int(kk)))=le9;
[Y,I]=sort(de,2);
Ihy_int_nb(kk,l)=I(2);
for tmp=3:n_hy_int
if the_hy_int(I(tmp))~=the_hy_int(I(2)) kk . . .
phi_hy_int(I(tmp))~=phi_hy_int(I(2))
if (the_hy_int(I(tmp))-the_hy_int(kk))*(the_hy_int(I(2)) ...
-the_hy_int(kk))<0 kk (phi_hy_int(l(tmp))-phi_hy_int(kk)) ...
*(phi_hy_int(I(2))-phi_hy_int(kk)) <0
Ihy_int_nb(kk,2)=I(tmp);
break
end
end
end
13=1(tmp);
de=(ihy_int-ihy_int(I(2)).*ones(1,n_hy_int)).~2+(jhy_int-jhy_int(I(2))
.*ones(1,n_hy_int))."2+(khy_int-khy_int(I(2)).*ones(1,n_hy_int))."2;
de(find(the_hy_int==the_hy_int(I(2))))=le9;
de(find(phi_hy_int==phi_hy_int(I(2))))=le9;
[Y_l,I_l]=sort(de,2);
Ihy_int_nb_l(kk,l)=I_l(2);
for tmp=3:n_hy_int
if theJhy_int(I_l(tmp))~=the_hy_int(I_l(2)) kk . . .
phi_hy_int(I_l(tmp))~=phi_hy_int(I_l(2))
if (the_hy_int(I_l(tmp))-the_hy_int(I(2)))*(the_hy_int(I_l(2)) ...
-the_hy_int(I(2)))<0 && (phi_hy_int(I_l(tmp))-phi_hy_int(1(2))) ...
*(phi_hy_int(I_l(2))-phi_hy_int(I(2)))<0
Ihy_int_nb_l(kk,2)=I_l(tmp);
break
end
end
end
de=(ihy_int-ihy_int(I3).*ones(l,n_hy_int)).~2+(jhy_int-jhy_int(I3) ...
.*ones(l,ii_hy_int)).~2+(khy_int-khy_int(I3).*ones(l,n_hy_int)).~2;
de(find(the_hy_int==the_hy_int(I3)))=le9;
de(find(phijiy_int==phi_hy_int(13)))=le9;
[de_2,1_2]=sort(de,2);
Ihy_int_nb_2(kk,l)=I_2(2);
for tmp=3:n_hy_int
if the_hy_int(I_2(tmp))~=the_hy_int(I_2(2)) kk . . .
phi_hy_int(I_2(tmp))~=phi_hy_int(I_2(2))
if (the_hy_int(I_2(tmp))-the_hy_int(I3))*(the_hy_int(I_2(2)) ...
-the_hy_int(I3))<0 && (phi_hy_int(I_2(tmp))-phi_hy_int(I3)) ...
*(phi_hy_int(I_2(2))-phi_hy_int(13))<0
Ihy_int_nb_2(kk,2)=I_2(tmp);
break
end
end
end
end
% for interior surface hz
;
n_hz_int=size(ihz_int,2);
Ihz_int_nb=[] ;
Ihz_int_nb_l=[] ;
Ihz_int_nb_2=[] ;
IpO=[];
Ippi=[];
for kk=l:n_hz_int
de=(ihz_int-ihz_int(kk).*ones(l,n_hz_int)).~2+(jhz_int-jhz_int(kk) ...
.*ones(l,n_hz_int)).~2+(khz_int-khz_int(kk).*ones(l,n_hz_int))."2;
if the_hz_int(kk)==0 I I the_hz_int(kk)==pi ^special location index
if the_hz_int(kk)==0
IpO=[IpO kk];
phi_hz_int(kk)=le9;
else
Ippi=[Ippi kk];
phi_hz_int(kk)=le9;
end
% avoid choose neighboring locations -> dthe=0 or dphi=0
de(find(the_hz_int==the_hz_int(kk)))=le9;
de(find(phi_hz_int==phi_hz_int(kk)))=le9;
de(find(the_hz_int==0))=le9; °/„ avoid these two locations
de(find(the_hz_int==pi))=le9;
[Y,I]=sort(de,2);
Ihz_int_nb(kk,l)=I(2);
for tmp=3:n_hz_int
if abs(phi_hz_int(I(tmp))-phi_hz_int(I(2)))>pi/2 kk
abs(phi_hz_int(I(tmp))-phi_hz_int(I(2)))<pi
Ihz_int_nb(kk,2)=1(tmp);
break
end
...
end
else
% avoid choose neighboring locations -> dthe=0 or dphi=0
de(find(the_hz_int==the_hz_int(kk)))=le9;
de(find(phi_hz_int==phi_hz_int(kk)))=le9;
de(find(the_hz_int==0))=le9; % avoid these two locations
de(find(the_hz_int==pi))=le9;
[Y,I]=sort(de,2);
Ihz_int_nb(kk,l)=K2);
for tmp=3:n_hz_int
if the_hz_int(I(tmp))~=the_hz_int(I(2)) kk . . .
phi_hz_int(I(tmp))~=phi_hz_int(I(2))
if (the_hz_int(I(tmp))-the_hz_int(kk))*(the_hz_int(1(2))
-the_hz_int(kk))<0 kk (phi_hz_int(I(tmp))-phi_hz_int(kk))
*(phi_hz_int(1(2))-phi_hz_int(kk))<0
Ihz_int_nb(kk,2)=I(tmp);
break
end
end
end
end
13=1(tmp);
de=(ihz_int-ihz_int(I(2)).*ones(l,n_hz_int)).~2+(jhz_int-jhz_int(1(2))
.*ones(l,n_hz_int)).~2+(khz_int-khz_int(1(2)).*ones(l,n_hz_int))."2;
de(find(the_hz_int==the_hz_int(I(2))))=le9;
de(find(phi_hz_int==phi_hz_int(I(2))))=le9;
de(find(the_hz_int==0))=le9; % avoid these two locations
de(find(the_hz_int==pi))=le9;
[Y_l,I_l]=sort(de,2);
Ihz_int_nb_l(kk,l)=I_l(2);
for tmp=3:n_hz_int
if the_hz_int(I_l(tmp))~=the_hz_int(I_l(2)) kk ...
phi_hz_int(I_l(tmp))~=phi_hz_int(I_l(2))
if (the_hz_int(I_l(tmp))-the_hz_int(I(2)))*(the_hz_int(I_l(2)) ...
-the_hz_int(I(2)))<0 && (phi_hz_int(I_l(tmp))-phi_hz_int(1(2))) ..
*(phi_hz_int(I_l(2))-phi_hz_int(I(2))) <0
Ihz_int_nb_l(kk,2)=I_l(tmp);
break
end
end
end
de=(ihz_int-ihz_int(I3).*ones(l,n_hz_int)).~2+(jhz_int-jhz_int(I3)
.*ones(1,n_hz_int))."2+(khz_int-khz_int(13).*ones(1,n_hz_int))."2;
de(f ind(the_hz_int==the_hz_int(13)))=le9;
de(find(phi_hz_int==phi_hz_int(I3)))=le9;
de(find(the_hz_int==0))=le9: % avoid these two locations
...
de(find(the_hz_int==pi))=le9;
[de_2,1_2] =sort(de,2);
Ihz_int_nb=2(kk,l)=I_2(2);
for tmp=3:n_hz_int
if theJiz_int(I_2(tmp))~=the_hz_int(I_2(2)) ...
&& phi_hz_int(I_2(tmp))~=phi_hz_int(I_2(2))
if (the_hz_int(I_2(tmp))-the_hz_int(I3))*(the_hz_int(I_2(2)) ...
-thejiz_int(13))<0 && (phi_hz_int(I_2(tmp))-phi_hz_int(I3)) ...
*(phi_hz_int(I_2(2))-phi_hz_int(13))<0
Ihz_int_nb_2(kk,2)=I_2(tmp);
break
end
end
end
end
% reference: spectial locations:
chz_IpO=chz_int_grid(:,IpO);
chz_Ippi=chz_int_grid(:,Ippi);
Ihz_int_nbpO=Ihz_int_nb(:,IpO);
Ihz_int_nbppi=Ihz_int_nb(:,Ippi);
Ihz_int_nbpO_l=Ihz_int_nb_l(:,IpO);
Ihz_int_nbppi_l=Ihz_int_nb_l(:,Ippi);
Ihz_int_nbpO_2=Ihz_int_nb_2(:,IpO);
Ihz_int_nbppi_2=Ihz_int_nb_2(:,Ippi);
C.6 Calculate the magnetic surface currents Ms
% the surface derivative of H field is done based on the
% 3-D FDTD code (Fortran) initially developed by Henri Tandradinata
ccc
ccc
ccc
ccc
local coordinates of the surface currents: iObsP,jObsP,kObsP
displacement imported from the acoustic-structural simulation: delta
object radius: a
index matrix of the neighboring H surface fields: Inb, Inbl, Inb2
ccc for x-component (interior)
do n2=l,obsInPx
i=i0bsPx(n2)
j=j0bsPx(n2)
k=k0bsPx(n2)
inbla=i0bsPx(Inblx(n2,D)
jnbla=j0bsPx(Inblx(n2,l))
knbla=k0bsPx(Inblx(n2,l))
inblb=i0bsPx(Inblx(n2,2))
jnblb=j0bsPx(Inblx(n2,2))
knblb=k0bsPx(Inblx(n2,2))
inb2a=i0bsPx(Inb2x(n2,l))
jnb2a=j0bsPx(Inb2x(n2,l))
knb2a=k0bsPx(Inb2x(n2,1))
inb2b=i0bsPx(Inb2x(n2,2))
jnb2b=j0bsPx(Inb2x(n2,2))
knb2b=k0bsPx(Inb2x(n2,2))
1
2
d_dsl=l/2*(Hx(inbla,jnbla,knbla)-Hx(inblb,jnblb)knblb))
/(l/2*a*(theMx(Inblx(n2,l))-theMx(Inblx(n2,2)))
*(l/2*a*(phiMx(Inblx(n2,l))-phiMx(Inblx(n2,2)))))
d_ds2=l/2*(Hx(inb2a,jnb2a,knb2a)-Hx(inb2b)jnb2b,knb2b))
1
/(l/2*a*(theMx(Inb2x(n2,l))-theMx(Inb2x(n2,2)))
2
*(l/2*a*(phiMx(Inb2x(n2,l))-phiMx(Inb2x(n2,2)))))
1
g_x=l/2*(d_dsl-d_ds2)/(l/2*a*(theMx(Inbx(n2,l))
-theMx(Inbx(n2,2)))*(l/2*a*(phiMx(Inbx(n2,l))
2
-phiMx(Inbx(n2,2)))))*deltax(n2)
end do
ccc for y-component (interior)
do n2=l,obsInPy
i=i0bsPy(n2)
j=j0bsPy(n2)
122
k=k0bsPy(n2)
inbla=i0bsPy(Inbly(n2,l))
jnbla=j0bsPy(Inbly(n2,l))
knbla=k0bsPy(Inbly(n2,l))
inblb=iObsPy(Inbly(n2,2))
jnblb=j0bsPy(Inbly(n2,2))
knblb=k0bsPy(Inbly(n2,2))
inb2a=i0bsPy(Inb2y(n2,l))
jnb2a=j0bsPy(Inb2y(n2,l))
knb2a=k0bsPy(Inb2y(n2,1))
inb2b=i0bsPy(Inb2y(n2,2))
jnb2b=j0bsPy(Inb2y(n2,2))
knb2b=k0bsPy(Inb2y(n2,2))
1
2
d_dsl=l/2*(Hy(inblaJnbla,knbla)-Hy(inblb,jnblb,knblb))
/(l/2*a*(theMy(Inbly(n2,l))-theMy(Inbly(n2,2)))
*(l/2*a*(phiMy(Inbly(n2,l))-phiMy(Inbly(n2,2)))))
1
2
d_ds2=l/2*(Hy(inb2a)jnb2a)knb2a)-Hy(inb2b)jnb2b)knb2b))
/(l/2*a*(theMy(Inb2y(n2,l))-theMy(Inb2y(n2,2)))
*(l/2*a*(phiMy(Inb2y(n2,1))-phiMy(Inb2y(n2,2)))))
1
2
g_y=l/2*(d_dsl-d_ds2)/(l/2*a*(theMy(Inby(n2,l))
-theMy(Inby(n2,2)))*(l/2*a*(phiMy(Inby(n2,1))
-phiMy(Inby(n2,2)))))*deltay(n2)
end do
ccc for z-component (interior)
c expect for special locations:
do n2=l,obsInPz
i=i0bsPz(n2)
j=j0bsPz(n2)
k=k0bsPz(n2)
inbla=i0bsPz(Inblz(n2,l)>
123
jnbla=j0bsPz(Inblz(n2,l))
knbia=k0bsPz(Inblz(n2,l))
inblb=i0bsPz(Inblz(n2,2))
jnblb=j0bsPz(Inblz(n2,2))
knblb=k0bsPz(Inblz(n2,2))
inb2a=i0bsPz(Inb2z(n2,D)
jnb2a=j0bsPz(Inb2z(n2,l))
knb2a=k0bsPz(Inb2z(n2,1))
inb2b=i0bsPz(Inb2z(n2,2))
jnb2b=j0bsPz(Inb2z(n2,2))
knb2b=k0bsPz(Inb2z(n2,2))
1
2
d_dsl=l/2*(Hz(inbla,jnbla,knbla)-Hz(inblb,jnblb,knblb))
/(l/2*a*(theMz(Inblz(n2,l))-theMz(Inblz(n2,2)))
*(l/2*a*(phiMz(Inblz(n2,l))-phiMz(Inblz(n2,2)))))
1
2
d_ds2=l/2*(Hz(inb2a)jnb2a)knb2a)-Hz(inb2b)jnb2b,knb2b))
/(l/2*a*(theMz(Inb2z(n2,l))-theMz(Inb2z(n2,2)))
*(l/2*a*(phiMz(Inb2z(n2,l))-phiMz(Inb2z(n2,2)))))
1
2
g_z=l/2*(d_dsl-d_ds2)/(l/2*a*(theMz(Inbz(n2,1))
-theMz(Inbz(n2,2)))*(l/2*a*(phiMz(Inbz(n2,l))
-phiMz(Inbz(n2,2)))))*deltaz(n2)
end do
c for special locations:
do n2=l,obsInPzpO
i=i0bsPzp0(n2)
j=j0bsPzp0(n2)
k=k0bsPzp0(n2)
inbla=iObsP(InblzpO(n2,1))
jnbla=j0bsP(Inblzp0(n2,l))
knbla=kObsP(InblzpO(n2,l))
inblb=iObsP(InblzpO(n2,2))
jnblb=jObsP(InblzpO(n2,2))
knblb=k0bsP(Inblzp0(n2,2))
inb2a=i0bsP(Inb2zp0(n2,1))
j nb2a=j ObsP(Inb2zp0(n2,1))
knb2a=k0bsP(Inb2zp0(n2,1))
inb2b=i0bsP(Inb2zp0(n2,2))
jnb2b=jObsP(Inb2zp0(n2,2))
knb2b=k0bsP(Inb2zp0(n2,2))
1
2
d_dsl=l/2*(Hz(inbla,jnbla,knbla)-Hz(inblb,jnblb,knblb))
/(l/2*a*(theMz(InblzpO(n2,l))-theMz(InblzpO(n2,2)))
*(l/2*a*(phiMz(InblzpO(n2,l))-phtMz(InblzpO(n2,2)))))
1
2
d_ds2=l/2*(Hz(inb2a,jnb2a,knb2a)-Hz(inb2b,j nb2b,knb2b))
/(l/2*a*(theMz(Inb2zp0(n2,l))-theMz(Inb2zp0(n2,-2)))
*(l/2*a*(phiMz(Inb2zp0(n2,l))-phiMz(Inb2zp0(n2,2)))))
g_z=l/2*(d_dsl-d_ds2)/(l/2*(l/2*a*theMz(InbzpO(n2,l))
1
*(l/2*a*theMz(InbzpO(n2,l))))*abs(sin(phiMz(InbzpO(n2,l))
2
-phiMz(InbzpO(n2,2)))))*deltaz(n2)
end do
do n2=l,obsInPzppi
i=i0bsPzppi(n2)
j=j0bsPzppi(n2)
k=kObsPzppi(n2)
inbla=iObsP(Inblzppi(n2,1))
jnbla=j ObsP(Inblzppi(n2,1))
knbla=kObsP(Inblzppi(n2,1))
inblb=iObsP(Inblzppi(n2,2))
jnblb=j ObsP(Inblzppi(n2,2))
knblb=kObsP(Inblzppi(n2,2))
inb2a=i0bsP(Inb2zppi(n2,1))
jnb2a=j ObsP(Inb2zppi(n2,1))
knb2a=k0bsP(Inb2zppi(n2,1))
inb2b=i0bsP(Inb2zppi(n2,2))
125
jnb2b=jObsP(Inb2zppi(n2,2))
knb2b=k0bsP(Inb2zppi(n2,2))
1
2
d_dsl=l/2*(Hz(inbla,jnbla)knbla)-Hz(inblb)jnblb,knblb))
/(l/2*a*(theMz(Inblzppi(n2,l))-theMz(Inblzppi(n2,2)))
*(l/2*a*(phiMz(Inblzppi(n2,1))-phiMz(Inblzppi(n2,2)))))
1
2
d_ds2=l/2*(Hz(inb2aJnb2a,knb2a)-Hz(inb2b,jnb2b,knb2b))
/(l/2*a*(theMz(Inb2zppi(n2,1))-theMz(Inb2zppi(n2,2)))
* (l/2*a* (phiMz (Inb2zppi (n2,1)) -phiMz (Inb2zppi (n2,2) ).)))
g_z=l/2*(d_dsl-d_ds2)/(1/2*(l/2*a*(pi-theMz(Inbzppi(n2,1)))
1
*(l/2*a*(pi-theMz(Inbzppi(n2,1)))))*abs(sin(phiMz(Inbzppi(n2,1))
2
-phiMz(Inbzppi(n2,2)))))*deltaz(n2)
end do
% the magnetic surface current is then calculated in matlab:
% parameters
load epsl.mat
load eps2.mat
load sigl.mat
load sig2.mat
load dx.mat
load dt.mat
load nL.mat
load nR.mat
% for x-component (interor)
for nn=nL:nR
temp=zeros(l,n_hx_int);
for mm=l:nn-l
f_x=((l/epsi)*exp(-sigl*(dt*(nn-l-mm))/epsl)-(l/eps2)
*exp(-sig2*(dt*(nn-l-mm))/eps2)).*ones(1,n_hx_int);
Mx(nn,:)=temp-dt.*(mx_l(mm,:)-mx_r(mm,:)).*f_x;
temp=Mx(nn,:);
end
end
...
keyboard
% for y-component (interor)
for nn=nL:nR
temp=zeros(l,n_hy_int);
for mm=l:nn-l
f_y=((l/epsl)*exp(-sigl*(dt*(nn-l-mm))/epsl)-(l/eps2)
*exp(-sig2*(dt*(nn-l-mm))/eps2)).*ones(1,n_hy_int);
My(nn,:)=temp-dt.*(my_1(mm,:)-my_r(mm,:)).*f_y;
temp=My(nn,:);
end
...
end
% for z-component (interor)
for nn=nL:nR
temp=zeros(l,n_hz_int);
for mm=l:nn-l
f_z=((l/epsl)*exp(-sigl*(dt*(nn-l-mm))/epsl)-(l/eps2)
*exp(-sig2*(dt*(nn-l-mm))/eps2)).*ones(l,n_hz_int);
Mz(nn,:)=temp-dt.*(mz_l(mm):)-mz_r(mm,:)).*f_z;
temp=Mz(nn,:);
end
end
...
Документ
Категория
Без категории
Просмотров
0
Размер файла
5 172 Кб
Теги
sdewsdweddes
1/--страниц
Пожаловаться на содержимое документа