PMCCPMCCPMCC

Search tips
Search criteria 

Advanced

 
Logo of nihpaAbout Author manuscriptsSubmit a manuscriptHHS Public Access; Author Manuscript; Accepted for publication in peer reviewed journal;
 
IEEE Trans Antennas Propag. Author manuscript; available in PMC 2010 April 22.
Published in final edited form as:
IEEE Trans Antennas Propag. 2010 January 1; 58(1): 145–154.
PMCID: PMC2858419
NIHMSID: NIHMS151248

A Sparsity Regularization Approach to the Electromagnetic Inverse Scattering Problem

David W. Winters, Barry D. Van Veen, Fellow, IEEE, and Susan C. Hagness, Fellow, IEEE

Abstract

We investigate solving the electromagnetic inverse scattering problem using the distorted Born iterative method (DBIM) in conjunction with a variable-selection approach known as the elastic net. The elastic net applies both 1 and 2 penalties to regularize the system of linear equations that result at each iteration of the DBIM. The elastic net thus incorporates both the stabilizing effect of the 2 penalty with the sparsity encouraging effect of the 1 penalty. The DBIM with the elastic net outperforms the commonly used 2 regularizer when the unknown distribution of dielectric properties is sparse in a known set of basis functions. We consider two very different 3-D examples to demonstrate the efficacy and applicability of our approach. For both examples, we use a scalar approximation in the inverse solution. In the first example the actual distribution of dielectric properties is exactly sparse in a set of 3-D wavelets. The performances of the elastic net and 2 approaches are compared to the ideal case where it is known a priori which wavelets are involved in the true solution. The second example comes from the area of microwave imaging for breast cancer detection. For a given set of 3-D Gaussian basis functions, we show that the elastic net approach can produce a more accurate estimate of the distribution of dielectric properties (in particular, the effective conductivity) within an anatomically realistic 3-D numerical breast phantom. In contrast, the DBIM with an 2 penalty produces an estimate which suffers from multiple artifacts.

Keywords: Breast cancer, electromagnetic tomography, FDTD methods, inverse problems, microwave imaging, regularization

I. Introduction

Solving the electromagnetic inverse scattering problem involves estimating the distribution of dielectric properties within a volume V based upon observations of the scattered electromagnetic field. It is well known that the inverse scattering problem is both nonlinear and ill-posed [1]. Given an infinite number of completely precise and noise-free measurements, the inverse scattering problem has a unique solution [1]. However, in the real world, measurements will be finite in number and will be limited both in terms of accuracy and precision. In addition, the dielectric properties are typically discretized in some manner to simplify computations on digital computers. It is for these reasons that the inverse scattering problem does not have a unique solution in practice.

There are a number of methods for solving the inverse scattering problem, including conjugate gradient methods (e.g., [2]–[6]), the distorted Born iterative method (DBIM) and equivalent Gauss-Newton methods (e.g. [7]–[12]), contrast source inversion methods (e.g., [13]–[16]), and various others (e.g., [17], [18]). For reasons that will soon become apparent, the results of this paper are formulated in terms of the DBIM. The DBIM replaces the nonlinear inverse scattering problem with a series of linear approximations of the form

Aob
(1)

where o is a K × 1 vector containing the discretized estimate of the dielectric properties contrast, b is an M × 1 vector containing the measurements of the scattered electromagnetic field, and A is an M × K matrix. In each iteration of the DBIM, A and b are functions of the properties contrast estimated in the previous iteration. The final solution is given by the summation of the estimated contrast vectors from each iteration. An important point is that for typical discretization schemes and realistic measurement systems, M << K, in which case (1) is an underdetermined set of equations. Hence, the systems of linear equations (1) are typically very ill-conditioned and directly applying the method of least-squares to (1) at each iteration of the DBIM results in a solution which bears little resemblance to the true distribution of dielectric properties. Regularization is necessary to stabilize the problem and to define a unique solution [1].

A common approach to regularization with the DBIM involves solving the set of linear equations at each iteration via penalized least-squares [9]. The penalty is chosen to favor a solution of a particular form, such as those that are continuous or smooth. The relative strength of the penalty is controlled by a regularization parameter. When the penalty involves the 2 norm of o, the approach is referred to as ridge regression [19] (also known as Tikhonov regularization [20]). Ridge regression achieves reduced overall mean square error through a bias-variance trade-off [21]. A problem with ridge regression is that every element of o in the estimated solution will generally be non-zero, even if the true solution only involves a subset of the elements. Consequently the ridge regression solution may contain artifacts that are not present in the true distribution of dielectric properties in V, and this can decrease the imaging accuracy. This solution stratedgy is referred to in this paper as DBIM-RR.

Sparse approximation methods have been widely applied in the context of linear inverse problems of the form (1) (e.g., [22]–[26]). The goal is to find a solution that accounts for the observed data using only a small subset of the elements of o. If the true distribution of dielectric properties contrast within V is sparse (that is, only a small number of the elements of o are non-zero), then applying sparse approximation methods at each iteration of the DBIM should result in a solution that is free of the artifacts associated with ridge regression. While the contrast o is not necessarily sparse in general, we assume there exists a set of basis functions that approximate o in a sparse manner. Let

oΨφ
(2)

where Ψ is a K-by-R real-valued matrix of basis vectors, [var phi] is an R-by-1 complex-valued vector of coefficients, and R is the number of basis functions in the expansion. If o is sparse with respect to Ψ, then a small number of the elements of [var phi] are nonzero. Such approximations are common in other fields. For example, it is well known that wavelets can represent images with a small number of coefficients [27]. While some bases may lead to more sparse representations than others, one does not need to know the optimally sparse basis to reap the benefits of sparse approximation.

A direct approach to the sparse approximation problem involves penalizing the number of elements of [var phi] involved in the solution to (1), which requires an exhaustive search among all combinations of elements and thus is not computationally practical [25]. An indirect approach known as convex relaxation involves instead penalizing the 1 norm of [var phi]. A growing body of research indicates that the minimum 1 solution is equal to the sparsest solution (or very close to it) under a variety of circumstances (e.g., [24]–[26], [28], [29]). There are several numerical approaches for solving the 1 penalized least-squares problem, including greedy algorithms [30]–[32], and convex programming [33]–[35]. However, in our experience, applying these purely 1 regularized methods within the DBIM leads to numerical difficulties due to the problem size, ill-conditioning, and iterative aspects of the DBIM.

