Superlens Enhanced 2-D Microwave Tomography With Contrast Source Inversion Method

The contrast source inversion (CSI) algorithm is one of the primary techniques used for solution of non-linear inverse problems in microwave tomography. In this paper, we describe a modification of the CSI method adapted to imaging of the 2-D objects in the presence of the focusing media under TM-polarization. The focusing media is presented in the form of the Veselago lens. The data domain and imaging domain are properly positioned with respect to the location of the lens. Specifically, the sensors are located at the focal points of the lens with respect to the location of the individual pixels discretizing the contrast source. Such positioning of the source and observation locations in the presence of the lens, eliminates rank deficiency in the formulation of the inverse problem and results in significant improvements to both convergence speed of underlying conjugate gradient iterations and the accuracy of the image reconstruction in the CSI method.


I. INTRODUCTION
M ICROWAVE Tomography (MWT) has various important applications in medical imaging [1], nondestructive testing [2], security screening [3], structural health monitoring [4], remote sensing [5], and other areas [6]. While qualitative imaging techniques [7], such as Ground Penetrating Radar [8] and Synthetic Aperture Radar [9] methods, have been well developed and can often provide useful information about the location and shape of the imaged objects, it is the more accurate quantitative methods which have a potential to greatly expand the areas of MWT applicability. The quantitative methods of MWT solve the non-linear inverse problem of reconstructing electric permittivity and magnetic permeability of the objects from the information on how the object scatters the electromagnetic fields. The two primary methods of quantitative MWT are the CSI method [10] and the Distorted Born Iterative Method (DBIM) [11] (a.k.a. Gauss-Newton Inversion (GNI) method [19]). Both methods are iterative in nature and provide about equal abilities in quality of reconstruction and the overall computational cost [12], [13].
The CSI method, however, is much easier to implement as it does not require inversion of either the forward or inverse scattering operators at each iteration of the underlying non-linear optimization scheme (Conjugate-Gradient (CG) method in CSI and Gauss-Newton method in DBIM). The CSI method also offers a great flexibility in accounting for the complex media in which the imaging experiments are conducted. This method has been used for imaging of the objects in both free-space environment [14] and metal enclosures [15]. The latter has been shown to improve the quality of the CSI images through the modification of the nonsymmetric enclosure positioning with respect to the imaged objects and hence producing a broader variety in the scattered field information.
The discretized inverse scattering operators, i.e., the operators mapping the scattered field to the contrast sources producing them, most commonly have incomplete rank [16]. Inversion of such operators in MWT presents an ill-posed problem [17] which does not have a unique solution unless the problem is regularized [18]. Typically, the regularization is done mathematically through augmentation of the inverse scattering operators with additional constraints of a non-Maxwellian nature (e.g., the boundedness of the solution in Tikhonov regularization [19]). In [20], [21], it was shown, however, that the rank-deficiency of the discretized inverse scattering operator and thus the ill-posedness of the inverse scattering problem can be eliminated in rigorous compliance with Maxwell's Equations. For that purpose, the kernel of the forward scattering operator must be the Green's function of a focusing media and the sensor positions must be properly situated with respect to the locations of the sources in the imaging domain. Under these conditions, the inverse scattering problem can be solved in the same rigorous manner as the forward scattering problem.
In this paper, we study how the elimination of rank deficiency in the discretized inverse scattering operator established through staging of the imaging experiment in appropriate focusing media impacts the performance of the CSI method. As the focusing media, we use the idealized model of the Veselago lens [22], [23]. The lens is formed by a dielectric slab of double negative parameters, finite thickness, and infinite extent. Such lens, upon both the values of relative permittivity and permeability approaching negative unity, has the ability to reconstruct at its focal point both evanescent and propagating waves contained in the spectrum of the point source [24]. Such superlenses [25] have been experimentally demonstrated both in microwave [26] and optical regimes [27]. As a result, the Green's function of the lens has singular behavior at the location of the source point and near singular behavior at the corresponding focal point. We position the pixelated imaging domain (imaging domain discretized with piece-wise basis functions) on one side of the lens, while sampling the scattered field on the other side of the lens at the locations symmetrically corresponding to the centers of the object pixels. Under these conditions and sufficient focusing established by the lens, the discretized inverse scattering operator can be shown to attain its full rank. We subsequently perform the full non-linear inversion using CSI algorithm that is appropriately modified to account for the presence of the lens. Since the functional which is minimized in the CSI method is constructed via Fréchet derivatives [28] taken with respect to the contrast source function, the rigorous modification of the CSI method accounting for the presence of focusing media is effected through the use of the lens Green's function in the forward and inverse scattering operators instead of the homogeneous background Green's function featured in the conventional CSI algorithm.
The numerical results show that performing non-linear inversion with the CSI method in the focusing media instead of the conventional non-focusing media [12], [14], [29], substantially improves the accuracy of the image reconstruction. The new CSI algorithm is shown to be more robust also. It exhibits consistent convergence to an accurate image of the object within a single frequency experiment and regardless of the object shape. To validate the CSI algorithm and study it's properties the 'E'-shaped and 2-layer square inhomogeneous objects are considered.

