A Multiplicative Regularizer Augmented With Spatial Priors for Microwave Imaging

The standard weighted ${L_{2}}$ norm total variation multiplicative regularization (MR) term originally developed for microwave imaging (MWI) algorithms is modified to take into account structural prior information, also known as spatial priors (SPs), about the object being imaged. This modification adds one extra term to the integrand of the standard MR, thus being referred to as an augmented MR (AMR). The main advantage of the proposed approach is that it requires a minimal change to the existing MWI algorithms that are already equipped with the MR. Using two experimental data sets, it is shown that the proposed AMR 1) can handle partial (incomplete) SP and 2) can, to some extent, enhance the quantitative accuracy achievable from MWI.


I. INTRODUCTION
Microwave imaging (MWI) is a noninvasive imaging method, in which the quantitative images of the relative complex permittivity profiles of the objects of interest (OI) can be created. This is often done by applying an electromagnetic inverse scattering algorithm (or simply an inversion algorithm hereafter) to the microwave scattered field data collected outside the OI. In some situations, the measured data might not have sufficient information content to reconstruct the OI to the desired accuracy. In such cases, one option is to provide the inversion algorithm with some prior information. This communication is focused on a method to inject prior structural information, also known as spatial priors (SPs), into the inversion algorithm. This SP may come from a different higher resolution imaging modality such as magnetic resonance imaging (MRI) [1] or ultrasound imaging [2], [3]. For example, in [4], an image registration approach was proposed to use the MRI images in MWI reconstruction. Besides, a combined imaging technology named magnetic resonance microwave tomography was proposed in [1]. Instead of using SP in MWI reconstruction, another approach is to perform joint inversion of electromagnetic data with another type of data such as acoustic data to improve the achievable image accuracy [5]. The joint inversion of such multiphysics data can be, for example, linked through appropriate similarity regularization terms such as cross-gradient functions [6]. The joint inversion is not within the scope of this communication.
An incorporation of SP into MWI inversion algorithms has already been considered by several authors using different approaches [1], [2], [7]- [15]. These methods can be classified (at least) under three categories: 1) those that favor similarity of the reconstructed permittivity values in each region within the SP, which we refer to as similaritybased approaches, e.g., [9]; 2) those that use an inhomogeneous background medium derived from the SP for the inversion algorithm, e.g., [13]; and 3) those that use special basis functions obtained based on the SP for the expansion of the unknown permittivity profile in the inversion process, e.g., [12]. The proposed algorithm in this communication enforces SP in a different manner: via providing edge information for the inversion algorithm based on the SP. (An edge represents the transition from one region to another in the complex permittivity image.) We refer to this type of SP algorithms as edgebased approaches. Noting that edges are associated with the gradient of the image, the edge-based approaches are somehow related to similarity-based approaches through the gradient operator. It should be noted that the proposed edge-based algorithms can retrieve the permittivity of those regions that are absent in the incomplete SP. (That is, in such regions, these SP algorithms switch back to the regular inversion.) This is an indication that the proposed edge-based approach belongs to the class of soft prior regularization as opposed to hard prior regularization [8].
The main feature of this algorithm is that it only requires a minimal change to the existing the Gauss-Newton inversion (GNI) and contrast source inversion (CSI) algorithms equipped with the standard multiplicative regularizer (MR), e.g., [16], [17]. (The term "standard" indicates the common MR used in the literature: the weighted L 2 norm total variation multiplicative regularizer, e.g., [17].) Using the proposed SP technique, the existing MR-GNI and MR-CSI algorithms can be easily changed to accommodate the available SP by augmenting the integrand of the standard MR by an extra term. Thus, we refer to it as an augmented MR (AMR). When the proposed AMR is used with GNI and CSI, we refer to the resulting algorithms as AMR-GNI and AMR-CSI, respectively. Herein, we first briefly review the standard MR and then describe the small change that needs to be done in order to incorporate SP into the existing MR. We then evaluate the performance of AMR-GNI and AMR-CSI against two experimental data sets for 2-D scalar MWI. Besides, we compare the performance of the proposed algorithms against a recently proposed SP algorithm that includes an MR term and a similarity-based SP term in conjunction with the GNI algorithm which we refer to as the MRSP-GNI algorithm [14]. The MRSP-GNI algorithm belongs to the class of similarity-based approaches. This is in contrast to the AMR-GNI and AMR-CSI algorithms, which are edge-based approaches.