In this paper we investigate solving the electromagnetic inverse scattering problem using the DBIM in conjunction with a variable-selection approach known as the elastic net [21]. We denote this approach as DBIM-EN. The elastic net applies both 1 and 2 penalties to (1) in order to regularize each iteration of the DBIM. The elastic net thus incorporates both the stabilizing effect of ridge regression and the sparsity encouraging effect of the 1 penalty. One advantage of the elastic net is that the solution is computable even for very small 1 regularization parameters. We follow [21] and use a modified version of the moderately greedy algorithm known as Least Angle Regression (LARS) [31] for implementing the elastic net at each iteration of the DBIM. LARS efficiently produces a sequence of solutions which simplifies the selection of an 1 regularization parameter. We consider two very different 3-D numerical testbeds - a cube with an interior dielectric properties distribution that is exactly represented by a linear combination of a small number of wavelets, and an anatomically realistic breast phantom with heterogeneous dielectric properties not exactly sparse in any set of basis functions. Examples given using these testbeds demonstrate the efficacy and broad applicability of our elastic net approach to the inverse scattering problem.

We note that Baussard et al. [36] proposed an inverse scattering algorithm which is capable of producing a sparse solution. However, the algorithm in [36] is presented in the context of a specific set of multiresolution spline basis functions, and sparsity is achieved via a heuristic refinement process. In this paper we propose searching for a sparse solution to the inverse scattering problem for an arbitrary set of basis functions, and we encourage sparsity via the principle of convex relaxation.

Throughout this paper, electromagnetic field vectors and dyads are denoted by upper case letters with an overline (e.g., Ē). Position vectors are shown as lower case letters with an arrow overline (e.g., r), while all other vector quantities are indicated by lowercase boldface type (e.g., v). Matrices are denoted by uppercase boldface type (e.g., M); the matrix transpose and complex-conjugate transpose operations are represented by superscripts T and H respectively. The function notation f(x|y) denotes x as the independent variable and y as a parameter.

II. Methods

We begin this section by introducing the DBIM [9]. The brief review is followed by a discussion of regularization and sparsity, including descriptions of ridge regression and the elastic net. We conclude this section by discussing the choice of regularization parameters.

A. Distorted Born Iterative Method

Suppose that data is acquired using an antenna array located outside of V, the volume throughout which the dielectric properties are to be estimated. In sequence, each antenna transmits an electromagnetic signal into V while the other antennas in the array act as receivers. Consider the case where the mth antenna is transmitting an electromagnetic signal at angular frequency ω. The nonlinear integral equation that relates the continuous spatial distribution of dielectric properties within V to the scattered electric field at the nth receiving antenna is given by [9]

Escat(rnrm)=E(rnrm)Einc(rnrm)=ω2μ0VGb(rn,r,)E(r,rm)((r,)b)dr,
(3)

In Equation (3), Ēscat is the mathematically defined scattered field, Ē is the total field (known at rn but unknown inside V), and Ēinc is the known incident field. The position vectors of the mth transmitting and the nth receiving antennas are given by rm and rn, respectively. Inside the integral, Gb is the dyadic Green's function of the homogeneous background, while [sm epsilon] and [sm epsilon]b are the spatially-varying complex relative permittivity of the object and the spatially invariant complex relative permittivity of the background, respectively. The difference between the object and background relative permittivity is known as the contrast function, which is denoted by o(r) [9].

We solve this nonlinear problem by using a series of simplifying assumptions. Under the Born approximation [9], the integral in (3) is linearized by replacing the total electric field Ē in the integral with the known incident field Ēinc. The scalar approximation [3] assumes that only the z-directed component of the incident field is non-zero (Einc=[EincxEincyEincz]T=[00Eincz]T) and only the z-directed component of the electric field is measured by the receiving antennas. In theory the scalar approximation results in a loss of information, but in practice it has been shown that it does not significantly impact imaging performance [37]. These approximations yield the following simplified integral equation:

Escatz(rnrm)ω2μ0VGbzz(rn,r,)Eincz(r,rm)o(r,)dr,
(4)

where Gbzz(rn,r,) represents the z-z component of the Green's function dyad.

Equation (4) can be discretized via the Riemann sum under the assumption that all quantities are constant over volume elements (voxels) of volume V. Applying this discretization scheme to the set of approximations (4) for all transmit-receive pairs results in the following set of linear equations:

A0o1b0
(5)

In (5), A0 is an M-by-K matrix, where M is the number of transmit-receive pairs in the antenna array and K denotes the number of voxels. The K-by-1 vector o1 contains the dielectric properties contrast for the K voxels in V, while b0 is an M-by-1 vector and has elements equal to Escatz(rnqrmq)=Ez(rnqrmq)Eincz(rnqrmq), for q = 1, . . . , M. We emphasize that b0 does not lie entirely in the span of the columns of A0, due to the linear (Born) approximation.

Solving (5) results in a discrete approximation ô1 of the true distribution of contrast o(r). A better approximation can be obtained by adding ô1 to the background and using a series of computational electromagnetics simulations to calculate the new incident electric field and inhomogeneous Green's function based on ô1 + [sm epsilon]b1 [11]. The symbol 1 is here used to represent a vector of all ones. The following new set of linear equations is obtained upon substituting these updated field quantities into (4):

A1o2b1
(6)

Equation (6) is solved, resulting in an estimate ô2, and the process is repeated for multiple iterations. Full-wave computational electromagnetics simulations are conducted at every iteration in order to calculate the updated incident electric field and Green's function. These simulations are often collectively referred to as the ”forward solver.” We note that the scalar approximation described above in the context of (4) only applies to the inverse solution; that is, the z-component of the electric field and the z-z component of the Green's function dyad are recorded from the full-wave forward solution and incorporated into the linear system of (6).

The vector bi contains the residual fields after the ith iteration. Once the norm of bi ceases to decrease significantly from iteration to iteration, the DBIM has converged and the estimated contrast is given by o^=io^i. The estimated distribution of dielectric properties is thus [sm epsilon]b1 + ô.

B. Basis Functions, Regularization and Sparsity

