Stroke Detection and Monitoring by Means of a Multifrequency Microwave Inversion Approach

Abstract

In the area of biomedical diagnostics, microwave imaging techniques have been recently proposed for performing brain stroke detection and monitoring. Indeed, theoretically, these techniques make it possible to meet the timeliness requirements of such a diagnosis with portable systems. Moreover, relying on the use of microwaves, they are noninvasive and allow continuous monitoring of critical patients. In this paper, the microwave imaging problem is solved by exploiting multifrequency data by an inexact-Newton method formulated in the framework of non-constant exponent Lebesgue spaces. First, the method is numerically validated with three-dimensional head models affected by anatomically-realistic strokes. Then, a further assessment through experimental data obtained with a cylindrical phantom is conducted. A quite accurate reconstruction of the variations of dielectric properties inside the patient’s head due to the insurgence of stroke is obtained in both numerical and experimental cases, showing the potentiality of the proposed approach.

1. Introduction

In stroke diagnostics, it is essential to intervene early and use non-invasive techniques to employ immediate and appropriate treatments. Indeed, it is necessary to quickly recognize which kind of stroke is involved, whether ischemic or hemorrhagic, and to act properly. Moreover, since a relevant number of cases of the disease have recurrences, continuous monitoring with non-invasive techniques is often required in stroke-affected patients [1,2,3].

In this context, although computerized tomography (CT) and magnetic resonance imaging (MRI) remain the gold diagnostic standard techniques, there is a relevant interest in the development of complementary techniques, to help the initial diagnosis and monitoring of strokes [4,5,6]. In particular, microwave diagnostic techniques offer a good opportunity to meet timeliness and portability requirements [7,8,9,10,11,12]. First of all, unlike the traditional CT and MRI procedures, which require apparatuses of a considerable size, cost, and power consumption, microwave techniques give the possibility of stroke diagnosis through the use of equipment that is both easily portable and reasonably priced [13]. In addition to the cost and time advantages, microwave methods for stroke detection have the benefit of being completely safe and non-invasive for patients, which is also beneficial for continuous bedside monitoring. Compared to CT, microwave radiation does not cause the ionization of human tissues, and the size of measurement apparatus is significantly smaller. Compared to MRI, microwave systems do not require intense magnetic fields and therefore can be safely applied in proximity to magnetic materials or medical devices. Moreover, the power consumption of microwave stroke imaging systems is very low, as each antenna requires an incident power in the order of milliwatts [14], which is typically smaller than the power radiated by mobile phones during standard voice communications. For the above reasons, in recent years, several prototypes of microwave stroke diagnostic systems have been proposed [15,16,17,18,19].

In terms of detection capabilities, compared to the gold standards for medical imaging, which can provide reconstructed images of the internal parts of the head with millimetric accuracy, microwave imaging methods have a much lower resolution, which is in the order of centimeters. This is the main limitation of such techniques. This disadvantage makes microwave imaging a complementary technique in clinical diagnosis, which cannot readily substitute the gold standards. Nevertheless, it can be very useful in emergency situations and for post-event continuous monitoring.

Since microwave stroke diagnostics implies the solution of a very challenging inverse problem, the core of all systems is constituted by the adopted processing algorithm. A very interesting review of the existing systems and algorithms is provided in [13]. Basically, there are three main classes of microwave techniques that allow stroke diagnosis based on electromagnetic field measurements collected by antennas placed around the patient’s head.

The first is based on classification methods, which after proper training are able to provide a binary or multi-class output which indicates the presence of stroke and/or its type [9,20]. Advantages of these techniques are speed, accuracy and the possible absence of electromagnetic modeling (in fully data-driven techniques); disadvantages are related to the limited information that can be retrieved and the need for suitable training.

The second class of algorithms is represented by qualitative imaging techniques, in which simplified information about the presence, shape and/or location of stroke are usually provided [21,22,23,24,25,26,27]. The main advantages of qualitative techniques are their speed, and their robustness to noise and to inaccuracies in the electromagnetic model. However, the disadvantages of these methods are related to the limited retrievable information and the absence of a quantitative meaning for the reconstructed images, where it can be difficult to discriminate between artifacts and useful data.

The third class of microwave stroke diagnostic techniques concerns quantitative imaging methods, which aim at fully reconstructing the dielectric properties of the internal parts of the head by solving the full nonlinear inverse scattering problem [28,29,30,31,32,33,34,35]. The main advantage of these methods is the richness of the information that can be acquired, which can be extremely useful for a precise diagnosis, as well as for monitoring the evolution of disease. The disadvantage is a complication of the electromagnetic model and therefore of the inverse problem that should be solved, which needs specific inversion techniques.

Despite the significant advances in microwave stroke diagnosis, most quantitative techniques make use of data acquired at a single frequency [14,36] or with frequency-hopping strategies, i.e., processing multiple frequencies one at a time [37,38]. However, the benefits of methods that simultaneously process data acquired at different frequencies (i.e., multifrequency approaches) are well known and adopted in several contexts [39,40,41].

