Theoretical Premises and Contemporary Optimizations of Microwave Tomography

Abstract

Microwave imaging has long been proposed as an effective means for biomedical applications—breast cancer detection and therapy monitoring being the most prominent because of the endogenous dielectric property contrast between malignant and normal breast tissue. While numerous numerical simulations have been presented demonstrating feasibility, translation to actual physical and clinical implementations have been lacking. In contrast, the Dartmouth team has taken somewhat counterintuitive but fundamentals-based approaches to the problem—primarily addressing the confounding multipath signal corruption problem and exploiting core concepts from the parameter estimation community. In so doing, we have configured a unique system design that is a synergism of both the hardware and software worlds. In this paper, we describe our approaches in the context of competing strategies and suggest rationales for why these techniques work—especially in 2D. Finally, we present data from actual neoadjuvant chemotherapy exams that confirm that our technique is capable of imaging the tumor and also visualizing its progression during treatment.

1. Introduction

The primary basis for microwave breast imaging is that there is considerable dielectric property contrast between malignant and normal tissue. This contrast is a complicated issue and its understanding has evolved substantially over time. At the most simplistic level, microwaves are especially useful as “water detectors.” At low microwave frequencies, water typically has relative permittivity values between 75 and 80 while that for fat ranges between 5 to 10. The earliest assumptions about breast tissue were that it was primarily fatty and tumors contained far more water on account of their rapid replication and proliferation [1]. However, recent studies have provided a more nuanced appreciation of normal breast tissue composition. In a study by Woodard and White [2], in the context of assessing the physical content of tissue for radiographic purposes, they found that the two constituents of breast tissue—adipose and mammary tissue—varied considerably from woman to woman, and concomitantly, so did their water content. In fact, the water content variations were: adipose—11.4–30.5%, and mammary—30.2–72.6%. More recent studies have tried to provide a more complex presentation of the properties. Studies by Lazebnik et al. [3], Sugitani et al. [4], Martellosio et al. [5], and Cheng and Fu et al. [6] have all presented values for all three tissue types referred to as adipose, fibroglandular and tumor, respectively. Relatively consistently, the studies have typically shown the adipose properties to be quite low, those for the fibroglandular to be substantially elevated and those for the tumors to be the highest.

However, there is considerable variation between results. As has been pointed out by Meaney et al. and [7, 8] Salahuddin et al. [9], there are important weaknesses and even flaws in the methodologies of these studies which could easily skew the desired results. Two of the more critical problems revolve around the use of the ubiquitous open-ended dielectric probe and the frequency sampling regimen used in algorithms for fitting the data to established Cole-Cole curves for broadband parameter estimation [3]. With regards to the former, it is well known that the penetration depth or sampling volume in front of the dielectric probe is on the order of 1/6th that of the probe diameter [10, 11]. Given that the probe diameters of the most common commercial probe (Keysight Technologies, Santa Clara, CA) are on the order of about 2 mm diameter, this suggests that the penetration depth is on the order of 0.3 mm. However, various reports explicitly state the assumed sampling volume to be 3 mm deep or more [12]. The significance is that it is relatively straightforward for a pathologist to perform a tissue analysis over a sample that is 3 mm thick. For a volume as small as 0.3 mm, the analysis is substantially less informative. In these cases, the probe measurements are really only assessing superficial dielectric properties. In addition to compositional issues near the surface, this is also where the greatest temperature (cooling after excision) and moisture (drying after excision) gradients appear especially during measurements of excised tissue. Both temperature and water content can have dramatic influence over the tissue properties. As a whole, these factors are not adequately addressed in recent reports and open the door for substantial variation.

With respect to the frequency sampling regimens, this is primarily a problem related to how the properties are fitted to the Cole-Cole equation to determine broadband coefficients. The norm for this is to acquire data in a logarithmic frequency fashion and was established by Cole and Cole [13]. The reason for this is that when the dielectric property data is plotted in the Cole-Cole plane, the sample points form nearly a circle and the coefficients are derived from the circle’s features such as, the radius and center location. When the data is collected in a logarithmic fashion, the points end up being spaced nearly evenly around the arc of the circle, while for the linearly spaced frequencies, the data points are highly concentrated in a localized corner of the plot. For fitting purposes, the evenly spread out data points produce more accurate findings. The problem becomes further compounded when the data spans multiple relaxation zones (the most prominent cross over point for tissue properties occurs within the range of about 1.5–2.5 GHz) [14]. When a linear sampling regimen is used, the lower relaxation zone is dramatically under-sampled compared to that for the higher zone such that the Cole-Cole fitting is grossly skewed to bias the results for the higher zone. While many of the reports properly sampled the data logarithmically with respect to frequency, some have not [3].

More recently there have been several clinical breast imaging studies [15, 16, 17]. The Preece study was based on the University of Bristol radar imaging technique which only produces intensity maps, where localized hot spots are indicative of tumors. Even so, the clinical data suggests that their system is capable of distinguishing tumor from normal tissue, even in dense breasts suggesting that there must be property contrast between the tumor and fibroglandular tissue. The reports by Poplack et al. [16] and Meaney et al. [17] suggest that there is contrast in both the permittivity and conductivity between malignant and normal tissue. The more interesting aspect is that the contrast is only statistically significant for the conductivity images. Given that there are questions related to the absolute dielectric property levels of the different tissue types for all ex vivo studies, we have attempted to develop a hypothesis based on previous literature. In particular, extensive studies by Foster et al. showed distinctly different tissue permittivity and conductivity relationships as functions of water content [18]. In these situations, the fat and fibroglandular tissue would have either no water or low levels of free water (i.e. the bound water in the fibroglandular tissue is primarily bound to long sugars and proteins). In these situations, the ionic content of the solutions is not able to exhibit higher conductivity because there is little to no free water necessary for conduction. However, the tumor has considerable free water and exhibits high conductivity. This mechanism is one means for explaining the relatively high conductivity contrast between the tumor and normal tissue. Conversely, the permittivity acts more as a measure of bulk water. In this situation, the fibroglandular tissue has relatively high water content compared with fat, but not as high as that for tumor. But the percentage difference between that for tumor and fibroglandular is typically in the range of 10–20% compared with what could amount to factors of 2× or more for the conductivity. This can explain the more subdued permittivity contrast we observe for the imaging studies between normal and malignant tumor. These features need to be explored further—preferably in clinical studies.

