Demonstration of Quantitative Microwave Imaging Using an Ideal Veselago Lens

Quantitative microwave imaging employing an ideal Veselago lens is demonstrated using synthetic 2D experiments. Reconstructions of the complex-valued dielectric constant of an arbitrary object-of-interest are obtained from noiseless data using a non-iterative technique. A closed-form Veselago lens Green’s function is utilized in the formulation of the problem. The contrast sources corresponding to a single illumination of the object-of-interest is recovered using total-field measurements at several receiver points located in the lens’ focusing region. Collecting the field in the presence of the ideal Veselago lens produces well-conditioned matrices for the discretized inverse source problem. The contrast is then directly obtained by first calculating the total-field within the object. Synthetic examples having features as small as $\lambda /25$ show that the resolution is only limited by the signal-to-noise ratio of data. Experiments with noisy data reveals that the inversion procedure directly mirrors the noise in the data. Two new filtering processes are introduced based on the truncated Singular Value Decomposition of the data operator and the domain operator involved in the inversion process. These successfully deal with noisy reconstructions. Results show that using the new filtering methods, even with a single illuminating source, outperforms the least-squares method with several illuminating sources.

Quantitative MWI aims to reconstruct the distribution of the electromagnetic properties of an object-of-interest (OI) using the non-invasive electromagnetic interrogation of the OI, that is, without probing the inside of the object. The data that is used to create quantitative images of the real and imaginary parts of the complex-valued permittivity consist of measurements that are made using receivers located outside of the imaging domain, wherein the OI is placed. Sets of data are collected for each position of the transmitter that illuminates the OI.
The reconstruction of electromagnetic properties, in the form of a multidimensional contrast function, requires that one solves, what is conventionally called, the Inverse Scattering Problem (ISTP). The ISTP is a non-linear illposed problem wherein the non-linearity arises because one This work is licensed under a Creative Commons Attribution 4.0 License. For more information, see https://creativecommons.org/licenses/by/4.0/ does not have access to either the total-field or the contrast within the imaging domain. A standard technique in MWI imaging is to first solve for the contrast sources, defined as the product of the total-field and contrast at any particular spatial location. This Inverse Source Problem (ISP) is linear but remains ill-posed [15], [16], [17]. In fact, much of the ill-posedness of the ISTP arises because of the ill-posedness of the ISP. If one could reconstruct the contrast source perfectly then one could easily determine the total field within the imaging domain and, thereby, the contrast [18], [19]. Thus, if an imaging procedure can produce an accurate representation of the contrast sources generated within the OI, for a given illuminating source, then the ISTP becomes well-posed. The instability and non-uniqueness of non-linear ill-posed problems make them difficult to solve.
Most ISTP research has focused on introducing some form of regularization in the mathematical formulation or algorithmic solution of the problem [1], [2], [20], [21], [22], [23], [24]. The non-linearity of the ISTP is mainly handled by utilizing iterative algorithms, although some non-iterative methods can be considered as typically lowresolution techniques, [25]. The ISTP iterative algorithms attempt to minimize some cost-functional that incorporates the measured data with a model that predicts contrast and/or contrast sources. These optimization type algorithms can get trapped in local minima, producing reconstruction artifacts. In addition, the instability of these problems means that small errors in the data measurements can lead to large changes in the solution. Thus, even small amounts of noise can make significant differences in the solution. In some recent contributions, the non-linearity of the ISTP is instead handled by properly recasting the original scattering experiments into new virtual ones, which are focused in some pivot points inside the unknown targets [26], [27], [28]. These virtual experiments allow to condition the spatial behavior of the internal field or of the contrast source, which are also unknowns, and somehow simplify the ISTP at hand.
In 2D, the time-harmonic MWI can be formulated in terms of Helmholtz boundary value problem (BVP), where the wavenumber is an unknown function of position that is to be determined. The ISP is formulated as the solution of this BVP taking the form of a first-kind Fredholm integral equation called the "data equation". The kernel of the data equation is the "Green's function" (GF) for the Helmholtz BVP. From a broader mathematical perspective, the integral operator maps the contrast-source function, i.e., the product of the contrast and the total field inside the imaging domain, to the scattered-field at the measurement points (receiver locations). When MWI is performed in free-space, using the free-space GF, the data operator becomes compact, and the inverse of a compact operator is well-known to be ill-posed [20]. From a qualitative point of view, it is the smoothing properties of the free-space GF that causes a loss of information in mapping contrast-sources to received data.

