# Microwave breast imaging with prior ultrasound information

Yingying Qin, Thomas Rodet, Marc Lambert, Dominique Lesselier

► To cite this version:

Yingying Qin, Thomas Rodet, Marc Lambert, Dominique Lesselier. Microwave breast imaging with prior ultrasound information. IEEE Open Journal of Antennas and Propagation, 2020, 1 (472-482), 10.1109/OJAP.2020.3019953. [hal-02905443](https://HAL Id: hal-02905443)

HAL Id: hal-02905443

https://centralesupelec.hal.science/hal-02905443

Submitted on 17 Sep 2020

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Microwave breast imaging with prior ultrasound information

Yingying Qin, Thomas Rodet, Marc Lambert, and Dominique Lesselier Senior Member IEEE

***Abstract***—Aiming at early detection of tumors, microwave breast imaging is investigated with a priori information on tissue boundaries yielded from ultrasound reflection data. A regularization term is to incorporate the information that two neighboring pixels should exhibit similar dielectric properties when not on a boundary while a jump would be allowed otherwise. This regularization is enforced in the distorted Born iterative method and in the contrast source inversion method. Comprehensive numerical experiments are carried out, involving a simple synthetic model and on two anatomically-realistic MRI-derived numerical breast phantoms. Imaging quality appears greatly improved with this regularization when tissue boundaries are indeed provided. Improvement is observed also when real parts and imaginary parts are retrieved in separate fashion.

***Index Terms***—breast imaging, microwave, ultrasound pri-ors, distorted Born iterative method, contrast source inversion method, tissue boundary regularization

## I. INTRODUCTION

Breast tumors are some of the most common tumors among women. Early detection is critical at an early stage of cancer progression [1]. Therefore, to develop technologies to image a small tumor at low cost and with low risk is an important issue. Currently, X-ray mammography is still the gold standard for this detection. Despite of the high-resolution of imaging result, X-ray mammography has a number of limitations including low sensitivity, ionizing radiation, discomfort from breast compression, and detection quite affected by breast density.

Microwave imaging has been investigated as an alternative or at least a complementary imaging modality. Several investigations on the electromagnetic (EM) properties in different types of tissue have been led, e.g., [2], [3]. Based on the difference in dielectric properties between tumorous and normal tissues, the contrast appears relatively higher than the one associated to X-ray. Also, microwave imaging is non-ionizing, low-cost, and enables easy examination. The techniques proposed for microwave imaging can be roughly divided into two groups: radar and tomography. Radar techniques [4]–[6] rely on ultrawide-band pulse illuminations to identify regions

with high contrast from backscattered signals. This is efficient and indeed indicates the location of strongly scattering parts. However, less detailed information about the breast is made available [7]. To achieve a possibly better retrieval of the distribution of EM parameters within the breast, tomographic microwave imaging is widely used, refer to [8]–[10] among many others. In this approach, several transmitting and receiving antennas are set around the breast. The transmitting ones illuminate it sequentially and scattered fields<fnref n="1" /> are acquired by the receiving ones. Simulation is usually needed to get the scattered field when breast geometry and EM parameters are known, which is a forward problem linear w.r.t. the incident field. To reconstruct the contrast given the scattered field is an inverse problem and it is as well-known nonlinear due to multiple scattering.

This inverse problem can be cast into an optimization problem where the misfit between measurements and simulation results from a numerical model is minimized. Under some conditions, it can be solved without iteration, e.g., a weak scatterer the size of which is no much larger than the wavelength, within the framework of the Born approximation (the total field being replaced by the incident one). In practice, these assumptions are usually not valid and imaging results remain unsatisfactory. This can be overcome by correcting the misfit iteratively. The Born Iterative Method [11], Distorted Iterative Method (DBIM) [12], Contrast Source Inversion (CSI) method [13], and Subspace-based Optimization Method (SOM) and its several variants [14]–[17], work this way. Recently, convolutional neural netwoks (CNN) has been investigated as a tool to solve the inverse scattering problem [18]–[20]. A well-trained network can provide the complex permittivity of the object given the scattered fields or the preliminary results of some traditional inversion algorithms.

It is also well-known that the inverse problem is ill-posed. To alleviate ill-posedness and stabilize the inversion, regularization is usually applied in additive or multiplicative fashion. In DBIM, Tikhonov regularization is standardly enforced to reach a robust estimation. Further, a priori information can be incorporated into the regularization term. Two-fold SOM (TSOM) confines the reconstruction within a low-dimension subspace. Huber regularization [8], [21] and weighted L2-norm total variation (TV) [22] smoothen small-scale noise while trying to preserve discontinuities. For piecewise constant profiles, value picking (VP) [23] regularization can be applied. Level set is also a regularization technique suitable for binary

Correspondance author: dominique.lesselier@12s.centralesupelec.fr

Yingying Qin is with Université Paris-Saclay, ENS Paris-Saclay, CNRS, Systèmes et Applications des Technologies de l'Information et de l'Energie, 94235, Cachan, France, and Université Paris-Saclay, CNRS, CentraleSupélec, Laboratoire des signaux et systèmes, 91190, Gif-sur-Yvette, France.

Thomas Rodet is with Université Paris-Saclay, ENS Paris-Saclay, CNRS, Systèmes et Applications des Technologies de l'Information et de l'Energie, 94235, Cachan, France.

Marc Lambert is with Université Paris-Saclay, CentraleSupélec, CNRS, Laboratoire de Génie Electrique et Electronique de Paris, 91192, Gif-sur-Yvette, France.

Dominique Lesselier is with Université Paris-Saclay, CNRS, Centrale-Supélec, Laboratoire des signaux et systèmes, 91190, Gif-sur-Yvette, France.

<footnote n="1">¹In practice, one should extract them from total fields as collected, and also field values are not observed but antenna-related S-parameters.</footnote>
cases [24] though it now adresses a broad range of cases [25].

The main drawback of microwave imaging is the relative low resolution due to the long wavelength. To achieve higher spatial resolution, a higher frequency is needed [26]. However, the dimension of the scatterer gets correspondingly larger compared to this wavelength and the inverse problem more difficult to solve while the penetration depth may be affected also. Besides, tissue heterogeneity and the fact that geometries and EM parameters vary from person to person may render the challenge quite complicated. Then, additional model-specific information from other imaging modalities can be very useful.

The structural region information can be extracted from images provided by the high-resolution modality as prior information. As an example, in [27], MRI-derived horizontal and vertical boundaries are incorporated into a Bayesian framework for functional image reconstruction. In [28], different tissue clusters are extracted from ultrasound (US) images by K-means to get a better distribution of dielectric properties which results into more accurate tissue-specific time-delays in Delay and Sum algorithm in EM reconstruction. In [29], the structural information is extracted from MRI and Diffuse Optical Tomography (DOT) involving a finite element method (FEM) is considered. In [30], a smaller regularization parameter is distributed in a Tikhonov regularization scheme to pixels identified as part of tumor from X-ray images in DOT breast imaging. In [31], MRI images are segmented into different regions to provide a FEM mesh and a Laplacian-type regularization follows to minimize variation in each region in near infrared (NIR) tomography. In [32], the structural information is extracted from US reconstruction with K-means clustering algorithm. Tissue permittivity values are assigned to these regions to form an inhomogeneous background and assist the EM reconstruction by the FEM-CSI algorithm. In [33], high-resolution images are segmented into different regions and pixels in the same region are constrained to have similar dielectric parameters in EM reconstruction.

Considering US imaging can offer high-resolution images with interior tissue boundaries from reflection algorithms when the travel time of the acoustic signal is recorded and an average sound speed is assumed, being emphasized that US data can be acquired simultaneously with EM ones so that no registration is needed (the hypothesis is the one of a pending breast). US imaging is chosen to offer the additional information in the present work. The tissue boundary information is incorporated into a traditional method for microwave breast imaging with a regularization term, which imposes on two adjacent pixels that the EM properties are the same when not on the boundaries and only undergo changes at interfaces of tissues. Then, one incorporates the US-information-guided regularization term into DBIM and CSI. Besides, CSI with separate constraints on real and imaginary parts is developed. For comparison, results of DBIM with Tikhonov regularization and CSI with Huber regularization are also shown. A synthetic breast model is used to validate the algorithm, then two anatomically-realistic MRI-derived numerical breast phantoms are considered.

The contribution is organized as follows. The forward problem is described in Section II. Regularization as well as inversion algorithms are considered in Section III. Experiments

on breast phantoms are discussed in Section IV. One concludes about present results and outlines ways ahead involving fusion procedures in Section V.

## II. FORWARD PROBLEM

One henceforth considers a two-dimensional non-magnetic case with transverse magnetic (TM) polarization. Time-harmonic waves are assumed with dependence $\exp(-i\omega t)$. The breast is located inside a domain of interest (DoI) $D$. The known background medium is characterized by a complex relative permittivity $\epsilon_b$, permeability $\mu_b$ and wavenumber $k_b = \omega\sqrt{\epsilon_0\epsilon_b\mu_b}$. The unknown scatterers (of same permeability $\mu_b$) have complex relative permittivity $\epsilon_r(\mathbf{r})$ and wavenumber $k(\mathbf{r})$ functions of position. $N_i$ transmitters illuminate the DoI successively and the scattered fields are collected by $N_r$ receivers evenly located on an exterior circle $S$.

The problem can be associated to two integral equations

$$E^t(\mathbf{r}) = E^i(\mathbf{r}) + \int_{D} (k^2(\mathbf{r}') - k_b^2)g(\mathbf{r}, \mathbf{r}') E^t(\mathbf{r}') d\mathbf{r}', r \in D \quad (1)$$

$$E^s(\mathbf{r}) = \int_{D} (k^2(\mathbf{r}') - k_b^2)g(\mathbf{r}, \mathbf{r}') E^t(\mathbf{r}') d\mathbf{r}', r \in S \quad (2)$$

where incident field $E^i(\mathbf{r})$ and total field $E^t(\mathbf{r})$ represent the electric field inside $D$, object absent or present, resp. $E^s$ is the scattered field collected by the receivers. The scalar Green's function is $g(\mathbf{r}, \mathbf{r}') = \frac{i}{4} H_0^{(1)}(k_b|\mathbf{r} - \mathbf{r}'|)$, $H_0^{(1)}$ 1st-kind 0th-order Hankel function. Denote the contrast as

$$\chi(\mathbf{r}) = \frac{k^2(\mathbf{r}) - k_b^2}{k_b^2} \quad (3)$$

and integral operators

$$G_d[x](\mathbf{r}) = k_b^2 \int_D g(\mathbf{r}, \mathbf{r}') x(\mathbf{r}') d\mathbf{r}', r \in D \quad (4)$$

$$G_s[x](\mathbf{r}) = k_b^2 \int_D g(\mathbf{r}, \mathbf{r}') x(\mathbf{r}') d\mathbf{r}', r \in S \quad (5)$$

The equations above simplify into

$$E^t(\mathbf{r}) = E^i(\mathbf{r}) + G_d[\chi E^t](\mathbf{r}), r \in D \quad (6)$$

$$E^s(\mathbf{r}) = G_s[\chi E^t](\mathbf{r}), r \in S \quad (7)$$

To handle the problem numerically, discrete forms of the equations have to be derived, here via a pulse-basis point-matching Method of Moments (MoM). The DoI is discretized into $M = N_x \times N_y$ subwavelength cells with centers at $\mathbf{r}_m$, $m = 1, 2, \dots, M$. The dielectric properties are considered homogeneous in each cell. Every square cell is approximated by a small disk with same area and with equivalent radius $R$ whether needed. The equations above become

$$\mathbf{E}^t = \mathbf{E}^i + \mathbf{G}_d \mathrm{diag}(\chi) \mathbf{E}^t \quad (8)$$

$$\mathbf{E}^s = \mathbf{G}_s \mathrm{diag}(\chi) \mathbf{E}^t \quad (9)$$

In this form, $\chi$ is a $M \times 1$ vector. The $M \times M$ matrix $\mathbf{G}_d$ is

$$\mathbf{G}_d(m, m') = \begin{cases} \frac{ik_b \pi R}{2} J_1(k_b R) H_0^{(1)}(k_b |\mathbf{r}_m - \mathbf{r}_{m'}|), & m \neq m' \\ \frac{ik_b \pi R}{2} H_1^{\{1\}}(k_b R) - 1, & \text{otherwise} \end{cases} \quad (10)$$
where $J_1$ is the 1st-kind Bessel function and $H_1^{(1)}$ the 1st-kind 1st-order Hankel function. The $N_r \times M$ matrix $G_s$ is

$$G_s(s, m) = \frac{ik_b \pi R}{2} J_1(k_b R) H_0^{(1)}(k_b |\mathbf{r}_s - \mathbf{r}_m|) \quad (11)$$

where $\mathbf{r}_s$ is the position of receiver s.

A source-type framework is widely used. Consider the contrast current $J(\mathbf{r}) = \chi(\mathbf{r}) E^t(\mathbf{r})$, the formulation becomes

$$J(\mathbf{r}) = \chi(E^i(\mathbf{r}) + G_d[J](\mathbf{r})), r \in \mathcal{D} \quad (12)$$

$$E^s(\mathbf{r}) = G_s[J](\mathbf{r}), r \in \mathcal{S} \quad (13)$$

The discrete form reads as

$$\mathbf{J} = \mathrm{diag}(\chi)\mathbf{E}^i + \mathrm{diag}(\chi)\mathbf{G}_d\mathbf{J} \quad (14)$$

$$\mathbf{E}^s = \mathbf{G}_s\mathbf{J} \quad (15)$$

### III. INVERSION ALGORITHMS

Due to the large wavelength, resolution in microwave imaging is expected to be poor. Ultrasound imaging with its high resolution is employed to assist it. An ultrasound-guided smoothness (UGS) regularization term is proposed to incorporate a priori information. This regularization is introduced first and then its implementation into DBIM and CSI is shown. DBIM with Tikhonov regularization and CSI with Huber regularization are presented in parallel.

#### A. Smoothness constraint

Assume that interior boundaries of the breast model follow from US imaging. For easier incorporation, the tissue boundaries are depicted via two images to show discontinuities in vertical and horizontal directions, resp., where dielectric properties also change. The principle is quite intuitive, two vertically or horizontally adjacent pixels should exhibit similar parameters when none lies at the boundaries. Otherwise, regularization is imposed. This constraint is expressed as

$$F_{\text{UGS}} = \sum_{i} \sum_{j} b_v(i, j) \|\tilde{\chi}(i, j) - \tilde{\chi}(i, j + 1)\|^2 + \sum_{i} \sum_{j} b_h(i, j) \|\tilde{\chi}(i, j) - \tilde{\chi}(i + 1, j)\|^2 \quad (16)$$

Notice that $\tilde{\chi}$ is a $N_x \times N_y$ matrix and $\chi = \mathrm{vec}(\tilde{\chi})$. $b_v$ and $b_h$ are included in the prior information and indicate if the pixel is at a boundary as

$$b_v(i, j) = \begin{cases} 0, & (i, j) \text{ on horizontal boundaries} \\ 1, & \text{else} \end{cases} \quad (17)$$

$$b_h(i, j) = \begin{cases} 0, & (i, j) \text{ on vertical boundaries} \\ 1, & \text{else} \end{cases} \quad (18)$$

Similarly with other smoothness constraints, the regularization term is based on the module of the gradient of the contrast, coefficients $b_v$ and $b_h$ guiding the smoothness while preserving edges with US information. Compared with other methods of incorporating US prior into EM inversion, the boundary is directly used but not the tissue region from segmentation.

This regularization term can be written in matrix form as

$$F_{\text{UGS}}(\chi) = \|\mathbf{D}_{\text{v}}\chi\|^2 + \|\mathbf{D}_{\text{h}}\chi\|^2 \quad (19)$$

with matrices $\mathbf{D}_v$ and $\mathbf{D}_h$ providing the difference of two adjacent pixels in vertical and horizontal directions. Take $\mathbf{D}_v$ as an example, it can be written as

$$\mathbf{D}_v = \begin{bmatrix} 1 & -1 & & & & & \\ & 1 & -1 & & & & \\ & & \ddots & \ddots & & & \\ & & & 0 & 0 & & \\ & & & & \ddots & \ddots & \\ & & & & & 1 & -1 \\ & & & & & & 1 \end{bmatrix} \quad (20)$$

Notice that the diagonal elements are 1 only when corresponding to pixels not at the boundaries, otherwise the value is zero, no constraint is imposed, and a jump is allowed.

This regularization term is quadratic and incorporated in an additive way, thus it will not introduce additional nonlinearity.

#### B. Distorted Born Iterative Method

Consider the contrast $\chi$ as a small perturbation $\delta\chi$ w.r.t. an inhomogeneous background $\chi_n$; one has $\chi = \chi_n + \delta\chi$ [12]. The secondary incident field by this inhomogeneous background is

$$\mathbf{E}^{\text{bac}} = (\mathbf{I}_{\text{M}} - \mathbf{G}_{\text{d}}\mathrm{diag}(\chi_n))^{-1}\mathbf{E}^i \quad (21)$$

where $\mathbf{I}_{\text{M}}$ is the $M$-dimensional identity matrix, with this background, the scattered field is

$$\mathbf{E}^s = \mathbf{G}_s\mathrm{diag}(\chi_n)\mathbf{E}^{\text{bac}} + \mathbf{G}_s^{\chi_n}\mathrm{diag}(\delta\chi)\mathbf{E}^t \quad (22)$$

Here, $\mathbf{G}_s^{\chi_n}$ is the inhomogeneous background Green's function. The cost functional is the sum of the difference between measured and calculated data and a regularization term

$$F(\delta\chi) = \sum_{p=1}^{N_i} \|E_p^s - \mathbf{G}_s\mathrm{diag}(\chi_n)E_p^{\text{bac}} - \mathbf{G}_s^{\chi_n}\mathrm{diag}(\delta\chi)E_p^{\text{bac}}\|^2 + \alpha F_r(\delta\chi) \quad (23)$$

within which $\mathbf{E}^t$ is replaced by $\mathbf{E}^{\text{bac}}$ with the Born approximation. The optimization procedure is summarized below.

At each iteration, update $\mathbf{G}_s^{\chi_n}$ by

$$\mathbf{G}_s^{\chi_n} = \mathbf{G}_s(\mathbf{I}_{\text{M}} - \mathrm{diag}(\chi_n)\mathbf{G}_d)^{-1} \quad (24)$$

Solve the forward problem

$$E_{p,n}^{\text{bac}} = E_p^i + \mathbf{G}_d\mathrm{diag}(\chi_n)E_{p,n}^{\text{bac}} \quad (25)$$

Calculate the scattered field

$$E_{p,n}^s = \mathbf{G}_s\mathrm{diag}(\chi_n)E_{p,n}^{\text{bac}} \quad (26)$$

Solve the optimization problem with 0th-order Tikhonov regularization

$$\min : F(\delta\chi) = \sum_{p=1}^{N_i} \|E_p^s - E_{p,n}^s - \mathbf{G}_s^{\chi_n}\mathrm{diag}(\delta\chi)E_{p,n}^{\text{bac}}\|^2 + \alpha \|\delta\chi\|^2 \quad (27)$$
or with UGS regularization

$$
\begin{equation}
\begin{split}
\min: F(\delta\chi) = \sum_{p=1}^{N_i} & \|E_p^s - E_{p,n}^s - \mathbf{G}_s^{\chi^n} \mathrm{diag}(\delta\chi) E_{p,n}^{\mathrm{bac}}\|^2 \\
& + \alpha (\|\mathbf{D}_v \delta\chi\|^2 + \|\mathbf{D}_h \delta\chi\|^2)
\end{split}
\tag{28}
\end{equation}
$$

It can be dealt with directly as

$$
\delta\chi = [\mathbf{K}^*\mathbf{K} + \alpha\mathbf{D}]^{-1}\mathbf{K}^*\mathbf{b}, \quad (29)
$$

D = I_M for Tikhonov regularization, D = D_h^*D_h + D_v^*D_v
for the UGS one. K is a (N_i N_r) x M matrix and element
K(i+N_r(p-1), j) = G_s^x(i, j)E_{p,n}^b(j) and b is a (N_i N_r) x 1
vector with b(i + N_r(p-1)) = E_p^s(i) - E_{p,n}^s(i). K^* is the
conjugate transpose of K. To conclude, update the contrast as

$$
\chi_{n+1} = \chi_n + \delta\chi \quad (30)
$$

Considering the imaginary part of the contrast contributes less due to the high difference in magnitude with the real part; consequently, one can separate them and the UGS regularization term becomes

$$
F_{\text{UGS}}(\delta\chi) = F_{\text{UGS}}(\Re\{\delta\chi\}) + \beta F_{\text{UGS}}(\Im\{\delta\chi\}) \quad (31)
$$

A larger regularization parameter (β > 1) can be assigned
to the imaginary part based on the prior information of the
magnitude of real and imaginary parts of contrast. To update
δχ, the linear equation

$$
\begin{bmatrix} \mathbf{K}_d + \alpha_r \mathbf{D} & \mathbf{K}_m \\ -\mathbf{K}_m & \mathbf{K}_d + \alpha_i \mathbf{D} \end{bmatrix} \begin{bmatrix} \Re\{\delta\chi\} \\ \Im\{\delta\chi\} \end{bmatrix} = \begin{bmatrix} y_r \\ y_i \end{bmatrix} \quad (32)
$$

needs to be solved, with $y_r = \Re\{\mathbf{K}^*\mathbf{b}\}$, $y_i = \Im\{\mathbf{K}^*\mathbf{b}\}$,
$\mathbf{K}_d = \mathbf{K}_r^*\mathbf{K}_r + \mathbf{K}_i^*\mathbf{K}_i$, $\mathbf{K}_m = \mathbf{K}_i^*\mathbf{K}_r - \mathbf{K}_r^*\mathbf{K}_i$, $\mathbf{K}_r = \Re\{\mathbf{K}\}$
and $\mathbf{K}_i = \Im\{\mathbf{K}\}$. The regularization parameters are $\alpha_r = \alpha$
and $\alpha_i = \alpha\beta$.

C. Contrast Source Inversion

The contrast source inversion method is based on the source-
type integral equations wherein the contrast source is regarded
as an independent parameter. The cost function is a sum of
normalized mismatches in data and state equations

$$
F(J_1, \dots, J_{N_i}, \chi) = \frac{\sum_{p=1}^{N_i} \|E_p^s - \mathbf{G}_s J_p\|^2}{\sum_{p=1}^{N_i} \|E_p^s\|^2} \\
+ \frac{\sum_{p=1}^{N_i} \|\operatorname{diag}(\chi) E_p^i + \operatorname{diag}(\chi) \mathbf{G}_d J_p - J_p\|^2}{\sum_{p=1}^{N_i} \|E_p^i\|^2} \tag{33}
$$

Notice that in the classical CSI [13] the second term normalization is $\sum_{p=1}^{N_i} \|\mathrm{diag}(\chi)E_p^i\|^2$, one simplifies it as $\sum_{p=1}^{N_i} \|E_p^i\|^2$ [16]. To improve the quality of the reconstruction and incorporate high frequency components of the image, one introduces prior information. Inspired by [8], [21], one incorporates Huber regularization into CSI.

The Huber function can be expressed as

$$
h(x) = \begin{cases} |x|^2, & |x| \le \gamma \\ 2\gamma|x| - \gamma^2, & \text{else} \end{cases} \quad (34)
$$

with $\gamma$ as the threshold. This function is used to estimate the
difference in dielectric properties between the pixel and its
neighborhood. The total difference is measured by

$$
F_{\text{HB}} = \sum_{v} \sum_{v' \in N_v} h(\chi_v - \chi_{v'}) \quad (35)
$$

where $N_v$ represents the neighborhood of $v$. Eight neighbors
are used for one pixel.<fnref n="2" /> The first derivative of Huber regular-
ization w.r.t. the contrast is

$$
g_{\text{HB}}^{\chi} = \sum_{v' \in N_v} \omega_{v'} \quad (36)
$$

with

$$
\omega_{v'} = \begin{cases} \chi_v - \chi_{v'}, & |\chi_v - \chi_{v'}| \le \gamma \\ \gamma(\chi_v - \chi_{v'}) / |\chi_v - \chi_{v'}|, & \text{else} \end{cases} \quad (37)
$$

Another regularization term is the smoothness constraint
guided by US information described in Section III-A. The first
derivative of UGS regularization is

$$
g_{\text{UGS}}^{\chi} = (\mathbf{D}_{\mathbf{v}}^{*} \mathbf{D}_{\mathbf{v}} + \mathbf{D}_{\mathbf{h}}^{*} \mathbf{D}_{\mathbf{h}}) \chi \quad (38)
$$

With these regularizations, the criterion becomes

$$
F(J_1, \dots, J_{N_i}, \chi) = \frac{\sum_{p=1}^{N_i} \|E_p^s - \mathbf{G}_s J_p\|^2}{\sum_{p=1}^{N_i} \|E_p^s\|^2} \\
+ \frac{\sum_{p=1}^{N_i} \|\mathrm{diag}(\chi)E_p^i + \mathrm{diag}(\chi)\mathbf{G}_d J_p - J_p\|^2 + \alpha F_r(\chi)}{\sum_{p=1}^{N_i} \|E_p^i\|^2} \tag{39}
$$

There is no need to solve the forward problem at each iter-
ation, so it is more efficient. With function F(J₁, ..., J<sub>N<sub>i</sub></sub>, χ),
it is difficult to optimize contrast current and contrast simul-
taneously due to their dependence. Here one follows [34] to
optimize them alternately by minimizing the cost functional
w.r.t. **J**, assuming known χ, and update χ, assuming known
**J**. The optimization procedure is in Algorithm.

Similarly, the regularization term can also be enforced on
the real and imaginary parts separately as

$$
F_r(\chi) = F_r(\Re\{\chi\}) + \beta F_r(\Im\{\chi\}) \quad (40)
$$

During update of the contrast, the gradients of the cost
functional w.r.t. to real and imaginary parts are

$$
\begin{align}
g_n^{\Re\{\chi\}} &= \sum_p \Re\{(\chi_{n-1} E_{p,n}^t - J_{p,n})(E_{p,n}^t)^*\} \nonumber \\
&\quad + \alpha_r g_{F_r}^{\Re\{\chi\}} (\Re\{\chi\}) \tag{41} \\
g_n^{\Im\{\chi\}} &= \sum_p \Im\{(\chi_{n-1} E_{p,n}^t - J_{p,n})(E_{p,n}^t)^*\} \nonumber \\
&\quad + \alpha_i g_{F_r}^{\Im\{\chi\}} (\Im\{\chi\}) \nonumber
\end{align}
$$

with $\alpha_r = \alpha$ and $\alpha_i = \alpha\beta$.

<footnote n="2"><sup>2</sup>In [8], there is a coefficient 0.5 since the difference between two neigh-
 boring pixels is calculated twice. Here it is in the regularization parameter.</footnote>
**Algorithm** CSI with HB/UGS regularization

**Input:** $E^s$, $E^i$, $G_s$, $G_d$

**Output:** $\chi$

Initialize $\chi_0$ and $J$ by back propagation [35] with

$$J_{p,0} = \frac{\|G_s^* E_p^s\|^2}{\|G_s G_s^* E_p^s\|^2} G_s^* E_p^s \text{ and } \chi_0 = \frac{\sum_{p=1}^{N_i} J_{p,0} E_{p,0}^*}{\sum_{p=1}^{N_i} \|E_{p,0}\|^2}$$

**for** $n = 1 : iter_{max}$ **do**

Calculate data equation error

$$\rho_{p,n-1} = E_p^s - G_s J_{p,n-1}$$

Calculate state equation error

$$r_{p,n-1} = \operatorname{diag}(\chi_{n-1}) E_{p,n-1}^t - J_{p,n-1}$$

Update contrast current:

$$\text{Gradient } g_{p,n}^J = -\frac{G_s^* \rho_{p,n-1} / \sum_p \|E_p^s\|^2}{r_{p,n-1} - G_d^* \operatorname{diag}(\chi_{n-1})^* r_{p,n-1}} - \frac{\sum_{p=1}^{N_i} \|E_p^i\|^2}{\sum_{p=1}^{N_i} \|E_p^i\|^2}$$

Polak-Ribière conjugate gradient search direction

$$v_{p,n} = g_{p,n}^J + \frac{\Re\langle g_{p,n}^J, g_{p,n}^J - g_{p,n-1}^J \rangle}{\langle g_{p,n-1}^J, g_{p,n-1}^J \rangle} v_{p,n-1}$$

Update $J$ by $J_{p,n} = J_{p,n-1} + \beta_{p,n} v_{p,n}$

$$\text{Update total field by } E_{p,n}^t = E_p^i + G_d J_{p,n}$$

Update contrast:

Gradient

$$g_n^\chi = \sum_p (\chi_{n-1} E_{p,n}^t - J_{p,n})(E_{p,n}^t)^* + \alpha g_{\text{HBF/UGS}}^\chi$$

Polak-Ribière conjugate gradient search direction

$$d_n = g_n^\chi + \frac{\Re\langle g_n^\chi, g_{n-1}^\chi - g_{n-1}^\chi \rangle}{\langle g_{n-1}^\chi, g_{n-1}^\chi \rangle} d_{n-1}$$

Contrast

$$\chi_n = \chi_{n-1} + \alpha_n d_n$$

**end for**

<figure><img src="image_2.png" /><figcaption>Fig. 1: Real (a) and imaginary (b) parts of synthetic Model 1.</figcaption></figure>

## IV. NUMERICAL SIMULATIONS

In this section, one considers numerical experiments on three breast models: a synthetic one, Model 1, to validate algorithms, and two realistic breast phantoms, Models 2 and 3. As mentioned, in US imaging, reflective boundaries where acoustic impedance changes can be detected by reflection algorithms when the travel time $t$ of acoustic signals is recorded and an average sound speed $c$ assumed. Here, tissue boundary information is not derived from carrying out these methods in full but from simpler simulations: three sources are put around the object; each one emits ultrasonic pulses in different directions (from -45° to 45° towards the center); the position of the boundary along each direction is calculated as $d = ct/2$; based on the angle of acoustic wave it is decided whether the boundary is vertical or horizontal and the images subsequently follow. The experiments are conducted at a single frequency of 1 GHz.

### A. Reconstruction of synthetic breast model

This simple synthetic breast model consists of skin, fatty, fibroconnective/glandular and tumorous tissues. Each tissue type is of regular shape and has uniform permittivity. The skin is 2 mm thick and the tumor of 6 mm-diameter. The immersion medium is chosen from [8], letting $\epsilon_b = 10 + 4i$. Relative

permittivities of tissues at 1 GHz are 39.89 + 15.64i, 4.80 + 0.82i, 48.82 + 15.89i, 56.27 + 17.96i, resp.. Figure 1 depicts the distribution of real and imaginary parts. The DoI is of size 0.1 m × 0.1 m. 40 antennas are distributed evenly on a circle of 0.057 m radius, operated as sources and receivers simultaneously. For the forward problem, the DoI is discretized into 200 × 200 pixels and the problem tackled by a conjugate-gradient fast Fourier transform (CG-FFT) algorithm. Additive Gaussian noise is added to the synthetic data with SNR = 30 dB.

For inversion, the domain is discretized into 80 × 80 pixels. One assumes that the breast is in a circle with 0.045 m radius (the radius of outer boundary of the breast model is 0.042 m). The reconstruction is confined within the region and pixels outside it restricted to same dielectric properties as the coupling medium. Retrieved real and imaginary parts of the relative permittivity are bounded by $1 \le \Re\{\epsilon_r\} \le 70$ and $0 \le \Im\{\epsilon_r\} \le 40$ at each iteration.

The regularization parameter $\alpha$ can be chosen by L-curve [36], generalized cross-validation [37], or other methods. When it is determined, it can be chosen as $\alpha_r$ in separate inversion. $\beta$, as the ratio of regularization parameter of imaginary part to real part, is chosen based on the relative magnitude of these two parts. $\gamma$ in Huber regularization is the threshold below which a quadratic cost is used to smooth noise and above which a linearly-varying cost is added to penalize it in a lesser extent for edge-preserving. It is selected as the smallest difference in the contrast (a margin can be left).

First, one experiments on DBIM with Tikhonov regularization (DBIM-TK) and CSI with separate Huber regularization (CSI-HB-S). The DBIM regularization parameter is $\alpha = 0.05$. With Huber, the threshold is $\gamma = 0.5$ and the regularization parameter $\alpha_i = 0.1$ for the imaginary part, $\alpha_r = 0.001$ for the real part. CSI with no regularization is run also in comparison and its results shown.

Results are in Figure 2. The glandular part is reconstructed as a ring with relatively higher relative permittivity but is still recognized. The contrast value of the center part is usually underestimated. Dielectric properties are not smooth in each region and the small tumor is not found. Now, one incorporates the prior information from US imaging by UGS regularization. The tissue boundary information incorporated is in Figure 3. As explained, these two images show the discontinuities in vertical and horizontal directions.

The regularization parameter is set to 1 in DBIM, 0.01 in CSI. Besides, one separates real and imaginary parts in DBIM-
<figure><img src="image_3.png" /><figcaption>Fig. 2: Retrieved real (left) and imaginary (right) parts of Model 1 with. (a),(b) DBIM-TK reconstruction; (c),(d) CSI reconstruction; (e),(f) CSI-HB-S reconstruction.</figcaption></figure>

<figure><img src="image_4.png" /><figcaption>Fig. 3: Boundary information in vertical (a) and horizontal (b) directions.</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>Fig. 4: Retrieved real (left) and imaginary (right) parts of Model 1. (a),(b) DBIM-UGS; (c),(d) DBIM-UGS-S; (e),(f) CSI-UGS; (g),(h) CSI-UGS-S.</figcaption></figure>

<table><caption>TABLE I: Relative error with Model 1 with SNR = 30 dB</caption><thead><tr><th>Method</th><th>DBIM-TK</th><th>DBIM-UGS</th><th>CSI</th><th>CSI-HB-S</th><th>CSI-UGS</th><th>CSI-UGS-S</th></tr></thead><tbody><tr><td>Err(ℝ{ε<sub>r</sub>})</td><td>0.385</td><td>0.335</td><td>0.370</td><td>0.345</td><td>0.342</td><td>0.322</td></tr><tr><td>Err(ℝ{ε<sub>r</sub>})</td><td>0.674</td><td>0.483</td><td>0.596</td><td>0.389</td><td>0.469</td><td>0.385</td></tr><tr><td>Err(ε<sub>r</sub>)</td><td>0.4245</td><td>0.353</td><td>0.400</td><td>0.350</td><td>0.358</td><td>0.329</td></tr></tbody></table>

UGS (DBIM-UGS-S) and CSI-UGS algorithm (CSI-UGS-S). The regularization parameters are $α_r = 1$ and $α_i = 10$ for DBIM-UGS-S and $α_r = 0.01$ and $α_i = 10$ for CSI-UGS-S. Results are in Figure 4.

With this US information guided regularization, the results are smoother, while edges are well preserved with both algorithms. Besides, the tumor is well distinguished within the glandular part in the reconstruction of the real part. Upon separation of real and imaginary part reconstructions and assigning a large regularization parameter to the imaginary part, one sees obvious improvement in its retrieval.

To evaluate the imaging results quantitatively, one computes the relative error of the permittivity as

$$Err(\zeta) = \frac{\|\zeta_{est} - \zeta_{true}\|_2}{\|\zeta_{true}\|_2} \quad (42)$$

where $ζ$ is the parameter to evaluate and subscripts "true" and "est" represent true value and estimation result, resp. The errors with the algorithms above are in Table I. Reconstructions are more accurate than those with the algorithms without prior information. Also, one sees a small decrease in the error result of the real part and a large one in the imaginary part when the two parts of contrast are retrieved separately.

The CPU time is also compared between the cases whether the real and imaginary parts are reconstructed separately. In average, it takes 13.89 seconds for one iteration in DBIM-UGS,14.38 in DBIM-UGS-S, 1.79 in CSI-UGS, and 0.73 in CSI-UGS-S, performed on Intel Core i7-8850H CPU (2.60 GHz) with 32 GByte memory. Considering both computation cost and imaging quality, DBIM-UGS and CSI-UGS-S are
<figure><img src="image_6.png" /><figcaption>Fig. 5: Inexact boundaries tested in experiments with Model 1: (a),(b) fake tumor; (c),(d) incomplete boundary.</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>Fig. 6: Retrieved real (left) and imaginary (right) parts of Model 1. (a),(b) case with a fake tumor; (c),(d) case with incomplete boundary.</figcaption></figure>

henceforth used to incorporate US information in the following experiments.

Robustness of the algorithm is tested based on DBIM-UGS. Here the reflective boundaries are assumed to have been obtained from some reflection mode imaging methods, however, in fact, those cannot be perfect. The acoustic signal may bounce between two interfaces so there will be artifacts and some boundaries may be missing due to a small variation in acoustic impedance. So, it is worthwhile to discuss the effects of inexact boundaries.

One has carried out several experiments with two kinds of inexact boundaries, which are shown in Figure 5. In the first kind, there is a fake tumor, and in the second one, interfaces between tumor and glandular part are incomplete. With these boundaries, retrievals are depicted in Figure 6. Since the regularization term is intended to suppress the discontinuity

<table><caption>TABLE II: Relative error with Model 1 with SNR = 10 dB</caption><thead><tr><th>Method</th><th>DBIM-TK</th><th>DBIM-UGS</th><th>CSI-HB-S</th><th>CSI-UGS-S</th></tr></thead><tbody><tr><td>Err(ℝ{ε<sub>r</sub>})</td><td>0.436</td><td>0.351</td><td>0.352</td><td>0.325</td></tr><tr><td>Err(ℝ<sup>3</sup>{ε<sub>r</sub>})</td><td>0.857</td><td>0.584</td><td>0.4015</td><td>0.396</td></tr><tr><td>Err(ε<sub>r</sub>)</td><td>0.497</td><td>0.3825</td><td>0.358</td><td>0.333</td></tr></tbody></table>

<figure><img src="image_8.png" /><figcaption>Fig. 7: Real (a) and imaginary (b) parts of Model 2.</figcaption></figure>

between two adjacent pixels not at a border, it has no constraint otherwise. That is more severe when boundaries are not complete. This is consistent with the experiments. If a fake tumor in the US information, it appears in the map with lower contrast, while the real one still has higher permittivity and one can detect it. When the boundaries are incomplete, the tumor has no evident border and is difficult to find. Thus, if the boundary information cannot yield the interface between tumor and other tissues, it is hard to image it with microwaves.

Now, a higher level of noise is considered with SNR = 10 dB. Experiments are conducted on DBIM-TK, DBIM-UGS, CSI-HB-S and CSI-UGS-S with the same regularization parameters. The relative errors are summarized in Table II. Due to the high level of noise, the quality of the imaging is degraded for all algorithms yet the one proposed is affected to a lesser extent.

## B. Reconstruction of realistic breast phantoms

To better validate the approach, one tests it on two more realistic breast phantoms as slices extracted from breast phantoms of the UWCEM repository [38], categorized in different classes according to radiographic density.

1) *Class 2 breast phantom*: Model 2 (ID 070604) has a scattered fibroglandular density. Debye parameters are from [9]. Figure 7 shows the phantom.

The cell size is 0.5 mm and one uses this grid to solve the forward problem. 40 antennas are set evenly on a circle of 0.08 m radius. Additive Gaussian noise of 30 dB is added to the data. For inversion, one adopts a 2 mm cell size, resulting in 88 × 77 pixels. Other configurations are the same as with experiments before. The reconstruction results of DBIM-TK and CSI-HB-S are in Figure 8. The regularization parameter is $α = 0.05$ for DBIM and one sets $α_r = 0.001$ and $α_i = 0.1$ in CSI for Huber regularization, its threshold being $γ = 0.5$.

The main glandular part is well retrieved by both algorithms, yet with a smaller size. Both fail in imaging the fine structure in the phantom. Yet, imaging is satisfactory as one can at least see the main structure.

Figure 9 shows the vertical and horizontal boundaries for the experiments next. As observed, discontinuities occur more
<figure><img src="image_9.png" /><figcaption>Fig. 8: Retrieved real (left) and imaginary (right) parts of Model 2. (a),(b) DBIM-TK; (c),(d) CSI-HB-S.</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>Fig. 9: Tissue boundaries of Model 2.</figcaption></figure>

<table><caption>TABLE III: Relative error with Model 2</caption><thead><tr><th>Method</th><th>DBIM-TK</th><th>DBIM-UGS</th><th>CSI-HB-S</th><th>CSI-UGS-S</th></tr></thead><tbody><tr><td>Err(ℝ{ε<sub>r</sub>})</td><td>0.465</td><td>0.385</td><td>0.443</td><td>0.383</td></tr><tr><td>Err(ℝ{ε<sub>r</sub>})</td><td>0.591</td><td>0.447</td><td>0.473</td><td>0.455</td></tr><tr><td>Err(ε<sub>r</sub>)</td><td>0.481</td><td>0.393</td><td>0.447</td><td>0.392</td></tr></tbody></table>

frequently in the breast phantom and several points presenting discontinuities are adjacent in the same direction with null corresponding diagonal elements in $D_h^*D_h + D_v^*D_v$. As a result, in DBIM, even with regularization term $F_{UGS} = α (||D_vδχ||^2 + ||D_hδχ||^2)$, the problem may remain ill-conditioned.

This can be tackled by adding a small penalty term on the points where discontinuity shows in prior information. One sets $b_{v/h}(i, j) = γ_t$ when the point is on a boundary and $D_v/h$ also changes accordingly. $γ_t$ is chosen as a small value meanwhile keeping matrix $K^*K + α(D_v^*D_v + D_h^*D_h)$ well-conditioned. Regularization parameters are chosen as $α = 10$ and $γ_t = 0.1$ for UGS regularization in DBIM-UGS. In CSI-UGS-S, they are set as $α_r = 0.01$ and $α_i = 1$ for the real and imaginary parts.

Results of DBIM-UGS and CSI-UGS-S are in Figure 10. The results improve and finer structures are imaged when the boundary information is incorporated. The relative errors of results above are proposed in Table III. From this error calculation, one observes that imaging quality is indeed enhanced when US high-resolution information is incorporated.

2) Class 3 breast phantom: Model 3 (ID 080304) is heterogeneously dense and its structure is quite complicated, see Figure 11. In Model 3, one inserts a synthetic 1 cm-diameter tumor at position (1.2 cm, -0.5 cm). Its relative permittivity is 59.98 + 19.83i, corresponding to the 75th percentile curve

at 1 GHz in [3]. Configurations are as with Model 2. Cell size remains 2 mm for inversion and the DoI has 83 × 51 pixels.

One chooses $α = 0.001$ for the Tikhonov regularization, $α_r = 0.001$ and $α_i = 0.1$ for the Huber regularization. Figure 12 displays the reconstruction results. Those by CSI-HB-S are smoother than those by DBIM-TK since the latter enforces no constraint on the spatial gradient of the contrast. Neither algorithm however can image the glandular part well and both fail in detecting the tumor.

Boundaries now used are shown in Figure 13. For DBIM, one sets $γ_t = 0.2$ and the other parameter values are like in the Model 2 experiment.

<figure><img src="image_11.png" /><figcaption>Fig. 11: Real (a) and imaginary (b) parts of Model 3.</figcaption></figure>
<figure><img src="image_12.png" /><figcaption>Fig. 12: Retrieved real (left) and imaginary (right) parts of Model 3. (a),(b) DBIM-TK; (c),(d) CSI-HB-S.</figcaption></figure>

<figure><img src="image_13.png" /><figcaption>Fig. 13: Tissue boundaries of Model 3</figcaption></figure>

<figure><img src="image_14.png" /><figcaption>Fig. 14: Retrieved real (left) and imaginary (right) parts of Model 3. (a),(b) DBIM-UGS; (c),(d) CSI-UGS-S.</figcaption></figure>

When boundary information is incorporated, results improve. Glandular and tumorous tissues are retrieved well. The small tumor can be distinguished from the glandular part, i.e., can be detected. Yet, singular points emerge with higher value. The regularization parameter $\gamma_t$ must be properly chosen: too small, one may see those singular pixels, too large, the result may be too smooth to detect the tumor. The relative error of Model 3 is in Table IV, showing improvement of imaging.

<table><caption>TABLE IV: Relative error with Model 3</caption><thead><tr><th>Method</th><th>DBIM-TK</th><th>DBIM-UGS</th><th>CSI-HB-S</th><th>CSI-UGS-S</th></tr></thead><tbody><tr><td>Err(ℝ{ε<sub>r</sub>})</td><td>0.486</td><td>0.427</td><td>0.4675</td><td>0.413</td></tr><tr><td>Err(ℝ<sup>3</sup>{ε<sub>r</sub>})</td><td>0.683</td><td>0.515</td><td>0.511</td><td>0.537</td></tr><tr><td>Err(ε<sub>r</sub>)</td><td>0.510</td><td>0.437</td><td>0.472</td><td>0.4275</td></tr></tbody></table>

## V. CONCLUSION

In the present work, one has proposed a regularization term to incorporate US information into microwave imaging. The regularizer tries to smoothen the dielectric properties between two adjacent pixels when none is on the tissue boundaries indicated by the US information, thus there is no need to segment US images to decide whether or not two pixels belong to the same tissue region, only tissue boundaries where reflection occurs in US imaging are used. Thus, it is easier to implement. One concludes that with microwave data only, it is quite difficult to image a small tumor inside the breast.

Adding high resolution information enables to well estimate its location and shape. One has also tested the robustness of the algorithm. Since the algorithm is imposed on pixels not on boundaries, this has more influence on the result when some boundaries are missing. If the interface of a small tumor and normal tissue is not detected by US imaging, it is hard to find it by microwave imaging. One has also shown that to retrieve real and imaginary parts in separate fashion can improve the results. Now, further attention should be on presentation of the boundary information from US imaging and attempting to find a way for a better exploitation.

As for joint (fused) inversion of electromagnetic and acoustic data for breast imaging, it is of good potential. Edge-preserving regularization [39] can be performed by introducing auxiliary variables indicating whether or not a pixel is on an edge. Edge markers could be obtained from the last parameter profile and guide the next optimization as regularization term. Alternate minimization would be used to update acoustic contrast, edge markers and dielectric contrast. Also, increasingly popular convolutional neural networks (CNN), now involving a two-stream CNN [40], combining feature maps at a certain level, should produce tissue types, the last CNN layer being a classifier (this is expected to be easier than outputting EM and US parameters at each pixel).

## ACKNOWLEDGEMENTS

The authors acknowledge the careful reading of Peipei Ran and her contribution to outlining the paths ahead.

## REFERENCES

[1] L. Tabar, M.-F. Yen, B. Vitak, H.-H. T. Chen, R. A. Smith, and S. W. Duffy, "Mammography service screening and mortality in breast cancer patients: 20-year follow-up before and after introduction of screening," *The Lancet*, vol. 361, no. 9367, pp. 1405-1410, 2003.

[2] M. Lazebnik, L. McCartney, D. Popovic, C. B. Watkins, M. J. Lindstrom, J. Harter, S. Sewall, A. Magliocco, J. H. Booske, M. Okoniewski, et al., "A large-scale study of the ultrawideband microwave dielectric properties of normal breast tissue obtained from reduction surgeries," *Physics in Medicine & Biology*, vol. 52, no. 10, p. 2637, 2007.

[3] M. Lazebnik, D. Popovic, L. McCartney, C. B. Watkins, M. J. Lindstrom, J. Harter, S. Sewall, T. Ogilvie, A. Magliocco, T. M. Breslin, et al., "A large-scale study of the ultrawideband microwave dielectric properties of normal, benign and malignant breast tissues obtained from cancer surgeries," *Physics in Medicine & Biology*, vol. 52, no. 20, p. 6093, 2007.

[4] S. K. Davis, B. D. Van Veen, S. C. Hagness, and F. Kelcz, "Breast tumor characterization based on ultrawideband microwave backscatter," *IEEE Transactions on Biomedical Engineering*, vol. 55, no. 1, pp. 237-246, 2007.
[5] W. C. Khor, M. E. Bialkowski, A. Abbosh, N. Seman, and S. Crozier, "An ultra wideband microwave imaging system for breast cancer detection," *IEICE Transactions on Communications*, vol. 90, no. 9, pp. 2376-2381, 2007.

[6] M. Klemm, J. A. Leendertz, D. Gibbins, I. J. Craddock, A. Preece, and R. Benjamin, "Microwave radar-based breast cancer detection: Imaging in inhomogeneous breast phantoms," *IEEE Antennas Wireless and Propagation Letters*, vol. 8, pp. 1349-1352, 2009.

[7] E. C. Fear, "Microwave imaging of the breast," *Technology in Cancer Research & Treatment*, vol. 4, no. 1, pp. 69-82, 2005.

[8] F. Bai, A. Franchois, and A. Pizurica, "3D microwave tomography with Huber regularization applied to realistic numerical breast phantoms," *Progress In Electromagnetics Research*, vol. 155, pp. 75-91, 2016.

[9] J. D. Shea, P. Kosmas, S. C. Hagness, and B. D. Van Veen, "Three-dimensional microwave imaging of realistic numerical breast phantoms via a multiple-frequency inverse scattering technique," *Medical Physics*, vol. 37, no. 8, pp. 4210-4226, 2010.

[10] S. H. Zainud-Deen, W. M. Hassen, E. Ali, K. H. Awadalla, and H. Sharshar, "Breast cancer detection using a hybrid finite difference frequency domain and particle swarm optimization techniques," in *2008 National Radio Science Conference*, pp. 1-8, IEEE, 2008.

[11] Y. Wang and W. C. Chew, "An iterative solution of the two-dimensional electromagnetic inverse scattering problem," *International Journal of Imaging Systems and Technology*, vol. 1, no. 1, pp. 100-108, 1989.

[12] W. C. Chew and Y.-M. Wang, "Reconstruction of two-dimensional permittivity distribution using the distorted Born iterative method," *IEEE Transactions on Medical Imaging*, vol. 9, no. 2, pp. 218-225, 1990.

[13] P. M. Van Den Berg and R. E. Kleinman, "A contrast source inversion method," *Inverse Problems*, vol. 13, no. 6, p. 1607, 1997.

[14] X. Chen, "Subspace-based optimization method for solving inverse-scattering problems," *IEEE Transactions on Geoscience and Remote Sensing*, vol. 48, no. 1, pp. 42-49, 2009.

[15] Y. Zhong and X. Chen, "Twofold subspace-based optimization method for solving inverse scattering problems," *Inverse Problems*, vol. 25, no. 8, p. 085003, 2009.

[16] Y. Zhong and X. Chen, "An FFT twofold subspace-based optimization method for solving electromagnetic inverse scattering problems," *IEEE Transactions on Antennas and Propagation*, vol. 59, no. 3, pp. 914-927, 2011.

[17] Y. Zhong, M. Lambert, D. Lesselier, and X. Chen, "A new integral equation method to solve highly nonlinear inverse scattering problems," *IEEE Transactions on Antennas and Propagation*, vol. 64, no. 5, pp. 1788-1799, 2016.

[18] Z. Wei and X. Chen, "Deep-learning schemes for full-wave nonlinear inverse scattering problems," *IEEE Transactions on Geoscience and Remote Sensing*, vol. 57, no. 4.

[19] "Deepnis: Deep neural network for nonlinear electromagnetic inverse scattering,"

[20] H. M. Yao, E. I. Wei, and L. Jiang, "Two-step enhanced deep learning approach for electromagnetic inverse scattering problems," *IEEE Antennas and Wireless Propagation Letters*, vol. 18, no. 11, pp. 2254-2258, 2019.

[21] F. Bai, A. Piżurica, A. Franchois, S. Van Loocke, D. De Zutter, and W. Philips, "Weakly convex discontinuity adaptive regularization for microwave imaging," *IEEE Transactions on Antennas and Propagation*, vol. 61, no. 12, pp. 6242-6246, 2013.

[22] P. Mojabi and J. LoVetri, "Microwave biomedical imaging using the multiplicative regularized Gauss-Newton inversion," *IEEE Antennas and Wireless Propagation Letters*, vol. 8, pp. 645-648, 2009.

[23] J. De Zaeytijd, A. Franchois, and J.-M. Geffrin, "A new value picking regularization strategy-application to the 3-D electromagnetic inverse scattering problem," *IEEE Transactions on Antennas and Propagation*, vol. 57, no. 4, pp. 1133-1149, 2009.

[24] A. Litman, D. Lesselier, and F. Santosa, "Reconstruction of a two-dimensional binary obstacle by controlled evolution of a level-set," *Inverse Problems*, vol. 14, no. 3, p. 685, 1998.

[25] O. Dorn and D. Lesselier, "Level set methods for inverse scattering," *Inverse Problems*, vol. 22, no. 4, p. R67, 2006.

[26] M. Slaney, A. C. Kak, and L. E. Larsen, "Limitations of imaging with first-order diffraction tomography," *IEEE Transactions on Microwave Theory and Techniques*, vol. 32, no. 8, pp. 860-874, 1984.

[27] G. Gindi, M. Lee, A. Rangarajan, and I. G. Zubal, "Bayesian reconstruction of functional images using anatomical information as priors," *IEEE Transactions on Medical Imaging*, vol. 12, no. 4, pp. 670-680, 1993.

[28] M. Omer, P. Mojabi, D. Kurrant, J. LoVetri, and E. Fear, "Proof-of-concept of the incorporation of ultrasound-derived structural information into microwave radar imaging," *IEEE Journal on Multiscale and Multiphysics Computational Techniques*, vol. 3, pp. 129-139, 2018.

[29] G. Gulsen, O. Birgul, M. B. Unlu, R. Shafiiha, and O. Nalcioğlu, "Combined diffuse optical tomography (DOT) and MRI system for cancer imaging in small animals," *Technology in Cancer Research & Treatment*, vol. 5, no. 4, pp. 351-363, 2006.

[30] A. Li, E. L. Miller, M. E. Kilmer, T. J. Brukilacchio, T. Chaves, J. Stott, Q. Zhang, T. Wu, M. Chorlton, R. H. Moore, et al., "Tomographic optical breast imaging guided by three-dimensional mammography," *Applied Optics*, vol. 42, no. 25, pp. 5181-5190, 2003.

[31] B. A. Brooksby, S. Jiang, H. Dehghani, B. W. Pogue, K. D. Paulsen, J. B. Weaver, C. Kogel, and S. P. Poplack, "Combining near-infrared tomography and magnetic resonance imaging to study in vivo breast tissue: implementation of a Laplacian-type regularization to incorporate magnetic resonance structure," *Journal of Biomedical Optics*, vol. 10, no. 5, p. 051504, 2005.

[32] N. Abdollahi, D. Kurrant, P. Mojabi, M. Omer, E. Fear, and J. LoVetri, "Incorporation of ultrasonic prior information for improving quantitative microwave imaging of breast," *IEEE Journal on Multiscale and Multiphysics Computational Techniques*, vol. 4, pp. 98-110, 2019.

[33] L. M. Neira, B. D. Van Veen, and S. C. Hagness, "High-resolution microwave breast imaging using a 3-d inverse scattering algorithm with a variable-strength spatial prior constraint," *IEEE Transactions on Antennas and Propagation*, vol. 65, no. 11, pp. 6002-6014, 2017.

[34] P. M. van den Berg, A. Van Broekhoven, and A. Abubakar, "Extended contrast source inversion," *Inverse Problems*, vol. 15, no. 5, p. 1325, 1999.

[35] K. Belkebir, P. C. Chaumet, and A. Sentenac, "Superresolution in total internal reflection tomography," *JOSA A*, vol. 22, no. 9, pp. 1889-1897, 2005.

[36] P. C. Hansen and D. P. O'Leary, "The use of the 1-curve in the regularization of discrete ill-posed problems," *SIAM Journal on Scientific Computing*, vol. 14, no. 6, pp. 1487-1503, 1993.

[37] G. H. Golub, M. Heath, and G. Wahba, "Generalized cross-validation as a method for choosing a good ridge parameter," *Technometrics*, vol. 21, no. 2, pp. 215-223, 1979.

[38] E. Zastrow, S. Davis, M. Lazebnik, F. Kelcz, B. Van Veen, and S. Hagness, "Database of 3D grid-based numerical breast phantoms for use in computational electromagnetics simulations," tech. rep., Department of Electrical and Computer Engineering University of Wisconsin-Madison, 2008.

[39] P. Charbonnier, L. Blanc-Feraud, G. Aubert, and M. Barlaud, "Deterministic edge-preserving regularization in computed imaging," *IEEE Transactions on Image Processing*, vol. 6, no. 2, pp. 289-311, 1997.

[40] K. Simonyan and Z. Andrew, "Two-stream convolutional networks for action recognition in videos," in *Proceedings of the 27th International Conference on Neural Information Processing Systems (NIPS'14)*, vol. 1, pp. 568-576, 2014.