Importance of phase unwrapping for the reconstruction of microwave tomographic images.

Microwave image reconstruction is typically based on a regularized least-square minimization of either the complex-valued field difference between recorded and modeled data or the logarithmic transformation of these field differences. Prior work has shown anecdotally that the latter outperforms the former in limited surveys of simulated and experimental phantom results. In this paper, we provide a theoretical explanation of these empirical findings by developing closed form solutions for the field and the inverted electromagnetic property parameters in one dimension which reveal the dependency of the estimated permittivity and conductivity on the absolute (unwrapped) phase of the measured signal at the receivers relative to the source transmission. The analysis predicts the poor performance of complex-valued field minimization as target size and/or frequency and electromagnetic contrast increase. Such poor performance is avoided by logarithmic transformation and preservation of absolute measured signal phase. Two-dimensional experiments based on both synthetic and clinical data are used to confirm these findings. Robustness of the logarithmic transformation to variation in the initial guess of the reconstructed target properties is also shown. The results are generalizable to three dimensions and indicate that the minimization form with logarithmic transformation offers image reconstruction performance characteristics that are much more desirable for medial microwave imaging applications relative to minimizing differences in complex-valued field quantities.

1. Introduction

The use of electromagnetic waves and their strong interaction with the media through which they propagate has yielded important promises for medical imaging applications over the past several decades [1-8]. The images typically represent maps of the two main electrical constitutive parameters: the permittivity ε(ω) and the conductivity σ(ω). These parameters are usually frequency dependent, indicated by their function of ω, where ω = 2πf and f is the frequency in Hertz. In the microwave frequency regime, around 1 GHz in our case, biological tissues are characterized by parameters that are highly frequency dependent, owing to their strong dependency on water content [9, 10]. Consequently, the presence of tumors, rich in water, produces significant permittivity and conductivity contrasts with respect to healthy tissues which can be captured by various methods.

The Born iterative method [11] and distorted Born iterative method [12] have been shown to successfully recover two-dimensional permittivity profiles, albeit usually applied to smooth functionals or low contrast objects [13] and still undesirably dependent on initial guess [14]. Furthermore, their converge rate can be degraded in the presence of noisy data. These pitfalls have led to the investigation of more numerically robust approaches where the minimization problem is solved in a least-square sense typically using Newton-like methods [15] with some level of regularization [16]. Alternatively, radar-based approaches have been reported to yield successful results in simulation [17], and has been adapted to a clinical setup in [18]. These approaches only produce maps of maximum signal strength and their effectiveness is limited in high contrast, heterogeneous targets such as is often the case for radiographically dense breasts [19]. In parallel to deterministic methods, stochastic approaches deserve to be mentioned as they often offer the possibility of spanning the search space efficiently without the drawback of remaining trapped in local minima, which are two important properties for microwave imaging. A review of some methods is proposed in [20] with a direct application to biomedical imaging in [21]. The major drawback at this point, however, seems to still be the computational time which, as reported, needs to be reduced in order for these methods to be fully applicable in a clinical setup.

In the present paper, we implement a deterministic Newton-like approach which produces two-dimensional images with sub-minute efficiency. The associated update equation is written in terms of electromagnetic properties (permittivity and conductivity) and is applied iteratively to improve initial estimates. The nonlinearities involved along with the finite and discrete sampling of the body’s response provided by the measurement data containing instrumentation noise, lead to solution instability and multiplicity that require regularization. Various formulations of the minimization problem have been proposed. One approach minimizes the complex-valued differences between the simulated fields at each receiver and their measured counterparts using a Gauss-Newton iterative approach [22]. The least square problem is formulated in terms of the normal equations (written here without regularization) as [23]