A. FORMULATION
Consider TM z polarized time-harmonic scattering (in e iωt time convention) in 2-D from a non-magnetic dielectric cylinder occupying cross-sectional area D, oriented along z-axis, and having arbitrary cross-section and spatial distribution of complex relative dielectric permittivity (ρ ). We will define the contrast χ(ρ ) and contrast source w(ρ ) functions, as follows: where ρ is the position vector, = ω (ρ )ε 0 μ 0 are the wavenumber of free-space, background media and dielectric cylinder, respectively. Note, that χ(ρ ) = 0 only when ρ ∈ D.
We can formulate the forward and inverse scattering problems using the following operator notation: where A is the domain where the position vector ρ resides, ρ and ρ are the observation and source points, respectively, and G(ρ, ρ ) is 2-D Green's function satisfying Thus, by virtue of (2) we can define operator G S : D → S producing the scattered electric field at sensor locations S due to contrast sources w inside the imaging domain D enclosing cylinder area D and operator G D : D → D producing the scattered electric field at the imaging domain D due to contrast sources w in it.
If the dielectric permittivity distribution (ρ ) is known, then the forward scattering problem can be formulated, as follows: where χ(ρ ) is a known contrast function and E(ρ ) is an unknown field due to this contrast in imaging domain ρ ∈ D.
Similarly, inverse scattering problem is given by where both contrast function χ(ρ ) and field inside the scatterer E(ρ ) are unknown, while E sc (ρ), ρ ∈ S is known. It is common to refer to the problem of finding the contrast function χ(ρ ) as inverse scattering problem and the problem of finding the contrast sources w(ρ ) as inverse source problem [19]. A typical set up for a solution of the inverse problem is depicted in Fig. 1. The scattered field data E sc is collected at the sensors residing in the domain S for J excitations that illuminate the object by different source configuration. In general, the imaging domain D is different from the object domain D since the exact size and shape of the imaging object are unknown and are subject of reconstruction.
Forward scattering problem (4) is in a form of a secondkind Fredholm [30] integral equation and, therefore, it is well-posed. However, if all sensors are located outside the imaging domain D S = ∅ (non-invasive imaging procedure), the inverse source problem (5) is in the form of the first-kind Fredholm integral equation with smooth kernel which makes it ill-posed.