II. THE STANDARD METHOD OF MICROWAVE IMAGING
In this work, we consider the case where MWI is conducted in the presence of an ideal Veselago Lens (VL). A VL is simply a slab made of double negative material (negative values of permittivity, and permeability, μ). Veselago postulated such a lens and characterized its unique focusing property by introducing a negative refraction index [29]. This focusing property is depicted in Figure 1, where a ray emanating from a point-source located in front of the VL slab, having thickness d, at location A produces focus-points at B and C. All the forward diverging rays generated by the point-source are refracted towards, and converge on, the singular point B inside the VL. The point B is situated at an equal distance from the first VL boundary as the source. The rays emanate from point B to the second boundary, where they are again re-focused to the outside focus point, C, a distance 2d away from the point source. All of the energy emanating from the point-source towards the lens reaches the focus point C. The GF for both 2D and 3D electromagnetic field cases has been provided in an efficient closed-form representation in [30]. This GF effectively translates the singularity of the free-space GF from the point A to point C with phase change.
The use of these unique focusing properties of the VL are well-known, and attempts have been made to create such a lens for imaging purposes [31], [32]. In those applications have focused on qualitative MWI. Some recent publications have dealt with quantitative MWI using an "almost" perfect lens containing some loss to allow the spectral-domain representation of the GF they used to converge (see, e.g., [33], [34]). The authors of [33] took advantage of the use of the spectral-domain representation of the GF in their iterative inversion technique that starts with a one-dimensional reconstruction based on the Born approximation followed by a full 3D holographic inversion. Dealing with the cumbersome spectral-domain representation of the GF does not fully reveal the potential of using an ideal VL for quantitative MWI. The recently developed closed-form representation of the VL GF allows us to directly investigate this potential by studying the matrix operators associated with the discretization of the ISP and ISTP integral operators for the VL imaging scenario [30]. We will show that one of the reasons for the less-than-perfect quantitative imaging demonstrated in [33], [34], even with noiseless data, is the use of the nonideal VL. In those works, loss was introduced into the VL to allow the spectral-domain representation of the VL GF to converge and not to simulate what might be achievable in practice, in terms of creating a VL.
In this research, we prescribe quantitative MWI using an ideal VL and formulate the associated 2D inverse problem. In the mathematical formulation, we use a closed-form definition of the VL Green's function, which was formerly introduced in [30]. We show that, through a direct non-iterative inversion method, the contrast sources of an arbitrary OI can be perfectly reconstructed. This reconstruction is obtained using one single illuminating source and the total-field is measured at a set of receivers located inside the focusing region of the lens. We then use the reconstructed contrast sources to find the contrast of the OI directly by calculating the total-field inside the OI. Furthermore, we show that the presence of the ideal VL in the field collection process produces well-conditioned matrices for the discretized inverse source problem. We define and use two synthetic test targets with feature as small λ/25 to show the impact of the VL in obtaining high-resolution numerical results that are limited only by the signal-to-noise ratio of the measurements. The high-quality reconstructed images verify the remarkable effect of the lens due to replacing the Green's function with a well-conditioned matrix in the discretized inverse problem. We show that using the lens and, consequently, the mirroring property of the new operator effectively alleviates the ill-posedness of the associated inverse scattering problem. The new process offers a direct non-iterative inversion algorithm leading to simple computations. Moreover, we have introduced two new noise filtering methods based on the truncated singular value decomposition of the data and domain operators, leading to a remarkable improvement of the reconstructed images from noisy data. The noise filtering method is not limited to the MWI approach introduced here, and it can be applied in other computational imaging techniques to mitigate the noise.
We first review the standard MWI setup in Section II. In Section III, we describe the VL ISTP setup and provide the mathematical formulation of the forward problem equation as well as the data and domain equations for the inverse problem using the closed-form representation of the VL GF. The matrix operators associated with the descretization of these operators are also described. In Section IV, we describe two simple filtering approaches aimed at suppressing direct noise translation, and in Section V, we provide results for the inversion of noisy synthetic data associated with two test targets. At the end of Section V we provide a numerical study that reveals how the VL mirrors the noise in the acquired data, and how the proposed filtering approaches mitigate this effect.
A schematic of the standard MWI ISTP is provided in Figure 2. The OI is assumed to reside inside the imaging domain D. The domain V with boundary = ∂V, is typically assumed to be unbounded consisting of free-space outside of D. The OI is represented as an inhomogeneous contrast function, χ(r, ω), that may, in general, depend on the positionr and the angular frequency ω. In what follows, we assume time-harmonic problems, having a e jωt time dependence, where ω is fixed, and therefore it will be dropped from the notation.
The OI is surrounded by the set of illuminating transmitters {T 1 , . . . T T } and receivers {R 1 , . . . R M } located outside of the imaging domain D with boundary ∂D. The receivers are used to collect both the incident-field and total-field data. The former is collected in the absence of the OI, while the latter is measured in the presence of the OI. This data can be collected with one or more of the transmitters illuminating the region, but typically, data is collected while only one transmitter is on. The difference between the total field and incident-field is called the scattered-field. This "data" contains information about the material properties of the OI and is "inverted" using various ISTP algorithms.
In this work, we limit ourselves to 2D formulations where the TM z polarization can be assumed. Under such an assumption, the ISTP is formulated in terms of a scalar field, u(r), representing the z-component of the electric field, E uâ z . The total field is decomposed into the incident and scattered-field according to the definition where u i and u s denote the incident and scattered-fields. The wavenumber at any location in space is written as k = ω √ r 0 μ r μ 0 where 0 and μ 0 are the permittivity and permeability of free-space, respectively, and, r , and μ r are the relative permittivity and permeability. In terms of the wavenumber, the contrast is defined as where k 0 is the free-space wavenumber. Note that, here, we are only interested in the contrast related to permittivity not permeability. The contrast source for a particular total field, u(r), is written as The data to be inverted is the scattered-field at each receiver location,r m , and is written as u s (r m ) for m ∈ {1, 2, . . . M}. Using the GF for the problem, G 0 (r|r ), representing the field due to a point source located atr , the expected scattered-field atr m can be written as With the data as the known quantities, this represents a nonlinear integral equation for the unknown contrast and total field within the imaging domain. The scattered and total fields in this equation correspond to a particular illuminating source. Equation (4) is sometimes referred to as the "data equation". This equation can also be thought of as a linear integral equation for the contrast sources, W(r) = χ(r)u(r), generated by each illuminating field. In this form, as mentioned in the introduction, where we are solving for W(r), it is called the ISP. The free-space GF, G 0 (r|r ), is the kernel of these integral equations. It is given analytically as where H (2) 0 denotes the Hankel function of the second kind and zero-th order [35]. The ill-posedness of the ISTP results from the fact that the receiver locations are necessarily outside of the imaging domain D and therefore, the GF operator which takes contrast sources to measured field values ouside D is compact. When the operator is discretized, for example using the Method of Moments (MoM), it results in an ill-conditioned least-squares matrix problem.