As was alluded to above, there are multiple near field microwave imaging techniques. These consist primarily of (a) radar techniques, (b) tomography or inverse problem techniques, (c) thermoacoustic imaging, and (d) holography. The radar techniques are generally a type of synthetic aperture radar (SAR) methods that typically utilize either backscatter or transmission data. For the backscatter approaches, a large amount of wideband data is acquired for many positions around the object surface and then time delays are synthetically added or subtracted from each measurement to focus sequentially at each pixel within the imaging domain. The contributions from each measurement are summed at each pixel and the resulting intensity maps are displayed as the images [19]. Multiple simulation efforts have been developed with the most advanced phantom and patient experiments performed by the Fear group at the University of Calgary [20]. Transmission techniques have also been developed—primarily by the group at the University of Bristol [21]. Theirs utilizes a fixed array of wideband antennas which directly contact the breast—albeit they utilize a coupling gel to enhance coupling and minimize unwanted contributions from multipath signals. This has advanced from simulation and phantom experiments to more extensive clinical trials [15].

Tomographic and inverse problems have been studied extensively in simulation with only very few translating to phantom and clinical work. These approaches typically utilize mostly transmission data and require nonlinear inverse algorithms to produce actual maps of the tissue permittivity and conductivity [22]. Simulation work includes studies by Rocca et al. [23], Fhager and Persson [24], Catapano et al. [25] and Shea et al. [26]. Extensive phantom and ex vivo animal studies have been performed by Semenov et al. [27]. The most comprehensive clinical work has been performed by Meaney et al. which includes studies of normal patients [28, 29], a diagnostic study comparing images of patients with and without tumors [16], and a study monitoring the progression of tumors during neoadjuvant chemotherapy [17]. For these studies, the imaging technique can distinguish tumors from normal tissue and benign lesions to a level of significance for lesions 1 cm and larger. For the therapy monitoring study, the technique was able to determine whether a tumor was responding to treatment or not within the first 30 days of the chemotherapy regimen.

The thermoacoustic techniques generally apply a low duty cycle, high power microwave pulse which is selectively absorbed by the malignant tissue and subsequently causes a mechanical vibration which can be detected by ultrasound transducers. The images are produced by synthetically combining the signals from the different ultrasound transducers. The technique has seen limited success in both phantom and clinical work [30]. The holography approach is being studied primarily at McMaster University by Dr. Nikolova. To date, they have developed an initial prototype which has produced promising phantom results [31].

For this chapter, we focus on tomographic or inverse problem approaches. Microwave tomography and/or microwave inverse problems have now been studied at great length for several decades [25, 26, 32]. The preponderance of efforts has been in simulation with very few advancing to actual implementation and/or clinical exams. While factors such as cost, exam time and image reconstruction complexity/costs are often cited as prime reasons for failure, our experience has led us to focus on two factors that inhibit progress. These are the need to explicitly contend with the problems of multi-path signal corruption [33] and the need for variance stabilizing transformations in the reconstruction process [34]. Viewing microwave imaging from the context of these two challenges clearly illustrates flaws in conventional approaches and inhibits overall progress.

Multi-path signal corruption has been acknowledged for decades. In the context of classic radar and telecommunications applications, signals such as ground clutter constructively and destructively add to the desired signals and excite unwanted artifacts such as ghosting [35]. For far-field applications such as radar, these artifacts can often be as little as a minor nuisance. However, for near field imaging situations, the corruption can be extreme—to the degree that the multi-path signal can easily completely overwhelm the desired one [33]. For far-field situations, the primary mechanism for multi-paths are reflections off of neighboring structures or surfaces and recombining with the original signal [36]. However, for near field cases, they very often propagate as surface waves along interfaces of support structures and the coupling medium or along the outside of antenna feedlines. One likely reason why this phenomenon goes unconsidered is the fact that these structures are simply not included in the models for numerical simulations [26, 37]. For many implementations, the computational costs are already enormous when simply including just the antennas, coupling medium and the target. Adding complexities such as the feedlines, support structure and the coupling medium tank would simply overwhelm the capabilities of modern computers. Consequently, for most simulation efforts, these structures are simply ignored—in fact, even neighboring antennas are usually also eliminated in the name of computational speed and cost. The unfortunate result is essentially a precise rendering of an unrealistic scenario.

There are few options in compensating for this challenge. The primary factor is that the multi-path signals originate from the same desired signal—i.e. it is the same frequency. Because of this, sophisticated filtering approaches are not effective. Techniques such as time-gating have been proposed in different implementations but have not resulted in any published results for microwave tomography [38]. One of the more challenging aspects of time gating is the need for a very broad band signal with fine sampling between frequencies. Given that the measured microwave data is most often acquired in the frequency domain, this can add dramatically to the acquisition time. In addition, many of the proposed antennas simply do not operate over sufficient bandwidth to make this possible.

