Abstract
Three contributions that can improve the performance of a Newton-type iterative quantitative microwave imaging algorithm in a biomedical context are proposed. (i) To speed up the iterative forward problem solution, we extrapolate the initial guess of the field from a few field solutions corresponding to previous source positions for the same complex permittivity (i.e., “marching on in source position”) as well as from a Born-type approximation that is computed from a field solution corresponding to one previous complex permittivity profile for the same source position. (ii) The regularized Gauss-Newton update system can be ill-conditioned; hence we propose to employ a two-level preconditioned iterative solution method. We apply the subspace preconditioned LSQR algorithm from Jacobsen et al. (2003) and we employ a 3D cosine basis. (iii) We propose a new constrained line search path in the Gauss-Newton optimization, which incorporates in a smooth manner lower and upper bounds on the object permittivity, such that these bounds never can be violated along the search path. Single-frequency reconstructions from bipolarized synthetic data are shown for various three-dimensional numerical biological phantoms, including a realistic breast phantom from the University of Wisconsin-Madison (UWCEM) online repository.
1. Introduction
During the last decade significant progress has been reported in fully 3D quantitative microwave imaging algorithms, for example, [1–11]. The aim of such algorithms is to reconstruct a 3D spatially varying complex permittivity profile from measured scattered field data. Consequently they consist of a numerical inversion scheme which employs electromagnetic field equations that exactly account for the vectorial nature and scattering of the electromagnetic waves. Research efforts have focused in particular on biomedical applications [3, 5, 11], since it has been pointed out that the electromagnetic properties at microwave frequencies are sensitive to tissue types and various physiological and pathological parameters [12–15]. The interest in these algorithms has been further stimulated by the possible benefits of microwave breast cancer screening as a supplemental diagnosis and/or therapy monitoring modality [4, 16–20]. Quantitative microwave imaging thus may open new perspectives as a moderate-cost, noninvasive, and nonionizing medical imaging technique.
It is well known that reconstructing the complex permittivity is a nonlinear and ill-posed inverse problem. The solution of this problem can be formulated as the minimizer of an appropriately regularized cost function that is composed of a data fit term and a regularization term or factor which incorporates a priori information about the object. Starting from an initial estimate, the permittivity profile is then updated progressively by means of an iterative optimization technique. Different approaches for choosing the cost function and the optimization technique have been proposed. Let us mention here the “modified gradient” type approach [3, 4, 6, 21, 22] and the “Newton-type” approach. In this paper we consider this last approach, where the cost function only depends on the complex permittivity and is optimized with a Newton-type technique. In this case a set of forward problems is solved in each iteration of the optimization. Whereas the “modified gradient” approach avoids these computationally demanding forward problem solutions, the “Newton-type” approach needs far less iterations to converge.
In previous work we have reported reconstructions from experimental data with the Gauss-Newton method applied to various types of cost functions for a variety of 3D piecewise-homogeneous objects [8, 23], including the 3D objects from the Fresnel database [24, 25]. In the present paper we deal with biological objects and we discuss three contributions that can improve the performance of a Newton-type algorithm with regard to the challenges posed by such objects. The reconstruction of a biological object is challenging when it is inhomogeneous with highly contrasted permittivities and when the number of reconstruction variables is large. Such conditions have an impact on the computational burden for the forward and update problems. In this paper we have chosen to implement the proposed contributions into the multiplicative smoothing (MS) regularized Gauss-Newton algorithm from De Zaeytijd et al. [8]. Let us mention that we also have implemented them into the stepwise relaxed value picking regularized Gauss-Newton algorithm for the reconstructions of the 3D Fresnel targets in [24], but their details were not reported there.
Our first contribution focuses on the forward problem solution and addresses the choice of the initial guess of the total field inside the reconstruction domain for the iterative forward solver. Given the fact that, in each iteration of the reconstruction algorithm, for all permittivity profiles in the line search the forward problem repeatedly needs to be solved for all illuminations, a reduction in the number of iterations of the forward solution leads to important overall computational savings. In Tijhuis et al. [26, 27] it was shown that a “marching on in anything” procedure, which extrapolates the initial guess from a few field solutions corresponding to previous values of a varying physical parameter, such as frequency and illumination angle, considerably enhances the speed of convergence. In [28] a “marching on in source position” extrapolation was applied in the first two iterations of a 2D distorted-wave iterative Born algorithm; the corresponding field solutions—for all sources and the two permittivity profiles—were stored and used to initialize a “marching on in profile” extrapolation for the remaining iterations of the reconstruction; a similar approach was employed in the line searches of a 2D quasi-Newton algorithm. In the 3D case a “marching on in profile” extrapolation is memory expensive; hence we applied the “marching on in source position” extrapolation for all iterations in [8]. In this paper we extrapolate the initial guess from a few solutions corresponding to previous source positions for the same permittivity, as we did in [8], but we also include a Born-type approximation that is computed from a field solution corresponding to one previous permittivity profile for the same source position. We employ this technique with the stabilized biconjugate gradient fast Fourier transform (BiCGSTAB-FFT) volume integral equation (VIE) forward solver reported in [29], but it can be applied with other iterative forward solvers as well, for example, [30, 31].
Our second contribution deals with the numerical solution of the update system in the Gauss-Newton optimization algorithm. In the early 2D microwave reconstruction algorithms, for example, the distorted Born iterative method (DBIM) [32] and the Newton-Kantorovich method [1], an update equation , with the Jacobian matrix, is obtained from linearizing the scattered field around the current permittivity iterate . The solution of this equation is formulated in the least squares sense, that is, the solution to the normal equations of the linearized (nonregularized) data error cost function—this corresponds to a Gauss-Newton update of the nonlinear data error cost function [33–35]. The resulting “least squares” update system with a dense (approximate) Hessian matrix is ill-conditioned and is regularized with a Tikhonov method—note that the regularization thus is applied to the permittivity update rather than to the permittivity. It is solved with a direct inversion method for the case of small objects in, for example, [1, 34]. More recently a Tikhonov regularized update system has been employed for breast imaging in 2D [36] and in 3D [17]; in the latter it is solved with a conjugate gradient technique. An equivalent approach consists in solving the nonregularized update system with the conjugate gradient for least squares (CGLS) algorithm, where regularization is achieved by terminating the CGLS algorithm after an adequate number of iterations. This approach has been tested in 2D with measured biomedical data [37, 38] and in 3D with synthetic data for realistic Magnetic Resonance Imaging (MRI-) derived numerical breast phantoms [18] and with clinical breast data [19].
In the present paper we start from a regularized data error cost function, for example, [8, 11, 23–25, 28, 39–42], in which the regularization is defined on the permittivity (see, e.g., [8, 43] for interpretations) and easily allows incorporating a priori information, for example, [42]. Application of the Gauss-Newton technique then yields an update system where the (approximate) Hessian matrix is composed of contributions from the data error and from the regularizing function. We observed that this update system still can be ill-conditioned, even in the presence of the positive definite matrix . In [8, 23] we solved the system with a BiCGSTAB routine. In [11] a CGLS iterative technique is combined with an implicit Jacobian calculation and a parallelized algorithm is tested with measured biomedical data. In this paper we employ a two-level preconditioned iterative method to obtain a computationally efficient and stable solution of the Gauss-Newton update. We apply the subspace preconditioned LSQR algorithm from Jacobsen et al. [44] with a 3D cosine basis. This algorithm extracts from the original update system a smaller system, which inherits its ill-conditioning but can be rapidly solved in a direct manner. The remainder of the system is better conditioned and can be solved iteratively with much less iterations than the original system.
Our third contribution is concerned with the use of permittivity constraints. Upper and lower bounds for the complex permittivity, which are based on a priori knowledge about the object, aim to improve the convergence of the algorithm. Often they are enforced, for example, by assigning after each update solution the constraint’s values to those reconstruction variables that violate the constraints [34]. In [45] such enforced update is reiterated a few times with the CGLS algorithm. In [40] constraints are incorporated in the Gauss-Newton framework by means of a nonlinear transformation that maps the constrained permittivity values on new, unconstrained optimization variables. In this paper we propose a technique that is inspired by the use of parameter transformations, but which is only applied in the line searches. We propose a new, constrained search path in the Gauss-Newton optimization, which incorporates in a smooth manner lower and upper bounds on the object permittivity, such that these bounds never can be violated along the search path. The technique can be implemented without modification of the original optimization scheme and line search algorithm.
The outline of the paper is as follows. In Section 2 the 3D regularized Gauss-Newton algorithm of De Zaeytijd et al. [8] is summarized, since the contributions of this paper are illustrated with this algorithm. In Section 3 the combined “marching on in source position-Born-type approximation” extrapolation procedure is detailed. The subspace preconditioned LSQR implementation is discussed in Section 4 and Section 5 treats the new constrained search path. In Section 6 single-frequency reconstructions from synthetic data are shown for various numerical biological phantoms, including a realistic 3D MRI-derived numerical breast phantom from the University of Wisconsin-Madison (UWCEM) online repository. Section 7 presents the conclusions.
2. Formulation
The object with unknown relative complex permittivityis immersed in a homogeneous background medium with relative complex permittivity , where is the relative permittivity, the conductivity, the permittivity of vacuum, the 3D position vector, and the angular frequency. In the following, the dependency of the time-harmonic fields is implicitly assumed and the -dependency of the complex permittivity is omitted, since a single frequency is adopted. In order to reconstruct the (complex) permittivity within a bounded investigation domain (see Figure 1(a)), scattering data are collected by successively illuminating the object with different known incident electric fields , which are radiated by elementary electrical dipoles positioned in points exterior to and with polarizations along unit vectors ; for each illumination the resulting scattered electric field components are measured in positions along unit vectors [8]. All these (complex scalar) scattered field components are stored in one -dimensional measured field vector , with . Note that (resp., ) refers to a unique combination of a position (resp., ) and a polarization (resp., an orientation ).
(a)
(b)
For the numerical reconstruction, we approximate the unknown relative complex permittivity with a piecewise constant function that assumes one value in each cell of a uniform cuboid grid , which comprises cubic cells with size (Figure 1(b)). These values are the reconstruction variables and they are collected in the -dimensional permittivity vector . The solution of the inverse scattering problem is then defined as the minimizer of a regularized cost function . In this paper we employ a MS regularized cost functionwhere is the (nonlinear) least squares data error, with a normalization constant and with the -dimensional simulated scattered field vector, which contains the scattered field component values as computed with the forward solver [29] for a given permittivity vector ; is a regularization parameter and is the smoothing regularizing function from [8]; see (A.1) in the Appendix.
For the minimization of we apply a Gauss-Newton technique with approximate line searches. In [8] the permittivity vector is updated from iterate to iterate aswhere is the (modified) Gauss-Newton search direction computed from update system (4) and is the descent step length, that is, a positive real number obtained by means of an inexact line search algorithm [46, pp. 34–38], such that is approximately minimized along the direction . In Section 5 of the present paper, we propose an alternative, constrained search path for (3). The (modified) Gauss-Newton search direction for cost function (2) is the solution of the systemwhere stands for conjugate transpose , where and where the subscript indicates quantities evaluated in . Furthermore, is the Jacobian matrix, with ; note that we employ the (independent) variables and (similar to [47]) and that we omit the subscript wherever this is more convenient; the vector and the (positive definite) matrix contain first-order derivatives () and second-order derivatives (), respectively, of the regularizing function ; see the Appendix for the explicit expressions. The expression between brackets in the RHS of (4) is the gradient with respect to the variables of cost function (2) and the factor in the LHS is obtained after neglecting in the approximate Hessian matrix (i.e., the submatrix that remains after neglecting the elements and ) a number of terms that are specific for the MS regularization (hence we call it a modified Gauss-Newton technique [8]). In [8] system (4) is solved iteratively with the BiCGSTAB routine [48].
Note that the search direction from (4) is also the least squares solution of the systemor the minimizer of the regularized linear least squares problemwhere and is the Cholesky factor of .
The computational effort in the reconstruction algorithm is mainly due to the multiple forward problem solutions and to the computation of the Gauss-Newton search directions. Indeed, for each iterate the forward problem is solved for illuminations, in order to compute the scattered field vector , which is needed to calculate the cost function —this is used in the stopping criterion and at the beginning of the line search—and to calculate the gradient in the RHS of (4), and in order to compute the total field on the grid, which is needed to calculate the Jacobian matrix (7) (this may require additional forward problem solutions); see further; the computation of the search direction involves the solution of a large and usually ill-conditioned linear system (4); furthermore an multi-illumination forward problem has to be solved for each trial value of in the line search algorithm. Note that the elements of the Jacobian matrix are given by [8]where is the background wave number, is a 3D unity pulse function with support cell , and is the total field on the computational grid resulting from dipole excitation in the transmitting position and oriented along (this quantity is available from the forward problem solution for ); is the total field on the computational grid resulting from a dipole excitation in the receiver position and oriented along ; hence additional forward problem solutions are needed for those receiver positions and orientations that do not coincide with transmitting positions and orientations. In the following sections we propose some ways to improve the computational efficiency.
3. Initial Guess for the Forward Solver
In this section we further reduce the computational effort for the forward problem solution by extending the “marching on in source position” extrapolation procedure [27, 28] with a Born-type approximation. We apply this technique to the BiCGSTAB-FFT iterative forward solver reported in [29, 49], but it can be implemented with other iterative forward solvers as well. Let us first remind the reader of some main principles of this solver. The 3D domain integral equation for the electrical field is expressed in terms of the unknown electric flux density and of a mixed potential formulation for the contrast current densities. This expression is discretized with a Galerkin Method of Moments, by employing a number of vectorial rooftop functions to test the equation and to expand the electric flux density. Note that the vector potential is not expanded directly as in [30] and that the -singularity of the Green function is handled by singularity subtraction. The rooftop functions are defined on a uniform cuboid grid , which is an integer subdivision of the inversion grid ; hence the forward problem cell size can be chosen as a fraction of the reconstruction cell size ; is usually around one-tenth of the wavelength in the considered medium. With using the normalized contrast function, defined as , the resulting linear system is written aswhere the -dimensional vectors , , and contain the tested incident field values, the unknown expansion coefficients of the electric flux density, and the products , respectively, on the grid for illumination and permittivity vector . The signs indicate one of both support cells of a rooftop function. The products and correspond to the tested total and scattered fields, respectively, on the grid. Note that the elements of the dense matrices do not depend on the permittivity; hence they need to be calculated only once for a series of scattering simulations with varying contrast. Equation (8) is solved for the coefficients with the BiCGSTAB method [48], starting from an initial guess . The FFT method [30, 50, 51] is used to speed up the matrix-vector multiplications in each iteration of BiCGSTAB while the multiplication is fast because is a sparse matrix. The resulting BiCGSTAB-FFT method has a memory requirement of and a computational effort of , with the number of iterations needed to reach the desired accuracy, expressed by a relative accuracy threshold .
In our inversion scheme the forward problem has to be solved several times: for a varying permittivity vector and with each of those for varying illuminations . In [27] it is shown that the number of iterations to solve (8) can be significantly reduced by means of a “marching on in anything” technique provided that is not much smaller than the relative error introduced by noise and the discretization. This technique proposes an adequate choice for the initial guess based on available solutions which correspond to slightly different illumination and/or object configurations. Suppose we have a few vectors , , that can be regarded as approximations for . The initial guess is then calculated as the linear combinationwhich minimizes the error between the LHS and RHS of (8). The coefficients thus are a solution of the linear systemwhich is a small system, since usually is sufficient.
In this paper we propose to use as the approximations to on the one hand a few solutions , which were computed for the same permittivity vector but for nearby transmitter positions—that is, “marching on in source position”—and on the other hand a Born-type approximation , which is calculated aswithwhere the vectors have elementsIn this expression, the vector contains the tested values of the total field , which is generated by the applied current density and the Born-type contrast current densitywhere is the new permittivity function (gathered in the new permittivity vector ) and is the total field corresponding to a previous permittivity function (gathered in the permittivity vector ) for the same illumination . The coefficients are available from this previous solution. Since the matrix is sparse and since the multiplications can be done with FFTs, the calculation of is fast. Note further that in (10) only the multiplication of with needs to be done, since the other products are the incident field vectors, which are available. The inclusion of the Born-type approximation in the “marching on” scheme thus allows for a simple extrapolation over the permittivity without having to store multiple solution vectors for a number of different permittivity profiles as would be the case in a “marching on in permittivity” scheme [28].
4. Subspace Preconditioned Gauss-Newton Update
The solution of system (4) for the Gauss-Newton search direction by means of a direct inversion method would require operations and operations would be needed to compute the product . Since the number of unknowns in a 3D inverse problem is large, it is more efficient to solve (4) iteratively. This requires per iteration a multiplication of an -dimensional vector with , followed by a multiplication of an -dimensional vector with . The computational complexity is much less than a direct inversion provided that the number of iterations can be kept small. However the condition number of —and therefore of —generally is very large. Even with the well-conditioned regularization term in (4) the conditioning of the system matrix remains problematic and worsens towards the end of the optimization in case of MS regularization, since is proportional to the least squares data error.
We therefore propose to solve (4) with the iterative subspace preconditioned LSQR algorithm (SPLSQR) of Jacobsen et al. [44], which is an implementation of a two-level iterative method [52] for the solution of large-scale regularized linear least squares problems such as (6). Since this algorithm is conceived for real system matrices, minimization problem (6) is reformulated aswithwhere and stand for the real and imaginary parts, respectively, and where is an matrix with , is a vector of length , and is a vector of length . As can be seen in Figure 2 for a generic (but relatively small) inverse problem, the singular value spectrum of the matrix with , that is, without regularization, gradually decreases over a large number of singular values without showing a clear threshold where this spectrum could be truncated. However, when , the regularization introduces a platform at the lower end of the spectrum of . Nonetheless the conditioning of still is not very good, because of the rapidly decreasing singular values in the first part of the spectrum.
The principle idea of the SPLSQR algorithm [44] is a splitting of the solution space into two orthogonally complementary subspaces and of dimensions and , respectively, with , spanned by the (orthogonal) columns of the matrices and , respectively. This means we look for a solution . It is furthermore desirable that has a small dimension and that the conditioning of is much better than the conditioning of . Indeed, introducing the QR factorizationwhere is an orthogonal matrix, is upper triangular, , and , it can be shown thatBecause has full rank, the same holds for and the first term in the RHS of (18) can be made zero for every . Therefore, (15) is equivalent toSince (20) is a small upper-triangular system, it is solved conveniently for by back substitution. The solution of (19) for is done with the LSQR algorithm [53], which although mathematically equivalent is numerically more stable than the CGLS method [44] and which can be formulated directly in terms of , thereby avoiding the construction of and multiplications with the large matrix . Calculating the QR factorization (17) and performing the multiplications with or are relatively cheap if done with Householder transformations, because is small.
With an appropriate choice of (and ) the number of LSQR iterations to solve (19) can be kept small. If, for instance, consists of the principal right singular vectors of , then the spectrum of , denoted as , coincides with up to , the th singular value of . Furthermore, since is orthogonal to , the subspace is spanned by the least significant right singular vectors of ; hence coincides with . It can be seen from Figure 2 that if is large enough, practically completely lies within the plateau in , introduced by the regularization. This means that the conditioning of in this situation is very good. It also implies that inherits the ill-conditioning of , but this does not pose a problem, since (20) is solved directly. However, it is computationally too expensive to construct this SVD subspace; hence a subspace that “resembles” it is chosen. It is well known that for most applications the right singular vectors of an ill-posed problem become more oscillatory as the corresponding singular values decrease. Therefore, we propose in this paper a truncated 3D discrete cosine basis (the DCT-2 basis from [54]), defined on the cubic grid , which is used for the real and imaginary part of the complex search direction vector separately; that is, the matrix now has the formwhere contains per column one of the DCT basis vectors with the lowest spatial frequencies in the -, -, and -directions. In this case the multiplication of with can be evaluated efficiently by performing 3D discrete cosine transforms [54] on the rows of (actually on the first and the second half of these rows separately) and retaining only the components with the lowest spatial frequencies in the -, -, and -directions, respectively. This subspace can be considered as a coarse grid approximation to the actual solution space. Figure 2 shows for this subspace and it can be seen that it indeed coincides with for large singular values.
5. A New Constrained Line Search Path
In order to improve the convergence of the optimization technique, but mainly to prevent the permittivity and conductivity values from becoming nonphysical or too high to handle with the chosen forward discretization cell size , including a priori knowledge concerning the expected upper and lower bounds on the complex permittivity in the biological object by means of constraints is recommended. Let us denote the upper and lower bounds by and , respectively, for the real part of the reconstruction variables and by and , respectively, for the imaginary part; that is,where and stand for the real and imaginary parts, respectively. Such constraints can be incorporated in the Gauss-Newton framework using a nonlinear transformation that maps the constrained permittivity values on new, unconstrained optimization variables [40]. Here, we follow an approach which is inspired by the use of parameter transformations but which is only applied in the line search. Furthermore our transformation differs from the ones reported in [40]. More specifically, we propose to replace the search path in line search (3) with a smooth, constrained path that entirely lies within the constraints, if the starting point does so, and which starts along a descent direction. With this constrained path, it is sure that the cost function will be reduced if is not a local minimizer and that the constraints will not be violated. This path is defined aswhere is the descent direction, solution of (4), and where we propose for the vector function the expressionPath (23) only considers the constraints that can be violated along path (3) with , hence the distinction between cases and in (24). For small (3) and (23) coincide, but in the vicinity of the constraints path (23) is bent away from (3) and it has a limit point on the constraints, as is illustrated in Figure 3(a) for a problem with only one complex optimization variable. Although theoretically the optimization variables can never reach their bounds with this procedure and the path always starts along a descent direction, it is possible that very little progress is made if one or more variables are very close to one or more of their bounds. Indeed, in such a situation the line search path deviates from its initial direction already for very small -values and starts to run along the projection of on the nearest boundaries of the constrained optimization domain as is illustrated in Figure 3(b). It is possible that although , calculated with (4), is a descent direction, its projection on these boundaries is not and the cost function starts to increase again for very small values of . In this case the line search is terminated after only a negligible reduction of the cost function. If is the steepest descent direction, then its projection on boundaries like the ones considered here (i.e., upper and lower bounds on the optimization variables), for example, , with a unit vector, remains a descent direction (or is zero), since the projection of on the gradient direction satisfies [46, 49]. Therefore, it is possible in our implementation to switch temporarily to the steepest descent direction when the abovementioned problems with the constraints are encountered.
(a)
(b)
Note that the use of the new search path implements the constrained optimization in a smooth way, without requiring any adjustments to the update system and the line search algorithm. Indeed, the same inexact line search algorithm as in [8] is used, comprising bracketing and sectioning phases and polynomial interpolation; see [46, pp. 34–38]. It also avoids sequences of time consuming minimizations of reduced problems, as is needed in active set methods. However, mapping (24) is highly nonlinear and a large step in may correspond to a negligible step in the actual optimization domain when many optimization variables are close to their bounds. The line search algorithm then may require many steps, hence many forward problem simulations, before a local minimum along the search path is reached. Although the Born-type extrapolation from Section 3 results in a rapid solution of additional forward problems when the permittivity profile changes little, a large number of forward problem solutions may increase the total computation time significantly. Therefore we recommend avoiding too stringent constraints. For example, instead of letting the maximum bound on the imaginary part be at most as required by the passivity condition, we allow for a (sufficiently small) positive value. Note finally that although our implementation uses constant values throughout the grid for the constraints , , , and , these constraints could vary locally in the grid.
6. Numerical Examples
In this section we discuss reconstructions from synthetic data of a number of simplified 3D biological models as well as of a realistic 3D MRI-derived breast phantom from the UWCEM online repository. All computations are performed on a 64-bit computer with 2 GHz Dual Core AMD Opteron processor and 8 GB RAM. FFTs are calculated using the FFTW library [55]. The LSQR solver is a self-written nonoptimized code.
In some of the examples we added white Gaussian noise to the simulated scattered fields. We define the signal-to-noise ratio (SNR) aswhere is the variance of the noise,with the exact permittivity profile. We also define the noise level as the least squares data error obtained for : .
In this paper we did not focus on optimizing the transmitter and receiver configurations for each example, by trying to use the smallest possible number of transmitters and receivers in the most appropriate locations for an adequately chosen number of permittivity unknowns. In [56] estimates are given for the number of degrees of freedom in each component of a single-view 3D scattered field; the corresponding information is accessible by uniformly positioning receivers over a measurement sphere around the scatterer. However it is stressed in [56] that their results do not apply to the case where transmitters and/or receivers are located in the scatterer’s reactive zone. In all our following examples transmitters and receivers are in the scatterer’s near field; they are positioned on aspect-limited surfaces and the numbers of receiver positions surpass estimates provided in [56]; two polarizations are employed for each transmitter and receiver. If the number of permittivity unknowns is larger than the information content in the data—which also can occur with large data sets—then it is the aim of the regularization to reduce the dimension of the solution space.
Let us mention that our fast forward solver, which employs a cuboid grid located completely inside the measurement surface, has a drawback with respect to the spherical and circular measurement surfaces and targets employed in the examples below, in the fact that the distance between these measurement surfaces and targets cannot be reduced beyond a limit determined by the “radius” of the cuboid; in our examples we even added some extra spacing; the circularly and elliptically cylindrical reconstruction domains employed with FEM and FDTD solvers in [17, 33] are more convenient in this respect.
6.1. Numerical Arm Phantoms
We consider two simplified arm phantoms: an adult’s arm and a child’s arm. The adult’s arm is from [3] and consists of muscle with permittivity and bone with permittivity and is immersed in water with permittivity . For the child’s arm we employ the same permittivities. The working frequency is 1 GHz, which yields a background wavelength of cm. The investigation domain is a cube with side 10 cm (). The initial guess for the reconstructions is the background medium.
6.1.1. Adult’s Arm Phantom
This phantom is an 8.7 cm wide finite cylinder of muscle containing a 3.3 cm wide finite cylinder of bone. As with the example in [3], which was inspired by the simulated antenna configurations in [57], we employ three parallel rings containing transmitter and receiver antennas and surrounding domain ; see Figure 4. We also use 18 transmitter and 90 receiver positions evenly distributed over the three rings, but in each transmitter position we apply two polarizations, along and the azimuthal direction , to an electrical elementary dipole and for each illumination, the scattered field is computed in all 90 receiver positions along these two directions. The data vector thus contains field components or complex numbers; note that due to reciprocity of these data values are redundant, which is the number of transmitter position pairs, , times the number of data values per transmitter position. The synthetic data are generated with a discretization of cells with size 3.33 mm () and Gaussian noise yielding a SNR of 30 dB is added (the noise level is ). The exact permittivity profile is depicted in Figures 5 and 6, where the permittivity cell is twice as large as that employed for the discretization of the fields.
(a)
(b)
(c)
(d)
(a)
(b)
(c)
(d)
For the inversion we use a grid with cells with size 6.67 mm (); hence the number of unknowns is . Constraints on the complex permittivity are necessary for the convergence of the reconstruction and are implemented in the line search as described in Section 5. The stopping criterion is set to , which accounts for the noise and the discretization error. From the reconstructions in Figures 5 and 6 we observe that the shape, dimensions, and position of the muscle and bone are correct. The permittivity values are rather close to the exact values, in particular for the real part in Figure 5. Some smoothing is visible due to the smoothing regularization. These results were obtained after 8 iterations (i.e., 8 solutions of update system (4)), including 21 multiview forward problem solutions, where the “marching on in source position” technique combined with a Born-type approximation from Section 3 was used for the initial guess of the BiCGSTAB-FFT iterative forward solver. The total execution time was 2 hours and 30 minutes.
6.1.2. Child’s Arm Phantom
This phantom is narrower than the previous phantom, it is positioned obliquely (Figure 7), and we employ larger sets of data and unknowns. It consists of a 6 cm diameter finite cylinder of muscle and a bone structure, which is modeled as a 2cm diameter finite cylinder and a 4 cm diameter sphere. We only consider the portion of the arm that falls within the investigation domain . The successive illuminations come from 120 electrical dipoles along two orientations, and the azimuthal direction , in 60 positions that are evenly distributed over 5 horizontal circles with radius 10.14 cm (=) and with vertical spacing 2.36 cm (=). For each illumination, the scattered field is calculated in all 60 positions along the same directions, which leads to complex numbers, nearly half of which are redundant due to reciprocity. The synthetic data are generated with a discretization of cells with size 3.33 mm (≈); no noise is added this time.
For the inversion we use a slightly coarser grid—in order to avoid committing an inverse crime—with cells with size 4 mm (≈), which yields complex permittivity unknowns. From —this inequality is even stronger if the redundant data due to reciprocity are removed—it follows that the rank of the Jacobian matrix is smaller than ; hence update system (4) would be singular if no regularization were applied. When the iterations proceed, parameter in (4) becomes smaller and the system evolves toward a singular system, which motivates the use of the subspace preconditioning of Section 4. Note that in practice a perfect data fit is not possible due to measurement and discretization errors or—in our simulation study—due to the misfit of the simulation grids for the data generation and the reconstruction; hence system (4) does not become singular. The SPLSQR algorithm is used with ; that is, only the discrete cosine basis vectors with the 8 lowest spatial frequencies in the -, -, and -directions are retained.
The result of Figures 8 and 9 is obtained in 7 iterations (i.e., 7 solutions of update system (4)) after a total execution time of about 5 hours and 45 minutes if the relative accuracy of the BiCGSTAB-FFT iterative forward solver is set to . The data error at this point is and is not significantly reduced by further proceeding with the iterations. The constraints that were imposed on the permittivity unknowns are and . The fact that we allow for slightly positive imaginary parts is motivated by the observation that too severe constraints can stall or even stop the convergence into a local minimum. The regularization parameter in (2) was set to —note that with MS regularization the choice of this parameter is not very critical [8]. From Figures 8 and 9 it can be concluded that a nice reconstruction is obtained, where the structures in the arm are clearly visible in the correct permittivity ranges. It deviates from the actual profile by its smoother appearance, which is due to the type of regularization we employed.
(a)
(b)
(c)
(d)
(a)
(b)
(c)
(d)
The “marching on” technique of Section 3 has been used with ; that is, the initial guess for a forward problem solution is obtained as a linear combination of Born-type approximation (11) and the solutions for three previous transmitter positions on the same horizontal circle in Figure 7. There is no extrapolation over different dipole circles. Without the “marching on” scheme the total computation time for the present example increases to about 8 hours, an increase of about 40%.
6.2. Simple Breast Phantoms
We consider two numerical malignant breast phantoms: one without and one with a skin layer. The breast is a hemisphere with a radius of 5 cm that is filled with a homogeneous, averaged breast tissue with permittivity . This value is obtained with a single-pole Debye dispersion model by adopting the Debye parameters from a numerical breast phantom in [58]. For the phantom with skin, the breast is covered with a 2.5 mm skin layer with permittivity (dry skin). We chose a lossless background medium with permittivity , which matches the real part of the fatty breast tissue . Some guidelines for selecting a matching medium are given in [22]. The operating frequency is 1 GHz, yielding a background wavelength cm. The tumor is a 2 cm diameter sphere with permittivity [58] at the position (−2 cm, 2 cm, 0 cm).
The antenna configuration is depicted in Figure 10—note that an experimental hemispherical radar antenna array is reported in [59]. There are 72 antenna positions on a hemisphere with radius 9 cm surrounding the front side of the breast, with 12 positions evenly spaced on each of the 6 meridians. Only the 48 central positions indicated in bold are used for the transmitting dipole, which is subsequently oriented along a meridional and an azimuthal polarization. The investigation domain , with dimensions , is also indicated. The scattered field is computed in all 72 antenna positions for both orientations and Gaussian noise yielding a SNR of 30 dB is added (the noise level is ). The resulting dimension of the data vector is ; note that due to reciprocity of these data values are redundant, which is the number of transmitter position pairs, , times the number of data values per transmitter position. The cell size in used for the data generation is 5 mm (=, ≈) for the phantom without skin (see Figure 11 (top)) and 2.5 mm (=, ≈), that is, the thickness of the skin, for the phantom with skin, in order to facilitate modeling of the skin layer (see Figure 12 (top)). We also employ these respective cell sizes for the forward problem solutions during the reconstructions.
(a)
(b)
(a)
(b)
(c)
(d)
(a)
(b)
(c)
(d)
6.2.1. Simple Breast Phantom without Skin
For the reconstruction of this phantom the cell size of the grid is 5 mm; hence the number of complex permittivity unknowns is . The relative accuracy threshold for the BiCGSTAB-FFT iterative forward solver is set to , since there is no point in much lowering it given the noisy data. The initial guess is the background medium everywhere in . The regularization parameter is and the permittivity constraints are and . The preconditioned LSQR algorithm with is used to solve for the Gauss-Newton update to an accuracy of .
After 6 iterations the least squares data error reaches the noise level, that is, , and the optimization is stopped. With further iterations, the reconstruction error—that is, the error between the reconstructed and exact permittivity profiles—increases again, while the cost function is only marginally reduced. The resulting images are shown in Figure 11 (bottom), where the tumor is clearly visible at the right location and with a higher relative permittivity and conductivity than the surrounding breast tissue. Considering the limited aperture data and the small dimensions of the tumor with respect to the background medium wavelength (<), this result is quite satisfactory. The total execution time for this reconstruction was 31 minutes. When the acceleration techniques proposed in this paper were not applied, that is, without the “marching on” technique and by using a nonpreconditioned BiCGSTAB solver for the Gauss-Newton update, the optimization still was not finished after 160 minutes and in that time only 3 iterations were completed. The main cause of this is the large number of iterations needed by the nonpreconditioned BiCGSTAB solver to converge to the desired accuracy, especially after the third iteration when the data error is already close to the noise level.
6.2.2. Simple Breast Phantom with Skin
The breast phantom with skin is more realistic but it renders the inverse problem more challenging: there is a nonnegligible scattering from the skin interface, which can obscure the signal from the tumor, and the small thickness of the skin layer results in a finer discretization . For this phantom we demonstrate possibilities and limitations of the reconstruction algorithm with regard to a partial reconstruction based on a priori information. We assume the skin layer to be known and we perform the optimization only for the permittivity cells inside the inner skin contour. Furthermore, we reduce the number of permittivity unknowns by optimizing for aggregates of permittivity cells that are cubes with side 5 mm, that is, 8 cells of the forward grid , or portions of cubes next to the skin. This yields complex permittivity unknowns. The permittivity constraints are set to and . The initial guess for this reconstruction is for all cells inside the breast, which is in between the breast permittivity and the tumor permittivity .
Note that the preconditioned LSQR algorithm with the 3D discrete cosine basis defined on the cuboid grid cannot be used for this partial optimization example. It may be possible to develop another subspace of breast specific basis functions for this problem. The nonpreconditioned BiCGSTAB solver thus is used instead, but with the relatively small number of reconstruction variables the iteration count remains acceptable. The accuracy for this solver is set to and the regularization parameter is . The resulting images after 5 iterations and about 2.5 hours are shown in Figure 12 (bottom). The tumor is again clearly visible at the right location and with a higher relative permittivity and conductivity than the surrounding breast tissue, but these values again are lower than the exact values. Some higher values also appear in a number of spots near the (high permittivity) skin interface, which might be a compensating behavior for reconstruction errors elsewhere in the image. The longer computation time results from the finer discretization in the forward problem.
6.2.3. Reconstruction of Effective Permittivities
Inspired by a homogenization strategy reported in [58] for 2D time domain ultrawideband detection of breast tumors, we use our inversion algorithm to estimate a homogenized effective breast permittivity from the scattering data of the malignant breast phantom without skin. This value could be used, for example, as an approximation for the healthy tissue permittivity in other—qualitative—inversion methods, such as the linear sampling algorithm [22, 60, 61], where anomalies are detected against a known background by solving only a linear problem. It also could be used as a guideline for adapting the matching medium or as an initial estimate for a more detailed quantitative reconstruction. We use only one (aggregated) reconstruction variable for this reconstruction—we keep the same discretization as before for the forward problem solutions—and its support coincides with the breast phantom. The reconstruction is obtained after 2 iterations in 5 minutes and is close to the healthy tissue permittivity .
For the malignant breast phantom with skin we estimate two effective permittivities, one for the skin layer and one for the breast tissues, in a similar way as presented for the first phantom. Starting from the background permittivity , the result after 3 iterations and 52 minutes is and , values that are within 6% of the actual complex permittivities and in the exact phantom.
6.3. Realistic 3D MRI-Derived Numerical Breast Phantom
The UWCEM Numerical Breast Phantom Repository contains a number of anatomically realistic MRI-derived numerical breast phantoms for breast cancer detection and treatment applications [62]. These phantoms capture the structural heterogeneity of normal breast tissue and incorporate the realistic dispersive dielectric properties of normal breast tissue from 0.5 to 20 GHz reported by [14, 15]. We selected for this paper Phantom 1 from ACR class 1 (ID 071904), which is a mostly fatty breast phantom with some glandular and fibroconnective inhomogeneities. The complex permittivity in a sagittal slice through this phantom at a frequency of 2 GHz is depicted in Figure 13 [63]. For the background medium we chose a material with permittivity (the corresponding conductivity is S/m), which differs from the background medium employed for the simple breast phantoms by the addition of some losses. The background wavelength is cm.
(a)
(b)
Since the cell size of this MRI-derived phantom is very small, only 0.5 mm, and since such a high resolution is not needed at the considered frequency nor desirable due to the high memory and computational costs, we derived a coarser permittivity model with cell size 2.5 mm (=) from this phantom by local averaging [61]. This coarser model is depicted in Figure 14 in the same sagittal slice. An artificial 2 cm diameter spherical tumor with permittivity is added rather close to the chest wall to make its detection more challenging. We also removed the muscle layer that was added at the base of the original phantom [63], since in our free-space measurement setup this thin high-contrast layer would cause too much scattering at its (nonrealistic) interfaces with the background medium.
(a)
(b)
The coarser permittivity model from Figure 14 is used to generate the data. The dipole configuration is depicted in Figure 15. It consists of 168 dipoles in 84 positions on an ellipsoidal surface around the front side of the breast with polarizations in two orthogonal directions tangential to this surface. All these dipoles are used to sample the field, but only 48 of them in 24 positions (indicated with the larger black dots) are used to illuminate the phantom because of memory limitations (in particular with respect to the size of the Jacobian matrix). This yields a total of complex numbers; note that due to reciprocity of these data values are redundant, which is the number of transmitter position pairs, , times the number of data values per transmitter position.
The permittivity grid and the forward field grid for the reconstruction contain cells with size 5 mm (≈), which yields complex permittivity unknowns. As with the child’s arm phantom the rank of the Jacobian matrix is smaller than ; hence regularization is indispensable. For this example the subspace dimension is , where the discrete cosine base with, respectively, the 8, 10, and 7 lowest spatial frequencies in the -, -, and -directions is employed, that is, a coarse grid approximation with roughly one-third of the resolution of the full permittivity grid . The regularization parameter is . To test the abilities of the reconstruction method we perform a blind reconstruction; that is, we do not perform reconstruction solely within a known breast contour as in [17, 18], and the initial estimate is the background medium everywhere in . The constraints on the permittivity are and . The “marching on” scheme is used with ; that is, the initial guess for the BiCGSTAB-FFT iterative forward solver is obtained as a linear combination of Born-type approximation (11) and the solutions for three previous transmitter positions. Figure 16 shows the result after 13 iterations (or 18 hours and 25 minutes); the changes in the permittivity profile then are less than 1 percent. The shape of the breast and the overall structure of the inhomogeneities are clearly visible in the reconstruction. The tumor is located correctly as well. Even the small lump of fibroglandular tissue near the nipple can be resolved. We observe some smoothing that might be due to the type of regularization. However, the reconstructed relative permittivity and conductivity values of the inhomogeneities are lower than the exact values: in the center of the tumor the reconstructed permittivity is instead of . In [17, 18] quantitative reconstructions of numerical breast phantoms from the UWCEM repository are presented for a smaller 1 cm diameter inclusion. The reconstructions are from multiple frequency data values and employ a finer discretization grid than we did. From a visual inspection of Figure 8 in [17] at 1.5 GHz, the inclusion hardly can be distinguished in the permittivity curves, but it clearly appears in the difference image of reconstructions from data with and without inclusion (Figure 9 in [17]), the quantitative values still being too low though. Similar observations hold for [18] at 2.5 GHz, where a contrast agent was applied to the tumor. In [19] a 1.2 cm size tumor is reconstructed in a mostly fatty breast from single-frequency clinical data at 1.3 GHz and the identification of this tumor is facilitated by comparison with a reconstruction of the other—healthy—breast where no tumor is present.
(a)
(b)
Finally we employ this example to illustrate the effect of the subspace preconditioning. We compare the efficiency of the SPLSQR algorithm to the efficiency of a conventional iterative solver without preconditioning in solving (4). The conventional solver is the BiCGSTAB routine [48]. For the first () and the final () Gauss-Newton iterations we let the SPLSQR algorithm with the same parameters as mentioned earlier solve (19) to a relative accuracy of and we calculate the resulting accuracy on original system (4). We then let BiCGSTAB solve (4) to that accuracy and compare the number of iterations and the time needed by both methods. The comparison is summarized in Table 1. In the beginning of the optimization the SPLSQR algorithm is more than 4 times faster than BiCGSTAB and towards the end of the optimization, when the preconditioning becomes more crucial due to the decrease of parameter in (4), the speedup factor approaches 20. The reduction in solution time is less than the reduction in the number of iterations in the case , since the computation of QR factorization (17) is also included in the former. Note that the time needed to solve (4) a single time without preconditioning at is 90% of the total reconstruction time with preconditioning. Note furthermore that we use the BiCGSTAB routine from the PIM library [48], which is an optimized Fortran library, while we implemented the SPLSQR ourselves in C; thus the speedup might even be more significant with an optimized implementation.
7. Conclusions
We have employed a three-dimensional quantitative reconstruction algorithm in the context of biomedical microwave imaging. The algorithm is based on a Gauss-Newton optimization scheme with line searches that is applied to a regularized nonlinear cost function. The forward solver has been improved by the combination of a “marching on in source position” and a “Born-type approximation” extrapolation procedure for choosing the initial guess of the total field on the computational grid. A two-level subspace preconditioned iterative method has been employed for the solution of the complex permittivity update system. In particular, we have applied the subspace preconditioned LSQR algorithm from Jacobsen et al. (2003) with a 3D cosine basis. A significant improvement in the computational efficiency with respect to using the BiCGSTAB solver is observed here, but it would be worthwhile to compare its performance also with a CGLS iterative scheme. A modification of the line search path based on a parameter transformation has been proposed to allow for a smooth incorporation of constraints in the optimization method. Numerical experiments have been conducted at one single frequency, thereby avoiding difficulties with the dispersive nature of body tissues, and with bipolarized antenna configurations. They have shown that the method yields promising results for biomedical imaging, in particular for the challenging problem of breast tumor detection. Multifrequency schemes, dedicated antenna positioning configurations, and the incorporation of more a priori information in the reconstruction strategy may further improve the resolution of the images and the estimation of the permittivity. Possible benefits of using bipolarized versus monopolarized transmitters and/or receivers could be explored further, although a bipolarized configuration may substantially affect the complexity of a measurement apparatus. Future developments also should include the use of a problem-specific regularization strategy to cope with the smearing effect shown by the smoothing regularization.
Appendix
In this paper we use the smoothing regularizing function from [8] (we set the normalization constant in equation of [8]):where the cell index is replaced with the triplet in the -, -, and -directions. Whenever the triplet in (A.1) indicates a cell outside (i.e., , , , , , and ; see Figure 1(b)), the corresponding value of is equal to the background permittivity . The gradient and Hessian matrix of are given bywith andwhere is a real and constant matrix with . The explicit expressions for the elements of and are
in these expressions represents the set of neighboring cells of cell , that is, the six cells that share a face with cell . These also include the virtual neighboring cells just outside when is on the border of .
Conflict of Interests
The authors declare that there is no conflict of interests regarding the publication of this paper.
Acknowledgments
Jürgen De Zaeytijd was a Research Assistant of the Fund for Scientific Research-Flanders (F.W.O.-Vlaanderen). The authors would like to thank the Computational Electromagnetics Laboratory of the University of Wisconsin-Madison for making their numerical breast phantoms available.