B. GREEN'S FUNCTION OF VESELAGO LENS
Veselago lens has a unique property of reassembling at the focal point not only the propagating waves as the conventional lenses do [30] but also the evanescent waves [22], [23]. The evanescent waves while traveling through the lens acquire the same gain as their attenuation in the process of propagation outside the lens, hence, arriving to the focal point with their initial magnitude with which they depart from the source point. Depiction of the field produced by a point source in the presence of the Veselago lens computed with OpenFMM software [31] is shown in Fig. 2. Besides the singular field behavior in the vicinity of the source point location, one can see two focal points corresponding to the points of ray intersection (Fig. 3). One focal point is located half way through the lens' thickness and the other one is on the opposite side of the lens symmetrically situated with respect to the center of the lens and the location of the source point.
The operator definition (2) can be rewritten for the imaging experiment scenario in the presence of the Veselago lens (or any other focusing media), as follows: where G foc (ρ, ρ ) is the Green's function of the focusing media [22], [32], and A is the domain where the position vector to the observation point ρ resides. Thus, we can define an operator F S : D → S producing the scattered electric field at sensor locations S due to contrast sources w inside the imaging domain D inclosing cylinder area D and an operator F D : D → D producing the scattered electric field in the imaging domain D due to contrast sources w. Sensors have to be located inside the focal point of the lens to allow for the efficient inversion. In this work, the Green's function of the Veselago lens was computed directly via error-controlled evaluation of Sommerfeld integrals [32].

C. DISCRETIZATION
The Method of Moments (MoM) discretization of the operators G S and G D given in (2) or their focusing media counterparts F S and F D given in (6) requires the discretization of the imaging domain D and sensor domain S. The domain S is naturally discrete as the scattered field is observed at a finite number of sensors K. The position of the kth sensor is defined by the radius vector κ k , k = 1, . . . , K.
The imaging domain D is discretized with M square cells and the position of the mth cell centroid is defined by the radius-vector ρ m . The unknown contrast sources w(ρ ) are discretized by using piece-wise basis functions [33] Thus, substitution of (7) into (2), followed by pointsampling of the scattered field in the ranges of the operators G D and G S casts (2) in the matrix form where w is an M × 1 vector of unknown contrast source expansion coefficients in (7), E sc,D is an M × 1 vector of values of the scattered field at M centroids of the cells discretizing the imaging domain D, E sc,S is an K × 1 vector of values of the scattered field at K sensor locations, and matrix elements of G D and G S are the inner products [33] where m, m = 1, . . . , M, and k = 1, . . . , K.
Similarly, point-sampling of the scattered field in the ranges of the focused media operators F D and F S results in their discretized forms F D and F S , and the matrix entries are defined as inner products where m, m = 1, . . . , M, and k = 1, . . . , K. A general example of the operators G S and F S in the inverse problem setup is depicted in Fig. 3.
The evaluation of the integrals in (10)-(13) can be done using standard Gauss-Legendre quadrature rules [34] after the singularity is extracted and evaluated analytically using the technique in [35].

D. G S AND F S OPERATORS
The nearly singular field at the focal point is the key to elimination of the rank deficiency in the operator F S producing the scattered field at sensor locations, provided each kth sensor position (defined by κ k ) is situated at the focal point of the mth pixel center ρ m with respect to the lens' center, as shown in Fig. 3. Under such particular selection of the sensor locations and the corresponding discretization of the imaging domain, the matrix F S of the discretized operator F S becomes diagonally dominant, as the contribution of any pixel to its corresponding sensor is substantially larger than to the other sensors due to the focusing properties of the lens (Fig. 2, Fig. 3). In other words, the lens performs the spatial separation of the scattered fields emanating from the individual pixels of the discretized contrast source w and routes these fields predominantly to their corresponding focal points.
The full-rank matrix F S indicates well-posedness of the underlying inverse scattering problem formulation. In [20], it was shown that such well-posed formulation of the inverse problem allows for its direct solution in analogous manner to the solution of the forward scattering problem. For comparison, in Fig. 4 we present the singular values plot of the matrix G S resulting from the discretization of the operator G S in the case when the imaging experiment is staged in the homogeneous background media and the matrix F S resulting from the discretization of the operator F S in the case when the imaging experiment is conducted in the presence of the focusing media (Fig. 3). The matrix G S is seen to be rank deficient and hence not being able to produce a unique solution. It is a manifestation of the ill-posedness in the pertinent formulation of the inverse scattering problem in the non-focused imaging setup.