One technique employed by the Dartmouth group is the use of a lossy coupling medium [39]. This poses unique challenges, but when considering the imaging problem in a comprehensive manner, there is considerable merit to it. The main drawback is that propagating across even a short span in a lossy imaging medium can easily exceed the measurement dynamic range of high quality commercial vector network analyzers (VNA). This is not a trivial concern, but when properly understood, it is possible to devise systems to adequately accommodate this. Dynamic range considerations are discussed in more detail in the Methods section. The primary benefit is that the highly attenuating medium dramatically suppresses the unwanted surface waves. Figure 1 shows a set of simulations for a monopole antenna radiating into a coupling medium where the active part of the antenna is positioned several centimeters above a Plexiglas plate [40]. For the low attenuating liquids (low conductivity), the surface waves easily reach the Plexiglas via coaxial modes traveling on the outsides of the feedlines. Once sufficiently powerful signals reach the low loss plate, planar modes are excited and the waves propagate unimpeded everywhere. However, as the conductivity of the liquid is increased, the coaxial modes are dramatically reduced to the point that for the last example, no signal reaches the Plexiglas. In this instance, the desired signal still propagates in a well-behaved beam pattern from the active part of the monopole antenna. In effect, we have traded a nearly impossible problem—i.e. uncontrolled propagation from damaging multi-path signals—for the need for a high dynamic range VNA which is just a difficult problem. It should be noted that multi-port VNA’s with the necessary dynamic range are now commercially available from some vendors—albeit at a significant cost.

Figure 1.

Simulations of the field patterns and associated surface waves for differing coupling bath conductivities (S/m) (Permission granted to re-print images by Human Press [40]).

As a side note, in explicitly dealing with multi-path issue, we have been able to realize several opportunities that would not have been possible had we taken a more conventionally intuitive approach. For instance, we have found that the monopole antennas are ideal in this setting. First, the naturally occurring resistive loading of the lossy bath dramatically improves the monopole antenna bandwidth [41]. These can be used from roughly 500 MHz to 3 GHz with a −10 dB return loss across the band—well in excess of that for most conventional antennas. While these antennas are essentially isotropic radiators, their low profile allows them to be packed tightly around the target which dramatically reduces the propagation distances compared with more conventional aperture antennas. This shortening of the propagation distance easily compensates for the signal loss suffered due to having lower directivity. Ultimately, the close packing feature is a substantial advantage over the loss of directivity. In addition, because of the lossy bath, antenna low physical profile and the isotropic radiation pattern, there is essentially no mutual coupling between antennas even when spaced as closely as 2 cm apart [42]. Mutual coupling can be significantly debilitating for antenna array performance. Finally, the combination of the low loss medium and the low profile antennas makes the broadcast waves appear to be propagating in a mostly dielectric medium. This allows us to exploit the discrete dipole approximation (DDA) as an efficient means for simulating the signal forward solution which is recognized as the largest time cost in the reconstruction process [43]. We are now able to recover 2D images in 6 s and estimate that fully 3D images can be reconstructed in a few minutes—all without the aid of parallel processors and graphical processing units (GPUs). Overall, the advantages of using a lossy coupling medium have led to dramatic innovations.

Finally, while much has been written about the mathematics and algorithms regarding microwave tomography and/or inverse problems, the most common approach is a non-linear, iterative one which broadly falls under the category of multi-parameter estimation problems. These problems have been studied extensively within the probability and statistics community for which a host of definitions and techniques have been developed to optimize algorithmic performance. In particular, work by Box and Cox in the 1950s and 1960s devised ways to assess the performance of different single-step and iterative algorithms and derived a wide range of suitable transformations to improve performance [44]. Their focus centered around characterizing problems where the data was inherently heteroscedastic and developing transformations to make the data more homoscedastic—i.e. amenable to standard, least squares multi-parameter estimation techniques. Two of the more ubiquitous examples include the log transforms used in X-ray CT and optical coherence tomography (OCT) [45, 46]. The basic assumptions are that the error function difference between the measured and computed field values should have a zero mean and a normal distribution. A convenient way to test this is to simply examine the residual vector after the reconstruction [47]. For the X-ray CT case, the image reconstructions are simply not possible without the transformation. This is also the case with OCT, for which the log transform is now widely adopted [48].

In earlier work, we demonstrated that the residual data for the microwave case was highly heteroscedastic when applied to an algorithm operating directly on the complex field data (Figure 2a) [34]. Our interpretation was that it was largely due to the wide signal strength dynamic range of the field values. However, once we applied the log transformation, the residuals were significantly more normal with a near zero mean (Figure 2b). The main challenge here is that when taking the log transform of a complex number, the result is the log magnitude and its phase [49]. The phase term immediately implies that there could be some form of unwrapping necessary. For the X-ray CT case, the detected signals are all real numbers, so there is no phase term. In OCT, the governing equation for the light is the transport equation and the phase is generated by harmonically modulating the light with a 100 MHz signal [46]. The major point here is that the wavelength associated with the 100 MHz is very large such that there is never enough scattering to generate phase changes greater than +/− 180 degrees. In effect, the data is always unwrapped. However, for the microwave breast imaging case, the dielectric scatterers are often physically on the same order size as the wavelengths and the scattered fields can quite frequently change phase by values well in excess of +/− 180 degrees. For the parameter estimation problems utilizing complex data, it is critical that all of the measured and computed phase pairs be on the same Riemann sheet for the algorithms to work properly. We briefly describe methods for unwrapping the phase in Section 2.2 which are covered in more detail in Meaney et al. [50].