III. MWI USING THE IDEAL VESELAGO LENS A. PROBLEM STATEMENT
Access to analytic expressions for the ideal VL allows one to easily incorporate it into the mathematical formulation of performing imaging with a VL. A scheme of MWI in the presence of a physical VL is shown in Figure 3. As shown, the OI is placed within the imaging domain, D, located to the left of a slab (the VL) of thickness d with surfaces located at x = d 1 and x = d 2 . The imaging domain can have any height, but its width is limited to the width of the VL because only contrast sources within a distance d from the lens are focused onto the opposite side of the lens.Here we choose a square imaging domain with height also equal to d.
As in the standard approach, two sets of measurements are performed: one with the OI within D, producing the totalfield, and one set of measurements without the OI in D, producing the incident-field. The VL slab is assumed to be present for both of these measurements.
The receivers are located within a measurement domain D R , on the opposite side of the VL from where the OI is located. The measurement domain is a translated version of the imaging domain, where the translation is exactly 2d to the right (as shown in Figure 3). For each centroid of a cell in a discretization of D a corresponding receiver point is introduced within D R . Thus, for example, in Figure 3 we see a vertical column of receiver points labeled R 1 , R 2 · · · R M located a distance of 2d from the corresponding column of the centroids of the cells inside the discretized OI located inside the imaging domain D denoted R 1 , R 2 · · · R M . In other words, corresponding to the i th column of cells of the discretized OI located a distance l i from the VL, there will be a column of receivers on the other side of the VL located a distance of d − l i from the VL.
The unique focusing properties of the VL allows one to non-invasively measure the scattered-field data inside the inspected domain. This represents a crucial difference with respect to the standard MWI ISTP, as a much larger amount of independent data is available and it also alleviates the ill-posedness of the problem, as discussed in the following. This phenomenon has some similarities with the magnetic resonance imaging based electrical properties tomography, which involves the solution of an inverse scattering problem wherein the RF magnetic field inside the human body is processed [36], [37], [38].

B. THE GREEN'S FUNCTION FOR THE VESELAGO LENS
The formulation of the ISP in the presence of the VL consists of the exact same data equation as in (4) but with the scalar GF being replaced with that of the VL. The new GF operator, G VL (r|r 0 ), represents the field at a receiver point due to a point-source located atr 0 in the presence of the VL slab.
An analytic expression for the VL GF was recently introduced in [30]. For a point source located atr 0 = (x 0 , y 0 ) and the VL described as in Section III-A the VL GF, G VL , is the solution to the Helmholtz partial differential equation (PDE) boundary value problem (BVP) given as The three distinct Regions I, II and III correspond to the regions to the left of, within, and to the right of the VL slab.
With the illuminating source located along the x = x 0 = 0 line, these regions are identified as: and x > d 2 , respectively. Notice that the same PDE operator is utilized in each region, because even within the VL slab the wavenumber is k 0 . The effect of the VL is taken into account by the boundary condition on the slab as described in [30]. The solution to the VL GF is given as where H (1) 0 denotes the Hankel function of the first kind of zeroth order, where R = r −r 0 , R 1 = r −r 1 , R 2 = r −r 2 ,r 1 = (2d 1 , y 0 ), andr 2 = (2d, y 0 ). Note that the coordinates r 1 and r 2 are the locations of the singularities produced by the ideal VL, and they both depend on the location of the point source, r 0 . Also, note that, x 0 , the xcoordinate of r 0 , is chosen so that Regions II and III are subdivided into two regions. This solution was obtained by "patching-together" five solutions of the PDE and taking into account power-flow considerations between the regions. The two planes passing through these singularities have been referred to as "power-transfer singularity planes", [30]. In the sequel, the abbreviated notation G VL (r|r 0 ) will be used, because the location of the ideal Veselago lens and consequently the associated parameters d 1 and d will remain fixed throughout this work.