where k and k are the real and imaginary parts of the wavenumber [k(ω)]2 = ω2ε(ω)μ0 + iωμ0σ(ω), respectively, Δ denotes the increment at each iteration, J and J are the Jacobian matrix and its transpose, respectively, and ℜ and ℑ represent the real and imaginary part operators, respectively. The electric field quantities of the right-hand side are obtained from measurement data (Ē) and computed with our forward Finite-Difference Time-Domain solver (Ē). We refer to the use of the complex-valued fields, directly in terms of their real and imaginary parts, as the complex algorithm, which is both intuitive and easy to formulate. This algorithm does not retain information on any phase wrapping of the signals since exp(iϕ + i2mπ) = exp(iϕ) when m is an integer. We show in this paper that the absolute phase contains crucial information for the proper estimation of the electric permittivity and conductivity profiles. This situation is analogous to height inversion in interferometric SAR, where the height value crucially depends on the ability to unwrap the phase at the receiver level [24, 25]. When the phase wrapping information is lost, proper reconstruction of the permittivity and conductivity profiles is possible only in simple cases: when the contrast in the field of view is low, when the reconstruction occurs at low frequencies (i.e. with long wavelength signals), or when a good initial guess is provided to the algorithm. These three conditions effectively ensure that the phases of the measured and computed electromagnetic signals in the imaging region remain within 2π of each other which corresponds to a constant m throughout the imaging region, thus avoiding jumps in phase to different Riemann sheets.

A second method used to formulate the minimization criterion is based on the functional relationship between the estimation parameters and the scattered field. In one dimension, this relationship is known to be exponential (a derivation is presented in Section 2) making it natural to linearize the problem using logarithmic transformation [26]. In this case, differences between the logarithmic amplitudes and the phases of the computed and measured fields are minimized [22], an approach we term the logarithmic algorithm, for which the normal equations become (again written without regularization)

where Γ denotes the logarithmic amplitude, Φ represents the phase, and J is the new Jacobian matrix. An immediate disadvantage of this method appears to be the complexity of the Jacobian matrix J. However, it has been shown that J can be constructed from the Jacobian of the complex algorithm directly [27], thus alleviating any additional computational burden.

In this paper, we provide a theoretical explanation of the dependency of the estimated permittivity and conductivity on the absolute (unwrapped) phase of the measured signals at the receivers relative to the source transmission by developing closed form solutions for the field and inverted electromagnetic property parameters in one dimension. The fundamental strength of the logarithmic algorithm is that it retains this absolute phase information of the electromagnetic fields, which requires an unwrapping procedure across multiple Riemann sheets [27, 28] without which phase information is lost and the reconstructed property profiles become incorrect. A two-dimensional configuration is also considered, identical to the one described in [29]. Comparison between the algorithms and results quantification is performed first on a series of phantom reconstructions for which the truth (position, shape, properties) is known. The analysis shows that the logarithmic algorithm yields improved images compared to those obtained with the complex algorithm from various perspectives. First, it quickly and consistently converges to the same and correct solution regardless of initial guess, whereas the images obtained from the complex algorithm may not only be inaccurate but also converge very slowly. Second, for a lossy background configuration, the complex algorithm suffers from disappearing targets at higher frequencies and inaccurate reconstruction of larger targets, both of which do not occur with the logarithmic algorithm. Finally, both algorithms are tested in a clinical situation where the truth is obtained from magnetic resonance images of the patient against which our microwave images are compared. Here again results show superiority of the logarithmic algorithm and illustrate how our one dimensional mathematical analysis has direct applicability to clinical two dimensional hardware systems.

2. Mathematical considerations in one dimension

This section presents a one dimensional simplification of the more typical two dimensional tomographic conditions considered in the next section. As such, we consider a simple dielectric and lossy slab upon which a plane wave impinges, creating measurable reflected and transmitted waves. The well-known field solutions corresponding to this configuration can be inverted analytically, yielding a closed-form solution for the unknown permittivity and conductivity. The resulting mathematical formulae show that knowledge of the absolute phase is essential for correct reconstruction of both parameters. A time dependence of e−i is assumed and subsequently omitted.

2.1. Formulation

The one dimensional configuration consists of a slab of thickness d and relative permittivity ε = ε′ + iε″, embedded in a background of permittivity ε and conductivity σ = 0, as depicted in Fig. 1. Both media are non-magnetic and the background is lossless, which does not reduce the generality of our conclusions. A plane wave of wavenumber (where k0 = ω/c is the wavenumber in free-space) impinges on the slab at normal incidence. The reflection and transmission coefficients are measured as S-parameters, defined by [30]

Fig. 1.