Here, a multifrequency microwave tomographic method built on a Newton-based scheme [42,43,44] is adopted to reconstruct the dielectric properties of a cross-section of the head. The proposed strategy, in particular, solves the inverse scattering problem under a two-dimensional (2D) scalar approximation by exploiting the benefits of regularization in Lebesgue spaces with variable exponents $L^{p(·)}$, where an adaptive procedure defines the values that the exponent function $p(·)$ assumes [14]. As will be shown, the possibility of simultaneously using multifrequency data allows a more accurate solution of the inverse problem, leading to more precise reconstructions of the dielectric properties of the tissues inside the patient’s head. The presented multifrequency approach is therefore a step forward with respect to the technique described in [12], where a single-frequency $L^p$-space algorithm has been presented, and compared to the imaging methods reviewed in [10,11]. Among them, it is worth noting that the quantitative approaches described in [7,8], already experimented in clinical tests, are also based on single-frequency and frequency-hopping schemes, respectively. Of course, our method is in principle able to provide more information compared to classification approaches [9].

A preliminary assessment of the proposed inversion method, in the presence of simplified numerical configurations, has been presented in [45,46]. In this paper, the analysis of the proposed multifrequency method is extended by considering for the first time the presence of anatomically realistic three-dimensional models of the stroke-affected region, whose properties have been defined based on a dataset of MRI images [47]. Furthermore, the methodology has been tested here for the first time by using experimental data involving a phantom composed of cylindrical targets filled with tissue-mimicking liquids.

The paper has the following organization. Section 2 introduces the mathematical formulation of the problem and the inversion procedure; Section 3 validates the method, first by numerical simulations and then through experimental data. Finally, Section 4 concludes the paper.

2. Formulation of the Imaging Method

The microwave imaging problem is addressed by considering a two-dimensional investigation domain $D$ placed in a horizontal transverse plane of the patient’s head (Figure 1) at a height $z=h_m$.

The adopted measurement configuration consists of a multi-view and multi-static setup with $H$ antennas arranged around the patient’s head that sequentially illuminate $D$ during the data acquisition phase. Antennas lie on the same horizontal plane of $D$ and are distributed with uniform angular spacing on the profile of an ellipse of axes $a_m$ (along $x$ axis) and $b_m$ (along $y$ axis). In the $h$-th view ($h=1,\dots ,H$) one antenna acts as source, whereas the remaining $H-1$ positions are occupied by receiving probes and constitute the measurement domain $M_h$.

A time-harmonic dependence of electromagnetic fields of type $e^{j2\pi ft}$ is also assumed and omitted in the following formulation, $f$being the frequency. In the inversion procedure, transmitting antennas are modelled as line-current sources, whereas the receiving antennas are assumed as ideal field probes. As in other related works, the mutual coupling between antennas is neglected [37,48]. The adopted assumptions of homogeneous properties of the target along the $z$ axis and of $z$-polarized transverse-magnetic (TM) incident electromagnetic fields allow consideration of the problem of microwave tomography through a two-dimensional and scalar formulation [10].

2.1. Definition of the Inverse Problem

In order to formulate the inverse problem, let us start by considering what happens at a fixed frequency $f$; later, we will combine the measurements at all the available frequencies $f_n$, $n=1, \dots , N$ to simultaneously use them for the inversion.

The patient’s head is surrounded by a homogeneous coupling medium with relative dielectric permittivity ${\epsilon }_r^b\left(f\right)$ and (equivalent) electric conductivity ${\sigma }^b\left(f\right)$. In a reference configuration (which is assumed to be available, at least approximately) the investigation domain is characterized by a relative dielectric permittivity ${\widehat{\epsilon }}_r\left(\mathbf{r},f\right)$ and an electric conductivity $\widehat{\sigma }\left(\mathbf{r},f\right)$, where $\mathbf{r}=\left(x,y\right)$ is the position vector in the cross-sectional plane where $D$ lies. In the configuration under test, the dielectric properties are assumed to be ${\epsilon }_r\left(\mathbf{r},f\right)={\widehat{\epsilon }}_r\left(\mathbf{r},f\right)+∆{\epsilon }_r\left(\mathbf{r}\right)$ and $\sigma \left(\mathbf{r},f\right)=\widehat{\sigma }\left(\mathbf{r},f\right)+∆\sigma \left(\mathbf{r}\right)$, $\mathbf{r}\in D$, i.e., the changes versus the reference dielectric profile are supposed as frequency-independent.

In this way, the problem of finding the unknown dielectric properties of $D$ is recast as finding $∆{\epsilon }_r\left(\mathbf{r}\right)$ and $∆\sigma \left(\mathbf{r}\right)$, which are grouped in the unknown vector

$$∆\mathbf{x}\left(\mathbf{r}\right)=\left[\begin{array}{cc}{\displaystyle \frac{∆{\epsilon }_r\left(\mathbf{r}\right)}{n_{ϵ}}}& {\displaystyle \frac{∆\sigma \left(\mathbf{r}\right)}{n_{\sigma }}}\end{array}\right],$$

(1)