C. THE VESELAGO LENS INVERSE PROBLEM 1) GENERATING SYNTHETIC DATA
To generate synthetic data for a given contrast, we must first solve for the total-field within the OI using the equation With the illuminating source at positionr 0 , the incident-field is just u i (r) = G VL (r|r 0 ). Note that, on the left side of the ideal VL, the GF is identical to that of free-space: Thus, within D the fields are the same as in the standard, free-space, problem.
Once the total-field within the OI is found, the scatteredfield data at each receiver point, u s (r m ),r m ∈ D R , can then be computed using the VL data equation written as Integral equation (8) is discretized using MoM in which a total-field vector, u ∈ C N , where the k th element represents the total-field at the k th centroid, is introduced. The discretized forward operator, in (8), consisting of G VL , χ and the scalar k 2 0 , becomes a square matrix denoted as G χ . This matrix maps total-fields over the cells to the scattered-fields at the centroids of each cell. The total-field inside the OI for the known contrast is then computed as where u i is the vector of incident-fields at each centroid. This results in the same total-field as in the standard MWI problem because G VL = G 0 . Details of computing G χ are given in the Appendix A. For creating synthetic data the total-field calculated via the MoM is used in the data equation (9) to compute the scattered-field at each of the receiver points,r m . The same discretization of the imaging domain, D, as in the MoM solution for the total-field is utilized. For a given 2D contrast function, the discretized VL data operator in equation (9), G χ VL ∈ C M×N , can be computed by integrating over each cell of the imaging domain. The rows of this matrix correspond with one of the M receiver points, while the columns correspond with the cells of the discretized D.
Performing this integration must take into account the special nature of the VL GF given in (7), in that, if an integration region within a cell, translated by a distance 2d to the right, lies to the left of the receiver point, then H (2) 0 must be used in the integral, whereas if the region is to the right of the receiver point, then H (1) 0 must be used. The resulting matrix equation for generating data, d ∈ C M , using the total-field, u ∈ C N , obtained by (10) is In a practical setup, measured data must be calibrated to remove differences between the model used in the ISTP and the actual measurement setup, [23], [39], [40], [41]. For our synthetic experiments, to simulate data that would be obtained by measurements, noise is added to the synthetic data. The noise is added as follows: where NP indicates noise percentage, RV is a random vector created by a uniform distribution of numbers belonging to the interval (−1, 1) and RMS is the root mean square of data values. This is a standard technique when creating synthetic data, [42].

2) GENERATING ISP MATRIX
As was previously mentioned, the receiver points constitute a square grid of M points. For the ISP, the imaging domain is re-discretized into M cells with the coordinates of the centroids of each of these cells,r i , i ∈ {1, 2, . . . , M}, located a distance 2d to the left of the receiver points. The contrast source is reconstructed at each of these points. Thus, a column vector of contrast sources, w, is introduced where the i th element of w is W i u(r i )χ (r i ). These are assumed to be constant in each cell and pulled out of the integration in (9). The integration is performed in a similar manner as in creating the forward data operator, G χ VL . The new operator is denoted G VL ∈ C M×M . This operator takes the discretized contrast-source vector to a corresponding scattered-field data vector. It will be inverted to obtain the reconstructed contrast source from a given synthetic noisy data vector of scatteredfields, d N . This is a well-conditioned square matrix because of the focusing nature of the VL GF (i.e., the kernel of the integral operator). In other words, the role of the VL is to produce an invertible one-to-one mapping between the scattered-field at the M receiver points and the contrast at the M centroids of the discretized imaging domain, χ. Symbolically, we have Once the contrast source vector is obtained, it is utilized to find the unknown contrast vector, χ, using the definition of contrast sources: where indicates the element-wise vector product. Thus, to find the contrast we need to first reconstruct the totalfield, u, within the OI using the w. At each of the centroids of the reconstruction grid,r m for m ∈ {1, 2, . . . , M} ∈ D R , the total-field is calculated as The incident-field in (15) is obtained by using definition of the VL GF in (7) as u(r i ) = G VL (r i |r 0 ), for the illuminating point source located atr 0 .Therefore, using the reconstructed contrast sources, found in (14), the reconstructed total-field inside the OI can be found by the matrix equation wherein G D VL ∈ C M×M . By virtue of the special nature of the VL GF, the fact that the reconstruction grid centroids are on the same side of the VL as the OI means that G VL = G 0 in equation (15). The reconstruction grid is not necessarily the same as that used in the forward problem, so the reconstructed total-field will not be identical to the synthetic total-field. As will be shown in the results section, when an "inverse crime" is committed wherein a reconstruction mesh that is identical to the forward mesh is used, and no noise is added, we obtain perfect reconstructions of both the contrast and the total-field. This perfect reconstruction is obtained even for the case of a single illumination source.
When noise is added and the inversion mesh is not identical to the forward mesh the reconstructions degrade, even with the well-posedness of the VL ISP. This is expected because although G VL is well conditioned, its condition number for the cases to be presented in the Numerical Results section is still on the order of 10 3 . Thus, the relative error in the reconstructed contrast sources that is due to the relative error in the scattered-field data, because of the added noise, can be amplified by as much as the condition number. This will result in error in the total-field obtained using (16), making the subsequent error in the reconstructed contrast, using (14), considerable.

