*Article* **Incoherent Radar Imaging for Breast Cancer Detection and Experimental Validation against 3D Multimodal Breast Phantoms**

**Antonio Cuccaro 1, Angela Dell'Aversano 1, Giuseppe Ruvio 2,3, Jacinta Browne 4,5 and Raffaele Solimene 1,6,7,\***

	- <sup>3</sup> Endowave Ltd., Galway 8, Ireland

**Abstract:** In this paper we consider radar approaches for breast cancer detection. The aim is to give a brief review of the main features of incoherent methods, based on beam-forming and Multiple SIgnal Classification (MUSIC) algorithms, that we have recently developed, and to compare them with classical coherent beam-forming. Those methods have the remarkable advantage of not requiring antenna characterization/compensation, which can be problematic in view of the close (to the breast) proximity set-up usually employed in breast imaging. Moreover, we proceed to an experimental validation of one of the incoherent methods, i.e., the I-MUSIC, using the multimodal breast phantom we have previously developed. While in a previous paper we focused on the phantom manufacture and characterization, here we are mainly concerned with providing the detail of the reconstruction algorithm, in particular for a new multi-step clutter rejection method that was employed and only barely described. In this regard, this contribution can be considered as a completion of our previous study. The experiments against the phantom show promising results and highlight the crucial role played by the clutter rejection procedure.

**Keywords:** microwave imaging; incoherent imaging; clutter rejection; breast cancer detection

## **1. Introduction**

Global statistics have demonstrated that breast cancer is the most frequently diagnosed invasive cancer and the leading cause of death due to cancer among female patients [1]. In recent years the incidence of breast cancer in developed countries has continued to rise; but at same time, the rate of mortality has undergone a substantial decline [2]. This is due to the improvements in medical cancer treatment and in the implementation of screening programs as well as to improved imaging techniques [3].

As shown in [4], and as can be naturally expected, the survival of a patient is strongly determined by the stage of the disease at the time the treatment starts. Therefore, early diagnostics is crucial. This requires further improvements of the capabilities of current diagnostic modalities. In addition, over the last few years, this has steered efforts towards the development of new imaging modalities with the aim of supplementing the ones currently employed in the clinical practice.

Current conventional imaging modalities are X-ray mammography, digital breast tomo-synthesis, ultrasound and magnetic resonance imaging (MRI), with mammography actually being the golden standard in breast cancer imaging [5]. Among new imaging

**Citation:** Cuccaro, A.; Dell'Aversano, A.; Ruvio, G.; Browne, J.; Solimene, R. Incoherent Radar Imaging for Breast Cancer Detection and Experimental Validation against 3D Multimodal Breast Phantoms. *J. Imaging* **2021**, *7*, 23. https://doi.org/10.3390/ jimaging7020023

Academic Editors: Leonardo Rundo, Carmelo Militello, Vincenzo Conti, Fulvio Zaccagna and Changhee Han Received: 30 December 2020 Accepted: 25 January 2021 Published: 1 February 2021

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2021 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

modalities, in this paper we focus on microwave breast imaging (MBI). Microwave imaging has triggered a great deal of research over the last decades because it offers a number of potential advantages related to the use of non-ionizing radiation, it does not require to compress the breast and requires a relatively cheap technology [6]. All these features along with the progress achieved in this field [7], show that MBI is actually a "promising imaging modality" [8,9].

Many algorithms for microwave imaging have been tailored for breast diagnostics [10]. Some of them reconstruct a 3D volume; others are based on a sliced approach and, for example, they reconstruct repeated coronal slices of the breast and thus reduce the imaging algorithm complexity and accelerate image reformatting [11]. In any case, microwave breast imaging entails solving a non-linear ill-posed inverse scattering problem since diffraction effects cannot be ignored as in X-ray tomography.

Microwave imaging algorithms can be coarsely grouped in two broad categories, depending on the way non-linearity is dealt with.

When the aim is to reconstruct the dielectric/conductivity profile of the breast tissue under examination, "quantitative" algorithms must be adopted. In these cases, the nonlinearity of the problem must be taken into account and the reconstructions are basically achieved by iterative optimization procedures that try to minimize some cost function of the misfit between the available data and the model ones [12]. As such, non-linear inversions are generally computationally very intensive [13] and can suffer from convergence and reliability problems due to false solutions [14]. However, we need to mention that some hybrid approaches, that exploit a priory information provided by other modality, can help mitigate these issues.

The imaging problem is drastically simplified if the imaging method is based on linearized scattering models [15]. In this case, imaging results into linear procedures that are robust and computationally effective. However, only qualitative information can be obtained. Indeed, the corresponding images are more like hot maps where strong inhomogeneities are highlighted. Therefore, linear methods can be conveniently employed if the main objective is to detect and localize targets with a significant dielectric contrast as compared to the surrounding background tissues.

Linear imaging methods are commonly addressed as radar approaches and are the ones we are concerned with in this contribution.

Among the linear methods, beam-forming (BF) is probably the most popular in MBI. Basically, it consists in time-shifting the signals received over the measurement aperture in order to isolate signals scattered from (and hence to focus at) a particular synthetic focal point belonging to the spatial area to be imaged [16]. The BF approach is attractive for the excellent compromise between the achievable performance and the procedure complexity. In [17] the classical delay and sum (DAS) beam-former is used for breast cancer imaging, but many different versions of DAS beam-former have been employed and proposed in literature. For example the delay multiply and sum (DMAS) beam-former is proposed in [18] and the enhanced DAS (EDAS) beam-former in [19]. Besides BF, many other linear inversion methods can be found in literature. For example, a number of linear inversion methods that rely on the spatial spectral representation of the solutions of the wave equation have been developed in different applicative contexts [20]. Among them we mention the range migration [21], the Stolt migration [22], the wave-interference migration [23] and the Holographic Imaging (HI) [24]. These methods are very appealing since their implementation requires computing a Fourier transformation [25,26] which can be effectively achieved via a Fast Fourier Transform (FFT) algorithm. While the latter typically requires the scattered field data to be collected over a planar (rectilinear for 2D cases) measurement aperture, since the Cartesian spatial coordinates naturally match the spatial Fourier transform setting, the extension to deal with circular configurations (more suitable for breast imaging) was previously pursued, for example, in [27].

A detailed analytical comparison between beam-forming and holographic methods has been carried out (for a scattering scenario pertinent to breast imaging) in [28], where the role of critical parameters, such as the operating frequency range, the number of scatterers and data discretization, was considered. Instead, in [29] it is shown that all these methods are variants of the so-called generalized holography.