Configuration of the problem: a wave impinges on a slab of infinite lateral dimensions, and creates a reflected and transmitted wave that can be measured.

where is the index of refraction of the slab and

is the reflection coefficient at the first interface [31]. The problem consists of reconstructing the complex permittivity ε knowing only S11 and S21 from measured data.

2.2. Inversion formulae for a single point

Henceforth, we suppose that S11 and S21 are known for a single realization, which typically represents measurements at one given frequency. We are, therefore, interested in reconstructing a scalar ε, as opposed to estimating a function of frequency (the generalization of which is straightforward [32]).

The mathematical formulation for estimating the constitutive parameters from reflection and transmission measurements was first described in [30], generalized recently to reconstruct the effective parameters of metamaterials [33, 34]. Based upon these works, we proceed by noting that the reconstruction of ε can be achieved by first inverting for n in Eqs. (3). After a few algebraic manipulations, it can be shown that

The choice of the sign in Eq. (5) is dictated by the physical requirement that the real part of the impedance and the imaginary part of the index of refraction n″ of the medium (reconstructed subsequently in Eq. (6)) both be positive (the condition n″ > 0 may yield numerical instabilities for lossless media and should therefore be used with caution. In the present case, however, we avoid this issue by only considering lossy media, motivated by the medical imaging application.) Once the sign is established, Eq. (5) essentially indicates that is a known quantity which can be evaluated from measurement data. We can therefore invert for n as

Note the presence of the term 2mπ where m is a positive or negative integer due to the fact that ei2m = 1. The inversion formula in Eq. (6) indicates that the reconstruction of the imaginary part of the index of refraction can be accomplished directly, whereas the inversion of its real part is complicated by the presence of the branch cut in the complex logarithm function, which appears as the unknown parameter m. It should be noted that the choice of m is irrelevant as far as the exponential function is concerned, i.e. it is not visible in the results from the forward solution.

The inversion of the complex permittivity is obtained directly from

and the conductivity is determined from σ = ωε0ε″. As a consequence, the parameter m influences both real and imaginary parts of the permittivity or, in other words, the value of the absolute phase influences the inversion of both the permittivity and the conductivity.

This result indicates that restricting the reconstruction of ε′ and σ to m = 0 yields erroneous results. Let us consider a simple situation relevant to medical imaging, with the following parameters: ε = 28.3, σ = 0, ε′ = 63.9, σ = 1.75 S/m. These property values correspond to a background composed of 80% glycerin and a phantom consisting of 40% glycerin at 1100 MHz. Note that σ = 0 in order to use Eqs. (3). In case of a lossy background, this assumption is tantamount to measuring the fields at the boundary of the slab, which can subsequently be related to the fields at any point if the background permittivity and conductivity are known. The assumption is therefore not directly relevant to the behavior of the inversion algorithm.

The reconstructed results are shown in Fig. 2 for both ε′ and σ and indicate that:

Fig. 2.

Reconstruction of the permittivity (ε′) and conductivity (σ) of the slab for various values of m and parameters k0d. The true values are ε′ = 63.9 and σ = 1.75 S/m.

As expected, the parameter m influences both the real and imaginary parts of ε, i.e. the permittivity and the conductivity.

The true value is obtained for different choices of m, as a function of the parameter k0d. Note that k0d represents both the frequency and the size of the slab. If it is understood that an increase in k0d is due to a variation in frequency, the frequency dependence of ε should be introduced which would quantitatively modify the curves in Fig. 2, but not their qualitative implication that a proper reconstruction of the constitutive parameters necessitates knowledge of the unwrapped phase. Moreover, the parameter m can be related to the number of periods of the electromagnetic wave inside the slab. Hence, for low values of k0d, is it not surprising to find that m = 0 gives the right solution. As k0d increases, however, more periods of the wave are included in the slab (either because the slab is thicker or because the wavelength is shorter) and the correct solution is obtained for higher values of m.

Consequently, the determination of m becomes the main purpose of an unwrapping procedure. For example, we can start with a low k0d value to ensure that m = 0, then progressively increase k0d while tracking the phase in order to capture the jumps in m. Interestingly, low values of k0d correspond to either thin samples, low contrasts, or low frequencies. It is therefore expected that the reconstruction performs best in these cases if determination of m is not possible, which is consistent with the two-dimensional imaging results commonly obtained by microwave tomographic systems.