We re-express the vector of dielectric properties contrast at any iteration of the DBIM as a linear combination of basis functions using (2). The DBIM for this basis function expansion consists of estimating the coefficients of the basis functions instead of the contrast at each voxel. The linear problem at each iteration of the DBIM becomes

AΨφb
(7)

While it is assumed that the number of measurements M is less than the number of discrete samples K, the number of basis functions R can be any size in theory. It is assumed that R is smaller than the single-precision numerical rank of A, which can be considered the effective number of measurements. In the examples considered in this paper, the rank of A is always around 100. Under these assumptions, (7) is an underdetermined system of equations. In addition, b does not lie entirely in the span of the columns of due to the linear approximation involved in obtaining (4). Solving the normal equations [38] results in the least-squares solution

φ^=(ΨHAHAΨ)1ΨHAHb
(8)

Owing to the ill-posed nature of the inverse scattering problem, is almost never full numerical rank in single precision arithmetic.

Tikhonov regularization and ridge regression stabilize (8) by adding in a fraction of the identity matrix [19], [20]:

φ^rr=(ΨHAHAΨ+γ2I)1ΨHAHb
(9)

The statistical interpretation of (9) is that the overall mean square error of φ^rr is reduced through a bias-variance trade-off [21]. In the context of the ill-posed inverse scattering problem, the addition of γ2I pads the smallest eigenvalues of ΨHAH and thus improves numerical stability. The disadvantage of ridge regression is that in general every element of φ^rr will be non-zero, even if the corresponding basis function plays no part in the true solution. This can introduce spurious artifacts in the estimated distribution of dielectric properties.

We address this issue by solving (7) via the elastic net [21]. This involves replacing (7) with the following optimization problem

minφAΨφb22+γ1φ1+γ2φ22
(10)

We denote the solution to (10) as φ^en. The elastic net can be thought of as a generalization of both ridge regression and of sparse approximation methods. When either γ1 or γ2 become very large, φ^en will tend to 0. When both γ1 and γ2 go to zero, then φ^en becomes (8) which is not computable. When just γ1 goes to zero, φ^en approaches (9). As γ2 goes to zero, then (10) becomes an 1-based sparse approximation problem. However, we note that the nonlinear and ill-posed nature of the inverse scattering problem can create numerical difficulties for greedy and convex optimization algorithms searching for a solution with γ2 = 0. For instance, convex optimization approaches typically need to be initialized with a solution that already satisfies the system of linear equations [35]. Finding such an initialization becomes difficult since (8) is not computable. Greedy algorithms that build up a solution one basis function at a time and use least-squares to determine the coefficients [31], [32] may be unstable. These difficulties can be eliminated by using a non-zero value for γ2 in (10).

We follow [21] and use a modified version of the moderately greedy algorithm known as Least Angle Regression (LARS) [31] for implementing the elastic net at each iteration of the DBIM. The modified LARS algorithm forms a sequence of solutions with monotonically increasing 1 norm, starting with the null solution and ending with φ^rr. This allows γ1 to be replaced by a new parameter s which is equal to the desired 1 norm of the solution [21]; obviously s is bounded from above and below by φ^rr1 and 0 respectively. A sparse solution is obtained by selecting an itermediate value for s. An advantage of using LARS is that its greedy nature simplifies the selection of a suitable value for s. Another advantage is that LARS allows a previously selected element from [var phi] to be discarded if it is determined not to be suitable for the solution. This attribute makes LARS more robust than other greedy algorithms such as orthogonal matching pursuits [31].

C. Selection of Regularization Parameters

We use the L-curve principle [39] to determine the regularization parameters for the DBIM-EN and DBIM-RR. The L-curve determines a suitable regularization parameter by balancing the norm of the linearized residual [var phi]b and the norm of the solution [var phi] [39]. The L-curve has been shown to be very robust in practice [39], although it is not optimal in any sense. The L-curve applied to ridge regression balances the 2 norms of the linearized residual and of the solution. There are two regularization parameters to determine for the elastic net: γ2 and s. We first determine γ2 by applying the L-curve to all of the basis functions at once. With this choice for γ2, we again use the L-curve to choose s by balancing the 1 norms of the linearized residual and of the solution. We justify the use of this two step procedure based upon the empirical observation that the shape of the 1 L-curves do not change appreciably for a wide range of values of γ2. We choose all regularization parameters using the point on the L-curve closest to the intersection of lines that are fit to points from the initial and tail regions of the L-curve.

III. Examples

We apply the DBIM-EN and the DBIM-RR to data simulated using two different 3-D computational testbeds. The first testbed consists of a generic object - a simple cube - whose distribution of dielectric properties is exactly represented by a linear combination of a small number of wavelet basis functions. We compare the DBIM-EN with those obtained using the DBIM-RR for ten different distributions of dielectric properties within the object. The second testbed consists of an anatomically realistic breast phantom with dielectric properties that correspond to the microwave frequency range. A previously reported inverse scattering algorithm also based upon the DBIM [40] is applied to the simulated data for the purpose of generating a low-resolution estimate of the distribution of dielectric properties within the breast phantom. The DBIM-EN and DBIM-RR estimate higher-resolution details within the breast phantom using this low-resolution estimate as an initial guess. For both examples we use a fast implementation of the modified LARS algorithm available at [41].

A. Microwave Imaging of a Heterogeneous Lossy Dielectric Cube

The purpose of the first testbed is to demonstrate the performance of our sparsity approach to inverse scattering for the ideal case where the distribution of dielectric properties is sparse in a known set of basis functions. We first describe the computational testbed, which is shown in Fig. 1. Then we present qualitative and quantitative results which indicate that the elastic net outperforms the ridge regression approach for this scenario.

Fig. 1
The first computational testbed consists of a 6.4 × 6.4 × 6.4 cm3 cubic object surrounded by a 40-element antenna array of 1.4-cm-long dipoles. The dielectric properties distribution within the cube is generated by a linear combination ...

1) Testbed