As remarked above, linear methods restrict the imaging stage to a mere detection and localization of strong in-homogeneities. However, even in light of this reduced task and under the simplified linear framework, the imaging problem is still extremely difficult for a number of reasons. One of the issue, is that all the previous methods "coherently" combine data collected at different frequencies. Therefore, the achievable performance is negatively affected by frequency dispersion of breast tissues (which are unknown or known with a considerable degree of uncertainty and vary from patient to patient) as well as by the antenna frequency response, which is hard to predict because it is in close proximity to breast. As shown in [30,31], this drawback can be mitigated by employing non-coherent imaging strategies. In particular in those papers we introduced and compared incoherent versions of beam-forming and MUSIC [32] (I-MUSIC) and showed that the performance remains stable by using different types of antennas although they were non-characterized, i.e., their frequency responses were not estimated nor enclosed in the model upon which the algorithms were based.

Another crucial aspect is the clutter that generally overwhelms the relatively weak signal coming from the cancer targets. Accordingly, before obtaining the image, data must be first processed in order to reduce the clutter due to the antenna's internal reflection, the skin interface and other non-tumor breast tissues. To this end, we employed a hybrid clutter removal method [33].

In this contribution we will give a quick review of the incoherent methods and detail the clutter mitigation procedure that was used, but actually not described, in [33]. In particular, the achievable cancer detection is checked by using experimental data collected by employing the multi-modal breast phantom developed in [33]. Accordingly, this paper focuses more on the image process and can be considered as a companion paper of [33] which, instead, mainly considered the development of the breast phantom and its tissue characterization for different imaging modalities.

## **2. Ideal Scattering Configuration and Beam-Forming**

In order to introduce the notation and to more easily describe the incoherent imaging methods we consider an idealized scattering scenario. More specifically, the scattering problem is considered for a two-dimensional scalar configuration (see Figure 1). Here, invariance is assumed along the *z*-axis and the electromagnetic incident field has a transverse magnetic TM polarization.

According to the measurement set-up commonly used for breast imaging, sensors are assumed to be located over a circle which is in close proximity and embodies the scattering region *D*. The position of each sensor is identified by the vector *ro* ∈ Γ*o*, Γ*<sup>o</sup>* being the circular measurement curve. The scattered signals are assumed to be collected only at the same position as the transmitter, while the latter can assume different positions over the circle in order to synthesize the measurement aperture. Hence, a multimonostatic configuration is considered. Note that, while the multimonostatic setting is by far the most common, more complex multiview/multistatic sensors' arrangement are often employed: these configurations are not considered herein.

As far as the background medium is concerned, it is assumed to be homogeneous and lossless. Of course, this does not match with realistic breast structures which consist of many layers of different materials that can have even articulated boundaries. Nonetheless, because the breast structure is actually unknown, and for the sake of simplicity, during the image stage an equivalent homogeneous background medium is usually considered.

According to the previous assumptions and under Born approximation [15], the scattered field in the frequency domain is given as:

$$E\_S(\omega, \underline{r}\_\sigma) = \left(\frac{\omega}{2\upsilon}\right)^2 P(\omega) \int\_D G^2(\omega \,/v, \underline{r}\_{0'} \underline{r}) \chi(\underline{r}) d\underline{r} \tag{1}$$

where *D* is the spatial region under investigation, *ω* is the angular frequency and *v* the background propagation speed. Moreover, *ES*(*ω*,*ro*) is the scattered field data, *P*(*ω*) is the temporal Fourier spectrum of the transmitted pulse, *G*(.) = 1/4*jH*(2) <sup>0</sup> (·) is the twodimensional scalar background Green function, with *H*(2) <sup>0</sup> (·) being the Hankel function of second kind and zero order. Finally, *χ*(*r*) is the so-called contrast function which describes the scatterers in terms of their shape and electromagnetic parameters. Note that in general *χ*(*r*) is also frequency dependent but here such a dependence has been neglected. In particular, exploiting the asymptotic expansion of the Hankel function, i.e., *H*(2) <sup>0</sup> (*x*) <sup>√</sup>2/*π<sup>x</sup>* exp [−*j*(*<sup>x</sup>* <sup>−</sup> *<sup>π</sup>*/4)]), Equation (1) can be recast as:

$$E\_S(\omega, \underline{r}\_o) = (j\omega / 2\pi v) P(\omega) \int\_D \frac{\exp\left(\frac{-j2\omega}{v} |\underline{r}\_o - \underline{r}|\right)}{|\underline{r}\_o - \underline{r}|} \chi(\underline{r}) d\underline{r} \tag{2}$$

**Figure 1.** Pictorial view of the scattering scene. Invariance is assumed along the *z*-axis.

In practice, the quantity that is actually measured is not the scattered field but rather the system scattering parameters. This basically entails taking in account the antenna frequency response. Accordingly, (2) modifies as:

$$S(\omega, \underline{\mathfrak{r}}\_o) = (j\omega/2\pi\upsilon)P(\omega) \int\_D \frac{\exp\left(\frac{-j2\omega}{\upsilon}|\underline{\mathfrak{r}}\_o - \underline{\mathfrak{r}}|\right)}{|\underline{\mathfrak{r}}\_o - \underline{\mathfrak{r}}|} \chi(\underline{\mathfrak{r}})d\underline{\mathfrak{r}} \tag{3}$$

where *P*˜(*ω*) = *H*2(*ω*)*P*(*ω*) now takes into account the squared (because the antenna acts as TX and RX) antenna response assumed to be solely dependent on the frequency *ω* and *S*(*ω*,*ro*) are the scattering measurements.

In order to introduce the beam-forming method, it is convenient to consider the time domain version of Equation (3). Hence, by Fourier transforming with respect to *ω*, we obtain the time domain scattering measurements as:

$$s(t, \underline{\mathbf{x}}\_{o}) = \int\_{D} \tilde{p}[t - \pi(\underline{\mathbf{x}}\_{o}, \underline{\mathbf{x}})] \chi(\underline{\mathbf{x}}) d\underline{\mathbf{x}} \tag{4}$$