In order to explore results that are more typical of the two dimensional medical imaging environment, we next consider a frequency dispersive case where the permittivity of the background corresponds to a liquid bath with 80% glycerin (we still assume the background to be lossless), while the parameters (permittivity and conductivity) of the slab correspond to a liquid of 50% glycerin. In addition, we fix the slab thickness at 1 cm, corresponding to a small tumor, and study the reconstruction of the slab’s parameters as function of frequency. The results, summarized in Fig. 3, indicate that the solution m = 0 is correct at lower frequencies, while the solution m = 1 is correct at higher frequencies. The cross-over between the two solutions occurs at about 1.8 GHz.

Fig. 3.

True and reconstructed frequency dependent parameters using Eq. (6) for m=0 and m = 1. The background permittivity (dashed curves) corresponds to an 80% glycerin bath while the slab permittivity and conductivity (solid curves) correspond to a 50% glycerin bath.

3. Two dimensional tomographic reconstruction

We subsequently consider the two-dimensional tomographic geometry described in [29]: the tomographic microwave setup is composed of 16 antennas equi-spaced in a circular fashion (with a diameter of 15 cm) where each antenna acts sequentially as a transmitter whose signal is captured by the 15 remaining receivers, yielding 15 ×16 = 240 measurement data at each excitation frequency (yielding effectively 120 independent observations due to the reciprocal character of the media). The antennas are submerged in a lossy liquid in order to minimize the multipath contributions at each receiver from reflections at the tank boundaries. A direct generalization of the study performed on the one-dimensional slab suggests that the complex algorithm should accurately reconstruct property profiles for objects that are either small, low contrast, or imaged at low frequencies, but poorly reconstruct objects that are large, higher contrast, and/or imaged at high frequencies. On the other hand, the logarithmic algorithm should accurately reconstruct all cases provided that the phases are properly unwrapped. We first investigate the effects of size and explore the influence of frequency in Section 3.2.

3.1. Reconstruction of large objects

The configuration under study consists of three circular phantoms of increasing sizes which are embedded in a background whose parameters are given in Table 1. The measured fields at each receiver are obtained numerically, using the forward model described in [35]. The data are then used in our iterative inversion algorithm, based on either the complex or the logarithmic forms. The results are shown in Figs. 4–5 and illustrate the higher quality of the images recovered using the logarithmic algorithm. In particular, the smallest inclusion is identified comparatively well with both algorithms even though results are still clearer with the logarithmic algorithm (see Fig. 4(a) and (d) for the permittivity profiles, and Fig. 5(a) and (d) for the conductivity profiles). For increasing sizes, the logarithmic algorithm computes a very good reconstruction, in both the contrast achieved and the shape of the target. The complex algorithm quickly fails: the circular shape cannot be identified and the levels in the permittivity profile are underestimated.

Table 1.

Phantoms and background properties for the reconstructed images in Figs. 4–5. The background permittivity corresponds to an 80% glycerin bath while the slab permittivity and conductivity correspond to a 50% glycerin bath

	Radius [cm]	ε at 1100 MHz	σ at 1100 MHz
Phantom 1	0.5	63.9	1.75
Phantom 2	1.73	63.9	1.75
Phantom 3	3.47	63.9	1.75
Background bath	7.5	28.3	1.26

Fig. 4.

Reconstructed permittivity profile for three phantoms of increasing size (see Table 1) using either the logarithmic or the complex algorithms. Note the scale difference in (a)-(d-f) without which the target is barely visible.

Fig. 5.

Reconstructed conductivity profile for three phantoms of increasing size (see Table 1) using either the logarithmic or the complex algorithms.

3.2. Robustness to increasing frequency

