close

Вход

Забыли?

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

?

Microwave imaging for two-dimensional electrical property distribution profiling

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