The testbed consists of a 6.4 × 6.4 × 6.4 cm3 lossy dielectric cube. The distribution of dielectric properties within the cube is an exact linear combination of 13 3-D Haar wavelets [42]. The cube is immersed in a coupling medium (relative permittivity: [sm epsilon]r = 2.6, conductivity: σ = 0.0 S/m) and is surrounded, as shown in Fig. 1, by a 40-element cylindrical antenna array consisting of five elliptical rings (10.4 × 11.6 cm) of eight electrically-small z-directed dipole antennas. The ring spacing in the z-dimension is 1 cm. Each dipole antenna is modeled by two segments of 6-mm-long copper wire separated by a 2 mm gap. Physical interaction between any of the 1.4-cm-long array elements is minimized by offsetting the placement of the dipoles in each ring by 22.5° from the placements in the neighboring rings.

A finite-difference time-domain (FDTD) computational electromagnetics simulation [43] is conducted to acquire microwave signals measured at all recording antennas for every transmitting antenna in the array. The spatial grid cell size in these simulations is 2 mm. In each simulation, a different antenna array element is excited with a modulated Gaussian pulse (−20 dB bandwidth: 500 MHz to 3.5 GHz) applied at the feed point. The bandwidth of the radiated signal 5 cm away in the coupling medium is 875 MHz to 3.75 GHz. The z-directed electric fields at the feed point of the other antennas in the array are recorded and transformed to phasors at 2.1 GHz.

The distribution of dielectric properties within the cube is given by a 3-D Haar wavelet expansion consisting of a scaling function and 3 levels of wavelets (R = 512 basis functions). The scaling function and the first level of wavelets have support over the entire cube, while the wavelets from the second and third levels have support over sub-volumes of 3.2 × 3.2 × 3.2 cm3 and 1.6 × 1.6 × 1.6 cm3, respectively. Each of the orthonormal basis functions in the expansion is normalized by the amplitude of the scaling function.

We conduct a study consisting of ten trials, where in each trial the scaling function and 12 of the Haar wavelets are randomly chosen to represent the true distribution of dielectric properties within the cube. The complex coefficient for the scaling function is chosen in each trial so that the baseline dielectric properties within the cube are [sm epsilon]r = 10.0 and σ = 0.125 S/m. Four wavelets are then randomly selected from each of the three levels in the Haar expansion. The coefficients for these randomly selected wavelets are ±(0.25 – j0.0535), where the sign is chosen randomly. For each of the ten trials, [sm epsilon]r ranges from 7.04 to 12.96, and σ ranges from 0.0511 to 0.199 S/m. An x-y cross section (z = 2.0 cm) of the distribution of [sm epsilon]r within the cube for one of the ten trials is shown in Fig. 2(a).

Fig. 2
Cross sections in the x-y plane (z = 2 cm) showing the results for one of the ten trials for the first computational testbed shown in Fig. 1. (a) The true distribution of relative permittivity; (b) Estimated permittivity using the DBIM in conjunction ...

2) Results

We use FDTD with a spatial grid cell size of 2 mm as the forward solver for the DBIM. Note that this is the same spatial grid cell size used for the “measurement” FDTD simulations. Committing this so-called inverse crime [1] allows the field residual (the measurement error, discussed later) to be driven to zero, which is desirable from the point of view of being able to set an upper bound on performance. This, together with the lack of measurement noise, allows the impact of ridge regression and the elastic net to be isolated. The example in Section III-B considers the case where the inverse crime is not committed. In these ten trials we make use of some a priori information about the cubic object. We use the DBIM to optimize just over the set of the 512 Haar wavelets; this implies knowledge of the location of the cube boundary. We also assume knowledge of the baseline dielectric properties of the cube, which are used as an initial guess by the DBIM to speed up convergence of the DBIM.

We consider the performances of the DBIM-EN and DBIM-RR and compare them to the ideal solution calculated with what is known as the oracle [42]. The oracle involves solving (7) by considering only the Haar wavelets which are known to be present in the true solution. Since this knowledge will never be available in practice, the oracle provides a lower bound on the error in the estimated coefficients (φφ^2). We simulate data acquisition with FDTD and then apply the DBIM to data collected for each of the ten trials using the oracle, ridge regression, and the elastic net. The DBIM converges in six iterations for all trials and approaches. The oracle approach produces solutions which involve only 13 of the 512 Haar wavelets. In contrast, the ridge regression solutions involve all 512 Haar wavelets. The elastic net solutions involve between 132 and 208 wavelets at each iteration with a mean of 162.

The results of the ten trials for all three approaches are summarized in Table I. The normalized field residual b62b021 indicates how well the three approaches fit the scattered electromagnetic field data. The normalized coefficient error φφ^2φ21 measures how close the estimates are to the true solution. Table I lists the mean values and standard deviations for these two metrics. All three approaches achieve similar levels of performance with regards to fitting the scattered electromagnetic field data. The slightly lower field residuals for the approaches using ridge regression and the elastic net can be explained by the extra degrees of freedom afforded those solution strategies and the ill-posed nature of the nonlinear inverse scattering problem. As expected, the oracle performs the best in terms of coefficient error since it only optimizes over the 13 Haar wavelets that are actually involved in the true solution. The mean coefficient error for the elastic net is approximately four standard deviations lower than the corresponding value for ridge regression. This indicates that the DBIM-EN results in a significant improvement over the DBIM-RR. The improved performance is evident in Fig. 2, which shows x-y cross sections of estimated [sm epsilon]r at 2100 MHz for one of the ten trials.

TABLE I
A summary of the results from the ten trials with the heterogeneous lossy dielectric cube testbed shown in Fig. 1. The performances of the DBIM using the oracle, ridge regression, and the elastic net are compared. The metrics shown are the normalized ...

B. Microwave Imaging of a Heterogeneous Breast Phantom

The purpose of the second testbed is to demonstrate the performance of our sparsity approach to the inverse scattering problem for a more realistic situation where the distribution of dielectric properties is not exactly sparse in a chosen set of basis functions. We consider an example from the field of microwave breast imaging [44]–[47]. Interest in this field has been fueled by data suggesting that the dielectric properties of breast tissue at microwave frequencies [48]–[52] are sensitive to certain physiological factors of clinical interest, such as water content, temperature, and vascularization. Consequently microwave tomography has the potential for characterizing normal breast tissue density - an important factor in assessing a patient's risk of breast cancer [53] - as well as detecting and monitoring malignancies - the focus of this example. We first describe the testbed containing the breast phantom shown in Fig. 3 and then present results which show that the elastic net approach again outperforms ridge regression.

Fig. 3
The computational testbed for the second example consists of an anatomically realistic numerical breast phantom surrounded by a 40-element antenna array of 1.4-cm-long dipoles.