IV. REGULARIZING THE CONTRAST INVERSION
Two standard techniques to mitigate the error in the recovered contrast are considered herein. The first, applied to the case of a single illuminating source, is to apply a Truncated Singular Value Decomposition (TSVD) to perform the contrast source recovery. The second method is to utilize more than a single illuminating source: each produces different contrast sources for the same contrast. As will be seen, reconstructions using these initial regularization techniques reveal some interesting features of inverting noisy data using G VL .

A. SINGULAR VALUE DECOMPOSITION
The G VL matrix includes small singular values which when inverted amplifies the noise. The small singular values can be removed from G VL when performing the inversion of the data. There are various ways to determine a suitable truncation parameter when performing TSVD [43] , here we simply consider the ratio of the singular values of G VL compared to the largest singular value. We choose a truncation parameter, τ svd , that eliminates singular values having a ratio to the largest value that is below this parameter.

B. LEAST SQUARES METHOD FOR SEVERAL ILLUMINATING SOURCES
The second method to regularize the problem is to utilize several illuminating sources. The set of all the noisy reconstructed contrast sources is used in a least-squared sense to determine the contrast. This is performed by minimizing the following sum of errors: where w T t is the contrast source obtained by solving (14) for each illuminating source T t , t = 1, 2, . . . , T.

A. IMAGING SETUP
All presented experiments are performed using a monochromatic illuminating source (i.e., single frequency), but all our dimensions are chosen with respect to the wavelength. The only essential dimension of the problem is that the left boundary of D should be within a distance d away from the VL, meaning that the largest total width of D will be equal to the width of the VL. In what follows, we provide all measurements in terms of wavelength, λ. For the specific results that follow, the illuminating source is positioned on the same side as the OI. When a single illuminating source is used it is placed a distance d/2 from the VL and one wavelength above the imaging domain. For the case of several illuminating sources, the sources are placed equidistant from each other on a semi-circle of radius 1.2λ from the center of the imaging domain (to the left of the VL).

B. TEST TARGET RECONSTRUCTIONS
As shown in Figure 4a, the first test target consists of two different 2D Gaussian functions inside the square λ × λ imaging domain. Note that for all the figures in this paper the color maps are chosen corresponding to the maximum and minimum values of the reconstructed contrast. The maximum values of permittivity are chosen to be 2 − 0.4j and 3 − 0.2j for the left and right Gaussian functions, respectively. This target was chosen because, although this lossy target possesses a continuous profile, it contains very small subwavelength features that are difficult to reconstruct for any imaging technique. Figure 4b shows the reconstructed permittivity obtained using a single illuminating source and the simple matrix inversion of (13) and (14). The imaging domain, D, is discretized with 25 cells in both the x-and y-directions (a total of 625 unknowns reconstructed from 625 scattered-field measurements). This perfect reconstruction is obtained when no noise is added to the generated synthetic data.
For the next example, the same target is used but noise is added to the synthetic data, using formula (12) with NP = 2%. The reconstructed permittivity in this case is shown in Figure 4c. The recovered contrast is considerably degraded, as expected. The two mitigation techniques previously described were also utilized with the noisy data. The contrast reconstruction using 40 and 100 illuminating sources and the Least Squares Method, described in Section IV, are shown in Figures 4d  and 4e, respectively. Random noise was added to each of the 40 and 100 scattered-field data sets corresponding to each of the 40 and 100 illuminating sources, respectively. Reconstruction results using TSVD with τ svd = 1% and 2% are shown in Figures 4f and 4g. As can be seen in the figures, the TSVD with τ svd = 2% provides the best reconstructions.
As shown in Figure 5a, the second test target is what we refer to as the "E-phantom." It has a permittivity of 3 − 0.2j and is placed inside a square λ × λ imaging domain with free-space as the background. As in the previous example the imaging domain was discretized into 25 cells in both the xand y-directions (a total of 625 unknowns reconstructed from 625 scattered-field measurements). This target was chosen to test the resolution of the VL imaging procedure as it has features as small as the width of one cell in the imaging domain discretization, λ/25. Figure 5b depicts the reconstructed permittivity obtained using a single illuminating source and the simple matrix inversion of (13) and (14). Again, here we have obtained perfect reconstructions when no noise is added to the synthetic data. When noise is added to the data, using formula (12) with NP = %2, the results degrade considerably as shown in Figure 5c. Using 40 and 80 illuminating sources with the Least Squares Method described in Section IV, the results shown in Figures 5d and 5e are obtained. In this case, for each data data set corresponding to a different illuminating source different random noise was added. Reconstruction results for a single illuminating source and using TSVD with τ svd = 1% and 1.2% are shown in Figures 5f and 5g. As can be seen in the figures, the multiple sources Least-Squares solution provides the best reconstructions, in the real-part, but at the cost of using 80 or even more sources. The TSVD with τ svd = 1% obtains good reconstructions of the real-part of the permittivity but struggles with the imaginary part.

