Abstract
The contrast source inversion (CSI) algorithm was introduced for microwave imaging in 1997 and has since proven to be one of the most successful algorithms for nonlinear microwave tomography. In the CSI algorithm, the nonlinear integral equation, which must be solved to extract the constitutive electromagnetic parameters of the object under test from the microwave measurements, is represented by two linear equations, known as the data and the object equations. In this paper, the data equation in the CSI algorithm is reformulated using the so-called log-phase formulation. In this formulation, the measured data is represented by the change in the logarithm of the amplitude and the change in the unwrapped phase. This formulation has previously been applied for nonlinear tomography within the framework of a Gauss-Newton based algorithm for detection of breast cancer. Here, significant improvements have been observed compared to the more commonly used real-imaginary formulation. The modified CSI algorithm is tested on both simulated data and on a measurement of a breast. It is shown that for imaging setups with large differences in the measured signals, the new formulation of the data equation significantly improves the performance of the CSI algorithm
1. Introduction
As far back as the late 1940s, researchers have reported on the contrast in constitutive electromagnetic parameters between healthy and cancerous breast tissue in the microwave region [1–5]. This contrast implies that the presence of a tumor will cause an incident electromagnetic field to scatter, thereby making it feasible to detect breast cancer by use of microwave imaging. Different approaches have been suggested for microwave imaging of the breast, with radar-based [6, 7] and tomographic approaches [8–10] being the most widespread. In the radar-based approaches, the position of point scatterers inside the breast is sought for by application of techniques inspired by those used in more classic radar applications. When using the tomographic approaches, one seeks to reconstruct the distribution of the constitutive electromagnetic parameters, that is, permittivity and conductivity, inside the breast.
Apart from the possible presence of a cancerous tumor in the breast, the different healthy tissue types in the breast will also cause scattering of the incident microwaves. This implies that the imaging problem is nonlinear and ill posed, making the use of microwave tomographic imaging a challenging task.
In 1997, the contrast source inversion (CSI) algorithm was introduced as a means for solving nonlinear tomographic microwave imaging problems [11]. In this algorithm, the electromagnetic inverse problem is formulated as a minimization problem using two coupled equations, known as the data and the object equations. In the data equation, the scattering objects are represented using the so-called contrast sources (from which the name contrast source inversion originates), while the object equation relates the physical contrast in constitutive parameters to the contrast sources. By iteratively solving the two equations, in turn updating the distribution of contrast sources and the distribution of constitutive parameters, the actual distribution of the constitutive parameters in the object under test can be reconstructed. Since its introduction, the algorithm has been successfully applied for microwave imaging in a number of different applications [12–15], including biomedical [16–18].
Apart from microwave imaging, the CSI algorithm has also been successfully applied in other areas, such as ultrasound, elastodynamics, and electrode logging [19–21].
One of the characteristic trades of the CSI algorithm is that it does not require the explicit solution of the forward scattering problem. This is different from the widely used Newton-based algorithms, in which each iteration of the algorithm requires the solution of the forward scattering problem [8, 9, 22, 23]. Since solving the forward scattering problem is most often a computationally expensive task, this feature is often considered one of the strongest arguments for using the CSI algorithm.
A number of modifications have previously been suggested for the CSI algorithm. These have most often been related to the regularization of the problem and include the use of a total variation term both as additive and multiplicative regularization [16, 24], the use of different algorithms for calculating the updates of the contrast sources and contrast [25], and other changes to the regularization scheme used in the algorithm [26]. In addition to this, multifrequency versions of the algorithm have been published [12].
While a large amount of literature has been published on different regularization schemes, the weighting of the measured data and reformulation of the scattering problem have received little attention. The authors of this paper have previously shown how the use of the so-called log-phase formulation, introduced in [22], improves the performance of Newton-based imaging algorithms [27–29], both in 2D and 3D. In the log-phase formulation, the cost function to be minimized in the imaging algorithm is formulated using the logarithm of the change in amplitude and the unwrapped change in the phase of the measured complex S-parameters. The improvement observed when using this formulation is in part caused by the fact that the log-phase formulation emphasizes relative changes, putting more weight on measurements made with antennas on opposite sides of the object under investigation, and in part by the use of the unwrapped phase, which implies that the algorithm is capable of handling phase changes of more than .
In this paper, a CSI algorithm, in which the data equation is reformulated using the log-phase formulation, is presented. In Section 2, the algorithm is presented, and in Section 3, the performance of the algorithm is illustrated by reconstructions of both simulated and measured data. The time notation , with being the imaginary unit, is assumed throughout this paper.
2. Algorithm
The algorithm presented in this paper is based on the CSI algorithm introduced in [24, Sec. 3], that is, the CSI algorithm using the Polak-Ribiére algorithm for updating the distribution of both the contrast sources and the contrast. Here, the algorithm will be derived for a two-dimensional, transverse magnetic imaging setup. This simplifies the expressions in that the electric fields can be expressed as scalars, and that the computational burden is considerably smaller than for a full three-dimensional inversion.
The CSI algorithm is based on the electric field integral equation which relates the known incident field transmitted by the imaging system and the unknown total field through the equation In this expression, is the Green's function for the background medium, and are position vectors, and the integration is performed over the imaging domain which is assumed to completely enclose the object under investigation. The subscript indicates that different antennas are used to irradiate the domain, yielding one distribution of the electric field for each transmitting antenna.
The background medium has the squared complex wave number given by with and being the permittivity and conductivity, respectively, of the background medium. The contrast in (1) is given by wherein is the squared complex wave number at the position . Hence, the contrast is zero if the permittivity and conductivity at a given position is equal to that of the background and nonzero at positions where an object with different constitutive parameters is present.
The expression in (1) contains the total field both on the left-hand side of the equation and in the integral on the right-hand side, leading to a nonlinear inversion problem with respect to the unknown distribution of the contrast . By introducing the contrast sources as the integral equation (1) can be rewritten as In this way, the nonlinear equation in (1) can be expressed as two coupled equations: one is the so-called data equation (5) which provides an expression for the linear relation between the contrast sources in the imaging domain and total field, , which can be measured by the receiving antennas of the imaging system.
The other equation is the so-called object equation which can be found by rewriting the expression in (4) as and provides a relation between the contrast sources and the actual contrast of the object in the imaging domain. This expression is seen to be nonlinear with respect to , which appears on both the left- and right-hand side of the equation.
In the CSI algorithm, these two equations are solved iteratively in turn, thereby minimizing the cost function [24, (1)] Here, is the measured data, and the subscripts and indicate the norms over the receiving antennas and the imaging domain, respectively. The operators and are the integral operators in (5) and (6), respectively, which map the influence of the contrast sources in the imaging domain to the measured data () and the influence from the contrast sources to the physical contrast in the imaging domain (). Finally, the summations are to be performed over all transmitters, and is a normalization factor given in [24, (11)].
This concludes the description of the basic CSI algorithm in this paper. The reader is referred to [24] and the paper in which the CSI algorithm was first introduced [11] for the complete details on the algorithm.
The new algorithm differs from that presented in [24] only in the way the measured data is represented. In the original algorithm, the data is used in its real-imaginary form, and the elements in the data equation are based on the difference between the signals measured for an empty measurement system and the signals measured with an object inserted. Expressed in terms of measured parameters, the data can be written as Herein and are the parameters measured with an empty system and the parameters measured with an object inserted in the imaging system, respectively.
In the new algorithm, the data equation is transformed from the real-imaginary formulation to the log-phase formulation [22]. This implies that the values of in the new algorithm are given as with indicating the unwrapped phase. When the natural logarithm is used, as is the case in this paper, the phase in this expression will be given in radians, while the real part will be the difference in the natural logarithm of the amplitude of the signals measured with and without the object in the imaging system.
The use of the log-phase formulation results in a nonlinear data equation (9) as opposed to the linear data equation (5) encountered when using the real-imaginary formulation. This implies that the input to the Polak-Ribiére algorithm, used to calculate the updates, must be changed slightly from those used in [24].
Following the notation in [24], the operator is introduced as the operator which maps the resulting field from the contrast sources to the receiving antennas in the log-phase formulation. This operator is most easily implemented by simply using the standard operator and then adding an additional layer in which the simple calculation in (9) is performed.
To determine the update directions and step length in the Polak-Ribiére algorithm, the complex conjugate of the derivative of the operator is needed. If the notation from [24] is used, the new updates can be determined by using the operator where is the operator used in [24, (17) and (19)], and is the complex conjugate of . Otherwise, the updates are determined in the same way as described in [24, Sec. 3].
In the simple implementation of the algorithm currently used by the authors, the computational demand of the log-phase formulation results in an increase of the computation time per iteration of just under . Hence, the added computational complexity should not be a limiting factor for the implementation of the algorithm.
To sum up, the algorithm presented here is identical to the algorithm presented in [24, Sec. 3] with the exception that the data equation is formulated using the log-phase formulation. This results in a slightly different calculation of the data error and requires that the adjoined operator used to determine the update direction of the contrast sources in [24, (17) and (19)] is replaced with the operator in (10).
3. Results
To illustrate the change in the performance of the algorithm caused by the reformulation of the data equation, three different imaging setups are presented. The first two are simulations of circular targets, and the third is a measurement of a breast with a relatively large tumor in it.
To allow for a comparison of the images, the RMS of the normalized error is introduced as where the complex permittivity is given by is the number of pixels in the imaging domain, and and are the true and reconstructed values of the complex permittivity in the individual pixels, respectively. Since the true distribution of the complex permittivity is not known for the patient measurement, the RMS of the normalized error will only be calculated for the simulated data sets.
Throughout this section, the CSI algorithm is set to terminate when the value of the cost function (7) reaches a value of 0.002. This value has been chosen ad hoc on the basis of a large number of image reconstructions performed with the imaging setup used in this paper. For other imaging systems, a different value may provide better results.
3.1. Imaging System
The imaging setup considered in this paper is described in [30] and consists of 16 monopole antennas positioned in a circular setup with a radius of . A schematic of the system is shown in Figure 1 and a photo is shown in Figure 2. The system is currently being used for breast imaging at Thayer School of Engineering at Dartmouth College and operates in the frequency domain, typically at a frequency between and . To maximize the coupling between the antennas and the interior of the breast, the imaging system is filled with a coupling liquid consisting of a mixture of glycerin and water.
When used for breast imaging, the patient is lying in a prone position with her breast suspended down into the center of the imaging system, and the antennas are moved vertically to seven different positions to image the entire length of the suspended breast. A thorough description of the system can be found in [30] and the references therein.
3.2. Simulation 1: Small Object
In the first of the simulated setups, the system is simulated at a frequency of , and the background has a relative permittivity of 20 and a conductivity of . These parameters of the permittivity and conductivity are similar to those found in a typical mixture of the glycerin-water coupling liquid used in the actual system.
A single cylindrical target is inserted into this background. The target has a radius of and a relative permittivity of 60 and a conductivity of and is located at . Gaussian noise simulating a noise floor below the level of the transmitted signal has been added to the simulated data. This corresponds to what is observed in the measurement system [30].
The inversion results are shown in Figure 3 wherein the dashed lines indicate the position of the scatterer. In all three examples presented in this paper, the imaging domain is divided into square pixels with a side length of . With the real-imaginary formulation, the CSI algorithm reaches the threshold value of the cost function after 421 iterations, while the algorithm using the log-phase formulation uses 954 iterations to reach the threshold. As mentioned above, the time used pr. iteration is approximately the same for the two different formulations, implying that the total calculation time for the log-phase formulation is approximately 2.3 times longer than that of the real-imaginary formulation.
(a) Perm., real-imag.
(b) Cond., real-imag.
(c) Perm., log-phase
(d) Cond., log-phase
The results obtained with the two different formulations are seen to be very similar, with the real-imaginary formulation reaching a slightly higher value of the permittivity in the center of the object. The log-phase formulation, on the other hand, yields slightly less artifacts in the background although this is not easy to see in the images.
In both of the formulations, the algorithm has trouble reconstructing the conductivity of the object. This is something which is often seen in microwave imaging, for example, [9, 17, 18, 23, 29].
The RMS of the normalized error for both the real-imaginary and the log-phase formulations is found to be 0.10. This corresponds well with the fact that the images reconstructed with the two different formulations look almost identical.
3.3. Simulation 2: Large Object
In the second setup, a second circular scatterer, surrounding the one modeled in the first setup, is included in the simulation. This scatterer has a radius of , is positioned with its center at , and has a relative permittivity of 10 and a conductivity of . Apart from the presence of this second scatterer, the simulation parameters are identical to those used in the first simulation.
The new object causes a large change in the measured signals. The phase of the signals measured with receivers on the opposite side of the imaging system from the transmitter changes more than 270 degrees when the object is inserted, and the amplitudes change more than . For the first simulated imaging case, described above, the corresponding values are a maximum change in phase of just over 50 degrees and a change in amplitude of approximately .
The results obtained using the two different formulations are shown in Figure 4. In these images, a clear difference is seen between the results obtained with the real-imaginary formulation and the results obtained with the log-phase formulation. In the permittivity images, the algorithm using the real-imaginary formulation is unable to reconstruct the position of the small object while it captures the circular shape of the large outer region. The log-phase formulation, on the other hand, captures the position and size of the small object much better and also have a permittivity value closer to that of the real value for the large circular scatterer. As with the first simulation, the conductivity images bear little resemblance with the actual object.
(a) Perm., real-imag.
(b) Cond., real-imag.
(c) Perm., log-phase
(d) Cond., log-phase
Using the real-imaginary formulation, the threshold value of the cost function was reached after 383 iterations while it took 994 iterations with the log-phase formulation. As a result, the log-phase formulation used approximately 2.6 times as long to reach a result as the real-imaginary formulation.
The RMS of the normalized error is 0.55 for the log-phase formulation and 1.24 for the real-imaginary formulation. This corresponds well with the fact that using the log-phase formulation results in images which, by visual inspection, seem to be in better agreement with the actual distributions.
3.4. Patient Measurement
To further investigate the performance of the algorithm using the log-phase formulation, it has been used for imaging the breast of a patient who had a tumor in her breast. The tumor was positioned relatively close to the chest wall of the patient at a 10 o'clock position when viewing en face. The measurement was performed at , and the coupling liquid (background of the reconstruction) had a relative permittivity of 13.8 and a conductivity of . To get a reference image, the Gauss-Newton-based algorithm described in [22] has also been applied for reconstructing images of the data.
In Figure 5, the images reconstructed from the data collected with the antennas positioned at the first of the seven imaging planes, closest to the chest wall, are shown. The imaging plane closest to the chest wall of the patient is usually the hardest to image since the cross-section of the breast is larger here than further away from the chest wall. Also, the effects of the presence of the chest wall implies that the assumption of the problem being two dimensional is compromised, leading to increased errors stemming from the two-dimensional modeling of the three-dimensional imaging problem. The measured signals change with as much as in amplitude and 280 degrees in phase when the breast is inserted in the imaging system compared to the measurement with the empty system, indicating that the imaging problem is highly nonlinear.
In Figure 5, the images reconstructed using the real-imaginary formulation are shown in (a) and (b), the images reconstructed using the log-phase formulation are shown in (c) and (d), and in (e) and (f), the results obtained using the Gauss-Newton algorithm are shown. For this imaging setup, the threshold value of the cost function was reached after 539 iterations with the algorithm using the log-phase formulation, while the algorithm using the real-imaginary formulation did not converge within the first 2000 iterations, at which point it was terminated.
In the images obtained using the log-phase formulation in the CSI algorithm and in the images obtained using the Gauss-Newton algorithm, the tumor is clearly visible in the permittivity images. Also, an increase in the reconstructed conductivity can be seen in the vicinity of the tumor, although the features in the conductivity images are more blurred than those in the permittivity images.
The images reconstructed with the real-imaginary formulation in the CSI algorithm show no clear contours. In the permittivity image, a ring of relatively low-valued artifacts shows the outline of the breast, but otherwise little useful information can be extracted from the image. In the conductivity image, high-valued artifacts can be seen along the rim of the imaging domain, close to the position of the antennas, but the image does not provide much information about the internal structure of the breast.
The results shown here and above for the second simulated data set clearly illustrate the improved performance of the CSI algorithm, when the data equation is reformulated using the log-phase formulation.
4. Discussion
In addition to the test cases shown here, the two different formulations of the CSI algorithm have been applied to a number of other simulations and patient measurements. For the patient measurements, the results are all similar to the one shown here, with the log-phase formulation yielding results very similar to those of the reference Gauss-Newton algorithm. The real-imaginary formulation, on the other hand, most often yields useless results, except for a few of the imaging planes positioned at the center of the breast. Here, the observed changes in amplitude and phase are low or moderate, and the two-dimensional modeling of the imaging problem is a better approximation than it is closer to the chest wall or nipple, where effects stemming from the three-dimensional nature of the imaging problem complicate the imaging problem.
Similar results are obtained for the simulations, where the use of the real-imaginary formulation yields good results in fewer iterations than the log-phase formulation for objects which results in only low or modest changes in the amplitude and phase of the measured signals. For imaging setups with larger changes in the measured signals, the log-phase formulation is capable of reconstructing satisfactory images of setups where the real-imaginary formulation fails to converge.
The reason for the improved performance of the CSI algorithm when using the log-phase formulation is believed to be the same as what has previously been observed when using the log-phase formulation in the context of a Newton-based algorithm [22, 29]. Here, it was found that the log-phase formulation outperforms the real-imaginary formulation in part because it puts more weight on those measurements where there is a large relative difference between the measurement of the object and the measurement of the empty system, whereas the real-imaginary formulation puts more weight on those measurements where the absolute change is large. This implies that the log-phase formulation compensates for any overall difference in the signal level, which might be present between two receive/transmit pairs. For example, a receive/transmit pair with the receiver positioned next to the transmitter will in general result in a higher signal level than and a receive/transmit pair with the receiver positioned on the opposite side of the imaging system from the transmitter, simply because of the distance between the antennas.
In the imaging setups presented in this paper, this implies that the measurements with the receivers positioned on the opposite side of the imaging system, which holds more information about the scatterer, are given more weight than those made with antennas positioned next to each other, thus leading to better images. A similar effect can be obtained by manually excluding antenna pairs, as illustrated in [31, 32].
Another reason for the improved performance of the log-phase formulation is that it utilizes the unwrapped phase of the signals. The phase can be unwrapped in different ways, as described in [22, 29], but the result is that the algorithm is better at handling imaging problems in which the phase change between the empty measurement and the object measurement is more than [28]. This corresponds well with the results in the previous section where the log-phase formulation has shown a clear improvement over the real-imaginary formulation in the reconstruction of the images of simulation 2 and the patient measurements in which the changes in the phase are also greater.
In this paper, the CSI algorithm has been set to terminate based on a threshold value of the cost function (7). A different approach would have been to simply terminate the algorithms after a fixed number of iterations, such as in [18, 24]. For comparison of the convergence rate between the two different formulations, however, the use of a threshold value is more convenient.
From the results in the previous section, it is noted that the CSI algorithm implemented with the log-phase algorithm requires a significant higher number of iterations to converge than the CSI algorithm implemented with the real-imaginary formulation. At least for those cases where the real-imaginary formulation actually does converge. Preliminary investigations indicate that, as a result of this, the Polak-Ribiére update step used in the current implementation of the algorithm is not as well suited for the log-phase formulation as it is for the real-imaginary formulation.
A contributing factor to the slower convergence of the algorithm using the log-phase formulation and the poorer performance of the Polak-Ribiére update algorithm could be that the data equation is nonlinear when the log-phase formulation is used. As mentioned in Section 2, this is different from the real-imaginary formulation, which results in a linear data equation.
To remedy the slower convergence of the log-phase formulation, an algorithm in which the update is calculated by solving a set of linear equations formulated as a matrix equation, using an overregularized solution, has been tested. This approach is similar to that presented in [25] and resulted in an algorithm in which the log-phase formulation generally converged in fewer iterations than the real-imaginary formulation and both formulations converged in fewer iterations than for the Polak-Ribiére algorithm. However, even though this algorithm resulted in both formulations converging in fewer iterations, the overall computation time increased compared to the simple update used in this paper due to the complexity of the update algorithm. Hence, the original Polak-Ribiére update algorithm was kept in this paper. The authors are, however, optimistic that by finding a more suitable update algorithm, this drawback of the log-phase formulation can be significantly reduced.
Finally, it should be noted that, in this paper, the most simple form of the CSI algorithm has been used in order to emphasize the change in the performance caused by the reformulation of the data equation. The new formulation of the algorithm can, of course, be combined with the extensions to the CSI algorithm which has already been published, such as the use of the multiplicative total-variation regularization [16] or multiple frequencies [12].
5. Conclusions
In this paper, a modified contrast source inversion (CSI) algorithm was introduced. The modification consisted of using the log-phase formulation for the data equation. This formulation has previously been applied in Gauss-Newton type algorithms and has been shown to improve the performance of the algorithm when reconstructing images of high-contrast targets.
For imaging setups with large differences in the measured signals, the test cases in this paper showed a clear improvement of the quality of the images when the formulation of the CSI algorithm was changed from the more commonly used real-imaginary formulation to the log-phase formulation.
Of special interest is the ability of the log-phase algorithm to reconstruct data from actual measurements. In this paper, it was shown that the CSI algorithm using the log-phase formulation is indeed capable of reconstructing images from measurements. And that the results are comparable to the images reconstructed with the Gauss-Newton-based imaging algorithm previously applied by the authors.
One drawback of the current implementation of the log-phase formulation in the CSI algorithm is that the algorithm is considerably slower to converge than when the real-imaginary formulation is used. To remedy this, alternatives to the Polak-Ribiére algorithm have been sought for. So far, a more efficient algorithm has not been found, but based on the initial work, the authors remain optimistic that a more suitable update algorithm can be implemented, thereby reducing the number of iterations needed to reach a solution with the log-phase formulation of the algorithm.