Figure 2.

Histograms of the residual error vector after reconstruction for the (a) non-transformed and (b) log transformed algorithms (Permission granted by Wiley Publishers to re-print graphs [34]).

Operationally, the key is that when starting at a baseline of a homogeneous medium, adding a contrasting target to the imaging zone can change the phases measured or computed for the different transmit/receive pairs substantially—often easily exceeding the single Riemann sheet bounds of −180 to +180 degrees. Depending on the antenna orientations, there are nominally some cases where the measurement and computed phases are already on the same Riemann sheet at the start of the algorithm, but there can also be many that are not. In cases where they are not, by mapping all of the measured and computed phases to the baseline Riemann sheet—i.e. −180 to +180 degrees—(which is effectively what happens when using the non-transformed algorithm), it essentially transforms that element into a “bad” data point. While not universally viewed in this context, the various solutions that have been proposed generally work with the net effect of forcing the starting measured and computed phase values to be on the same Riemann sheet. These include: (1) introducing a priori information [24], (2) frequency hopping [51], and (3) simply adding more data [25, 27]. A priori information works because it effectively assumes a sample target at the start of the image reconstruction process that has similar characteristics to the actual target. The desired result is that it generally positions the measured and computed phases for all or most transmit/receive pairs on the same Riemann sheet. There are several reports utilizing this technique in simulation that appear promising; however, these techniques have generally not advanced beyond simulation studies. These results unfortunately end up being biased based on the quality of the initial guess. Frequency hopping acts in a similar manner as for applying a priori information in that images are reconstructed at progressively increasing frequencies with the results from lower frequencies used as starting guesses for the subsequent higher frequency reconstructions. At the lowest frequency, assumptions can be made that the phase changes are modest and it can be assumed that all of the data is unwrapped. While the lower frequency images can be quite blurred because of the associated larger wavelengths, the algorithms converge more reliably. The property images that are transposed to the next higher frequency can be close enough to the actual images at the increased frequency such that it essentially positions the phases for all transmit/receive pairs on the appropriate Riemann sheet. The notion is that the resolution of each successively higher frequency increases until one reaches the highest frequency. Finally, a number of groups have advocated utilizing substantially more measurement data than that prescribed by the Dartmouth team [25, 27]. While these often advance the notion that there needs to be as much measurement data as number of unknowns, our own experience based on a wealth of literature suggests that this may not be necessary [52]. An equally valid interpretation in the context of this discussion is that they hope to increase the amount of “good” data such that the amount of “good” data simply outweighs the amount of “bad” data for which a least squares image process can achieve a reasonable image. For many real situations, each of these solutions may be unrealistic. Accurate a priori information may be difficult to generate at the time of or before the actual imaging session. Frequency hopping may be unrealistic because it is difficult to devise antennas with sufficient bandwidth to accommodate the algorithm. Finally, adding more data inevitably implies that the measurement system will require more channels which inherently leads to increases in the algorithm complexity and hardware costs.

The Dartmouth team has devised robust unwrapping techniques for both the measured and computed phases which are briefly summarized in Section 2.2 [50]. These generally exploit the wide operating bandwidth of the monopole antennas and even the nature of the algorithm convergence. One important consequence of these developments is that the amount of measurement data can be kept to a minimum (typically 16 antennas configured in a circle for 2D images) which dramatically reduces overall system cost and complexity [52]. The algorithm is fast and does not suffer from convergence to unwanted local minima even when starting from an initial estimate of the coupling bath properties. These approaches have been developed in the context of utilizing actual measurement data to maximize the benefit from each piece of measurement data while not imposing Riemann sheet criteria that essentially transform “good” data into “bad”. This approach is summarized in Section 2.2 and demonstrates excellent convergence behavior using real measurement data.

2. Methods

2.1 Hardware

The heart of the system concept is the monopole antenna array submerged in a lossy glycerin: water coupling bath. As discussed in the Introduction, the monopole antennas offer important advantages such as small size, ease of modeling and broad bandwidth while also presenting different opportunities. The major trade-off is that the signal attenuation is great, even for a typical physical distance of 15 cm across the lossy medium in the illumination chamber, such that most commercial vector network analyzers are inadequate. This has changed somewhat in that there are now commercial 24 channel systems from Rohde & Schwartz (Munich, Germany) and Keysight Technologies (Santa Clara, CA); however, their costs can easily exceed $100 K depending on the options necessary. For realizing an actual system, the most critical trade-offs are cost versus performance. In our case, this discussion involved numerical modeling concerns along with critical aspects of the microwave technology. For instance, as mentioned above, there is considerable debate regarding the required number of antennas. Our approach has been to configure just enough antennas for a single-plane array and then exploit mechanical motion to move the array vertically to achieve complete 3D coverage of the target—in this case the breast. We have demonstrated that 16 antennas is sufficient for the single array [52], which has kept our costs to modest levels. Figure 3 shows the antennas mounted to a horizontal plate which is supported by opposing motors that can bring the antennas to different levels surrounding the breast to enable full volumetric coverage of the breast. Below we summarize the more significant issues for the microwave electronics.

Figure 3.

Photographs of the antenna array, support structure and motors used for vertical, mechanical array motion: (a) close-up of mounting plate, (b) side view of tank, antennas and support structure, and (c) top view of the antennas in the tank.

2.1.1 Large dynamic range