C. EFFECT OF NOISY DATA ON VL INVERSIONS
The results that have just been presented reveal that the VL does indeed create a well-posed ISP by mirroring the contrast sources into the measurement domain. Unfortunately, during the direct-inversion process the noise that is present in the collected scattered-field data is directly mirrored back into the reconstructed contrast-sources. Even using multiple illumination sources is not always effective in removing this noise, and it may not be realistic to use such a large number of sources, e.g., 80 illuminating sources.
To track the effect of the noise on the VL inversion process, we explicitly introduce a variable, n, such that the noisy synthetic data is represented as The reconstructed contrast sources, w by (13) are now where w True = G −1 VL d is the complex vector of true values of the contrast sources and w = G −1 VL n is the error in the reconstruction due to the added noise.
The added synthetic noise that was added to the vector of synthetic data can be considered as high-frequency spatial noise when that vector is construed as a 2D image. This high frequency spatial noise is amplified by the inverse of G VL resulting in poor reconstructions. The noisy contrast sources for The true values in Figures 7a and 8a along with the previous targets are shown in Figures 7b and 8b. When we reconstruct the scattered-field within the imaging domain by operating on these noisy contrast sources with G D VL and add the known incident-field we get the reconstructed total-field. Surprisingly, the reconstructed total-fields are in good agreement with the true fields: noise is considerably dampened by the G D VL operator. The true and reconstructed total-fields associated with the targets shown in Figures (4a) and (5a), are displayed in Figures 7c and 8c, and Figures 7d and 8d, respectively.
Thus, we see that the G D VL operator has a low-pass filtering effect on the high spatial frequency noise that is in the reconstructed contrast sources. With this observation, introducing U = diag(u) and w = U χ , we note that we can write the scattered-field within D as Consequently, the contrast can be obtained as Algebraically, this is equivalent to χ = U −1 w or the elementwise division of the contrast sources by the total-field: w./u. But equation (21) gives us the opportunity to use the filtering effect of G D VL on w before multiplying by the inverse of G D VL U. We notice that applying the (G D VL U) −1 brings back the noise because of the poorly conditioned matrix G D VL . So instead we perform a TSVD of G D VL before the inversion: to get the final equation for the filtered inversion of the contrast: We call this TSVD domain-operator filtering. Note that this is independent of the ISP, or contrast-source recovery step, where we applied the inverse of the TSVD of G VL . The use of the spectral property of the mapping from the induced contrast sources to the fields inside the OI has also been exploited in [44], wherein, in the iterative procedure underlying the subspace optimization method, the ambiguous part of the induced sources are constructed by the singular vectors of domain operator, chosen according to their influence on the fields inside the domain of interest.
The complete inversion process, starting from the data, can be written as Results of adding this TSVD domain operator filtering step to the previous reconstructions are provided in Figures 9a-10c for the same examples previously presented. The first row in each figure shows the reconstructed contrast truncating only the G VL , operator with a truncation parameter of τ svd = 2% for the first target and τ svd = 1.2% for the E-phantom. The second row in each figure shows reconstructions using only the truncated domain operator filter with a τ svd = 1.5% for the first target and τ svd = 0.5% for the Ephantom. The last row shows results when both techniques are utilized. In this case, the same truncation parameter was utilized for both inversions, but these are different for the two targets as shown in the caption.
The errors in the reconstructions can be quantified by using the Root Mean Square Error (RMSE) of the images normalized by the maximum value in the true image, denoted by e χ , is calculated as where N is the number of cells inside the imaging domain. This error is calculated for the real, imaginary and absolute values of the reconstructed contrast and are shown in   Tables 1 and 2. The normalized RMSE for the different regularization techniques are provided in each row of the tables.
Comparing the normalized RMSEs in Table 1, the best result for the first target is obtained using TSVD on both G VL and G D VL with τ svd = 2%. Table 2 shows that the best results for the E-phantom are obtained using Least-Squares with 80 illuminating sources. However, the next best result for the real part of the reconstructed contrast of the E-phantom is obtained by applying TSVD on both G VL and G D VL with τ svd = 0.51%, which seems an acceptable error given that using 80 illuminating sources would be difficult in practice. The effect of truncating the SVD to obtain the reconstruction is clearly revealed by looking more closely at how the noise and data are mapped to the imaging domain. The SVD of the VL operator is written as G VL = P S Q * , where G VL is the VL Green's operator, S = diag(σ i ) ∈ C N×N is the diagonal matrix of singular values, P and Q are the unitary matrices in C N×N , and * indicates the Hermitian. Consequently, the pseudo-inverse of the VL Green's operator can be found as [1]. It is worthy of note that G VL does not produce zero singular values. However, in order to filter the mirrored noise, the very small singular values still need to be truncated. Using this pseudo-inverse and relation (18) the contrast sources can be found as where p i is the i th eigenvector. Setting up the ratios for i = 1, . . . , N, and plotting the obtained values corresponding to each column of the eigenvectors shows an increasing behaviour. This behaviour demonstrates how the type of noise used herein overtakes the reconstruction as its projection onto the eigenvectors corresponding to smaller singular values is larger than the projection of the data onto the same eigenvectors. A similar ratio can be formed for G D VL . The logarithmic graphs of the singular values of G + VL and the ratios obtained by (27) for the targets presented in this paper are shown in Figures 6a-6c. The figures also show where the SVDs have been truncated for each target.