1) Testbed

The numerical breast phantom shown in Fig. 3 is derived using a 3-D MRI dataset from a patient with “scattered fibroglandular” breast tissue, based on the American College of Radiology's BI-RAD system [54]. The scattered fibroglandular breast tissue is evident in the three cross sections of Fig. 4(a), (c), and (e). These orthogonal cross sections pass through the center of a 1-cm-diameter inclusion that has been added to the phantom to represent a malignant lesion.

Fig. 4
Three orthogonal cross sections of relative permittivity at 2100 MHz for the breast phantom testbed. The left column [(a),(c),(e)] shows the true distribution of relative permittivity at a spatial resolution of 0.5 mm. The right column [(b),(d),(f)] shows ...

Data is acquired using FDTD as with the testbed of Fig. 1, save for a few differences. The spatial grid cell size in this testbed is 0.5 mm; this smaller grid cell size is required to resolve the fine geometrical features of the breast phantom. The phantom testbed uses an antenna array whose elliptical rings have dimensions 9.6 × 12.4 cm. These array dimensions ensure that no antenna is closer than 1 cm to the surface of the breast phantom.

The numerical breast phantom is created following a procedure similar to those reported in [55], [56]. The intensity of the voxels in the MRI dataset is converted to dielectric properties via a piecewise linear mapping [6], [57]. The interior of the breast phantom is segmented into three distinct regions: adipose, fibroglandular, and transition. We adopt the dielectric properties reported in a recent large scale dielectric spectroscopy study [52] for the adipose and fibroglandular regions in the breast phantom.

Lazebnik et al. [52] reported the microwave-frequency dielectric properties for three breast tissue groups (Groups 1, 2, and 3) defined by their adipose content. Samples from Group 3 were composed primarily of adipose tissue with relatively low dielectric properties, and so we assign the dielectric properties reported for Group 3 to the adipose regions of the phantom. Groups 1 and 2 were comprised of samples with smaller amounts of adipose tissue and correspondingly higher dielectric properties, with samples from Group 1 having the least amount of adipose tissue and the highest average dielectric properties [52]. We choose to assign the properties from Group 2 to the fibroglandular region of the phantom, although using the properties from Group 1 would constitute a more challenging problem. We note that the combined use of Group 2 and Group 3 properties is still much more realistic than those considered in any other previous theoretical study of microwave imaging for breast cancer detection.

MRI voxel intensities in the adipose and fibroglandular regions are mapped to +/− 40% ranges about the mean properties assigned to each tissue type. Voxels in the transition region are mapped to the range spanning the maximum of the adipose range to the minimum of the fibroglandular range. The 2-mm-thick skin layer is modeled using the dielectric properties for dry skin [58]. The dielectric properties assigned to the spherical inclusion are adapted from [59] and are representative of malignant breast tissue properties in our frequency range of interest. Table II gives the dielectric properties used in the FDTD model at 2.1 GHz (the frequency at which the scattered signals are recorded).

TABLE II
Dielectric properties (at 2.1 GHz) of the various media present in the heterogeneous breast phantom testbed of Fig. 3.

2) Results

Note that in this example the spatial resolutions of the data-generation and the DBIM forward-solver FDTD simulations are not the same. The grid cell size is 0.5 mm for the data-generation FDTD simulations and 2 mm for the forward solver. This results in an inherent mismatch between the two sets of simulations and produces as much as 15% difference even when the exact same distribution of dielectric properties is simulated. The mismatch is partially corrected through the use of the calibration procedure proposed by Meaney et al. [44], but because this calibration procedure is not perfect, the inverse crime is not being committed in this example. The data-generation FDTD simulations are run in parallel on a computing cluster, while the forward solver is run on a desktop computer utilizing Acceleware hardware acceleration technology [60]. The forward solver simulations require about 30 seconds per antenna, or approximately 20 minutes per DBIM iteration.

Figure 4(b),(d),(f) shows the low-resolution estimate of the relative permittivity distribution at 2.1 GHz that we use as the initial guess in the DBIM-RR and DBIM-EN. The initial estimate of the conductivity is similar in appearance, except the grayscale spans 0 - 1.9 S/m. These estimates are obtained using the low-resolution inverse scattering algorithm reported in [40]. This algorithm, which is also based upon the DBIM, estimates the properties using a smooth set of basis functions with a nominal resolution of about 1 cm. The presence of the spherical inclusion is apparent in Fig. 4, although we note that the estimated contrast of the scatterer is less than the true contrast.

We investigate the feasibility of generating a higher-resolution image of the interior of the breast phantom using the DBIM-RR and DBIM-EN with a set of 3-D Gaussian basis functions, constructed as follows. A cuboidal volume (7.2 × 11.2 × 6.8 cm3) enclosing the breast phantom is first defined. Fifteen 1-D Gaussians are defined along each of the axes of the cuboidal volume. The standard deviations for these 1-D Gaussians are 3.6, 5.6, and 3.4 mm for the x, y, and z axes, respectively; the 1-D Gaussians are spaced along each axis of the cuboidal volume by 4/3 standard deviations. Three-dimensional Gaussian basis functions are defined using all 3375 combinations of the 1-D Gaussians and the Kronecker product [42]. For simplicity, we only include basis functions which have at least 95% of their support within the breast phantom interior; these chosen basis functions are truncated so that they are entirely supported within the breast phantom interior. This reduces the number of 3-D Gaussian basis functions from 3375 to 737. These 3-D Gaussian basis functions are not optimal for this example in any sense; they are selected for simplicity and to demonstrate that our sparsity approach to inverse scattering can be successful even when the true distribution of dielectric properties is not exactly sparse in the chosen set of basis functions.

Figures 5 and and66 show the results from applying the DBIM-RR and DBIM-EN using these 737 3-D Gaussian basis functions and the low-resolution initial guess of Fig. 4(b),(d),(f). The elastic net approach takes six iterations to converge, while the ridge regression approach converges after three iterations. The two sets of images of estimated relative permittivity at 2100 MHz (Fig. 5) appear similar, although we note that the contrast of the spherical inclusion is slightly higher in the elastic net estimate. The two sets of images of estimated effective conductivity (Fig. 6) are very different from each other. The elastic net estimate is sharper and has higher contrast, while the ridge regression estimate falsely indicates the presence of additional high contrast scatterers. These artifacts are suppressed in the elastic net estimate since on average only 234 of the 737 3-D Gaussian basis functions are involved in the solution at each iteration of the DBIM-EN.