where $n_{ϵ}={\epsilon }_r^b\left(f_1\right)$ and $n_{\sigma }={{2\pi f}_1ϵ}_0{\epsilon }_r^b\left(f_1\right)$ are the two constant normalization coefficients. It is worth noting that $∆\mathbf{x}\left(\mathbf{r}\right)=0$ for $\mathbf{r}\notin D$ is assumed.

Let us define as ${\widehat{E}}_{TOT}^{f,h}(\mathbf{r})$ and $E_{TOT}^{f,h}\left(\mathbf{r}\right)$, with $\mathbf{r}\in M_h$, the $z$-components of the total electric field in the reference configuration and in the configuration under test, respectively, at the frequency $f$, when the $h$-th antenna is in transmitting mode. To describe the difference in the measured fields $\Delta E_s^{f,h}\left(\mathbf{r}\right)=E_{TOT}^{f,h}\left(\mathbf{r}\right)-{\widehat{E}}_{TOT}^{f,h}\left(\mathbf{r}\right)$ as a function of the unknown $∆\mathbf{x}\left(\mathbf{r}\right)$, let us introduce the operators

$${\mathcal{A}}_{\left\{st,da\right\}}u\left(\mathbf{r}\right)=-k_b^2{\int }_{D}u\left({\mathbf{r}}^{\prime }\right){\widehat{g}}^f\left(\mathbf{r},{\mathbf{r}}^{\prime }\right)d{\mathbf{r}}^{\prime }, \mathbf{r}\in \left\{D,M_h\right\},$$

(2)

whose kernel is the Green’s function ${\widehat{g}}^f\left(\mathbf{r},\mathbf{r}\prime \right)$ associated to the reference dielectric configuration [10]. In (2), $k_b^2={\omega }^2{\mu }_0{\epsilon }_b^{*}\left(f\right)$, with ${\epsilon }_b^{*}\left(f\right)={\epsilon }_0{\epsilon }_r^b\left(f\right)-j{\sigma }^b\left(f\right)/2\pi f$ and ${\mu }_0$ being the magnetic permeability of the vacuum (a non-magnetic background medium is assumed). Thanks to these definitions, we have

$$\Delta E_s^{f,h}\left(\mathbf{r}\right)={\mathcal{A}}_{da}\left[\mathbf{K}\left(f\right)\cdot ∆\mathbf{x}\right]E_{TOT}^{f,h}\left(\mathbf{r}\right), \mathbf{r}\in M_h,$$

(3)

which is known as the data equation, where $\mathbf{K}\left(f\right)={\epsilon }_0/{\epsilon }_b^{*}\left(f\right)\left[\begin{array}{cc}{n}_{ϵ}& n_{\sigma }/j2\pi f{\epsilon }_0\end{array}\right]$, and

$$E_{TOT}^{f,h}\left(\mathbf{r}\right)={\widehat{E}}_{TOT}^{f,h}\left(\mathbf{r}\right)-{\mathcal{A}}_{st}\left[\mathbf{K}\left(f\right)\cdot ∆\mathbf{x}\right]E_{TOT}^{f,h}\left(\mathbf{r}\right), \mathbf{r}\in D,$$

(4)

called the state equation. Both are necessary to define the problem, and the combination of (3) and (4) results in the data and the unknown being linked by a nonlinear operator $F^{f,h}$

$$\Delta E_s^{f,h}=F^{f,h}\left(∆\mathbf{x}\right).$$

(5)

Measurements obtained from different views $h=1, \dots , H$ and frequencies are simultaneously exploited in the inversion. To this end, let us introduce the array $\Delta {\mathbf{E}}_s={\left[\Delta E_s^{f_{1,}1}, \dots , \Delta E_s^{f_{N,}H}\right]}^T$ and the operator $\mathbf{F}\left(∆\mathbf{x}\right)={\left[F^{f_{1,}1}\left(∆\mathbf{x}\right), \dots , F^{f_{N,}H}\left(∆\mathbf{x}\right)\right]}^T$. Applying (5) at each different frequency and view, the system below is attained:

$$\Delta {\mathbf{E}}_s=\mathbf{F}\left(∆\mathbf{x}\right),$$

(6)

where $∆\mathbf{x}\in X$, $\Delta {\mathbf{E}}_s\in Y$ and $\mathbf{F}: X\to Y$, with $X$, $Y$ being Lebesgue spaces with variable exponent functions, as detailed in the following paragraphs. The nonlinear Equation (6), which connects the available electric field data gathered at the various frequencies with the unknown $∆\mathbf{x}$, is inverted to find $\Delta \mathbf{x}$ with the procedure outlined in the next Section.

In the multifrequency solution of the inverse problem, an important parameter is the number of frequencies to work on (i.e., $N$). In order to study the influence of this parameter, a numerical analysis has been performed, assessing the behavior of the inversion algorithm for different numbers of considered frequencies in the selected frequency band.

2.2. Inversion Procedure