VI. CONCLUSION
We have demonstrated quantitative MWI incorporating an ideal VL and using a simple two-step non-iterative inversion technique. The first step reconstructs the contrast sources by directly inverting the ISP, and the second step reconstructs the contrast by first reconstructing the total-field within the imaging domain using those contrast sources and then simply dividing the contrast sources by the recovered total-field. Such a two-step procedure is not available with the traditional ISTP in free-space because of the ill-posedness of the free-space ISP. Even with regularization, the free-space ISP is difficult to solve due to non-radiating sources; the free-space ISTP is typically regularized by utilizing several illuminating sources where the contrast remains the same for each illumination. The 2D demonstration was performed by applying the recently introduced closed-form GF for the lossless VL. The VL ISP is easily set-up by replacing the free-space GF with that of the VL in the free-space ISP data equation. Discretization of the VL data-equation produces well-conditioned matrices that are easily inverted. When no noise is added to the data perfect quantitative images are obtained, at least for the several examples of OIs that were tried, including the two OI examples considered herein. Features as small as λ/25 are easily reconstructed, leading us to believe that the VL ISP is well posed with no nullspace. To be clear, the VL ISP data-operator formulated herein utilizes a measurement domain that is a 2d translation (i.e., twice the thickness of the VL) of the imaging domain. Within the measurement domain we take as many measurement points as there are cells within the imaging domain. This creates a one-to-one correspondence between the focus points of the cell centroids and the measurement points. With this set-up, the discretized VL matrix problem becomes well-conditioned and we believe that this is because the VL data-operator becomes well-posed (a formal proof would be required to make this statement definitive). It is difficult to say whether a well-posed ISP could be achieved with a different domain and/or range of the operator.
Although the VL data-operator is well-conditioned (having a condition number on the order of 10 2 ) we have shown that the small singular values and their associated eigenvectors will "mirror" the noise present in the collected data. To mitigate the noise mirroring issue, several approaches were introduced. We first considered simply using TSVD when inverting the VL data operator. Secondly we considered using several illuminating sources and a Least Squares technique for contrast recovery. The former approach shows better results in both of the numerical examples considered and has the advantage of using a single illuminating source. Finally, a new domain-operator filtering approach was introduced that can be used with either of the two previous methods to further improve the noise reduction.
The methods described in this paper are easily applied to the 3D vectorial case; the results are dependent on having an accurately computable Green's function for the ideal VL. Although it may not be achievable to actually create such a lens, it is important to have, at hand, the expected imaging improvements were such a lens created. Research into the potential to incorporate a virtualized version of a VL into the free-space MWI procedure is ongoing.