Fig. 5
Three cross sections of the estimated relative permittivity at 2100 MHz obtained using the DBIM-RR [(a),(c),(e)] and the DBIM-EN [(b),(d),(f)]. The output of a low-resolution inverse scattering algorithm [40], shown in Fig. 4(b),(d),(f), is used as an ...
Fig. 6
Three cross sections of the estimated effective conductivity at 2100 MHz obtained using the DBIM-RR [(a),(c),(e)] and the DBIM-EN [(b),(d),(f)]. The output of a low-resolution inverse scattering algorithm [40] is used as an initial guess for both algorithms. ...

IV. Summary and Conclusion

We demonstrated the feasibility of solving the electromagnetic inverse scattering problem using a basis function formulation of the DBIM in conjunction with a variable-selection approach known as the elastic net. The elastic net applies both 1 and 2 penalties to the system of linear equations that result at each iteration of the DBIM. The combined 1 and 2 regularizations stabilize the inverse problem and promote sparsity in the solution. A more typical approach known as ridge regression only involves the 2 penalty. We applied the DBIM-EN to data simulated using two different 3-D computational testbeds, and compared results with those obtained using the DBIM-RR.

First we presented reconstructions of a geometrically simple object whose heterogeneous dieletric properties distribution is exactly represented by a linear combination of a small number of wavelet basis functions (13 out of a possible 512). The DBIM-RR solution made use of all 512 basis functions, while the DBIM-EN solution made use of a subset of the basis functions. We used a representative set of reconstruction cross-sections to illustrate that the DBIM-EN performance is qualitatively closer than the DBIM-RR to that of the ideal approach, namely the DBIM with an oracle. We also conducted quantitative comparisons and demonstrated that the DBIM-EN consistently performed better than the DBIM-RR across ten trials involving different dielectric properties distributions within the object.

Second we presented reconstructions of a geometrically complex object - an anatomically realistic breast phantom - to demonstrate the feasibility of using the DBIM-EN for microwave breast imaging. We initialized the DBIM-EN as well as the DBIM-RR with a low-resolution estimate of the dielectric properties distribution within the breast phantom obtained using the inverse scattering algorithm reported in [40]. The DBIM-RR solution made use of the entire set of 737 3-D Gaussian basis functions chosen for this formulation, while the DBIM-EN solution made use of less than one third of the basis functions. The DBIM-EN produced an enhanced contrast image of a 1-cm-diameter inclusion, while DBIM-RR reconstructions falsely indicated the presence of additional scatterers.

Acknowledgments

The authors would like to thank K. Sjöstrand for permission to use their fast implementation of the modified LARS algorithm [41].

This work was supported by the Department of Defense Breast Cancer Research Program under award W81XWH-05-1-0365, the National Institutes of Health under grant R01 CA112398 awarded by the National Cancer Institute, and the National Science Foundation under grant BES 0201880.

Biography

• 

David W. Winters was born in Albany, NY in 1979. He received the B.S. degree from Rensselaer Polytechnic Institute in 2002 and the M.S. and the Ph.D. degrees from the University of Wisconsin-Madison in 2004 and 2007, all in electrical engineering. In 2004, Dr. Winters was awarded a three-year Department of Defense Breast Cancer Research Program predoctoral traineeship. Dr. Winters is currently employed by The MITRE Corporation

Barry D. Van Veen Barry D. Van Veen (S’81-M’86-SM’97-F’02) was born in Green Bay, WI. He received the B.S. degree from Michigan Technological University in 1983 and the Ph.D. degree from the University of Colorado in 1986, both in electrical engineering. He was an ONR Fellow while working on the Ph.D. degree.

In the spring of 1987 he was with the Department of Electrical and Computer Engineering at the University of Colorado-Boulder. Since August of 1987 he has been with the Department of Electrical and Computer Engineering at the University of Wisconsin-Madison and currently holds the rank of Professor. His research interests include signal processing for sensor arrays, magneto- and electroencephalography, and biomedical applications of signal processing.

Dr. Van Veen was a recipient of a 1989 Presidential Young Investigator Award from the National Science Foundation and a 1990 IEEE Signal Processing Society Paper Award. He served as an associate editor for the IEEE Transactions on Signal Processing and on the IEEE Signal Processing Society's Statistical Signal and Array Processing Technical Committee and the Sensor Array and Multichannel Technical Committee. He is a Fellow of the IEEE and received the Holdridge Teaching Excellence Award from the ECE Department at the University of Wisconsin in 1997. He coauthored ”Signals and Systems,” (1st Ed. 1999, 2nd Ed., 2003 Wiley) with Simon Haykin.

Susan C. Hagness Susan C. Hagness (S’91-M’98-SM04-F’09) received the B.S. degree with highest honors and the Ph.D. degree in electrical engineering from Northwestern University, Evanston, IL in 1993 and 1998, respectively. While pursuing the Ph.D. degree, she was a National Science Foundation (NSF) Graduate Fellow and a Tau Beta Pi Spencer Fellow.

Since August 1998, she has been with the Department of Electrical and Computer Engineering at the University of Wisconsin-Madison, where she currently holds the position of Philip D. Reed Professor. She is also a faculty affiliate of the Department of Biomedical Engineering. Her current research interests include microwave imaging, sensing, and thermal therapy techniques including breast cancer detection and treatment, electromagnetic inverse scattering, ultrawideband radar, dielectric spectroscopy, bioelectromagnetics, FDTD theory and applications in biology and medicine, and global modeling of carrier-field dynamics.