The inversion of (6) is addressed through an inexact-Newton approach [13], which employs an iterative process to find a regularized solution of the problem. This technique is based on two nested iterative cycles schematized in Algorithm 1 (main loop) and Algorithm 2 (internal loop), respectively. In more detail, Algorithm 1 is an iterative inexact-Newton method that performs a series of successive linearizations of the nonlinear Equation (6). In each iteration of Algorithm 1, the solution of the linearized equation is found by Algorithm 2, which is an iterative truncated Landweber-like regularization method formulated in Lebesgue spaces with variable exponents [49].

Algorithm 1. Inexact-Newton method.
Inputs: $\Delta {\mathbf{E}}_s\in Y$, $\delta p$, $p_{st}$,$J,K$, ${\Delta }_{th}$
1.1:	$∆{\mathbf{x}}_0\left(\mathbf{r}\right)\leftarrow 0$, $p_0\left(\mathbf{r}\right)\leftarrow p_{st}$
1.2:	for  $j=0,\dots ,J \mathbf{d}\mathbf{o}$
1.3:		${\mathbf{F}}_j^{\prime }\mathsf{\delta }\left(\mathbf{r}\right)=\mathbf{F}\left(∆{\mathbf{x}}_j\right)\left(\mathbf{r}\right)-\Delta {\mathbf{E}}_s\left(\mathbf{r}\right)$
1.4:		$\mathsf{\delta }\leftarrow$Algorithm 2 ($\Delta {\mathbf{E}}_s,{\mathbf{F}}_j^{\prime },K,p_j,{\Delta }_{th})$
1.5:		$∆{\mathbf{x}}_{j+1}\left(\mathbf{r}\right)\leftarrow ∆{\mathbf{x}}_j\left(\mathbf{r}\right)+\mathsf{\delta }\left(\mathbf{r}\right)$
1.6:		$p_{j+1}\left(\mathbf{r}\right)=2-\delta p\left(1-{\displaystyle \frac{\left|∆{\mathbf{x}}_{j+1}\left(\mathbf{r}\right)\right|}{\underset{\mathbf{r}\in D}{\mathrm{m}\mathrm{a}\mathrm{x}}\left|∆{\mathbf{x}}_{j+1}\left(\mathbf{r}\right)\right|}}\right)$
1.7:		if $($stop criterion is satisfied$)$ then
1.8:			Break
1.9:		End
1.10:	end
Output:   $∆\mathbf{x}\in X$

Algorithm 2. Landweber-like regularization method.
Inputs: $\Delta {\mathbf{E}}_s\in Y, {\mathbf{F}}_j^{\prime },K$, $p_j$, ${\Delta }_{th}$
2.1:	$\beta \leftarrow {1/‖{\mathbf{F}}_j^{\prime }‖}^2$
2.2:	$\mathsf{\delta }\leftarrow 0$
2.3:	$\mathbf{for} k=0,\dots ,K \mathbf{d}\mathbf{o}$
2.4:		$\mathsf{\delta }\leftarrow {\mathbf{J}}_{X^{*}}\left\{{\mathbf{J}}_X\left(\mathit{\delta }\right)-\beta {\mathbf{F}}_j^{\prime \mathbf{*}}{\mathbf{J}}_Y\left({\mathbf{F}}_j^{\prime }\mathsf{\delta }- \left[\Delta {\mathbf{E}}_s-\mathbf{F}\left(\mathsf{\delta }\right)\right]\right)\right\}$
2.5:		$\mathbf{if} ($$\mathrm{stop} \mathrm{criterion} \mathrm{is} \mathrm{satisfied})$ then
2.6:			break
2.7:		end
2.8:	End
$\mathbf{Output}: \mathsf{\delta }\in X$

Algorithm 1 takes in input the array $\Delta {\mathbf{E}}_s$, which contains the scattered fields measured for each frequency and view, as well as some parameters of the inexact-Newton method. Among them: $\delta p$, the exponent function range of variation; $p_{st}$, the initial value of the exponent function; $J$, the maximum number of main loop iterations; $K$, the maximum number of internal loop iterations; ${\Delta }_{th}$, the threshold on the relative variation of the data residual in two subsequent iterations.

At first, in Algorithm 1, assuming no a-priori information about the target, the solution method is initialized with an empty guess ${∆\mathbf{x}}_0=0$ and a fixed exponent function $p_0\left(\mathbf{r}\right)=p_{st}$ (line 1.1) in both the unknown space, $X$, and data space,$Y$. After that, a set of solution estimates $\Delta {\mathbf{x}}_j$ ($j=0, \dots , J$) is constructed by means of a series of successive linearizations of (6): line 1.3 represents the $j$th linearized equation, where ${\mathbf{F}}_j^{\prime }$ is the Fréchet derivative of $\mathbf{F}$ around $\Delta {\mathbf{x}}_j$, $\mathsf{\delta }\left(\mathbf{r}\right)$, $\mathbf{r}\in D$ is the unknown of the linearized equation, and $\mathbf{F}\left(∆{\mathbf{x}}_j\right)$ is the value of the operator $\mathbf{F}$ applied to $∆{\mathbf{x}}_j$. The equation of line 1.3 is solved in the unknown $\mathsf{\delta }$ by applying Algorithm 2, i.e., the truncated Landweber-like regularization approach formulated in $L^{p\left(\cdot \right)}$ spaces.