In general, the dynamic range is limited on the low end by the noise floor which is defined as $N=\mathit{kTB}$ where k is Boltzman’s constant, T is the temperature and B is the sampling bandwidth. The primary option for decreasing the noise floor is increasing the sampling time (i.e. additional averaging) which is the inverse of the bandwidth. Theoretically the noise floor could be made infinitely low if one allowed for the data acquisition to take forever. However, for practical systems, the lower end is also limited by the discretization resolution of the A/D converter. The most common way to address both is to add increased amplification in front of the receiver while also increasing the sampling time. The typical VNA is configured so that the lowest possible noise floor is roughly −100 dBm. Given that the transmit signals typically can range on the order of 1–10 mW, the maximum dynamic range is roughly 102 dB (assuming that we can realistically pick out a signal that is roughly 8 dB greater than the noise floor). In building our custom system, we have employed the Ettus B210 agile transceiver which by themselves can only see down to −110 dBm. By adding 20 dB gain, increasing the sampling time, and adding a low noise amplifier, we have simultaneously improved the overall noise figure of the system while also lowering the noise floor and boosting the signal. In this way, we are able to achieve substantially higher dynamic range capabilities in a low cost configuration. The Ettus boards cost a fraction of that of typical VNA’s and are programmable using general purpose instrumentation software including Matlab (Mathworks, Natick, MA) and LabView (National Instruments, Austin, TX) [53].

2.1.2 Channel-to-channel isolation

One hypothetical way to achieve a large dynamic range would be to use a 2-channel VNA coupled to a double-pole, 16 throw (DP16T) switch matrix to achieve coverage for all 16 channels. We have already discussed the challenges in achieving the dynamic range for a 2-channel system. Unfortunately, the multi-pole switch matrix introduces its own set of problems. In cases where the 2-port VNA already has a dynamic range approaching 140 dB, it is necessary for the switch matrix to have channel-to-channel isolation on the order of 140 dB. In practice it is especially difficult to achieve this and very rarely do switch manufacturers specify isolation greater than roughly 80–100 dB. There are ways to achieve this, but it generally requires substantial shielding and the addition of extra single-pole, single-throw (SPST) switches which introduces substantial insertion loss and inevitably detracts from the overall dynamic range.

Our latest concept is to exploit new software defined radio (SDR) technology as the key building block of our system [53]. In this case, we utilize the Ettus Research B210 (Austin, TX) which has two channels, each having a transmit/receive port and a dedicated receive port. Our design involves a physical separation of the boards used for transmit and receive (extra isolation), along with a set of switch/amplifier modules which add to the overall channel-to-channel isolation along with improving the dynamic range with a 20 dB gain amplifier in front of the receiver channel. One of the signals from the Transmit B210 is used as a reference signal to synchronize the transmit oscillators with those of the receive modules. Figure 4 shows a schematic of a representative 4-channel system along with an initial prototype showing the separate transmit and receive modules along with the switch and amplifiers. This design can be readily scaled up to a full 16 channel system.

Figure 4.

(a) Schematic diagram of a 4-channel measurement system and (b) a photograph of a first prototype (a is reproduced from “Meaney PM, Hartov A, Bulumulla S, Raynolds T, Davis C, Schoenberger F, Richter S, Paulsen KD, “4-channel, vector network analyzer microwave imaging prototype based on software defined radio technology,” Review of Scientific Instruments, vol. 90, pp. 044708-1—044708-14, 2019,” with permission of AIP Publishing [53]).

With respect to isolation, Figure 5a shows the shielded housings for the B210 boards. The primary feature is the ridge in the cover which lightly rests on a gold strap running across the board that physically separates the RF portion of the board from the lower frequency electronics and power supply circuitry. Tests have shown that the “one-way” leakage from these housings is on the order of −110 dB (from a nominally transmitted signal of 1 mW). While a signal of this amplitude may propagate along the outer surface of the electronics, because all of the other components of the system have comparable isolation specification, the greatest signal level that could leak back into the system would be on the order of −220 dB because of the “two-trip” attenuation in escaping the isolated board and then re-entering a companion board. Figure 5b shows the internal features of an earlier, custom transceiver module which included several amplifiers and switches along with a downconverting mixer. All feedthrus for the bias and control lines were coaxial low pass filters to restrict RF signals from escaping outside of the housings. Each component within the housing was sequestered to its own channelized compartment to limit internal leakage and the covers were configured with raised surfaces for tight fits of the two housing portions to further limit leakage to the outside. While somewhat cumbersome, in large quantity the housings can be fabricated at relatively low cost and have proven to be effective in limiting leakage, even to the extreme levels of −140 dBm.

Figure 5.

Photograph of the (a) shielded housings for the Ettus B210 SDR circuit boards, and (b) the shielded housing for an earlier transceiver module (b is reproduced from “Li D, Meaney PM, Raynolds T, Pendergrass SA, Fanning MW, Paulsen KD, “A parallel-detection microwave spectroscopy system for breast imaging,” Review of Scientific Instruments, vol. 75, pp. 2305–2313, 2004,” with permission of AIP Publishing [54]).

2.1.3 Synchronization

Commercial VNA’s are attractive because they produce measurements that guarantee built-in synchronization of the transmit and receive signals. However, as was discussed earlier and also in the subsequent sections, these network analyzers are unattractive because of their costs. There are a host of means for achieving synchronization of signals in different hardware configurations. For our SDR-based system, we exploit the phase locking of the transmit and receive oscillators with respect to a 10 MHz reference signal from their Octoclock (Ettus Research, Austin, TX). In addition, accurate and consistent triggering of the measurements is achieved via the PPS signal which is also supplied by the Octoclock. In the case with separate transmit and receive modules, we exploit the fact that the second output signal from the dedicated transmit board is synchronized to the output of the channel that feeds the antennas. This is subsequently used as a reference signal and is fed into an unused receive port of the receive channels for sampling and synchronization with the receive boards [53, 54].