In order to avoid spurious reflections from the tank boundaries, which would introduce signal errors at the receivers, we have opted to perform our measurements (and hence the simulations reported in this paper) in a lossy liquid, typically composed of some percentage or glycerin and water. The losses of the bath increase with frequency, e.g. from σ = 1.2595 S/m at 1100 MHz to σ = 1.7197 S/m at 1700 MHz for a bath composed of 80% glycerin. As a consequence, the signals from one transmitter to a receiver are significantly weaker at higher frequencies. For the sake of illustration, Fig. 6 shows the electric field in decibels radiated by one transmitter with the position of all the receivers indicated as black circles. The level at the receiver directly across the transmitters drops from about −90 dB at 1100 MHz to −124 dB at 1700 MHz. In addition to introducing hardware challenges at higher frequencies, this significant drop in field amplitude is a challenge for the reconstruction software. Typically, the complex algorithm emphasizes the fields that are much larger at the receivers closer to the transmitter. However, the important part of the information is often contained in the radiation propagating through the target [22]. When the complex formulation is used, this contribution is proportionally much smaller than the one corresponding to receivers in close proximity to the transmitter, which effectively carry very little information about the target. The information is used more effectively in the case of the logarithmic algorithm, which weights the fields transmitted through the targets more than those reflected off the target, collected in proximity to the transmitter. This weighting partially explains why the logarithmic algorithm systematically produces more accurate images, especially at higher frequencies.

Fig. 6.

Comparison between the field amplitudes in dB generated from the same transmitter at 1100 MHz and 1700 MHz. The black dots organized in a circular pattern denote the location of the 16 antennas.

The effect is illustrated in Fig. 7, where the permittivity profiles of a phantom (see the caption of the figures and Table 2 for details) are reconstructed using the logarithmic and the complex algorithms at frequencies ranging from 900 MHz to 1700 MHz, in 200 MHz increments. The logarithmic algorithm is able to retain the information about the target despite decreasing signal levels, and produces consistent images throughout the frequency range. The complex algorithm, however, faces the difficult task of reconstructing a permittivity profile with sequentially less information, and indeed, fails in reproducing the target which vanishes with increasing frequency.

Fig. 7.

Reconstructed permittivity profile for various frequencies of a phantom corresponding to a liquid of 90% of glycerin (and of radius radius=1.73 cm) in a background of properties corresponding to a bath of 80% of glycerin. Notice the disappearance of the target at higher frequencies when the complex algorithm is used. The properties of the two media are given in Table 2. Similar results have been observed on the conductivity images (results not shown).

Table 2.

Permittivity (ε) and conductivity (σ) of two different media as function of frequency

Frequency [MHz]	90% glycerin		80% glycerin
	ε	σ	ε	σ
900	18.54	0.775	32.01	1.053
1100	16.27	0.866	28.27	1.260
1300	14.78	0.941	25.43	1.435
1500	13.62	1.003	23.17	1.587
1700	12.71	1.062	21.32	1.720

3.3. Convergence rate

The inadequacy of the complex algorithm to produce accurate images in cases others than the simplest ones (small targets, low contrast, and low frequency) can be alleviated to some extent by providing good initial guesses, as shown for the example in [16]. The latter should be close enough to the true distribution, both in shape and level, so as to eliminate phase jumps between iterations. In medical imaging, however, a priori information on the property distributions of tissues is often not available. Moreover, access to prior information is typically possible only if the patient underwent tests using other methods such as magnetic resonance imaging. Hence, requiring a priori information increases the complexity, time, and cost of a diagnostic microwave imaging system and may only be warranted in certain circumstances. With its ability to fully utilize the information contained in the electromagnetic fields at each receiver, and in particular the absolute phase of the fields, the logarithmic algorithm is much less sensitive to initial guess inadequacies than the complex algorithm, which not only often converges to inaccurate solutions but also does so with a very low convergence rate.

We tested this effect by providing various initial guesses to the algorithms, in the form of a homogeneous centered permittivity and conductivity distribution over a circular region of radius 6 cm. Note that the position of the phantom is offset so that neither the boundary nor the center of the initial guess region overlap with the phantom.

Three initial guesses are studied. The first consists of assigning the region properties to be the same as those of the background:

This situation can be viewed as “no initial guess” and has been systematically used to generate the other results presented in this paper. The second initial guess assigns the region parameters as

while the third assigns the parameters as