E. CSI ALGORITHM IN FOCUSING MEDIA
The inverse source problem (5) can be solved iteratively as a non-linear optimization problem using the CSI algorithm [10], [14]. The cost-functional F combines two types of errors: mismatch in the scattered field at sensor locations ρ for the current guess of the contrast sources w given by and mismatch in the total field inside the imaging domain D for the current guess of the contrast sources w and contrast function χ given by where j is the excitation number, j = 1, . . . , J. The mismatches ρ j and r j in (14) and (15) are often called errors in the data equation and domain equation [15] and have the meaning of how well the current guess of the contrast function χ and contrast sources w satisfy (4) and (5), respectively. By introducing the following normalizations,

Algorithm 1 Contrast Source Inversion (CSI) Method [10] in Focusing Media
Inputs: True (measured) field at sensor locations E sc , E in , incident field inside the imaging domain E in,D , discretized operators F S , and F D , required relative tolerance tol Output: reconstructed contrast function χ , normalized mismatch at sensors δ S and in the domain equation δ D 1: Calculate initial guess w (0) and χ (0) 2: Calculate η S and η D(0) according to (16) (14) 13:

15:
if (δ S < tol) AND (δ D < tol) then 16: return χ (n) , δ S , δ D 17: else 18: (15) 19: Update η D(n) according to (16) 20: end if 21: end for the cost-functional F [14] can be written in a compact form, as follows: Minimization of the cost-functional F in (17) via CG method requires the computation of Fréchet derivatives [28] and leads to the following expression for the gradient [14] where¯ χ denotes the complex-conjugate, F S * and F D * are the adjoint operators to F S and F D defined in (6). In (18),

FIGURE 5. (a) Experiment setup for the CSI inversion in the presence of the Veselago lens. The observation points (Rx) are positioned in the lens focal points of the pixels discretizing the imaging domain D for the 'E'-shaped object. The interrogating field is produced by point sources (Tx) uniformly distributed around the imaging domain. (b) Depiction of the 2-layer square object used in the validation of the modified CSI algorithm.
F S is the component of the functional responsible for the mismatch in the scattered field at sensor locations ρ (14), and F D component is the mismatch in the domain equation (4) for the current guess of the contrast χ . The CSI algorithm in the homogeneous background media with wavenumber k b can be formulated by changing operators F S , F D , F S * , and F D * with scattering operators G S , G D , G S * and G D * defined in (2). The detailed derivation of the CSI algorithm and update parameters for CG method are presented in [10]. In this paper, we give the general pseudo-code for the CSI (Algorithm 1) using Polak-Ribiere search directions [36] for J = 1 excitation. Our algorithm formulation is slightly different from [14], as operators are generalized to allow for the inversion in the focusing media, and the stopping criteria is based not only on the mismatch at the sensor locations ρ (14), but also on the mismatch in the domain equation r (15). In Algorithm 1, * is used to denote the conjugate transpose and is the entrywise Hadamard product [37] and is the entrywise division.

F. INITIAL GUESS
The presence of the focusing media and its proper positioning with respect to the imaging domain D and sensor locations S has been observed to have a profound positive impact not only on the spectral properties of the discretized operators, but also on the quality of the initial guess in the solution.
The CSI algorithm is initiated by estimation of the contrast source value function w (0) using the back-propagation relationship between the true scattered field measured at the sensor locations E sc and the contrast sources [19] followed by the estimation of the total field E D(0) in the imaging domain using discretized operator F D acting on the estimated contrast source w (0) Fig. 10 and Fig. 11 show the initial guess for the contrast χ 0 obtained in the homogeneous background media ( b = 4ε 0 ) and in the presence of the lens for the 'E'shaped inhomogeneous object. It's not difficult to see that a substantially higher accuracy estimate of the imaged object is obtained when the initial guess search is aided by the focusing environment of the lens.