2.2 Software

While much has been made of a range of non-linear iterative schemes and various stochastic approaches for image reconstruction approaches, we have opted to use the robust and well-regarded Gauss-Newton algorithm. Figure 6 shows examples of the magnitude contour plots for 2D wave propagation from different antennas surrounding the medium with a high contrast object. The differences between the values measured at the receive antenna locations for the situation where an object is present minus that for the case for the homogeneous bath becomes the measurement data for the algorithm. Note that the techniques are completely translatable to 3D but are more easily understood in the context of 2D.

Figure 6.

2D magnitude contour plots within an imaging zone containing a high contrast scatterer for point source illuminations from different directions.

For the untransformed algorithm, the minimization statement examines the differences between the measured field values and the associated computed values at each receive antenna and for all views of the object from each iteration:

$$min{‖{\mathrm{E}}^{\mathrm{m}}-{\mathrm{E}}^{\mathrm{c}}\left({\mathrm{k}}^2\right)‖}^2$$E1

where ${\mathrm{E}}^{\mathrm{m}}$ and ${\mathrm{E}}^{\mathrm{c}}$ are the complex measured and computed electric field vectors, respectively, and ${\mathrm{k}}^2$ is the wavenumber squared which is essentially the image and embeds both the permittivity and conductivity, respectively. The associated log-transformed equation is [49]:

$$min{‖{\mathrm{\Gamma }}^{\mathrm{m}}-{\mathrm{\Gamma }}^{\mathrm{c}}\left({\mathrm{k}}^2\right)‖}^2+{‖{\mathrm{\Phi }}^{\mathrm{m}}-{\mathrm{\Phi }}^{\mathrm{c}}\left({\mathrm{k}}^2\right)‖}^2$$E2

where ${\mathrm{\Gamma }}^{\mathrm{m}}$, ${\mathrm{\Gamma }}^{\mathrm{c}}$, ${\mathrm{\Phi }}^{\mathrm{m}}$, and ${\mathrm{\Phi }}^{\mathrm{c}}$ are the measured and computed log-magnitude and phase vectors, respectively. The two equations are essentially the same except for the fact that the phase and magnitude need to be accounted for explicitly. In addition, both the measured and computed phases need to be properly unwrapped. The image is then updated at each iteration by:

$${\mathrm{k}}_{\mathrm{i}}^2={\mathrm{k}}_{\mathrm{i}-1}^2+\mathrm{S}∆{\mathrm{k}}^2$$E3

where ${\mathrm{k}}_{\mathrm{i}}^2$ and ${\mathrm{k}}_{\mathrm{i}-1}^2$ are the wavenumbers squared vectors at iterations i and i + 1, $∆{\mathrm{k}}^2$ is the image update at each iteration and S is the iteration step size. In all cases, the initial property estimate of the image is that of the homogeneous coupling bath which is the least biased estimate possible. The stopping criteria is generally met when the changes in the image ($∆{\mathrm{k}}^2$) are small—determined empirically. For all cases, regularization is necessary to stabilize the solution. For this situation, we have found both Tikhonov and Marquardt–Levenberg techniques to be suitable [55, 56].

For applications such as MR, phase unwrapping is normally performed as a function of position—i.e. the phase at a specific location is set as the reference phase and the phases at neighboring locations are compared with it [57]. If the differences between the values at new locations differ from the reference ones by more than 180°, the values for the new ones are adjusted by multiples of 360° until the differences are minimized. This process is continued until the entire spatial domain has been unwrapped. For the microwave measured values, this is impossible since the phases can only be known at the receive antenna locations—i.e. the phase can vary more than 180° between antennas. Our alternative is to unwrap the phases as a function of frequency. A baseline assumption implies that the scattered phases should be within the bounds of +/− 180° for the lowest frequency. The measured phases are forced to be within the baseline Riemann sheet, i.e. −180 to +180°. The values at the next frequency are unwrapped against these, followed by those at the next frequency being unwrapped with respect to the second set and so on. This is described in more detail in Meaney et al. [50]. For the computed phases, it would be possible to unwrap the values as a function of position. However, this is tedious, complicated and slow. We previously introduced a novel technique to unwrap the phases as a function of image reconstruction iteration [50]. As long as the phases do not change dramatically between iterations of the reconstruction process, it can be assumed that their differences will be less than 180°. To ensure this criterion, we deliberately force the step size (Eq. (3)) to be small during the earliest iterations so that the image, and hence the phases, differ only slightly between iterations. In this way, the phases only need to be computed at the antenna locations, and the values for the previous iteration need to be stored for comparison purposes. This technique is described in more detail in Meaney et al. [50].