Dr. Hagness served as an elected member of the IEEE Antennas and Propagation Society (AP-S) Administrative Committee from 2003 to 2005 and as an Associate Editor for the IEEE Antennas and Wireless Propagation Letters from 2002 to 2007. She is currently serving as the Chair of Commission K of the United States National Committee (USNC) of the International Union of Radio Science (URSI), and the Chair of the IEEE AP-S New Technologies Committee. She was the recipient of the Presidential Early Career Award for Scientists and Engineers presented by the White House in 2000. In 2002, she was named one of the 100 top young innovators in science and engineering in the world by the Massachusetts Institute of Technology (MIT) Technology Review magazine. She is also the recipient of theUniversity of Wisconsin Emil Steiger Distinguished Teaching Award (2003), the IEEE Engineering in Medicine and Biology Society Early Career Achievement Award (2004), the URSI Isaac Koga Gold Medal (2005), the IEEE Transactions on Biomedical Engineering Outstanding Paper Award (2007), the IEEE Education Societys Mac E. Van Valkenburg Early Career Teaching Award (2007), and the University of Wisconsin System Alliant Energy Underkofler Excellence in Teaching Award (2009).

References

1. Colton D, Kress R. Inverse Acoustic and Electromagnetic Scattering Theory. 2nd ed. Springer-Verlag Berlin Heidelberg; New York: 1998.
2. Kleinman RE, van den Berg PM. Iterative methods for solving integral equations. Radio Science. 1991;26(1):175–181.
3. Bulyshev AE, Souvorov AE, Semenov SY, Svenson RH, Nazarov AG, Sizov YE, Tatsis GP. Three-dimensional microwave tomography. Theory and computer experiments in scalar approximation. Inverse Problems. 2000;16:863–875.
4. Rekanos IT. Inverse scattering in the time domain: An iterative method using an FDTD sensitivity analysis scheme. IEEE Trans. Magn. 2002;38(2):1117–1120.
5. Takenaka T, Zhou H, Tanaka T. Inverse scattering for a three-dimensional object in the time domain. Journal of the Optical Society of America: A. 2003;20(10):1867–1874. [PubMed]
6. Winters DW, Bond EJ, Veen BDV, Hagness SC. Estimation of the frequency-dependent average dielectric properties of breast tissue using a time-domain inverse scattering technique. IEEE Trans. Antennas Propag. 2006;54(11)(2):3517–3528.
7. Chew WC, Wang YM. Reconstruction of 2-dimensional permittivity distribution using the distorted Born iterative method. IEEE Trans. Med. Imag. 1990;9(2):218–225. [PubMed]
8. Paulsen KD, Meaney PM, Moskowitz MJ, Sullivan JM. A dual mesh scheme for finite-element based reconstruction algorithms. IEEE Trans. Med. Imag. 1995 Sep;14(3):504–514. [PubMed]
9. Chew W. Waves and Fields in Inhomogeneous Media. IEEE Press; Piscataway, NJ: 1995.
10. Remis RF, van den Berg PM. On the equivalence of the Newton-Kantorovich and distorted Born methods. Inverse Problems. 2000;16:855–870.
11. Cui TJ, Chew WC. Inverse scattering of two-dimensional dielectric objects buried in a lossy earth using the distorted Born iterative method. IEEE Trans. Geosci. Remote Sens. 2001;39(2):339–346.
12. Zaeytijd JD, Franchois A, Eyraud C, Geffrin JM. Full-wave three-dimensional microwave imaging with a regularized Gauss-Newton method - theory and experiment. IEEE Trans. Antennas Propag. 2007 Nov;55(11):3279–3292.
13. van den Berg PM, Kleinman RE. A contrast source inversion method. Inverse Problems. 1997 Dec;13(6):1607–1620.
14. Abubakar A, van den Berg PM, Mallorqui JJ. Imaging of biological data using a multiplicative regularized contrast source inversion method. IEEE Trans. Microw. Theory Tech. 2002;50(7):1761–1771.
15. Zhang ZQ, Liu QH. Three-dimensional nonlinear image reconstruction for microwave biomedical imaging. IEEE Trans. Biomed. Eng. 2004;51(3):544–548. [PubMed]
16. Kooij BJ. Non-linear microwave tumor imaging in a heterogeneous female breast. Topical Meeting on Biomedical EMC, 16th International Zurich Symposium on Electromagnetic Compatibility; Zurich, Switzerland. Feb. 2005.pp. 53–58.
17. Caorsi S, Gragnani GL, Medicina S, Pastorino M, Zunino G. Microwave imaging based on a Markov random-field model. IEEE Trans. Antennas Propag. 1994 Mar;42(3):293–303.
18. Nakhkash M, et al. Application of the multilevel single-linkage method to one-dimensional electromagnetic inverse scattering problem. IEEE Trans. Antennas Propag. 1999;47(11):1658–1668.
19. Hoerl AE, Kennard RW. Ridge regression: Biased estimation for nonorthogonal problems. Technometrics. 1970;12(3):55–67.
20. Tikhonov AN, Arsenin VT. Solutions of Ill-Posed Problems. Winston; Washington, D.C.: 1977.
21. Zou H, Hastie T. Regularization and variable selection via the elastic net. J. R. Statist. Soc. B. 2005;67:301–320.
22. Levy S, Fullagar PK. Reconstruction of a sparse spike train from a portion of its spectrum and application to high-resolution deconvolution. Geophysics. 1981;46(9):1235–1243.
23. Tibshirani R. Regression shrinkage and selection via the lasso. J. R. Statist. Soc. B. 1996;58:267–288.
24. Donoho DL. Compressed sensing. IEEE Trans. Inf. Theory. 2006 Apr;52(4):1289–1306.
25. Tropp JA. Just relax: Convex programming methods for identifying sparse signals in noise. IEEE Trans. Inf. Theory. 2006 Mar;52(3):1030–1051.
26. Candes E, Romberg J. Sparsity and incoherence in compressive sampling. Inverse Problems. 2007 June;23(3):969–985.
27. DeVore RA, Jawerth B, Lucier BJ. Image compression through wavelet transform coding. IEEE Trans. Inf. Theory. 1992;38(2):719–746.
28. Donoho DL. For most large underdetermined systems of linear equations the minimal l(1)-norm solution is also the sparsest solution. Comm. Pure Applied Math. 2006 Jun;59(6):797–829.
29. Donoho DL. For most large underdetermined systems of equations, the minimal l(1)-norm near-solution approximates the sparsest near-solution. Comm. Pure Applied Math. 2006 Jul;59(7):907–934.
30. Davis G, Mallat S, Avellandeda M. Greedy adaptive approximation. J. Constr. Approx. 1997;13:57–98.
31. Efron B, Hastie T, Johnstone I, Tibshirani R. Least angle regression. Annals of Statistics. 2004;32(2):407–499.
32. Tropp JA. Greed is good: Algorithmic results for sparse approximation. IEEE Trans. Inf. Theory. 2004 Oct;50(10):2231–2242.
33. Chen SS, Donoho DL, Saunders MA. Atomic decomposition by basis pursuit. SIAM J. Sci. Comput. 1999;20:33–61.
34. Boyd S, Vanderberghe L. Convex Optimization. Cambridge University Press; Cambridge, UK: 2004.
35. Candes E, Romberg J. “L1-MAGIC: Recovery of sparse signals via convex programming. Applied and Computational Mathematics. California Institute of Technology; [Online]. Available: http://www.acm.caltech.edu/l1magic/
36. Baussard A, Miller EL, Lesselier D. Adaptive multiscale reconstruction of buried objects. Inverse Problems. 2004;20:S1–S15.
37. Shea JD, Kosmas P, Hagness SC, Veen BDV. Three-dimensional microwave imaging of realistic numerical breast phantoms via a multiple-frequency inverse scattering technique. Medical Physics. submitted. [PubMed]
38. Golub GH, Loan CFV. Matrix Computations. 3rd ed. John Hopkins University Press; Baltimore, MD: 1996.
39. Hansen PC, O'Leary DP. The use of the L-curve in the regularization of discrete ill-posed problems. SIAM J. Sci. Comput. 1993;14(6):1487–1503.
40. Winters DW, Shea JD, Kosmas P, Van Veen BD, Hagness SC. Three-dimensional microwave breast imaging: Dispersive dielectric properties estimation using patient-specific basis functions. IEEE Trans. Med. Imag. 2009 (in press) [PMC free article] [PubMed]
41. Sjöstrand K. Informatics and Mathematical Modelling. Technical University of Denmark, DTU; Matlab implementation of LASSO, LARS, the elastic net and SPCA. version 2.0. [Online]. Available: http://www2.imm.dtu.dk/pubdb/p.php?3897.
42. Mallat S. A Wavelet Tour of Signal Processing. 2nd ed. Academic Press; San Diego, CA: 1999.
43. Taflove A, Hagness SC. Computational Electrodynamics: The Finite-Difference Time-Domain Method. 3rd ed. Artech House; Norwood, MA: 2005.
44. Meaney PM, Fanning MW, Li D, Poplack SP, Paulsen KD. A clinical prototype for active microwave imaging of the breast. IEEE Trans. Microw. Theory Tech. 2000;48(11):1841–1853.
45. Bulyshev AE, Semenov SY, Souvorov AE, Svenson RH, Nazarov AG, Sizov YE, Tatsis GP. Computational modeling of three-dimensional microwave tomography of breast cancer. IEEE Trans. Biomed. Eng. 2001 Sep;48(9):1053–1056. [PubMed]
46. Zhang ZQ, Liu QH, Xiao CJ, Ward E, Ybarra G, Joines WT. Microwave breast imaging: 3-D forward scattering simulation. IEEE Trans. Biomed. Eng. 2003;50(10):1180–1189. [PubMed]
47. Meaney PM, Fanning MW, Raynolds T, Fox CJ, Fang QQ, Kogel CA, Poplack SP, Paulsen KD. Initial clinical experience with microwave breast imaging in women with normal mammography. Academic Radiology. 2007 Feb;14(2):207–218. [PMC free article] [PubMed]
48. Chaudhary SS, Mishra RK, Swarup A, Thomas JM. Dielectric properties of normal and malignant human breast tissues at radiowave and microwave frequencies. Indian J. Biochem. and Biophys. 1984;21:76–79. [PubMed]
49. Surowiec AJ, Stuchly SS, Barr JR, Swarup A. Dielectric properties of breast carcinoma and the surrounding tissues. IEEE Trans. Biomed. Eng. 1988;35:257–263. [PubMed]
50. Campbell AM, Land DV. Dielectric properties of female human breast tissue measured in vitro at 3.2 ghz. Phys. Med. Biol. 1992;37(1):193–210. [PubMed]
51. Joines WT, Dhenxing YZ, Jirtle RL. The measured electrical properties of normal and malignant human tissues from 50 to 900 MHz. Med. Phys. 1994;21:547–550. [PubMed]
52. Lazebnik M, McCartney L, Popovic D, Watkins CB, Lindstrom MJ, Harter J, Sewall S, Magliocco A, Booske JH, Okoniewski M, Hagness SC. A large-scale study of the ultrawideband microwave dielectric properties of normal breast tissue obtained from reduction surgeries. Phys. Med. Biol. 2007;52:2637–2656. [PubMed]
53. Kerlikowske K. The mammogram that cried Wolfe. N. Engl. J. Med. 2007 Jan;356(3):297–300. [PubMed]
54. American College of Radiology [Online]. Available: http://www.acr.org.
55. Davis SK. Ph.D. dissertation. University of Wisconsin - Madison; Madison, WI: Aug, 2006. Ultrawideband radar-based detection and classification of breast tumors.
56. Zastrow E, Davis SK, Lazebnik M, Kelcz F, Veen BDV, Hagness SC. Development of anatomically realistic numerical breast phantoms with accurate dielectric properties for modeling microwave interactions with the human breast. USNC/URSI Radio Science Meeting; July 2007. [PMC free article] [PubMed]
57. Converse M, Bond EJ, Van Veen BD, Hagness SC. A computational study of ultrawideband versus narrowband microwave hyperthermia for breast cancer treatment. IEEE Trans. Microw. Theory Tech. 2006;54(5):2169–2180.
58. Gabriel S, Lau RW, Gabriel C. The dielectric properties of biological tissues: III. Parametric models for the dielectric spectrum of tissues. Phys. Med. Biol. 1996;41:2271–2293. [PubMed]
59. Lazebnik M, Popovic D, McCartney L, Watkins CB, Lindstrom MJ, Harter J, Sewall S, Ogilvie T, Magliocco A, Breslin TM, Temple W, Mew D, Booske JH, Okoniewski M, Hagness SC. A large-scale study of the ultrawideband microwave dielectric properties of normal, benign, and malignant breast tissues obtained from cancer surgeries. Phys. Med. Biol. submitted. [PubMed]
60. Acceleware Corp. [Online]. Available: http://www.acceleware.com.