APPENDIX COMPUTING THE DESCRETIZED VESELAGO LENS OPERATORS
Computing the matrices of the discretized VL operators requires special attention. In the forward problem, solved using MoM, for each observation point defined as the centroid of each cell in the discretized imaging domain, D, for a given 2D contrast function χ(r), evaluating the elements of G χ requires the computation of the double integrals G χ m,n = − jk 2 0 4 y n + y /2 y n − y /2 where x and y are the width and height of the equal-sized cells inside the discretized D, R m,n = (x m − ξ) 2 + (y m − ζ ) 2 for x n − x /2 ≤ ξ ≤ x n + x /2 and y n − y /2 ≤ ζ ≤ y n + y /2. The variablesr n = (x n , y n ) andr m = (x m , y m ) are the coordinates of the source and observation points corresponding to the centroids of the n th and m th cells, for 1 ≤ m, n ≤ N. These integrals were computed using the MATLAB built-in function integral2 with a relative tolerance in the range of 10 −12 . This handles the singularity of the Hankel function and allows the accurate computation of the operator. Thus, the smooth contrast functions used as examples in this paper are not approximated using a pulse basis, as they are in other computational approaches such as Richmond's method [45].
The synthetic data consists of the scattered-field at each of the receiver points. The same discretization of the imaging domain, D, as in the MoM solution for the total-field is utilized. For a given 2D contrast function, the discretized VL data operator in (9), G χ VL ∈ C M×N , can be computed by integrating over each cell of the imaging domain. The rows of this matrix correspond with one of the M receiver points, while the columns correspond with the cells of the discretized D in the forward problem. The receiver points lie on the opposite side of the VL where the OI is located. The inversion mesh is discretized with cell sizes given by x and y (potentially different than x and y ). In computing the integrals, there may be a switch between the Hankel function of the first and second kind. The general integral for an observation pointr m = (x m , y m ), can be written as G χ VLm,n = y n + y /2 y n − y /2 where R 2m,n = ((x m − 2d) − ξ) 2 + (y m − ζ ) 2 , for 1 ≤ m ≤ M, and 1 ≤ n ≤ N. Based on the value of R 2m,n within the integration region, computing the double integrals can be classified into three cases: 1) When x n − x /2 < x m − 2d < x n + x /2, the singularity occurs within the integration region, and the integral must be split into two regions based on the GF: 0 k 0 R 2m,n dξ dζ, + jk 2 0 4 y n + y /2 y n − y /2 To ensure that neither of the two integration regions are too narrow in the x-direction, a minimum width of 10 −15 is used; that is, if the width of one of the integration regions becomes less than 10 −15 , the region is ignored for this case. 2) When x m −2d ≤ x n − x/2, the singularity occurs either on the left boundary or outside on the left side of the cell. In this case, the integral must be evaluated as 3) When x m −2d ≥ x n + x/2, the singularity occurs either on the right boundary or outside on the right side of the cell. In this case, the integral must be evaluated as For creating the square inverse operator, the integration to calculate G VL ∈ C M×M is performed in a similar manner as in creating the forward data operator, G χ VL . The m, n th element is obtained by evaluating the following integral: G VLm,n = y n + y /2 y n − y /2 (33) Here again, the value of R 2m,n produces three cases: 1) When x n − x/2 ≤ x m −2d ≤ x n + x/2, the singularity occurs within the integration region, and the integral must be split into two regions based on the GF: 0 k 0 R 2m,n dξ dζ, For this case there is no need to consider a tolerance because if the singularity occurs within the cell, it is always at center of the cell, and so the two sub-regions before and after the singularity will be of the same size having width x /2 (will not become narrow). 2) When x j − 2d < x i − x/2, we have G VLm,n = jk 2 0 4 y n + y /2 y n − y /2 3) When x m − 2d > x n + x/2, we have She received the M.S. Laurea degree (summa cum laude) in electronic engineering and the Ph.D. degree in information engineering from the Università Mediterranea of Reggio Calabria, Italy, in July 2012 and May 2016, respectively, where she is currently working as an Assistant Professor within a tenure track position. Her research activity mainly concerns electromagnetic inverse problems, with particular interest in: 1) inverse scattering problems from both a theoretical and applicative point of view and 2) field intensity focusing and shaping in non-homogeneous and unknown scenario, in the framework of hyperthermia treatment planning, wireless power transfer, and MRI shimming. She was a recipient of the "G. Barzilai" Award from the Italian Electromagnetics Society in 2014, while in March 2016, she received the Honorable Mention from IEEE-Antennas and Propagation Society (Central and Southern Italy Chapter) in the Best Student Member Paper Competition. Moreover, she was the recipient of the URSI Young Scientist Award in 2018.
TOMMASO ISERNIA (Fellow, IEEE) received the Laurea degree (summa cum laude) and the Ph.D. degree from the University of Naples Federico II, Naples, Italy. He is currently a Full Professor of electromagnetic fields with the Università Mediterranea of Reggio Calabria, Reggio Calabria, Italy, where he serves as the Supervisor of the LEMMA Research Group and as the Director of the Department of Information Engineering, Infrastructure and Sustainable Energy. His current research interests include inverse problems in electromagnetics, with particularly emphasis on phase retrieval, inverse scattering, and inverse source problems, as well as their applications to antenna synthesis and e.m. fields shaping for biomedical therapeutic applications. He was a recipient of the "G. Barzilai" Award from the Italian Electromagnetics Society in 1994. He has served as a member of the Board of Administrators of Consorzio Nazionale Italiano per le Telecomunicazioni from 2013 to 2019, and he is currently a member of "Senato Accademico" of Università Mediterranea.