In line 1.5, the solution of the linearized equation found by Algorithm 2 is used to update the global solution, which becomes $∆{\mathbf{x}}_{j+1}\left(\mathbf{r}\right)$. Then, the exponent functions are tuned according to the updated solution. In more detail, $X$ is considered as a variable exponent space where the exponent $p_j\left(\mathbf{r}\right)$ (at the $j$th iteration) depends on the position $\mathbf{r}$ in $D$, while $Y$ is a constant-exponent space with exponent computed as the spatial average value of $p_j\left(\mathbf{r}\right), \mathbf{r}\in D$. The variable exponent function for the next iteration, $p_{j+1}\left(\mathbf{r}\right)$, is chosen as reported in line 1.6, with $0\le \delta p<1$. In this way, lower values of the exponent are given to the regions of $D$ where no targets are detected, and higher values otherwise [49].

The exponent functions defined above are exploited by Algorithm 2, thanks to its $L^{p\left(\cdot \right)}$-space formulation. In particular, this algorithm takes in input the scattered field array $\Delta {\mathbf{E}}_s$, the Fréchet derivative of the scattering operator ${\mathbf{F}}_j^{\prime }$, the maximum inner iterations $K$, the exponent function $p_j$ and the threshold on the relative variation of the data residual ${\Delta }_{th}$. At the beginning, the value of the relaxation coefficient $\beta$ is set as ${1/‖{\mathbf{F}}_j^{\prime }‖}^2$, where $‖{\mathbf{F}}_j^{\prime }‖$ is the spectral norm of ${\mathbf{F}}_j^{\prime }$ (2.1), and the solution is initialized with an empty guess (2.2).

The main iteration of the Landweber-like algorithm is described in line 2.4, where ${\mathbf{J}}_{X^{*}}$, ${\mathbf{J}}_X$ and ${\mathbf{J}}_Y$ are the duality maps of spaces $X^{*}$, $X$ and $Y$, respectively ($X^{*}$ being the dual space of $X$). Detailed definition on the variable-exponent duality maps can be found in [37]. Moreover, ${\mathbf{F}}_j^{\prime \mathit{*}}$ is the adjoint of the linearized operator ${\mathbf{F}}_j^{\prime }$.

The algorithms terminate when the maximum numbers of iterations $J$ and $K$ are reached (lines 1.2 and 2.3, respectively) or when a stop criterion based on data residuals is satisfied (lines 1.8 and 2.6). This criterion terminates the iterations when the relative variation on data residuals in two successive steps is below a threshold value ${\Delta }_{th}$.

3. Results

The method is first validated using synthetic data obtained by three-dimensional numerical models of the head affected by stroke, as described in Section 3.1. Finally, the validation through experimental data concerning head models made with cylindrical targets is presented in Section 3.2.

3.1. Numerical Results with Realistic Stroke Models

As to the numerical case studies, the finite-difference time-domain (FDTD) open-source software gprMax version 3.1.5 [50] has been used for the three-dimensional solution of the forward electromagnetic problem to generate synthetic field data. The simulations are based on the AustinWoman voxel head model [51], where the dielectric parameters of the various head tissues have been derived from Debye models of the third and fourth order (when available), according to [52,53]. A coupling medium composed by a mixture of water and glycerin (70% volumetric) surrounds the head phantom, with the dielectric characteristics given in [37]. Regarding the imaging configuration, $H=16$ antenna positions have been placed around the head. Each one, in turn, is occupied by a source modelled by means of a Hertzian dipole, characterized by $z$-polarized radiation on the plane perpendicular to the antenna and passing through its center. Moreover, the excitation waveform feeding the sources in time domain is a Gaussian derivative pulse whose spectrum is centered at $f_0=$1 GHz, with −3 dB bandwidth between 480 MHz and 1.63 GHz. An FDTD simulation domain of dimensions $28.4$ cm $\times$ $32.0$ cm $\times$ 28.4 cm is assumed, and $N_d=\mathrm{3,408,000}$ cubic subdomains of side $e=2\mathrm{mm}$ have been considered for its discretization. Absorbing boundary conditions have been enforced with a perfectly matched layer of width equal to 10 subdomains on all sides of the simulation region. A time window of length $T=0.03\mathsf{\mu }\mathrm{s}$ has been analyzed in the FDTD simulation, considering a temporal step $\Delta t=3.85$ ps.