Finally, the most substantial computational cost for these types of algorithms are for the forward field solutions which need to be computed for each transmit antenna at each image reconstruction iteration. Different techniques have been proposed to reduce this time, most notably finite element (FEM) and finite difference time domain (FDTD) modeling [26, 58]. However, these can be quite slow, especially as the forward problem grows in size, and often requires multi-processor computers and GPU processors to keep the computation times modest. While our 2D FEM-based approach can usually produce images in roughly 5 minutes, we have explored ways to further reduce the time. In particular, we have been experimenting with the discrete dipole approximation (DDA) to improve efficiency [59]. This technique is generally not suitable when there are large, high contrast scatterers in the field of view—especially other bulky antennas; however, because our monopole antennas only slightly perturb the fields, and almost not at all when submerged in a lossy bath, this algorithm is ideal for our approach. The primary notion is that it utilizes Gauss-iterative based forward solver techniques whose most substantial time cost is a matrix–vector multiplication that needs to be performed repeatedly (O(N²)) where N is the vector length). However, it can be easily shown that the associated matrix can be formulated to be a symmetric, block Toeplitz matrix. Each block matrix can then be easily converted to a circulant matrix. Circulant matrix–vector multiplication can be performed efficiently (O(NlogN)) by use of the convolution theorem and the use of the fast Fourier transform (FFT). We have shown that this can reduce the forward solution time by a factor of 20–50 times with respect to efficient commercial software—i.e. full reconstructions in 6 seconds or less. This technique is described in detail in Hosseinzadegan et al. [59, 60].

3. Results

3.1 Algorithm convergence

While the final desired result of the iterative image reconstruction algorithm is the associated property images, the best way to assess the quality of the process is by analyzing the computed electric fields as compared with the measurement data [47]. As can be seen from Eq. (2), the minimization statement should converge to a point where measured and computed magnitudes and phases achieve a reasonable match. In essence, the final image can really be considered a bi-product of this process. A key reason why this distinction is important is that there are ample situations where researchers employ a range of a priori information to achieve reasonable results. In a laboratory or experimental setting this might appear reasonable; however, these ad hoc measures often introduce biases which can artificially skew the final results. The best unbiased way to examine the process is to focus on the actual terms being compared.

Figure 7 shows the normalized (a) magnitude and (b) phase values for the receiver positions associated with a single transmitter for both the measured values (heavy black lines) and computed values (narrower colored lines) at selected iterations [50]. The data was collected at 1300 GHz and both the magnitude and phase values have been normalized to the cases where there is no target in the bath—i.e. that the reconstruction process starts with an assumption of an empty bath which is reasonable since it is possible to exactly know the bath properties in all situations. Note that the measured phases exceed the +/− 180 degree bounds in some cases, implying that our algorithm is robust even for large, high contrast targets which often generate some measurements outside of the primary Riemann sheet. As a useful analogy, we refer to the curves for each measurement set as a projection, much like the nomenclature used in X-ray CT. It should be noted that in this example, a two-step algorithm was used whereby the first 50 iterations used a standard Levenberg–Marquardt technique to produce a good, albeit smoothed image. After which, a Tikhonov algorithm was used for 20 iterations to more carefully refine the images. The intermediate projection is a black, dashed line for both the magnitude and phase.

Figure 7.

1300 MHz (a) Magnitude and (b) phase projections for the measured data and multiple computed values at representative iterations. Note that the dashed black lines represent the projections at the end of the first step (Permission granted by Wiley Publishers to re-print graphs [50]).

For both the magnitude and phase, the projections monotonically approach those for the measurement data. The projections for the intermediate image have not yet matched the measurement data, but their curves retain many of the characteristics of the measurements but with a reduced overall strength and are also somewhat smoothed. Once the second step is employed, the computed magnitudes and phases rapidly approach the measured values. The convergence of the phase values is almost exact, while those for the amplitude are good but not as precise. The discussion in Section 4 provides a possible rationale for this mismatch. Regardless, the match is quite good and gives confidence that the process works well. It should be noted that this is for an actual breast exam where the measurement data inherently has 3D characteristics while the algorithm is a 2D one. In addition, these plots are for only one of the 16 antennas but is reasonably representative of all sets.

3.2 Therapy monitoring

While microwave imaging is poised to play a substantial role in the detection and diagnosis of breast cancer, it may be that its most important role will be in the area of neoadjuvant chemotherapy monitoring [17]. In this setting, it is important to be able to tell if the tumor is responding adequately to treatment at a sufficiently early time point so that treatment can be altered appropriately. Measures such as a physical exam present morphological information regarding the tumor response—i.e. size—however, these changes are generally quite slow to manifest compared to physiological changes. Early results suggest that the microwave images may be more physiological in nature [17].

Figure 8a shows representative sagittal MR scans of the ipsilateral breast for a woman with heterogeneously dense breasts and a 6.5 x 3.7 x 7.1 cm³ tumor in the right breast just before treatment started: (i) T2, (ii) T1 with gadolinium, and (iii) the T1 subtraction image for the cases with and without gadolinium. (i) Shows a large thickening of the skin surrounding most of the breast, while (ii) and (iii) show the extent of the tumor. Figure 8b shows the corresponding images after 85 days of treatment. The skin thickening had diminished substantially and the subtraction image indicates only small, minor enhancements.

Figure 8.