II. MULTIPLICATIVE REGULARIZER (MR) -REVIEW
The MR-GNI and MR-CSI algorithms use the following multiplicative regularization term [17]: to regularize their associated cost functionals at the nth iteration of the algorithms. In (1), A is the area of the imaging domain D, χ is the unknown relative complex permittivity contrast profile, χ n is the known estimate of χ at the nth iteration, r is the position vector spanning the imaging domain, ∇ is the gradient operator, and |·| is 0018-926X © 2020 IEEE. Personal use is permitted, but republication/redistribution requires IEEE permission. See https://www.ieee.org/publications/rights/index.html for more information.
the magnitude operator. The steering parameter δ 2 n , which is a real number, gets smaller as the inversion algorithm gets closer to the final solution. (For the expressions of δ 2 n in the MR-GNI and MR-CSI algorithms, see [16] and [17], respectively.) The operation of MR has been described in the previous works such as [17]- [19]; herein, we only provide a quick overview. To this end, let us consider the gradient operator of the MR at the nth iteration of the algorithm that operates on a given χ [20, Appendix D.3] where "∇" is the divergence operator, and b 2 n (r) is an inhomogeneous weight given as ( In the first few iterations, δ 2 n is dominant as compared to |∇χ n (r)| 2 at most of the positions r. Therefore, in the early iterations, b 2 n (r) can be approximated by the constant number 1/(Aδ 2 n ). Consequently, in the early iterations, we have where "∇ 2 " is the Laplacian operator. That is, in the early iterations, MR works similar to the Laplacian regularization which aims to smooth the reconstruction. After early iterations, the relative magnitude of |∇χ n (r)| 2 and δ 2 n should be considered locally in the imaging domain. If in a region within the imaging domain, δ 2 n is dominant compared to |∇χ n (r)| 2 , the above Laplacian approximation is still valid locally, and the regularization operator then attempts to smooth out that region. Otherwise, it favors to reconstruct an edge in that region.

III. AUGMENTED MULTIPLICATIVE REGULARIZER (AMR)
Let us now consider that we have some available SP about the OI. Given the SP, we can then guide the inversion algorithm regarding the locations, at which some edges (i.e., the boundaries between different regions) are to be expected. To this end, we not only rely on the relative magnitude of |∇χ n (r)| 2 and δ 2 n for edge detection as in (1) but also we introduce an extra term based on the available SP to further guide the inversion algorithm in edge detection. Thus, we modify the MR given in (1) to an augmented MR (or, simply AMR) as follows: As can be seen from (5), we have now augmented the integrand of the MR with an extra term, Q 2 |∇ P(r)| 2 , where P incorporates SP. In the implementation, P is a matrix containing some numbers, and each of them represents an expected region within the imaging domain D.
For example, if the SP contains four regions, we assign the values of 1, 0.75, 0.50, and 0.25 to each of these regions. 1 Besides, Q is 1 In this communication, assuming N regions, we start by assigning 1 to a given region, then 1 − (1/N ) is assigned to another region, and finally, 1 − ((N − 1)/N ) is assigned to the last region. The assignment of these numbers to different regions is ad hoc. Based on the numerical experience, the performance of the proposed algorithm is not very sensitive to the order, by which we assign these numbers to different regions. Furthermore, the reason behind choosing the maximum value of P to be one is based on assuming that a typical maximum contrast value in MWI (assuming the presence of a high permittivity matching fluid) is about one. Thus, throughout this communication, the maximum value of P was set to one in an attempt to make |∇χ n (r)| 2 and |∇ P(r)| 2 relatively balanced in (5). (This is despite the fact that the maximum contrast value in the imaging domain is 2 in the first example considered in this communication.) a scaling factor that controls the relative weight of |∇ P| 2 . In all the examples shown in this communication, this scaling factor has been set to Q = 1. Furthermore, in (5), "Edge I" refers to the ability of this regularizer to retrieve the edges from the reconstructed contrast at the previous iteration, and "Edge II" refers to the ability of this regularizer to retrieve the edges from the provided SP. One question that may arise is as follows: why do we use both |∇χ n (r)| 2 and |∇ P(r)| 2 in the AMR? The main reason lies in the fact that if the SP only provides partial (incomplete) structural information, the presence of |∇χ n (r)| 2 may still enable the inversion algorithm to retrieve the edges that are not present in the SP. Besides, if some of the edges provided by the SP are wrong, the presence of |∇χ n (r)| 2 may alleviate the effects of the wrong edges to some extent. Similar to the MR, the gradient of the AMR at the nth iteration will be where The operator L AMR n works similarly as L MR n with only one difference: besides the relative magnitude of |∇χ n (r)| 2 and δ 2 n , we now have an extra term, Q 2 |∇ P(r)| 2 , that also plays a role in edge reconstruction. Furthermore, note the simplicity of incorporating SP within the existing MR-GNI and MR-CSI algorithms. Herein, the only modification performed was to replace b 2 n given in (3) by a 2 n given in (7). On this replacement, the two algorithms are referred to as AMR-GNI and AMR-CSI.
Finally, let us address this question: what happens if we remove |∇χ n | 2 from the denominator of (5)? That is, what happens when we convert AMR to the so-called AMR (Type II) as As expected, the edges are now merely determined by the relative magnitude of Q 2 |∇ P(r)| 2 and δ 2 n . Let us assume that there is no available prior edge information in certain regions within the imaging domain; then, |∇ P(r)| will be zero in those regions. Thus, similar to (4), the regularization operator switches to the Laplacian operator in such regions. In summary, if we use (5), edge preserving might still be applied to the regions that are absent in the incomplete SP. However, if we use (8), the Laplacian operator will be applied to such regions. We refer to the use of (8) as AMR (Type II). 2

IV. EXPERIMENTAL RESULTS
We consider two previously reported experimental data sets to evaluate the performance of the AMR-GNI and AMR-CSI algorithms. The first data set is the so-called FoamTwinDielTM data set [22] collected by the Fresnel Institute in an anechoic chamber. The second data set was collected by the Electromagnetic Imaging Laboratory at the University of Manitoba from a human forearm [16] in the matching fluid of salt water. In all the cases presented below, all the inversion algorithms start with a trivial initial guess, which is the zero contrast for the AMR-GNI algorithm, and the back-propagated solution for the AMR-CSI algorithm.

A. FoamTwinDielTM Data Set
This target consists of three dielectric cylinders: two of which have a diameter of 31 mm, and a relative permittivity of r = 3 ± 0.3 and the other one has a diameter of 80 mm and a relative permittivity of 1.45 ± 0.15. The geometry of this target can be found in [22]. Herein, we have set the size of the imaging domain to 16 × 16 cm 2 discretized into 70 × 70 square cells. This target was illuminated from 18 different angles, and the resulting scattered field is collected at 241 data points per transmitter. (This procedure has been repeated at nine different operational frequencies ranging from 2 to 10 GHz with 1 GHz increment.) Similar to the procedure outlined in [23], we have created the SP shown in Fig. 1(a) for this target. However, note that the utilized SP in the proposed AMR approach is in the form of preferred edges through |∇ P| 2 as shown in Fig. 1(b). The scale of |∇ P| 2 can be understood by noting that the gradient operator involves division by x and y, and each of which represents the size of a discretized cell within the imaging domain along the x-and y-directions, respectively.
Let us now consider the inversion at the lowest frequency of operation, i.e., 2 GHz. The reason for this intentional choice of this frequency is for the blind 3 inversion using MR-GNI and MR-CSI to not yield sufficient resolution [ Fig. 1(c) and (d)]. (These blind reconstructions have already been shown by several groups including the authors, e.g., [23], and are only presented herein for the sake of side-by-side comparison with the inversion using the SP.) Incorporating the SP through the similarity-based approach using the MRSP-GNI algorithm [14] is shown in Fig. 1(e), which shows a high quantitative accuracy. Incorporating the SP in the form of preferred edges using the AMR-GNI and AMR-CSI algorithms results in the reconstructions as shown in Fig. 1(f) and (g), respectively. As can be seen, the AMR-GNI almost perfectly reconstructs the relative permittivities of the dielectric objects. Besides, the AMR-CSI reconstructed values, although not perfect, are still better than those obtained by MR-CSI. From the numerical experience, the use of AMR seems to work better with the GNI algorithm as compared to the CSI algorithm. By noting the presence of spurious edges in the MR-CSI (blind) reconstruction, see Fig. 1(d), the relatively poor performance of the AMR-CSI algorithm could be due to the common factor of MR-CSI and AMR-CSI, i.e., |∇χ n | 2 . To investigate this, we have tried inverting this data set with AMR-CSI (Type II) which removes |∇χ n | 2 from its regularizer. As shown in Fig. 1(h), the AMR-CSI (Type II) reconstruction is almost perfect, which is aligned with the speculation. 4 In the next step, we only provide partial (incomplete) SP to the inversion algorithms. The partial SP is shown in Fig. 2(a), where it is clear that the SP has no knowledge of the inner dielectric cylinder. The SP in the form of preferred edges is shown in Fig. 2(b), in which the edges associated with the inner cylinder are missing. The reconstruction results of the MRSP-GNI, AMR-GNI, AMR-CSI, and AMR-CSI (Type II) are shown in Fig. 2(c)-(f), respectively. As can be seen, the MRSP-GNI algorithm (a similarity-based approach) was not able to fully retrieve the inner dielectric cylinder as it was misguided by the partial SP. 5 However, both AMR-GNI and AMR-CSI are successful in the detection of the smaller dielectric cylinder within the larger one. Therefore, for this example, the edge-based SP approach outperforms the similarity-based SP approach. Finally, as it was the case with the previous example, the AMR-CSI still suffers from having some spurious edges. To alleviate this, we have also tried AMR-CSI (Type II), whose reconstruction is shown in Fig. 2(f). As can be seen, we can fix the issue of spurious edges; however, due to the absence of |∇χ n | 2 in the denominator of (8), we do not get reconstruction with sharp edges within the larger cylinder, and the reconstructed inner cylinder permittivity undershoots its expected permittivity.

B. Forearm Data Set
This data set [16], [24] was collected from a human forearm at 0.8 GHz using 24 coresident dipole antennas immersed in a matching fluid (salt water) with the relative complex permittivity of about 77 − j 17. The blind inversion of this data set, which was 3 The term blind inversion indicates that no prior information has been used in the inversion process. 4 We have also tried AMR-GNI (Type II); it provided similar reconstruction as in AMR-GNI, but it required tuning Q parameter in AMR (Type II) to be 10. For brevity, this is not shown here. 5 It should be noted that MRSP-GNI is, in fact, able to retrieve the inner cylinder if it is told that the available SP is only partial and is not accurate within the larger cylinder. As discussed in [14], the mechanism with which we can tell the MRSP-GNI algorithm that a given SP is partial is via the use of a probability vector. Herein, we have not used this probability vector. Subsequently, MRSP-GNI favors similar permittivity values within Region II of Fig. 2(a), thus missing the inner cylinder. On the other hand, in [23], the MRSP-GNI algorithm was told that the SP within Region II is not accurate, and thus, we were able to reconstruct the inner cylinder using MRSP-GNI.  originally reported in [16] and [24], is shown in Fig. 3 for the sake of comparison with inversion using SP. 6 (The imaging domain has the size of 10 × 9 cm 2 discretized to 100 × 100 cells.) One item to note in the blind inversions is the fact that although the presence of the two bones is clear, their reconstructed relative permittivities are overestimated. (According to [16], the relative complex permittivities of the bone and muscle are expected to be about 13− j 3 and 56− j 20 at the frequency of operation, respectively.) In the next step, we utilize the image from the MRI of the forearm shown in Fig. 4(a) to create the SP shown in Fig. 4(b) as described in the Appendix. The SP in the form of preferred edges is shown in Fig. 4(c). The reconstruction of the real and imaginary parts of the relative complex permittivity using the SP by the MRSP-GNI, AMR-GNI, and AMR-CSI algorithms is shown in Fig. 4(d)-(i), respectively. (Since the reconstruction using AMR-CSI (Type II) was similar to AMR-CSI reconstruction, it was not shown here for brevity.) As can be seen, the main improvement over the blind inversion is the enhanced quantitative retrieval of the bone's realpart permittivity. Other than this, the rest of the reconstruction does not show any significant improvements over the blind inversion. Furthermore, note that the MRSP-GNI reconstruction outperforms the AMR-GNI and AMR-CSI reconstrctions. This can be attributed to the fact that MRSP-GNI, due to its similarity-based approach, has smoothed out the permittivity values within each region of the imaging domain. Finally, we consider to use the partial and wrong SP shown in Fig. 5(a) and (b). The SP is partial as it misses one bone region, and it is wrong as it mistakenly assumes a square-shaped region [Region IV in Fig. 5(a)]. Furthermore, note the missing edges for one bone region as well as the wrong square edges in Fig. 5(b). The MRSP-GNI, AMR-GNI, and AMR-CSI reconstructions using these SP are shown in Fig. 5(c)-(h), respectively. As can be seen, all of these algorithms can detect the presence of the second bone which was missing in the partial SP. Furthermore, they have all alleviated the presence of the wrong SP regarding the square-shaped region. Note that although the square-shaped region is still visible in the reconstructions, its quantitative value is similar to the muscle region as it should.

V. CONCLUSION
A simple and easy-to-implement approach to modify the existing MWI algorithms employing MR was proposed to incorporate the complete or partial available SPs into the inversion process. This method adds one extra term to the integrand of the standard multiplicative regularizer. Based on the results presented, the incorporation of SPs seems to enhance the quantitative accuracy achievable from MWI.

APPENDIX
Image registration needs to be performed so that we can use the MRI image for inversion. Herein, we have used a preliminary and ad hoc method to extract the SP from the MRI image. First, after applying some thresholding on the MRI image, its resolution was decreased using the MATLAB command interpn to make its discretization the same as that of the microwave inversion. Then, we have compared this postprocessed MRI image with the blind reconstruction shown in Fig. 3. Based on the visual inspection, we have performed some translational movements (i.e., movements along thex-andŷ-directions) on the postprocessed MRI image to merely adjust its overall external contour with respect to that of the blind reconstructions. We have then used this translated postprocessed MRI image as the SP.