For the generation of simulated data, consideration of hemorrhagic strokes has been chosen. This kind of stroke was modeled by introducing into the brain of the considered head model an irregularly-shaped inclusion, characterized by the dielectric properties of blood, following the procedure described in [48]. Specifically, the adopted stroke models have been developed on the basis of the MRI images from the dataset reported in [47], in order to consider realistic stroke shapes. In particular, two head phantoms affected by stroke have been taken into account for numerical tests. The first (Stroke #1) has a lesion placed inside the left hemisphere, with center at $\left(12.5, 16.9, 18.7\right) \mathrm{cm}$ and main axes of length of $\left(2, 6, 4\right) \mathrm{cm}$ approximately. In the second case (Stroke #2), the affected region is in the right hemisphere centered at $\left(17.5, 18.9, 17.7\right) \mathrm{cm}$, and its main axes are $\left(4.2, 12, 6\right) \mathrm{cm}$ long. The second stroke presents a more complex shape than the first case, with jagged edges and inside cavities. More details on the definition of such realistic phantoms can be found in [48]. In Figure 2a, the 3-D model of the head is shown for the first case, while in Figure 2d the model for the second case is reported. In Figure 2b,c and in Figure 2e,f, the dielectric properties in the transverse plane located at $z=18.7$ cm at 600 MHz are reported, respectively, for both considered strokes.

As is well known, in brain stroke imaging, the selection of the operating frequency band is a compromise between the penetration of the electromagnetic waves inside the head and the achievable resolution of the reconstructed image. In this particular case, based on the results obtained in [12,37,48], the band between 600 MHz and 900 MHz has been proven to provide good reconstruction capabilities. Moreover, a further increase in the frequency does not yield significant improvements. Therefore, this frequency band has also been adopted in this work.

Figure 3 shows the behavior of the relative dielectric permittivity and electrical conductivity in the frequency band of interest for some of the main tissues characterizing the head and stroke models (i.e., blood, white and grey matter). It can be seen that, in the considered frequency band, the differences between the values of the dielectric properties are almost constant, thus confirming that the assumptions of a frequency-independent differential unknown $∆\mathbf{x}$ is suitable in this case.

Once FDTD simulations have been executed (of both stroke-affected models and a healthy head model, used as reference configuration), by means of Fast Fourier Transform (FFT) frequency domain field data ${\mathrm{E}}_{TOT}^{f_n,h}\left(\mathbf{r}\right)$ and ${\widehat{\mathrm{E}}}_{TOT}^{f_n,h}\left(\mathbf{r}\right)$ have been computed with $h=1,\dots ,H$ for $n=1,\dots , N$ equally spaced values, with a frequency step $∆f$ between $\left[600, 900\right] \mathrm{MHz}$. Finally, an additive white noise with zero-mean Gaussian distribution and $SNR=25\mathrm{dB}$ has been added to the data. For the inversion procedure, $D$ lies on the $xy$plane placed at the same height as the center of the antennas, i.e., at $z=18.7$cm.

As a first analysis, the behavior of the inversion method by varying the number of frequencies considered within the band of interest has been studied. In particular, four cases have been compared: single-frequency reconstruction at 600 MHz and multifrequency reconstruction for three different values of $∆f$: $∆f=150\mathrm{MHz}$ ($N=3$ frequencies), $∆f=100\mathrm{MHz} (N=4$frequencies) and $∆f=50\mathrm{MHz}$ ($N=7$ frequencies). The parameters $p_{st}=1.4$, $\delta p=0.6$, $J=50$, $K=100$, and ${\Delta }_{th}=0.001$ have been selected for inversion on the basis of previous heuristic analyses and results, such as those provided in [48,49]. Moreover, the investigation domain is composed of $N_i=1388$ square subdomains with $e=4$mm side.

Figure 4a,b shows the reconstructed distributions of the relative dielectric permittivity and electric conductivity in the single-frequency case for Stroke #1, while Figure 4i,j reports the same quantities for Stroke #2. Moreover, Figure 4c–h,k–p illustrate the reconstructed dielectric properties for Stroke #1 and Stroke #2 when $\Delta f$ is equal to $150 \mathrm{MHz}$, $100 \mathrm{MHz}$ and $50\mathrm{MHz}$, i.e., for all the considered frequency steps. For both stroke types, via a qualitative analysis, with the single-frequency inversion, the method produces a less accurate image, especially regarding the electric conductivity, which presents a significant artefact close to the targets. Instead, considering the multifrequency approach (e.g., for $∆f=$50 MHz, $N=7$), the performance improves significantly: the dielectric inhomogeneity due to the strokes becomes clearly evident, even though some background artifacts are still visible. The presence of artifacts, even with the proposed multifrequency approach, can be caused by multiple factors. The main reasons are related to the strong ill-posed nature of the inverse scattering problem, and by the fact that the scattered field due to the stroke is very weak. Therefore, artifacts may be ascribed to the presence of noise, as well as to the mismatch between the electromagnetic model used for the inversion (which is two-dimensional and simplified) and the actual measurement configuration. However, despite the possible presence of artifacts, thanks to the quantitative reconstruction of the dielectric properties, it is possible to differentiate between the stroke region and the background, since they exhibit different values of such parameters.

For comparison purposes, the following relative reconstruction errors have been evaluated considering the value of the contrast function at 600 MHz:

$$e_{obj}={\displaystyle \frac{1}{S_{obj}}}{\int }_{D_{obj}}{\displaystyle \frac{\left|{\chi }_R\left(\mathbf{r}\right)-{\chi }_A\left(\mathbf{r}\right)\right|}{\left|{\chi }_A\left(\mathbf{r}\right)+1\right|}}d\mathbf{r},$$

(7)

$$e_{bkg}={\displaystyle \frac{1}{S_{bkg}}}{\int }_{D_{bkg}}{\displaystyle \frac{\left|{\chi }_R\left(\mathbf{r}\right)-{\chi }_A\left(\mathbf{r}\right)\right|}{\left|{\chi }_A\left(\mathbf{r}\right)+1\right|}}d\mathbf{r},$$

(8)

$$e_{inv}={\displaystyle \frac{1}{S_{inv}}}{\int }_{D_{inv}}{\displaystyle \frac{\left|{\chi }_R\left(\mathbf{r}\right)-{\chi }_A\left(\mathbf{r}\right)\right|}{\left|{\chi }_A\left(\mathbf{r}\right)+1\right|}}d\mathbf{r},$$

(9)

where $D_{obj}$,$D_{bkg}$ are the regions of the stroke and the of the background, $S_{obj}$ and $S_{bkg}$ being their dimensions, respectively. Moreover, $S_{inv}$ is the dimension of the investigation domain, ${\chi }_R$ indicates the reconstructed contrast function and ${\chi }_A$ represents the actual contrast function.

The normalized mean square error in all the investigation domain $D$, $e_{NMSE}$, has also been computed as follows:

$$e_{NMSE}={\displaystyle \frac{{‖{\chi }_R-{\chi }_A‖}_2}{{‖{\chi }_A+1‖}_2}}.$$

(10)

Errors $e_{NMSE}$ and $e_{obj}$ are reported in Figure 5 for the reconstruction of Stroke #1 and #2 by the single- and multifrequency approach, varying $\Delta f$.

It can be noticed that both $e_{obj}$ and $e_{NMSE}$ significantly decrease by adopting the multifrequency approach as compared to the single-frequency case. On the one hand, the inaccuracy in the reconstruction of the stroke region measured by $e_{obj}$ significantly reduces in all cases when $\Delta f$ decreases, i.e., more frequencies are taken into account. On the other hand, the normalized mean square error for the contrast function, $e_{NMSE}$, shows a significant variation in both strokes between single-frequency approach and multifrequency inversion with $∆f=150 \mathrm{MHz}$. However, for Stroke #1, a significant variation can again be observed by further reducing $\Delta f$, whereas for Stroke #2 the error remains stable. In

Table 1. Numerical results. Percentage relative reconstruction error for the whole investigation domain ($e_{inv}$) and the background ($e_{bkg}$) versus the operating frequencies.

Operating Frequencies	Stroke #1		Stroke #2
	${\mathit{e}}_{\mathit{i}\mathit{n}\mathit{v}}$	${\mathit{e}}_{\mathit{b}\mathit{k}\mathit{g}}$	${\mathit{e}}_{\mathit{i}\mathit{n}\mathit{v}}$	${\mathit{e}}_{\mathit{b}\mathit{k}\mathit{g}}$
600 MHz	4.26%	3.02%	7.34%	5.74%
{600, 750, 900} MHz	3.37%	2.36%	6.90%	5.65%
{600, 700, 800, 900} MHz	3.12%	2.12%	6.10%	4.87%
{600, 650, 700, 750, 800, 850, 900} MHz	3.73%	2.87%	6.52%	5.42%

, the percentage relative reconstruction error for the whole investigation domain ($e_{inv}$) and on the background ($e_{bkg}$) versus the operating frequencies is reported, for both the cases presented in the numerical analysis. In both stroke models, a reduction in the global percentage error can be noticed by passing from the single-frequency reconstruction at 600 MHz to the multifrequency case with $N=3$ and $N=4$ frequencies. A further increase in the number of frequencies considered leads to a slight increase in the percentage error. This fact is due to the slight degradation of the background reconstruction with the highest number of frequencies, although accompanied by improvement in the reconstruction of the stroke region (as is visible in Figure 5).

In summary, the numerical results provided highlight that the use of the presented multifrequency method is particularly useful to address the ill-posed nature and nonlinearity of the inverse scattering problem at hand, leading to more accurate dielectric reconstructions.

3.2. Experimental Results

After the numerical validation, the second goal of the present work has been to assess the method through experimental data. The experimental setup is visible in Figure 6 and consists of a set of $H=16$ bowtie-type antennas placed around a cylindrical target, a vector network analyzer (VNA), a switch matrix that links the antennas and the VNA to collect data for all views and, finally, a PC that handles the data acquisition. More details about the setup can be found in [37]. In particular, the structure of the adopted antennas, shown in Figure 6c, has been introduced in [12]. Each antenna has an overall size of $4.5\times 2.5\times 0.7 \mathrm{cm}$ and is printed on two superimposed FR4 substrates. In the frontal face, it has a slotted bowtie structure. The rear metallized substrate, combined with two copper strips used to connect the substrates, act as a cavity backing. Such antennas have been successfully used, starting from 600 MHz, and the operating band extends up to 2.5 GHz [12,37,48].

In order to mimic the dielectric properties of the human head, glycerin-based mixtures have been adopted several times in the scientific literature [12,14,28,37,48]. As highlighted in [54], mixing glycerin with different volumetric percentages of water allows a wide range of dielectric properties, suitable for approximation of those of some biological tissues. In particular, in [12,37,48], a mixture with 70% (volumetric) of glycerin has been adopted to simulate the average properties of the head, and in [12,37,48] saline solutions have also been used to simulate the presence of a hemorrhagic brain stroke, due to the higher dielectric properties compared to the surrounding glycerin-based mixture. In the present case, as a reference configuration, a polypropylene graduated beaker filled with a 70% glycerin/water mixture [54] with radius $r_o=9$ cm, mimicking the average properties of the human head, has been adopted. The mixture has been carefully prepared and homogenized by using an immersion sonicator. Then, a second cylinder of radius $r_i=1$ cm at ($4.5, 0)$ cm has been placed inside the outer target and filled with 0.9% saline solution, to model the presence of stroke. A layer composed by a 70% glycerin/water mixture has been introduced between the outer cylinder and the antennas, using polyethylene bags to enclose the matching medium. The dielectric properties of the adopted mixtures have been verified just after data collection by means of short-circuited coaxial line measurements, as detailed in [37]. In particular, the dielectric properties of the 70% glycerin/water mixture in the band of interest are ${ϵ}_r=47.05\pm 2.63$ and $\sigma =0.79\pm 0.22 \mathrm{S}/\mathrm{m}$, while the dielectric properties of the 0.9% saline solution are ${ϵ}_r=75.57\pm 3.09$ and $\sigma =2.15\pm 0.05 \mathrm{S}/\mathrm{m}$. After data have been collected with the above-mentioned setup, the method has been run with an investigation domain divided into $N_i=1264$ square cells of dimensions equal to those considered in Section 3.1, with all ${\Delta }_{th}=0.03$. The single-frequency ($f=600\mathrm{MHz}$) and the multifrequency reconstruction with $∆f=50\mathrm{MHz}$ $(N=3$frequencies), $∆f=100 \mathrm{M}\mathrm{H}\mathrm{z} (N=4$frequencies), and $∆f=50\mathrm{MHz}$ ($N=7$ frequencies) have been obtained. Figure 7a,b shows the reconstructed relative dielectric permittivity and electric conductivity for the single-frequency case, while Figure 7c–h shows the retrieved dielectric properties when $\Delta f$ is equal to $150 \mathrm{MHz}$, $100 \mathrm{MHz}$ and $50\mathrm{MHz}$.

As can be observed in first figure, the dielectric permittivity is overestimated (the average value of $∆{ϵ}_r=28.49$ and $∆\sigma =1.36$ S/m) and the target is not correctly localized in the conductivity reconstruction. Conversely, in the set of figures related to the multifrequency reconstruction, the permittivity is close to the actual values, the conductivity better localizes the target and most of background artifacts disappear. In Figure 8, the relative error on the contrast function inside the target region and the NMSE in the investigation domain are reported. As can be noticed, $e_{obj}$ decreases significantly between single-frequency case and multifrequency reconstruction with $∆f=10$0 MHz $(N=4$), and then the error tends to flatten out. As for $e_{NMSE}$, there is a noticeable reduction when $∆f=$ 150 MHz ($N=3$) with respect to single-frequency reconstruction, while a smaller decrease can be observed by further reducing $\Delta f$. The percentage relative reconstruction errors for the whole investigation domain ($e_{inv}$) and on the background ($e_{bkg}$) are given in

Table 2. Experimental results. Percentage relative reconstruction error on the whole investigation domain (e_inv) and on the background (e_bkg) versus the operating frequencies.

Operating Frequencies	${\mathit{e}}_{\mathit{i}\mathit{n}\mathit{v}}$	${\mathit{e}}_{\mathit{b}\mathit{k}\mathit{g}}$
600 MHz	3.55%	3.09%
{600, 750, 900} MHz	1.78%	1.37%
{600, 700, 800, 900} MHz	2.31%	1.97%
{600, 650, 700, 750, 800, 850, 900} MHz	2.06%	1.73%

.

The reported results confirm that the multifrequency reconstruction, which is the main contribution of this work, is able to provide better results compared to the single-frequency method even in the presence of real measured data.

4. Conclusions

This paper focused on a multifrequency diagnostic method operating in variable exponent Lebesgue spaces. First, the method has been validated by simulating anatomically realistic numerical head models containing different types of strokes, whose shape has been designed on the basis of MRI images of stroke-affected patients. In particular, a validation study has been performed in which the single-frequency approach and multifrequency method with different numbers of input frequencies have been compared. Then, an experimental assessment has been carried out with the use of simplified cylindrical targets. It has been shown that the variable exponent multifrequency reconstruction, which represents the contribution of this work, led to improved quantitative imaging results compared to the single-frequency approach in all the studied cases, better addressing the ill-posed nature and nonlinearity of the inverse problem. In order to apply the developed method to real-world scenarios, one of the first steps will be the validation of the approach with experimental data collected with more realistic 3D-printed phantoms of stroke-affected heads instead of the cylindrical targets used in this paper. In addition, the inclusion of detailed models of the antennas, which take into account the mutual interactions between them, will be considered. Finally, the extension of the method to three-dimensional imaging will be studied.