Sagittal MR images of the right breast of a patient with a complete pathologic response: (a) images prior to therapy and (b) at day 85. For each set, (i) is a T2-weighted image, (ii) is a contrast-enhanced image from a spoiled gradient recalled (SPGR) sequence, and (iii) is a subtraction image between (ii) and a pre-contrast baseline, respectively (Reproduced from [17], with the permission of Springer Nature—License link: https://creativecommons.org/licenses/by/4.0/).

Figure 9 shows the 1300 MHz permittivity and conductivity images for planes 5, 6, and 7, respectively, where plane 1 is closest to the chestwall and plane 7 is closest to the nipple. Set (a) is for the contralateral breast (normal) and sets (b), (c), and (d) are for the ipsilateral breast at different time points during the treatment—prior to treatment, after 44 days of treatment, and just prior to surgery, respectively. In all cases, the outline of the breast is readily visible separating the tissue from the surrounding bath. For the contralateral breast, the majority of the tissue exhibits quite low permittivity and conductivity which is consistent with adipose tissue, while there are also centralized zones of slightly elevated tissue corresponding to the fibroglandular tissue. For the diseased breast at the first two time points, there is a pronounced ring of highly elevated tissue properties—primarily to the lower left quadrant. This feature is consistent with the edema that would be present from the skin thickening observed in the MR images. For these same time points, the localized, elevated properties of the tumor are also visible, albeit slightly diminished at the 44 day mark. By the final set, both the microwave observable skin thickening and tumor have essentially disappeared in agreement with the MR images shown in Figure 8b. It is noteworthy that the tumor has shrunk considerably with respect to size and property intensity by the midway point suggesting that the microwave images are capable of recognizing early tumor response.

Figure 9.

1300 MHz microwave tomographic images. Imaging planes five to seven are shown corresponding to the three closest planes to the nipple with the permittivity on the top row and conductivity on the bottom row. (a) Left (contralateral) breast prior to treatment, (b) right (ipsilateral) breast prior to treatment, (c) right breast 44 days into treatment, and (d) right breast immediately prior to surgery (Reproduced from [17], with the permission of Springer Nature—License link: https://creativecommons.org/licenses/by/4.0/).

4. Discussion

One of the more interesting aspects of our imaging system revolves around the notion of why a 2D algorithm works in the first place. While there have been efforts by others exploring 2D imaging, most have reverted to 3D implementations when 2D efforts failed. Notably, Semenov et al. [61] concluded that the reason their 2D system did not work was because of the inherent mismatch between the actual 3D wave propagation and the assumed 2D propagation in their algorithm. However, we contend that there are scientifically valid reasons why our 2D algorithm works well. The most important is that by choosing the monopole antennas, we are able to position them on a quite small diameter circle surrounding the target. Given the vertical orientation of the antennas, in the limit of a decreasing array circle diameter, the aspect ratio of the antenna length to the diameter of the circle naturally increases. This can be theoretically extrapolated to the point where the domain closely mimics that of cylindrical geometry. For competing systems, the motivation has often centered about acquiring as much measurement data as possible. This often implies that the antenna array must be configured on a larger diameter circle for the simple reason relating to space limitations associated with the physical size and the increased number of radiating elements. In Semenov et al. [62], thousands of antennas (effective number based partially on mechanical motion) were positioned on a 60 cm diameter circle for a roughly 8 hour, 3D data acquisition, single frequency exam time, while the diameter for the Dartmouth system is only on a 15.2 cm diameter where the broadband data acquisition takes roughly 10 minutes. For the latter, the field patterns have simply not devolved into their spherical radiation behavior within the small zone while they are considerably further along that process for the larger diameter system.

In Meaney et al. [63], we compared the magnitude and phase values for the scattering from a simple phantom illuminated in our imaging chamber. It is noteworthy that there was very little difference between the phase calculations using both 2D and 3D models. There was some deviation for the magnitude, but the differences never exceeded 0.5 dB. It is worth examining the far field behavior of the 2D and 3D waves, in spite of the fact that it is a simple estimation in this situation. For the phase in the main beam, the phases propagate as a function of R (distance from the antenna) for both the 2D and 3D cases. This presumably accounts for the close match of the previous calculations. However, for the magnitude, they decay as a function of 1/R and 1/R² for the 2D and 3D cases. Likewise, this likely partially accounts for the previously mentioned differences. As hypothesized above, these differences are likely mitigated because we operate within such a physically constrained space. As a final note, the lossy coupling bath most likely assists for this comparison by substantially attenuating (and almost eliminating) signals that might propagate out of the plane of interest only to be reflected back into it during the propagation process. Regardless, our 2D configuration in conjunction with our 2D algorithm has provided good images in a variety of settings from simulations, phantom experiments, animal experiments and even clinical trials.

5. Conclusions

The Dartmouth microwave imaging system is quite robust because of its synergism of hardware and numerical modeling expertise. Competing implementations of microwave imaging approaches generally ignore multipath signals. In fact, most numerical models preclude their existence by not incorporating structures which would excite surface wave propagation. The Dartmouth approach directly addresses this phenomenon which has led to a range of hardware and software innovations. Most notably, by being able to pack the antennas in such a small circle about the target, various factors have aligned to make 2D imaging a reality.

In addition, by incorporating fundamental concepts from parameter estimation theory, we have been able to implement a version of the image reconstruction process that essentially eliminates the a priori information requirement. While the log transform is most often used for purely real mathematical functions, there is precedent for its use with complex functions—notably for optical coherence tomography. The concept takes on an added level of complexity for microwave imaging because of the phase unwrapping criteria. By developing simple and effective ways to unwrap the phases, we now have an algorithm which does not require a priori information and does not converge to local minima or unwanted solutions. This type of robustness is essential when working in actual clinical situations.

Finally, the results bear out that microwave images can be diagnostically useful in the breast neoadjuvant chemotherapy monitoring setting. Initial pilot results suggest that the property changes observed as early as 30 days after treatment correlate well with final treatment response. Essentially the microwave properties provide a physiological assessment of the response which is known to occur earlier than measures based on morphological phenomena. While physiological models have been proposed that may explain the behavior observed in the imaging trials, more data is required to confirm this phenomenon. Regardless, these results suggest that microwave imaging may be a viable approach and may provide additional tools for clinicians in their fight to treat breast cancer.

Acknowledgments

This work was supported by NIH/NCI grant # R01-CA240760.

Conflict of interest

Drs. Meaney and Paulsen are co-inventors on several patents related to microwave imaging.