where *p*˜(.) is related to the transmitted pulse and is the Fourier transform of (*jω*/2)*P*˜(*ω*)/ |*r*<sup>0</sup> − *r*| and *τ*(*ro*,*r*) = 2/*v*|*r*<sup>0</sup> − *r*| is the round-trip delay.

Generally, the image obtained by the DAS beam-forming is given by:

$$I\_{BF}(\underline{\mathbf{r}}) = \int \mathcal{W}(t) \left[ \int\_{\Gamma\_{\mathcal{S}}} \mathbf{s}[t - T(\underline{\mathbf{r}}\_{\mathcal{U}}, \underline{\mathbf{r}}), \underline{\mathbf{r}}\_{\mathcal{O}}] d\underline{\mathbf{r}}\_{\mathcal{O}} \right]^2 dt \tag{5}$$

where *W*(.) is a suitable time window and *T*(*ro*,*r*) = *TW* − *τ*(*ro*,*r*), with *TW* = *maxro*,*r*{*τ*(*ro*, *r*)}. Accordingly, the received signals are "aligned" at the time instant *TW* and then summed. In particular, by setting *W*(*t*) = *δ*(*t* − *TW*) the reconstruction *IBF* becomes:

$$I\_{BF}(\underline{r}) = \left[ \int\_{\Gamma\_{\vartheta}} s[T\_W - T(\underline{r}\_{o'} \underline{r}), \underline{r}\_o] d\underline{r}\_o \right]^2 \tag{6}$$

and returning back to the frequency domain (details can be found in [28]):

$$I\_{BF}(\underline{r}) = \left| \int\_{\Omega} \int\_{\underline{\Gamma}\_{\rho}} S(\omega\_{\prime} \underline{r}\_{o}) \exp\left[j\omega \tau(\underline{r}\_{o}, \underline{r})\right] d\underline{r}\_{o} d\omega \right|^{2} \tag{7}$$

Equation (7) is functional to appreciate the difference with the incoherent approach to be shown in the sequel. In addition, it allowed a closed-form derivation of the point-spread function, that is the reconstruction of a point-like target *χ*(*rp*) = *δ*(*r* − *rp*), that permits to evaluate the achievable resolution in terms of the configuration parameters, including the frequency range and the data discretization [28]. In particular, it was shown that the common belief that in order to achieve a finer resolution a wider frequency band is required does not necessarily hold true. Indeed, while this statement is correct for aspect-limited configurations, for the case at hand, where measurements can be taken all around the scattering region (i.e., non-aspect limited configuration, see Figure 1), finer resolution can be obtained by moving a fixed frequency band towards high frequencies. This is an important result which has practical implications since it promotes the use of cheaper hardware and simplifies the antenna design, which does not necessarily have to work on an ultra-wide band. All details can be found in [28].

## **3. Incoherent Image Procedures**

As highlighted in (5), the measured scattering parameters depend on the antenna response. Indeed, this enters in shaping the frequency behaviour of *P*˜(*ω*) and in general introduces a frequency dependent propagation delay. The latter must be considered while setting the time window *W*(*t*) and the alignment time *Tw*. This requires near-field antenna characterization/equalization that can be pursued by a suitable set of measurements or numerical simulations. However, as the breast properties change from patient to patient, residual errors still remain. With uncertainty levels as high as the magnitude of the tumor scattered field, the imaging procedure's robustness is dramatically endangered. This is particularly true for dense breasts as they present lower tumor/healthy-tissue contrast. It can be noted that this drawback arises because frequency data are coherently summed. Therefore, a viable way to mitigate this problem is to devise imaging schemes which do rely on such a coherence and process each frequency data separately. This is the topic addressed in this section.

## *3.1. Incoherent Beam-Forming*

Basically, incoherent beam-forming (IBF) is achieved as follows:

$$I\_{IRF}(\underline{\mathbf{r}}) = \int\_{\Omega} \left| \int\_{\Gamma\_{\sigma}} S(\omega, \underline{\mathbf{r}}\_{o}) \exp\left[j\omega \tau(\underline{\mathbf{r}}\_{o}, \underline{\mathbf{r}})\right] d\underline{r}\_{o} \right|^{2} d\omega \tag{8}$$

where the basic difference, with respect to (5), is clearly that data are summed in amplitude along the frequency domain. Of course, it is interesting to elucidate how (8) relates physically (meaning) and in terms of the achievable performance to (5). As shown in [28], the time domain counterpart of (8) is

$$I\_{IBF}(\underline{r}) = \int \left[ \int\_{\Gamma\_{\sigma}} \mathbf{s} [t - T(\underline{r}\_{\sigma}, \underline{r}), \underline{r}\_{o}] d\underline{r}\_{o} \right]^2 dt \tag{9}$$

where basically the window function has been removed. Form the achievable performance point of view, in [28] the point-spread function was also analytically derived for the incoherent case and it was found that the main difference with respect to the coherent case is that side-lobes are slightly higher. However, the point-spread function main beams (and hence the resolution) are practically the same. Therefore, the cost to pay while using (8) in place of (5) is that side-lobe reconstruction increases a little bit (of course, the actual increase depends on the configuration parameters, especially by the frequency band) but this is largely rewarded since the need to estimate/compensate the antenna response is avoided.

## *3.2. Discrete Data Setting*

In the previous section we implicitly considered the situation where data are collected continuously all around the scattering scene. In practice, the number of data samples must be finite. Accordingly, in this section we recast the previous argument within a discrete data setting which, in turn, is also necessary to introduce the I-MUSIC, as shown below.

Therefore, say *ro*1,*ro*2, ··· ,*roNo No* measurement points taken uniformly over the measurement circle Γ*<sup>o</sup>* and *ω*1, *ω*2, ··· , *ωN f* the employed frequencies. In addition, denote as *r*1,*r*2, ··· ,*rNs* the coordinates of the pixels that divide the spatial region under test *D*. The finite dimensional (discrete) counterpart of (1) can then be written as:

$$\mathbf{S}^{\bar{i}}(\omega\_{\bar{i}}) = (\omega\_{\bar{i}}/2v)^2 \mathcal{P}(\omega\_{\bar{i}}) \mathbf{A}(\omega\_{\bar{i}}) \mathbf{b}(\omega\_{\bar{i}}) \tag{10}$$

where

$$\mathbf{S}^{l}(\omega\_{i}) = [\mathbf{S}(\underline{\mathbf{z}}\_{\sigma1}, \omega\_{i}), \mathbf{S}(\underline{\mathbf{z}}\_{\sigma2}, \omega\_{i}), \dots, \mathbf{S}(\underline{\mathbf{z}}\_{\sigma N\_{0}}, \omega\_{i})]^{\mathrm{T}} \in \mathbb{C}^{N\_{0}} \tag{11}$$

is the column vector of the scattering data collected at frequency *ωi*,

$$\mathbf{b}^{l}(\omega\_{i}) = [b\_{1}(\omega\_{i}), b\_{2}(\omega\_{i}), \dots, b\_{N\_{s}}(\omega\_{i})]^{T} \in \mathbb{C}^{N\_{s}} \tag{12}$$

is the vector of the pixel scattering coefficients), (·)*<sup>T</sup>* denoting the transpose, and **<sup>A</sup>**(*ωi*) ∈ *<sup>C</sup>No*×*Ns* is the *No* × *Ns* matrix propagator (indeed a discrete version of Equation (1)) whose *n*-th column has the form:

$$\mathbf{A}^{n}(\underline{\mathbf{r}}\_{n},\omega\_{i}) = \left[\mathbf{G}^{2}(\omega\_{i},\underline{\mathbf{r}}\_{\omega1},\underline{\mathbf{r}}\_{n}), \mathbf{G}^{2}(\omega\_{i},\underline{\mathbf{r}}\_{\omega2},\underline{\mathbf{r}}\_{n}), \dots, \mathbf{G}^{2}(\omega\_{i},\underline{\mathbf{r}}\_{\omega N\_{\vartheta}},\underline{\mathbf{r}}\_{n})\right]^{T} \tag{13}$$

where the Green function is the same as in Equation (2). Accordingly, the overall data scattering matrix is:

$$\mathbf{S} = [\mathbf{S}^1(\omega\_1)\,\mathbf{S}^2(\omega\_2)\,\cdots \,\mathbf{S}^{N\_f}(\omega\_{N\_f})] \in \mathbb{C}^{N\_0 \times N\_f} \tag{14}$$

Due to this discrete setting, Equation (8) can be particularized as:

$$I\_{IRF}(\underline{r}) = \sum\_{m=1}^{N\_f} \left| \sum\_{l=1}^{N\_o} S(\omega\_{m\prime} \underline{r}\_{ol}) \exp\left[j\omega\_m \underline{r}(\underline{r}\_{on\prime} \underline{r})\right] \right|^2 \tag{15}$$

where *r* ∈ *r*1,*r*2, ··· ,*rNs* . A crucial question to be addressed within the discrete setting is the choice of the minimum number of sensor positions that should be deployed around the scattering scene in order to obtain the same results as the ideal case (i.e., data collected continuously) or at least to avoid aliasing effects that can result in reconstruction crowed by spurious artifacts that can be mistaken for actual targets. In particular, it is shown that to avoid aliasing a sufficient condition is that the number of measurement points be:

$$N\_0 \ge 4k\_{\max} R\_c \tag{16}$$

where *kmax* is the wave number corresponding to the highest adopted frequencies and *Rc* < *R* the radius of the circular investigation domain. Basically, Equation (16) guarantees that data are properly "spatially" sampled for each adopted frequency. However, because of the multifrequency data, and the related mutifrequency reconstruction process, some degree of under-sampling can be tolerable for part of the frequency band. This is because aliasing spurious artifacts are frequency dependent. Thus, their positions change with the frequency. By contrast, the main contribution of the reconstruction always peaks at the actual scatterer's location. Therefore, even if condition (16) is not satisfied, while summing up different frequency contributions in (15), artifacts tend to be averaged out whereas the main beam (due to scatterer) is not.

## *3.3. Incoherent MUSIC*

The starting point is the construction of the correlation matrix for each frequency, that is:

$$\mathbf{R}(\omega\_i) = \mathbf{S}^i(\omega\_i)\mathbf{S}^{iH}(\omega\_i) = \mathbf{A}(\omega\_i)\mathbf{B}(\omega\_i)\mathbf{A}^H(\omega\_i) \tag{17}$$

where **b***H*(*ωi*) and **A***H*(*ωi*) are the Hermitian vector and matrix of **b**(*ωi*) and **A**(*ωi*), respectively, and **B**(*ωi*) = **b**(*ωi*)**b***H*(*ωi*). According to [32], scatterers can be localized by finding the steering vectors which are orthogonal to the so called noise subspace. This requires computing the eigenspectrum of **R**(*ωi*) and the steering vectors which basically consists of the normalized columns of the propagator **<sup>A</sup>**(*ωi*), that is **Sv**(*rn*, *<sup>ω</sup>i*) = **<sup>A</sup>***n*(*rn*, *<sup>ω</sup>i*)/**A***n*(*rn*, *<sup>ω</sup>i*) being evaluated in correspondence to the trial position *rn* within the spatial domain *D*. Hence, scatterers' positions are identified where the pseudospectrum

$$\phi(\underline{\mathbf{r}}\_{n'}\omega\_{i}) = \frac{1}{||\mathcal{P}\_{\mathcal{N}}[\mathbf{S}\mathbf{v}^{n}(\omega\_{i})]||^{2}}\tag{18}$$

peaks, with PN [·] being the projection operator onto the noise subspace. However, as shown in [31,34], the correlation matrix is rank deficiency with rank one. Therefore, the scheme to identify the scatterers' location can be modified defining PN = (I−PS ), with PS [·] being the projector onto the signal space. Hence, the detection is achievable by adopting the only significant singular vector associated to the signal subspace.

Note that Equation (18) refers to single frequency data. Multiple frequencies can be incoherently combined [35] giving rise to

$$\Phi\_{I-MLSIC}(\underline{r}\_{\rm u}) = \prod\_{i=1}^{N\_f} \phi(\underline{r}\_{\rm u}, \omega\_i) \tag{19}$$

Eventually, (19) is the proposed as algorithm for cancer detection.

## *3.4. Numerical Comparison*

In this section a numerical comparison between the I-MUSIC and the beam-forming strategies is shown. Initially, single frequency (*f* = 3 GHz) data are considered and a background medium with *<sup>r</sup>* = 9. The investigation domain *D* is assumed to be a circle of radius *Rc* = 6 cm whereas measurements are taken over a concentric circle of radius *R* slightly greater than *Rc*. A point-like target is located in the centre of *D*. For the case at hand, Equation (16) suggests *No* > 45 to avoid artefacts. The reconstructions corresponding to this case are shown in panels (a) and (b) of Figure 2. Note that at single frequency there is no difference between BF and IBF. Accordingly, in Figure 2 we just refer to beam-forming. When the number of points is lowered, spurious artefacts actually corrupt the reconstructions. In particular, the bottom panels of the same figure depict the

reconstructions when the number of sensors is reduced by seven times. These results are perfectly consistent with the theory developed in [28,34].

**Figure 2.** Comparing I-MUSIC and beam-forming (BF) for single frequency data. The left column refers to I-MUSIC; the right one to BF. In panels (**a**,**b**) *No* = 49 whereas in panels (**c**,**d**) *No* = 7.

As mentioned above, frequencies are a good ally to mitigate artifacts when the data are under-sampled. This can be observed in Figure 3 where the three reconstruction schemes, i.e., I-MUSIC, BF and IBF, are compared for the same cases as in panels (c) and (d) of Figure 2 but by considering two different frequency bands. As can be seen, the frequency band greatly helps in reducing artefacts. In addition, as expected and according to previous discussion, IBF presents higher side-lobe levels than BF. Furthermore, I-MUSIC outperforms BF schemes since it allows for a more sharper target localization and better resilience to aliasing. Therefore, I-MUSIC is the method that has been selected for undergoing the experimental validation reported in the sequel.

**Figure 3.** Illustrating the role of the frequency band. In all the reconstructions only *No* = 7 sensors are considered. In the top panels the frequency band is [1, 3] GHz, in the bottom panels the frequency band is reduced to [2, 3] GHz. Finally, (**a**,**d**) refer to I-MUSIC, (**b**,**e**) to BF and (**c**,**f**) to IBF.

## **4. Experimental Analysis**

As mentioned in the introduction, this paper can be regarded as a companion paper of [33]. In that paper, we mainly focused on the design, construction and characterization of the breast phantom; microwave imaging algorithms were not described at all. While the detection algorithm was actually the I-MUSIC that we have already described in previous contributions (and whose main ingredients have been briefly recalled above in conjunction to the comparison with more classical BF methods) the clutter rejection algorithm deserves a more in-depth description. Therefore, the program for this section is to first briefly report about the measurement set-up and the breast phantom and then to move on to a detailed description of the clutter rejection method. Finally, a few experimental reconstructions are used to show the effectiveness of the I-MUSIC + de-cluttering procedure.

## *4.1. Measurement Set-Up*

The pictorial view of the measurement set-up is shown in Figure 4 and basically coincides with the measurement scheme adopted in [36]. A breast phantom was scanned by an antipodal Vivaldi antenna in the frequency range [0.5–5] GHz connected to a VNA. In particular, at a given height the antenna rotated around the phantom (with a 5◦ angular step) in order to synthesize a multimonostatic configuration (i.e., *TX* and *RX* were colocated) for a total 72 scanning positions. In general, data collected at different heights can be simultaneously employed to get a 3D reconstruction. However, here we exploited the sliced approach. The phantom and the antenna were immersed within a coupling medium with relative dielectric permittivity equal to 12. This was done for antenna miniaturization purposes and to reduce the dielectric discontinuity from the antenna side to the breast, which can hinder microwave energy penetration [35]. Accordingly, such a value of the dielectric permittivity was used to define the equivalent homogeneous reference background medium which was used to build up the scattering model upon which the detection algorithm was based. No information concerning the phantom nor the antenna response (which was not estimated or compensated) was exploited in the following image stage.

**Figure 4.** Schematic diagram showing the MBI scanning setup. The system antenna + phantom is immersed in a coupling medium. The antenna is connected to a Vector Network Analyzer (VNA) scanning the phantom at a fixed height in multimonostatic configuration. This allows collecting data for a single coronal slice.

## *4.2. Breast Phantom*

In [33] we developed multimodal anthropomorphic breast phantoms suitable for evaluating the imaging performance of microwave imaging in comparison to the established diagnostic imaging modalities of Magnetic Resonance Imaging, Ultrasound, Mammography and Computed Tomography. In that study, the aim was to build a bridge between the numerical simulation environment and a more realistic diagnostic scenario. To this end, the constructed anthropomorphic phantoms mimic breast tissues in terms of their heterogeneity, anatomy, morphology, and mechanical and dielectric characteristics and reproduce different healthy and pathologic tissue types for each of the modalities, taking into consideration the differing imaging and contrast mechanisms for each modality. In that study, two phantoms were developed: the phantom (named as 'Phantom A') had a simple and less morphologically accurate interface between mammary fat and fibroglandular tissue; the second ('Phantom B') had a more relevant complex fat and fibroglandular interface. Both were extracted from real patient MRI datasets. Apart from the different morphological structure, the phantoms had the same five different tissue-mimicking materials: skin, subcutaneous fat, fibroglandular tissue, pectoral muscle and tumor. The phantoms' construction used non-toxic materials, and they were inexpensive and relatively easy to manufacture. Both phantoms were characterized and scanned using conventional modalities (MRI, US, mammography and CT). The details concerning all the steps required for their manufacturing, characterization and imaging can be found in [33]. Their MRI

coronal slices are reported below and highlight the different tissue mimicking layers as well as the tumor. In particular, it can be seen that in both cases the tumor was located inside the fibroglandural structure. The tissue dielectric permittivity and the conductivity of the different tissues are reported in Figure 1 of [33] which shows that the dielectric contrast between tumor and fibroglandural tissue was at best 1.2:1 (hence extremely low) within the frequency band [0.5, 4] GHz.

## *4.3. Clutter Rejection Algorithm*

In order to obtain the reconstruction, in our case a single coronal slice, before the imaging stage, data had to be properly processed in order to reduce clutter disturbances that arose from the antenna internal reflections, the skin interface, and other breast tissues. As the clutter tends to mask the informative signal, it needed to be reduced before the image construction procedure was run. Different clutter rejection methods have been proposed in the literature. For example, in [36] a hybrid artefact removal algorithm for microwave imaging is used, while in [37] some of the most common algorithms used in Through Wall Imaging (TWI) applications were compared, including the simplest average trace subtraction strategy. In this paper, a new method for "extracting" the useful signal is proposed. It was based on a two-step entropy computation and a subspace projection stage. The first entropy step was used to set a time-gating procedure in order to remove the strong antenna's internal reflections and skin contribution; the second one was instead used to select the subset of sensors' positions where tumor contribution is stronger. Finally, a subspace projection procedure [38] was aimed at mitigating contributions due to breast inhomogeneities. After mitigating the clutter, I-MUSIC was employed to obtain the image.

For convenience we rewrote the scattering data matrix as follows:

$$\mathbf{S} = \begin{bmatrix} \mathbf{S}\_1 \\ \mathbf{S}\_2 \\ \vdots \\ \mathbf{S}\_{N\_g} \end{bmatrix} \tag{20}$$

where **S***<sup>i</sup>* are the rows of **S** and hence vectors whose indexes range over the frequencies, i.e., **<sup>S</sup>***<sup>i</sup>* <sup>∈</sup> *<sup>C</sup>Nf* . In order to compute the time gate, the first step is to transform the rows of **<sup>S</sup>** in time domain. Accordingly, upon applying an IDFT routine, we get:

$$\mathbf{s} = \begin{bmatrix} \mathbf{s}\_1 \\ \mathbf{s}\_2 \\ \vdots \\ \mathbf{s}\_{N\_0} \end{bmatrix} \tag{21}$$

with **<sup>s</sup>** ∈ *<sup>C</sup>No*×*Nt* being the time domain version of **<sup>S</sup>** and *Nt* is the number of retained time domain samples. Hence, the rows of **s** are basically the time-traces (A-scan in usual radar literature) collected over the different antenna positions. Note that if data had already been collected in time domain, this step would have not be necessary.

In Reference [39], an entropy-based metric was used to discriminate between clutter and target signals. The same idea was adopted here to seek a suitable time-gating. Accordingly, normalized time traces, **s**˜*n*, were constructed whose entries are given by:

$$\tilde{s}\_{\rm tr}(t\_{\rm m}) = \frac{|s\_{\rm tr}(t\_{\rm m})|}{\sum\_{l=1}^{N\_{\rm o}} |s\_{l}(t\_{\rm m})|} \quad \forall t\_{\rm m} = t\_{1\prime} \cdot \cdots \cdot t\_{N\_{\rm l}} \tag{22}$$

where *sn*(*tm*) is just the m-th entry of **s***n*, i.e., the n-th time trace. Now, *s*˜*n*(*tm*) ≥ 0 and ∑*No <sup>n</sup>*=<sup>1</sup> *s*˜*n*(*tm*) = 1, ∀*tm*. Therefore, for each instant of time, the vector of the normalized data could be assimilated to a probability density function. This observation suggested introducing the entropy measure as

$$\mathfrak{e}\_{\mathfrak{s}}(t\_{\mathfrak{m}}) = -\sum\_{n=1}^{N\_{\mathfrak{o}}} \mathfrak{s}\_{\mathfrak{n}}(t\_{\mathfrak{m}}) \log[\mathfrak{s}\_{\mathfrak{n}}(t\_{\mathfrak{m}})] \tag{2.3}$$

It was expected that *<sup>s</sup>* was high for those instants of time for which the received signals were similar across the different spatial acquisitions. Of course, this occurred when the antenna was receiving its internal reflections or the skin contribution. Figure 5 (left panel) shows a typical entropy behavior obtained for the collected data. As can be observed, *<sup>s</sup>* was nearly constant and high until the time *tmmin* (marked by the dashed red circle), where the entropy attained its first abrupt change compared to its maximum value log(*No*). According to previous discussion, signals coming from the phantom should start to be received only beyond *tmmin* . Signals before such an instant must be discarded. This can be enforced by adopting a time-gating with a time-windowing that removes signals for *tm* < *tmmin* , that is:

$$s\_{\mathcal{W}n}(t\_m) = \mathcal{W}(t\_m)s\_n(t\_m) \tag{24}$$

with

$$\mathcal{W}(t\_{\rm mf}) = \begin{cases} \begin{array}{c} 0 \quad \text{if } t\_{\rm mf} \le \; t\_{m\_{\rm min}} \\ 1 \quad \text{elsewhere} \end{array} \tag{25}$$

After time gating, the scattering data matrix was denoted as:

$$\mathbf{s}\_{W} = \begin{bmatrix} \mathbf{s}\_{W1} \\ \mathbf{s}\_{W2} \\ \vdots \\ \mathbf{s}\_{WN\_o} \end{bmatrix} \tag{26}$$

and a further entropy-based windowing was applied. More in detail, **s***Wn* were further normalized as follows:

$$\mathfrak{s}\_{\mathsf{Wm}}(t\_i) = \frac{|s\_{\mathsf{Wm}}(t\_i)|}{\sum\_{I=1}^{N\_{\mathsf{I}}} |s\_{\mathsf{Wm}}(t\_I)|} \tag{27}$$

**Figure 5.** Illustrating entropy behaviour for the case of phantom B. The (**left**) panel shows the entropy *s*(*tm*) and the red dashed circle identifies the time-gating value (3 ns). The (**right**) panel shows ˆ*s*, the orange and green shaded regions highlights the set of sensors' positions whose data can be retained.

Note that now each time trace underwent a different normalization with the normalizing factor being provided by the summation of the magnitude of its time samples. Once again, it follows that *<sup>s</sup>*ˆ*Wn*(*ti*) > 0 and <sup>∑</sup>*Nt <sup>i</sup>*=<sup>1</sup> *s*ˆ*Wn*(*ti*) = 1, for each sensor's position. Hence, the *Nt*-dimensional vectors of the normalized windowed A-scans could be assimilated as

above to a probability density function and the corresponding entropy can be computed as follows:

$$\mathfrak{E}\_s(\mathfrak{r}\_{on}) = -\sum\_{l=1}^{N\_l} \mathfrak{s}\_{Wn}(t\_l) \log[\mathfrak{s}\_{Wn}(t\_l)] \tag{28}$$

where *ron* is the *n*-th sensor's position index. The rationale behind this further entropy step is the following. If data at a given position are mainly contributed by clutter and noise then a relatively high level of entropy is expected. This is because the signal magnitude along time should be nearly the same. Instead, when the target significantly contributes then the entropy should decrease. Accordingly, the subset of measurements that effectively "see" the target can be roughly identified by looking for where ˆ*s*(*ron*) is below a threshold value. Say *ni* and *ns* indicate the positions in between ˆ*s*(*ron*) is below the given threshold, then a windowing is applied to keep only the time traces collected over the positions indexed between *ni* and *ns*. The entropy behaviour reported in Figure 5 (right panel) illustrates the previous discussion. In particular, only data whose positions belonged to the orange and green shaded regions should be retained during the image formation (the threshold was chosen heuristically). In particular, in the following reconstructions only sensors relative to the green zone were retained. Eventually, the second entropy step resulted in a selection of some of the rows of **sW** reported in (26), so that the data to be used were:

$$\mathbf{\dot{s}w} = \begin{bmatrix} \mathbf{s}\_{W\_{W\_i}} \\ \mathbf{s}\_{W\_{W\_{i+1}}} \\ \vdots \\ \mathbf{s}\_{W\_{W\_{N\_y}}} \end{bmatrix} \tag{29}$$

These two-step entropy strategies allowed us to select the time-gating to apply to each time trace and to select the sub-set of data (across the different sensor positions) to employ in the reconstruction stage. However, this did not yet ensure that in the remained traces there was only the tumor signal. On the contrary, the latter could still be overshadowed by clutter due to the internal inhomogeneity of the breast. To mitigate this clutter residue, it was convenient to return back into the frequency domain (even because the detection algorithm worked in such a domain) by a DFT routine. Hence, (30) becomes:

$$\mathbf{\hat{S}}\_{\mathbf{W}} = \begin{bmatrix} \mathbf{S}\_{\mathbf{W}n\_i} \\ \mathbf{S}\_{\mathbf{W}n\_{i+1}} \\ \vdots \\ \mathbf{S}\_{\mathbf{W}n\_{n\_y}} \end{bmatrix} \tag{30}$$

which is a *<sup>N</sup>*¯ *<sup>o</sup>* <sup>×</sup> *Nf* matrix, with *<sup>N</sup>*¯ *<sup>o</sup>* being the actual number of measurement positions. For the case at hand *<sup>N</sup>*¯ *<sup>o</sup>* <sup>=</sup> *ns* <sup>−</sup> *ni*. Then, it was reasonable to assume that clutter magnitude was higher than tumor signals. Accordingly, a clutter-rejection subspace-based technique was adopted. In particular, the retained scattering matrix was first expressed in terms of its singular value decomposition (SVD):

$$\mathbf{\hat{S}}\_{\mathbf{W}} = \mathbf{U} \mathbf{A} \mathbf{V}^{\mathbf{H}} \tag{31}$$

where **U** and **V** are unitary matrices containing the left and the right singular functions, respectively, and **Λ** is a diagonal matrix containing the singular values *λ*1, *λ*2, ··· , *λP*, in decreasing order, with *P* = min[*N*¯ *<sup>o</sup>*, *Nf* ]. Clutter could then be mitigated by disregarding the projection of the scattering matrix **Sˆ <sup>W</sup>** onto the singular functions corresponding to the highest singular values. The number of projections to discard generally required a priori information on the clutter, which were in general not available. However, as shown in [35], a conservative choice is to discard the projections of the scattering matrix over the singular

functions corresponding to the first or the first two highest singular values. Accordingly, the final de-cluttered data matrix was obtained as:

$$\mathbf{S\_d} = \sum\_{\substack{l=2 \ \text{or} \ l=3}}^p \lambda^l \mathbf{u}^l (\mathbf{v}^l)^H \tag{32}$$

with **ul** and **vl** being the *l*-th column vectors of **U** and **V**, respectively. Eventually, **Sd** is the data matrix passed to the I-MUSIC stage.

## *4.4. Reconstruction Results*

According to the sliced approach mentioned above, data collected at different heights were singularly processed to get the corresponding coronal slice reconstructions. In the sequel we just show only those ones obtained from data collected at the height corresponding to the centre of the tumor. In particular, although data were collected within the frequency band [0.5, 5] GHz, in the following reconstructions only the band *Bw* = [1, 3] GHz was exploited.

The rationale under the following examples is to appreciate the role played by the various steps the clutter rejection method consists of as well as the number of frequencies to be employed in the reconstructions.

The first example is shown in Figure 6 and refers to Phanton A whose MRI is reported in panel (f) to appreciate the breast internal morphology and for comparison purposes with respect to the microwave imaging. In that figure a blue dashed circle is also reported which identifies the spatial region used to perform the reconstructions. As can be seen, such a spatial region was larger than the the phantom coronal slice. In panels (a) to (e) the I-MUSIC indicator is reported. In particular, in panel (a) to (c) only 20 frequencies, uniformly taken within *Bw*, were exploited. More in detail, panel (a) shows the image obtained by pre-processing the data through only the time-gating procedure, by setting the time-gating at *tmmin* = 3 ns according to the first entropy step described above. As can be seen, the reconstruction just returned a hot spot roughly located at the centre of the imaging area. Since I-MUSIC tends to peak at the centre of targets, this means that only time-gating data were not enough to detect the tumor as data still appeared as if produced by a target whose equivalent centre was roughly in the centre of the scene. In panel (b), in order to improve clutter rejection, the subspace approach (achieved by discarding just the first projection of the scattering data matrix) was added to the time-gating. Once again, reconstruction peaks did not match with the expected tumor location, meaning that data were still dominated by a strong clutter contribution. Finally, in panel (c) we also enclosed the sensor selection procedure according to the second entropy computation described in the previous section. In particular, the 2D image was obtained by exploiting only the sensors whose index ranges from *ni* = 28 to *ns* = 56, which correspond to an angular coverage between 135◦ and 275◦. For the case at hand, hence, the entropy procedure selected the part of the measurement circular line that is closer to the tumor. This is consistent with the adopted multimonostatic configuration since the scattered field data were collected only in reflection mode. As can be seen, now the tumor was clearly detected and this highlighted that the sensor selection (often overlooked in literature) was a crucial step in addressing imaging in a highly cluttered scenarios. Panel (c) also shows that the image was strongly populated by secondary lobe contributions. This might be due to the employed reduced number of sensors (arising from the second entropy step) which mainly impacted the side lobe structure of the point-spread function. According to the discussion reported above concerning the role of the number of frequencies in mitigating aliasing artefacts, and in general side lobe structure, the quality of reconstruction can be improved by using more frequencies. This was actually the case as can be appreciated looking at panels (d) and (e), where the number of frequencies was increased to *Nf* = 40 and *Nf* = 100, respectively. In particular, panel (e) shows an extremely clear tumor detection. Additionally, in that panel we marked through a yellow circle the actual tumor position.

A few comments are in order concerning the obtained reconstruction matching to what I-MUSIC is expected to return. As explained above, I-MUSIC is a radar approach and hence, as such, it is mainly asked to provide tumor detection and rough location. Therefore, it does not aim at reproducing the breast tissue profile as MRI does. At the microwave regime, this task can be attempted by exploiting more sophisticated reconstruction methods that perform the non-linear inversion. Ideally, Φ*I*−*MUSIC* allows for very sharp tumor location (as compared to BF), as shown in the illustrative numerical examples reported above. However, because of uncertainties, clutter residues, and especially owing to the simplified model used while computing the steering vectors (we just used an equivalent homogeneous medium since the breast features are in general unknown), Φ*I*−*MUSIC* results smeared and delocalized from the actual tumor position. Indeed, this is a drawback common to any radar approach that relies on an assumed scattering model. Nonetheless, this method is extremely quick and does not require a priori information about the used antenna, which was completely ignored in the imaging procedure [33,35]. Finally, we once again remark that the circular boundaries appearing in the reconstructions (i.e., panels (a) to (e)) just delimited the spatial region within which the reconstruction was carried out. The actual boundary of the breast was removed by the clutter rejection procedure. The relative (with respect to the image area) size of phantom coronal slice can be appreciated in panel (f) of Figure 6, where the boundary of the image spatial region has been overlapped to the phantom MRI. Eventually, the spots highlighted by the I-MUSIC does not in general indicate the actual tumor positions. However, they allows a clear tumor detection and to highlight in which quadrant (of the coronal slice) it appears.

The second reconstruction example refers to phantom B and is reported in Figure 7. In this case, the entropy procedure returned the same time-gating as above and almost the same observation angular coverage. Indeed, by considering four different slice heights, *ni* ranged from 28 to 34 while *ns* remained 56. The same discussion as above applies here, with the tumor being very well detected. However, we remark that since phantom B is more complex (from a morphological point of view) than phantom A (see panel (f)) and exhibits less "circular" symmetry, it is reasonable that clutter space dimension was increased. To this end, during the subspace projection clutter reduction stage the first two (instead of only the first one) projections of the scattering matrix were discarded in Equation (32). In general, it difficult to a priori set the number of projections to be discarded. Here, we used the conservative and heuristic approach to discard at most two projections. Rejecting more projections can help in further reducing clutter but the risk is that the tumor signal can be discarded as well.

**Figure 6.** Reconstructions and MR image for phantom A. (**a**) Reconstruction with only time-gating, *Nf* = 20. (**b**) Reconstruction with time-gating + rejection of the first SVD projection of the scattering matrix, *Nf* = 20 (**c**) Reconstruction with time-gating + sensor selection + rejection of the first SVD projection of the scattering matrix, *Nf* = 20. (**d**) Reconstruction with time-gating + sensor selection + rejection of the first SVD projection of the scattering matrix, *Nf* = 40. (**e**) Reconstruction with time-gating + sensor selection + rejection of the first SVD projection of the scattering matrix, *Nf* = 100. (**f**) MR coronal slice image of Phantom A. In particular, in panel (**f**) the blue dashed circle indicates the circular boundary of the spatial region within which the reconstructions reported in the other panels have been achieved. This is highlighted even in panel (**e**). Moreover, in the latter, the yellow circle denotes the tumor location and size.

**Figure 7.** Reconstructions and MR image for phantom A. (**a**) Reconstruction with only time-gating, *Nf* = 20. (**b**) Reconstruction with time-gating + rejection of the first two SVD projections of the scattering matrix, *Nf* = 20 (**c**) Reconstruction with time-gating + sensor selection + rejection of the first two SVD projections of the scattering matrix, *Nf* = 20. (**d**) Reconstruction with time-gating + sensor selection + rejection of the first two SVD projections of the scattering matrix, *Nf* = 40. (**e**) Reconstruction with time-gating + sensor selection + rejection of the first two SVD projections of the scattering matrix, *Nf* = 100. (**f**) MR coronal slice image of Phantom B. In particular, in panel (**f**) the blue dashed circle indicates the circular boundary of the spatial region within which the reconstructions reported in the other panels have been achieved. This is highlighted even in panel (**e**). Moreover, in the latter, the yellow circle denotes the tumor location and size.

## **5. Conclusions**

Microwave breast imaging requires to deal with a number of issues which go far beyond the need do devise suitable inversion algorithms. Indeed, under the simplified linear framework subtended by the so-called radar approaches, which aim at a mere detection and localization of tumors, data must be properly pre-processed (before the imaging stage) to make sure that the resulting signals are actually useful to pursue the objective.

In this regard, one of the problems to be faced is the need to estimate or compensate for the antenna frequency response, especially when coherent wide band radar imaging methods are employed to obtain the scene image. Indeed, on one hand, antenna frequency response shapes the actually received pulse signals and modifies the overall round-trip delay; both these effects must be taken into account while implementing the beam-forming image procedure. On the other hand, because of the close proximity set-up usually adopted in breast imaging, the antenna couples with the unknown breast and its response becomes different from the free-space case. To overcome this drawback, it is shown that incoherent methods, that do not simultaneously use the frequency data but rather process each frequency separately and then combine the outcomes, can be employed with a minor reduction of the performance, especially if incoherence is used in conjunction to a MUSIC like algorithm (I-MUSIC).

Another crucial aspect is the clutter that overwhelms the target signals and can impair imaging. In this paper we have introduced a new multi-step clutter rejection method that is based on two entropy computations for time-gating setting and the selection of the sensors whose signal are less corrupted by clutter, followed by a standard subspace rejection procedure based on the SVD computation of the scattering data matrix.

The effectiveness of the de-clutter plus I-MUSIC has demonstrated against experimental data collected by using a multimodal phantom we previously developed and characterized in [33]. The results show that, for the considered phantoms, the proposed method very well succeed in detecting and localizing the tumor, though the dielectric contrast with respect to the surrounding fibroglandural tissue was only 1.2:1. This contribution can be considered as completing [33], where we mainly focused on the phantom manufacturing and characterization and only barely described the microwave imaging procedure.

As a concluding remark, we would like to remark that microwave breast imaging is a very broad research field and by this paper we did not intend to give a comprehensive account of the huge available literature. We have just focused the spot on our specific perspective.

**Author Contributions:** Conceptualization, R.S. and G.R.; methodology, R.S., A.C. and A.D.; software, A.C. and A.D.; validation, A.C., G.R. and J.B.; formal analysis, A.C. and A.D.; data curation, G.R., A.C. and J.E.B.; writing—original draft preparation, R.S. and A.C.; supervision, R.S. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research received no external funding.

**Institutional Review Board Statement:** Not applicable.

**Data Availability Statement:** No new data were created or analyzed in this study. Data sharing is not applicable to this article.

**Conflicts of Interest:** The authors declare no conflict of interest.

## **References**