III. NUMERICAL RESULTS
To validate a modification of the CSI algorithm accommodating for the presence of the focusing media, we perform numerical reconstruction of the relative electric permittivity distribution (1) for two different objects. The first one is a 2-layer square object consisting of the denser core scatterer of a square cross-section centered at the origin and having side length of 0.9λ 0 and the layer of a square cross-section tightly surrounding the core scatterer, centered at the origin and having the length of its side of 0.5λ 0 (Fig. 5). The relative permittivity of the core scatterer is 1 , while the relative permittivity of the surrounding square shell is 2 . The second imaged object is a lossless 'E'-shaped object' composed of two different materials with relative permittivity ε 1 = 4.2 and ε 2 = 4.5 depicted in Fig. 5. The loss (imaginary part of 1 and 2 ) for the 2-layer square object is specified for individual simulations in the caption for Fig. 10, while for lossless simulations, 1 = ε 1 = 4.5 and 2 = ε 2 = 4.2.
In both experiments, the imaging domain of λ 0 ×λ 0 at the 6.77 GHz frequency of the interrogating field is discretized  with 144 elements situated on equidistant 12 × 12 grid of square elements (pixels). To ensure proper focusing properties of the Veselago lens for the particular imaging domain and its discretization, the background material relative electric permittivity and magnetic permeability are taken to be 4 and 1. The relative permittivity and permeability of the Veselago lens are −4 and −0.999999, respectively. Note that the closer the lens relative permeability gets to −1, the more ideal are its focusing properties, and the more challenging becomes the numerical evaluation of its Green's function [24], [32]. The thickness of the lens is 2 m. The interrogating field is produced by a point source located 2λ 0 away from the origin. The properties of the background material with refraction index of 2 as well as material of the superlens with refraction index of −2 are taken to be close to the experimental values of the background medium index of refraction 1.925 and the lens index of refraction −1.925 reported in the pioneering experimental demonstration of microwave superlens [26]. In order to image a 2-D object we take a substantially thicker lens, however, than the one reported in [26] where the image of a single point source is reconstructed. The larger thickness of the lens allows for large distance between the lens and the data domain S and reduces the effect of fringing fields present between the lens and focal spot (see Fig. 2). In practical experimental set ups, the thickness of the lens could be substantially reduced albeit at some loss of image quality.
In both numerical experiments, J = 16 different locations of the point source are used to increase the reconstruction quality. The scattered field data used in the proposed CSI method was generated synthetically via solution of the forward scattering problem for the same object. To avoid the inverse crime and ensure accurate computation of the scattered field at the receiver points, in the forward scattering solution the imaging domain is discretized with 400 pixels.   Fig. 11 show the evolution of the reconstructed permittivity of the 2-layer square object and 'E'-shaped object, respectively, over the iteration of the CSI algorithms in the case when the imaging experiment is conducted inside the homogeneous background ( b = 4ε 0 and in the case when the focusing media of Veselago lens is present. All the other parameters in the numerical image reconstruction with the CSI algorithm were otherwise the same. The experiments in Fig. 6-9, 11 were conducted in the absence of either noise or loss in the materials (except for non-idealities of the lens). To demonstrate effect of the noise and ability of modified CSI to reconstruct imaginary part of the permittivity additive noise with different signal-to-noise-ratio (SNR) and non-zero conductivity were considered in the cases of 2-layer square object reconstruction in Fig. 10 (c)-(e). Specifically, the sensor locations (i.e., the observation points) remained at the same exactly positions in both experiments. They were positioned to be in a proper correspondence to the locations of the pixels discretizing the imaging domain with respect to the lens' plane of symmetry as depicted in Fig. 5. In the homogeneous background experiment, the lens was removed, however. One can see that while the CSI algorithm in the focusing media consistently improves the image quality as the CG iterations progress, the experiment staged in the homogeneous background did not produce a recognizable image of the object as a result of the CG search for either the 2-layer square or 'E'-shaped scatterer.
The mismatches between the true scattered field and the scattered field produced by the contrast source approximation in the conventional CSI algorithm and its proposed modification due to the presence of the focusing media are shown in Fig. 6 for the case of the 2-layer square object and in Fig. 7 for the case of the 'E'-shaped object. The CG iterations in the homogeneous background CSI algorithm showed similar convergence, prior to flattening out of the sensor mismatch as a function of the iteration count and its stagnation. However, as it can be seen from Fig. 10 and Fig. 11 the reconstructed permittivity remained far in value from the true one at every iteration of the CG search process. This is due to manifestation of the ill-posedness associated with the formulation of the inverse source problem in the absence of the focusing media. The non-uniqueness of the solution results in the existence of multiple contrast source functions w which can produce the same scattered field values at the sensor locations. In the CSI imaging experiment conducted in the focusing media, the ill-posedness of the inverse problem is eliminated. As a result, only one unique and true distribution of the contrast source w can produce the correct values of the scattered field at the sensor locations. Thus, the contrast w during the CG iterative process is seen to converge to its true values when the focusing media is present, and the underlying formulation of the inverse problem is well-posed.
The behavior of the error in the reconstructed contrast is shown in Fig. 8 for the 2-layer square object and in Fig. 9 for the 'E'-shaped object. The error behavior in the sought contrast over iterations is very different for the case of the imaging experiment staged in the homogeneous background and in the presence of the Veselago lens. In the case of the 'E'-shaped object, the error in the reconstructed contrast is growing in the process of CG iterative search when the focusing media of the lens is absent. While such behavior should be possible to avoid in principle by increasing the contribution into the functional due to the mismatch in the contrast source by reducing the contribution of data error η S ρ (n) 2 into the functional (17), we have found that even in the case of substantial dominance of the contrast source term η D r (n) 2 in the functional the error in the contrast still grows over CG iterations as demonstrated in Fig. 9.
The CSI algorithm has been known to successfully produce images in the absence of the focusing media [14], [15]. The reason the reconstruction results of poor quality are obtained in homogeneous background imaging set up are due to the following two factors. First, the location of the sensors in the homogeneous background imaging set-up is sub-optimal for the image reconstruction. Specifically, the observation points are located at the points where they would be best positioned if the lens was present. Optimal positioning of the observation points in homogeneous background CSI algorithm, however, is all around the object as there is no preferential direction to where the field is scattered by the object. Note, however, that in our experiments the positioning of the receive points was intentionally preserved to be the same as in the imaging experiments involving the lens in order to identify the difference in the quality of image reconstruction achieved due to the presence of the focusing media. Second, the imaging experiments in homogeneous media are typically conducted at multiple frequencies [14], [15] in order to obtain recognizable images in the CSI algorithm. The results of the image reconstruction in homogeneous background in this work are conducted at one frequency of the interrogating field which is the same frequency as the frequency used in the CSI image reconstruction when the lens was present.

IV. CONCLUSION
The paper describes a modification of the Contrast Source Inversion (CSI) algorithm which conducts image reconstruction of an object from its scattered electromagnetic field in the presence of a focusing media. The focusing media is formed by the Veselago lens which is properly positioned with respect to the object and the sensor locations at which the scattered field data is collected. Such set up of the imaging experiment in the presence of the focusing media is shown to eliminate the ill-posedness of the inverse scattering problem formulation and to greatly improve the accuracy and robustness of the CSI algorithm. The improvements of the proposed modified CSI algorithm are demonstrated in the numerical image reconstruction from the synthetic scattered field data produced by inhomogeneous objects of complex shape.