In all three cases the solutions were identical to Fig. 4(c) and Fig. 5(c) for the permittivity and conductivity reconstructed with the logarithmic algorithm, respectively, and identical to Fig. 4(f) and Fig. 5(f) for the parameters reconstructed with the complex algorithm. It should be noted that even though the third initial guess is closer to the real values of the phantom (ε = 63.9 and σ = 1.75 S/m) than the other estimates, the complex algorithm is still not able to reconstruct a good image but converges to the same solution as obtained with the other two initial guesses.

The convergence rate of the two algorithms has been found to be very different: fast for the logarithmic algorithm to very slow for the complex algorithm. It typically took roughly 15 iterations for the logarithmic algorithm to converge, regardless of initial guess, which highlights the ability of this method to quickly compensate for possible erroneous biases away from correct values. On the other hand, the complex algorithm was found to converge very slowly, requiring about 100 to 200 iterations to reach the same final (and yet often inaccurate) solution. This difference in convergence rate is intrinsically related to our Newton-based method: the logarithmic transformation linearizes the problem and translates the first derivative of the functional (or Jacobian in our case) to a good estimate of the convergence direction toward the global minimum. Such is not the case when the complex algorithm is used since the functional remains highly nonlinear and multi-valued. As a consequence, the logarithmic algorithm is not only more accurate but also more efficient at reconstructing the correct image which are two important aspects for proper inversion of clinical data.

4. Clinical case

The purpose of the previous section was to illustrate how our conclusions based on the one dimensional analysis are directly relevant to a two-dimensional situation where locations, shapes, and properties of the inclusions were known. We now test these algorithms on a clinical case using data acquired by the system described in [29, 36] in a configuration whereby a patient was lying on top of the microwave system with her breast buoyant in the liquid bath within the imaging region (i.e. surrounded by the monopole antennas). A set of 240 measurements (16 transmitters × 15 receivers per transmitter) were collected at multiple planes, from the chest wall to the nipple region, and multiple frequencies, from 900 MHz to 1500 MHz. The patient was 36 years old with radiographically heterogeneously dense breasts and presented with a widely dispersed tumor spreading from near the nipple to the chestwall and occupying space in the lower to lower right quadrant when viewed en face with overall dimensions of 6.5×3.7×7.1 cm. Interestingly, and not too uncommonly for patients with larger tumors, the right breast exhibited a large region of skin thickening surrounding much of the breast. This patient ultimately underwent neoadjuvant chemotherapy; however these images are pre-treatment. Fig. 8 shows a sequence of three MR sagittal images for this case: (a) T2, (b) T1 with fat suppression after injection of gadolinium contrast agent, and (c) the subtraction between (b) and the same prior to injection of gadolinium. The T2 image is not very interesting in that it essentially highlights areas of higher water content including the edema under the skin, the tumor and other high-water features in the breast. The T1 image with fat suppression and gadolinium enhancement is interesting because it primarily highlights the tumor and the edema below the skin-thickened surfaces. The subtraction image is much more interesting with respect to isolating just the tumor (red arrows) since the subsurface edema (blue arrows) is the same for both images and is eliminated as part of the subtraction process.

Fig. 8.

Sagittal MR images of the right breast of the test patient prior to therapy.

For the microwave images, we show only the coronal 900 MHz thru 1500 MHz images for a single plane about 5 cm fromthe chestwall reconstructed with both the complex and logarithmic algorithms (see Figs. 9 and 10). The logarithmic algorithm permittivity images are intriguing because they nominally show the breast perimeter along with a large elevated internal zone corresponding to the tumor which appeared to enhance with increasing frequency. In addition, there is a pronounced elevated permittivity and conductivity zone along the breast perimeter and along the field of view boundary. Initially we thought this was a type of artifact we previously observed in isolated exams when the breast tissue actually contacted the antennas. However, in this case, it is clear that this feature represents the high dielectric and conductive properties of the edema (blue arrows in Fig. 8). This was confirmed at later imaging sessions (not shown) where, as the tumor responded to therapy, the thickened skin and associated edema degenerated to the point of completely disappearing. These results illustrate that the logarithmic algorithm is quite capable of recovering clinically relevant images of a very complicated target.

Fig. 9.

Permittivity images of clinical breast data for various frequencies.

Fig. 10.

Conductivity images of clinical breast data for various frequencies (same clinical case as in Fig. 9).

In contrast, the complex algorithm produces permittivity images of much lower resolution, albeit identifying the tumor as well as the skin thickening effect at the proper locations. These reasonable results, however, have been obtained at the expense of a convergence rate about five times slower than with the logarithmic algorithm (100 iterations as opposed to 20). In the corresponding conductivity images, the tumor only vaguely appears at lower frequencies and is barely identifiable at 1500 MHz. This phenomenon is similar to the one illustrated in Fig. 7 whereby the phantom disappears at higher frequencies due to the inability of the complex algorithm to lock onto the proper Riemann sheet (i.e. the proper value of m in our one dimensional analysis). On the other hand, the logarithmic algorithm not only retains a good reconstruction as frequency increases but also exhibits improved resolution as desired.

These clinical results are but one illustration of the importance of properly unwrapping the phase information, as demonstrated by our one-dimensional analysis. Analysis of additional cases (not shown) consistently drove toward the same conclusion and led us to conclude that the logarithmic algorithm is indeed superior to the complex one in reconstructing biomedical images from our hardware system.

5. Conclusion

Several advantages of the logarithmic algorithm over the complex algorithm have been noted previously [22, 27]. In the present paper, we have exposed an underlying reason for these advantages by studying a simplified one-dimensional configuration in which the solutions for the fields as well as the estimation parameters can be expressed in closed form. Specifically, we have shown that both the permittivity and the conductivity depend on the absolute phase of the fields, reflected by their dependency on the parameter m which is related to the number of periods of the electromagnetic wave that are included inside the medium through which the wave is propagating. The straightforward solution, m = 0, appears to be valid only at low normalized k0d parameters, i.e. at low frequencies and small sample sizes. This result explains why these cases do not present major imaging challenges. For higher frequencies or larger sample sizes, and equivalently for targets with high contrast relative to their background, the solution m = 0 is often not valid and a proper image can only be obtained if the phase of the signals at each receiver is unwrapped. The one dimensional study presented here illustrates how increasingly large m’s are needed to reconstruct the constitutive parameters at increasingly large k0d. In addition, the proper use of the phase information requires a logarithmic transformation of the fields (thus converting the phase into an additive term), which in turn requires identification of the Riemann sheet of the complex logarithm function. The resulting two dimensional algorithm is shown to be robust to target size, to increasing frequency, and to initial guess, in addition to having a much faster convergence rate, which are important advantages in a clinical system. As with both one and two dimensional implementations, this concept is fully generalizable to three dimensional configurations. While we have not tested other implementations beyond our FDTD/Gauss-Newton approach, we expect that stochastic approaches might also benefit from our logarithmic transformation. As shown by our one-dimensional analysis, unwrapped phase is a critical piece of information in extreme cases (high contrast, high frequency, large scatterers, very inhomogeneous media) that goes beyond the minimization method itself.

While the logarithmic algorithm is advantageous, it does impose an important practical requirement to unwrap both the measured and computed phase values. The unwrapping process for the computed values is relatively straight forward as long as the forward solution is computed at all points in the domain, as is the case with finite element and FDTD approaches. Unwrapping the measurement data is a more significant challenge depending on the system hardware implementation. For our system, we utilize monopole antennas that can operate from 500 MHz to 3 GHz when submerged in our lossy coupling liquid. Data acquired in small frequency intervals can be used to exploit the assumption that the associated phase changes between frequencies are small, from which a robust unwrapping strategy can be devised [27]. Use of the logarithmic algorithm may be more difficult when broadband data are not available.

The logarithmic algorithm raises a theoretical challenge that still needs to be resolved before truly high frequency images can be confidently reconstructed. The issue has been discussed in [28] and concerns the appearance of field nulls, i.e. points where the field amplitude drops dramatically. If such nulls appear within the field of view, ambiguities in phase unwrapping can occur since tracking the phase along two different paths may yield two values of m, in turn resulting in non-unique images. This problem is more prominent at higher frequencies where these nulls form in close proximity to the target. At lower frequencies, the theoretical possibility of nulls remains but they typically appear outside the field of view and do not disrupt the imaging process. Strategies to systematically overcome this problem are under investigation.