Three-Dimensional Microwave Head Imaging with GPU-Based FDTD and the DBIM Method

We present a preliminary study of microwave head imaging using a three-dimensional (3-D) implementation of the distorted Born iterative method (DBIM). Our aim is to examine the benefits of using the more computationally intensive 3-D implementation in scenarios where limited prior information is available, or when the target occupies an area that is not covered by the imaging array’s transverse planes. We show that, in some cases, the 3-D implementation outperforms its two-dimensional (2-D) counterpart despite the increased number of unknowns for the linear problem at each DBIM iteration. We also discuss how the 3-D algorithm can be implemented efficiently using graphic processing units (GPUs) and validate this implementation with experimental data from a simplified brain phantom. In this work, we have implemented a non-linear microwave imaging approach using DBIM with GPU-accelerated FDTD. Moreover, the paper offers a direct comparison of 2-D and 3-D microwave tomography implementations for head imaging and stroke detection in inhomogenous anatomically complex numerical head phantoms.


Introduction
Microwave medical imaging (MWI) is gaining increased research interest [1][2][3][4][5] for its potential to offer low-cost, portable solutions to important healthcare needs such as early screening for breast cancer or detecting acute stroke inside an ambulance. Microwave tomography (MWT) is an important MWI technique which estimates the spatial distribution of dielectric properties in a reconstruction region by solving the electromagnetic (EM) inverse scattering problem [1]. EM inverse scattering processes scattered signals from the imaging domain to estimate the spatial distribution of its complex permittivity (i.e., dielectric constant and conductivity).
MWT relies on the difference in the dielectric properties of the different human tissues [6,7]. In breast imaging, for example, the dielectric contrast of cancerous vs. healthy tissue depends on the density of the tissue surrounding the tumor, and several studies have been published recently to assess this contrast with different methodologies [8][9][10][11][12][13]. High-frequency EM waves increase resolution but attenuate quickly, while low frequencies result in good penetration depth but low resolution. This trade-off dictates the working frequency range (e.g., 0.5-2.0 GHz for head imaging systems [14]).
It is well known that the inverse scattering problem of MWI is intrinsically nonlinear, and this may lead to false solutions due to the existence of local minima [15,16]. The non-linearity of the problem is combined with limited knowledge of prior information and a limited number of observation data points compared to the number of unknowns. Moreover, measurement noise and a mismatch between the model used in the imaging algorithm and the true experiment present additional challenges towards an accurate solution. Non-linearity, ill-posedness, and limited prior information in MWT require advanced optimization and regularization techniques with a careful selection of parameters that are problem-specific.
A few MWT experimental prototypes have been developed in the last twenty years [17][18][19][20][21][22], and various methods have been applied to the resulting EM inverse scattering problem [23]. Most of the work, however, has focused on two-dimensional (2-D) algorithms to avoid the added complexity of three-dimensional (3-D) implementations.
MWT 3-D implementations for medical imaging have been proposed in the literature for a few years now, but to the best of the authors' knowledge, they have been tested primarily with simulation data [24][25][26][27][28]. A notable exception is the 3-D MWT algorithm in [2], which produced clinical 3-D microwave tomographic images of the breast by combining the discrete dipole approximation (DDA) as the forward solver with the Gauss-Newton (GN) algorithm.
MWT 2-D problems are easier to model and solve numerically, and they can often produce reconstructions of sufficient quality. In contrast, 3-D algorithms can be computationally expensive and require solving a more severely ill-posed inverse problem due to the much higher ratio of unknowns relative to data points. In iterative GN implementations such as the distorted Born iterative method (DBIM), a 3-D forward scattering problem must be solved at each iteration. This can be computationally prohibitive without the use of semi-analytical approaches [2] or numerical methods implemented in high performance hardware such as graphic processing units (GPUs) [25,29]. CPU parallelization can also model 3-D MWI forward solvers efficiently with CG type inverse solvers as in [30], where the finite element method (FEM) was combined with the contrast source inversion (CSI) method to image an experimental breast phantom, and in [31], which uses FDTD as a forward solver for reconstructing a weak scatterer. The latter is implemented on a parallel Linux cluster with a 20-core CPU and a running time of 1 h for each iteration is required.
For microwave head imaging, studies have analysed scattering from 3-D head models using FDTD [32] or FEM [33], and the FFT-based volume integral equation (VIE) method [29]. Most imaging algorithms are 2-D with only a few papers presenting 3-D reconstructions. For example, a 2-D Newton-type method using S-parameters for inversion has been proposed and tested with a experimental SAM phantom [34]. A frequency-domain beamforming imaging algorithm with Bessel function has been applied to qualitative imaging of a realistic human head in 2-D [35]. This proposed method has been further applied to 3-D head imaging using a wearable EM cap with 16 planar antennas [36] and compared with 2-D results obtained by the polar sensitivity encoding method [37] that works on the principle of encoding S-parameters. A GN-type algorithm using the open source FreeFEM++ solver has been implemented for imaging a numerical head phantom, which solves the inverse problem by considering five "sub-problems" of cross-sections in parallel across five rings of the antenna array [38]. Recently, an FEM-based and FDTD-based iterative algorithm have been implemented on a single GPU card with Tikhonov regularization and a gradient-based inverse solver [4]. This work has performed 2-D and 3-D reconstructions with both a numerical and an anthropomorphic mannequin head phantom, using a matching medium with permittivity comparable to the brain material. 3-D imaging has also been proposed in [39], which combined FEM with TSVD for linear inversion.
It is important to note that the complex structure of the head cannot be known a priori, resulting in a highly non-linear and ill-posed EM inverse scattering problem. Assuming that some prior information may be available, it is therefore important to investigate performance under different imaging scenarios. To this end, we have recently applied the DBIM with the fast iterative shrinkage/thresholding algorithm (FISTA) [40] to 2-D numerical brain phantoms with limited prior information. Results showed that it is possible to reconstruct the target inside a head phantom of known shape with several unknown tissue layers [41].
This work investigates microwave head imaging further with the help a 3-D DBIM algorithm. This allows to test performance for more realistic scenarios, where weak target responses are captured by a limited number of antennas surrounding the head in a 3-D array. To this end, we have implemented an in-house GPU-based 3-D FDTD forward solver, which enables a computationally efficient 3-D implementation of our recently proposed DBIM-FISTA. By optimizing the FDTD code specifically for our application, we achieve a more computationally efficient implementation than previous work [3,25] which used GPU-based platform Acceleware [42] to perform MWT reconstructions of numerical breast phantoms of varying tissue density.
We also simplify the 3-D inversion further by considering a scalar Green's function, which assumes a single electric field component transmission by our linearly polarized antenna. We compare results from this approximation with the vectorial Green's function which considers all three electric field components inside the reconstruction domain. Importantly, our results show, for the first time in microwave head imaging, that the scalar 3-D approximation leads to almost identical results with the vectorial implementation. Our DBIM algorithm is also more efficient at each DBIM iteration, as the employed FISTA algorithm uses an accelerated momentum which leads to faster convergence than CG-type inverse solvers [43]. Moreover, as a shirnkage/thresholding algorithm, it helps to reduce the chance of getting stuck into local minima [44].
The algorithm is tested on different reconstruction scenarios with the MRI-based Zubal head phantom [45] and is validated with an experimental brain phantom. Beyond validating the algorithm, our study examines the benefits of imaging in 3-D by comparing with 2-D implementations for challenging target location and limited prior information of the head's structure. To this end, we use a cross-sectional axial slice from the 3-D brain model and compare reconstructed images for this slice using the 2-D and 3-D algorithm.
The remainder of the paper is organized as follows: In Section 2, we review our recently proposed DBIM-FISTA algorithm [40] for reconstructing the spatial distribution of dielectric properties inside the region of interest. We also presents the algorithm's implementation in 3-D with our in-house, GPU-based FDTD forward solver. It also details the experimental head phantom and the different head models based on the Zubal head phantom used to test our algorithm. Reconstruction results for the numerical models and experimental phantom are presented in Section 3, followed by discussion in Section 4.

The DBIM-FISTA Algorithm
The DBIM can solve nonlinear EM inverse scattering problems iteratively to reconstruct the spatial distribution of dielectric properties within a region V [25]. It is based on approximating the non-linear integral equation which describes the relationship of the electric field with the continuous spatial distribution of dielectric properties via the Born approximation. The non-linear integral equation of the fields scattered by the object to be imaged for each transmitter-receiver (TR) pair can be written as, whereĒ t ,Ē s ,Ē b are the total, scattered and background vector electric fields, respectively. The total field is measured at each antenna, but is unknown inside region V, and is approximated by E b under the Born approximation. The vectors r m and r n denote the transmitting and receiving antenna locations, ω is the angular frequency, µ 0 and ε 0 are the permeability and permittivity of free space, andḠ b is the dyadic Green's function for the background medium. The difference δ ε between the relative complex permittivities of reconstruction ε(r) and background ε b (r) is defined as Assuming an equal discretization step ∆h for each dimension and constant permittivity inside each voxel, (1) can be written as a discrete equation as, 1 Note that the current terms in (7) correspond to excitation currents for the transmitting antennas. Hence, only the last row of the matrix is non-zero for a z-polarized antenna, leading to a Green's function that is reduced to [25,47], Using (8), the vector product ofḠ b andĒ b in (3) is computed as, As our antenna records the z-component of the electric field, a scalar approximation of the scattered fieldĒ s =Ē s ·ẑ = E s z (10) in the LHS of (3) can be used [25,47,48]. Thus (3) is further reduced to a scalar equation as, where the z-component of product (9) is computed as, If we assume that the cross-products of xand y-components are negligible [25] when a z-polarized antenna is used, we can simplify the Green's function further to account only for the z-directed components as We must emphasize that the simplification of assuming a z-polarized source is only done by our imaging algorithm, and is not used to generate the data in CST, which models the exact antenna used in our experiments. Moreover, the assumption that only the z-component of the Green's function is non-zero is employed only in the inversion of the linear matrix at each DBIM iteration. The 3D FDTD forward solver used by the imaging algorithm can calculate all three field components. This allows comparing the z-only approximation of the Green's function with the vector formulation in (13). We have compared results from these two formulations for all our imaging scenarios and concluded that the scalar approximation produced results of similar quality. More details on this important contribution of this work are presented in Section 3.
By substituting the Green's function in (3), the ill-posed linear system is built as, where , with M transmit-receive pairs and N voxels of the reconstruction region V, p is the p-th TR pair (r n , r m ) p , q is the q-th voxel inside region V, and b = {E s (p)} is a M × 1 vector of the scattered fields. The computational complexity of (16) using the dyadic Green's function is O(M × 9N), while the computational complexity further reduced to O(M × 3N) and O(M × N) respectively when using (13) and (14). At each iteration K, the forward solver (FDTD for our implementation) is run first to obtain the background data and then used to build (15), which is solved by the inverse solver FISTA. Finally, the background profile is updated by (17) and the DBIM continues to next iteration K + 1.
We use a single-pole Debye model to simulate the frequency-dependent materials in FDTD as, where ε ∞ , ∆ε, τ and σ s are the parameters of the Debye model. Thus (15) is constructed as, where and represents the real and imaginary part respectively and As we only consider real valued computations, the following linear system is built as, wherē We also use the convolutional perfectly matched layer (CPML) absorbing boundary condition to terminate the FDTD computational domain.
A minimization problem with a regularization term is considered by FISTA as the following steps when l 1 norm is chosen, where λ is a regularization parameter. Thus the structure of FISTA is constructed as [43], where t 0 = 1 and t k+1 = ψ is the soft thresholding function y k is the solution at the k-th iteration. L k is a non-negative parameter which is selected based on the following strategy: Find the smallest non-negative integers i k with such that where η > 1 and Q L (x, y) = 1 2 The computational complexity for FISTA is O(kMN) with a higher order convergence rate O( 1 k 2 ) compared with the traditional CG-type algorithm. The improvement is caused by the use of the accelerated momentum t k , which also reduces the chance of getting stuck into a local minimum when combined with a regularization term [44]. Our implementation selects the initial value of x to be 0 and the regularization parameter as where δ is a factor with 0 < δ < 1.
The stopping criterion is based on the relative error between current F(x k ) and previous F(x k−1 ), defined as FISTA stops when the relative error e opt becomes smaller than a preset value, usually chosen between 10 −4 and 10 −2 .

Implementation of the DBIM-FISTA with a GPU-Based FDTD Forward Solver
We have used CUDA toolkit to accelerate our FDTD forward solver in GPU similar to [49], which focused however on modeling RF wave interactions in high-field MRI. The FDTD algorithm consists of electric and magnetic field updates, and each can be viewed as a kernel on GPU. The kernel, executed as a grid in GPU, can be divided into multiple blocks of threads, whose number can be adjusted based on the GPU platform. We can therefore increase the computational efficiency significantly by computing the blocks of threads simultaneously. Our algorithm is designed to use the high performance of GPU without changing the MATLAB based DBIM code, to take advantage of MATLAB's capability for matrix computing. To this end, we have implemented the 3-D FDTD algorithm in C++ with CUDA and then incorporated it with MEX functions in MATLAB, which is also used to implement the FISTA inverse solver. The environment of our implementation is MATLAB 2020 (Mathworks Inc., Natick, MA, USA), CUDA 10.2 (Nvidia Co., Santa Clara, CA, USA) and Visual Studio 2015 (Microsoft Co., Redmond, WA, USA) with Intel Xeon CPU E5-2640 v3 @ 2.60 GHz and GPU Tesla K20c with 5 GB memory.
The maximum number of threads for each block is 1024 and we use a 2-D block with size 32 × 32 = 1024 to make full use of each block. The number of threads for each block should be set to a multiple of 32 as the working unit in the GPU is a warp, which is a set of 32 threads [50]. Similarly, we use a 2-D grid GPU with size Imax 32 × Jmax 32 for each kernel, E i (i = x, y or z) for example, and the loop inside each kernel for updating the fields is in the following structure shown in Algorithm 1, where Imax, Jmax and Kmax represent number of voxels for the x, y, and z dimension, respectively.

Algorithm 1 Kernel of updating electric field
The GPU Tesla K20c has 13 streaming processors (SMs) and for each SM there are at most 2048 threads. Thus 13 × 2048 threads can be run simultaneously. The flowchart of the 3-D FDTD algorithm is shown in Figure 1. To reduce computational times, most of the FDTD functions are run on GPU (kernels), including the source update. At each DBIM iteration, the FDTD variables are first allocated space on GPU after initialization by MEX functions. Thus all FDTD-related variables are on GPU memory. Then the FDTD part on GPU starts to simulate the wave propagation for each transmitter antenna, and the calculated data to be used in DBIM is transferred back to MATLAB via MEX functions, which can be used directly by the MATLAB code. Finally the ill-conditioned linear system is built by (3) and solved by FISTA to obtain the reconstructed values.
The discretization size for FDTD is defined as where λ f is the wavelength and n d is the discretization step. To ensure accuracy and reduce numerical dispersion, n d > 10 is required. While ∆h is 30 mm in air, its value inside the brain is defined by a much smaller wavelength (λ f = 300 √ ε r ), thus requiring a greater n d and a smaller ∆h. We have therefore selected a cubic grid voxel size of 2 mm for each dimension in all our 2-D and 3-D FDTD simulations. This value can balance accuracy and computational burden. For example, using a 2 mm vs. a 1 mm cubic voxel side results in field values with a mean relative difference of 10 −3 but requires an 8 times bigger grid. The time step ∆t defined as where c is the speed of EM wave, will also increase for a higher resolution model and thus lead to more iterations. In our simulations, the iteration number ranges between 1000-2000.
To compare computational burden vs. grid resolution, we calculated run times for a two antenna system using 2 mm, 1.5 mm and 1 mm size cubic voxels. These are acceptable resolutions for a grid with physical volume of about 300 mm × 300 mm × 200 mm, which is required for our inverse problem with the numerical head phantom. Resolutions of 0.5 mm or lower exceed the GPU's memory requirements, and they also increase errors and instabilities for the inverse model [51]. The size and run times corresponding to these different grids when 16 antennas are used as transmitters and receivers are shown in Table 1. As this implementation is specific for the inverse problem at hand, it is more than 30% faster than our previously used codes implemented with commercial software package Acceleware, which is designed for general GPU-based FDTD simulations. Moreover, the mean absolute error of the fields received by each antenna calculated by Acceleware and our new code is less than 10 −8 , which shows our implementation is as accurate as Acceleware's. The running time for FISTA for the above-mentioned case with 2 mm resolution where the size of A is 240 × 449,949 at each DBIM iteration, is around 5-15 s depending on the number of FISTA iterations. This number typically ranges from 20-100, i.e., for each iteration the average running time is around 0.2 s.
As noted earlier, we have also employed the 2-D DBIM-FISTA algorithm [40] in this paper to compare performance with the 3-D implementation. The transverse magnetic (TM) mode of the EM wave is considered for the 2-D problem. The 2-D code is implemented in MATLAB without the need of GPU acceleration for the FDTD forward model, which requires approximately 1 s for each of the eight transmitting antennas and a total of 10-20 s for each DBIM iteration.

Numerical and Experimental Phantoms
The original MRI-derived Zubal head phantom [45] comprises 256 × 256 × 128 voxels with size of 1.1 mm × 1.1 mm × 1.4 mm. First, we imported the phantom into CST Microwave Studio to obtain the numerical "measured" data (S-parameters) of the full-wave 3-D interaction of the experimental antenna array with the phantom, which is used to test the 3-D algorithm. We then resized the model to fit the FDTD's grid of 2 mm for each side. We have also transformed the original Zubal head model from dozens of materials into an eight-material head model. The permittivity of the materials used in CST was obtained from the IT'IS foundation database [52]. These data were also used to develop single-pole Debye models for our FDTD code by curve fitting. As these two approaches are not identical, there are discrepancies between the CST and FDTD Debye models.
The 3-D structure of the Zubal head phantoms in CST is shown in Figure 2. As for the 2-D cases, a slice of the Zubal head phantom inside the brain is used in FDTD as shown in Figure 2d. The model includes eight layers with tissue types, color codes and respective Debye parameters (for our FDTD models) shown in Table 2. The head models are surrounded by 90-10% glycerol-water mixture, which we have used previously as immersion liquid in our recent experiment work [53,54]. The permittivity of the glycerol-water mixture is ε r ≈ 15.9 − 14.2j at 1.0 GHz. The Debye parameters ε ∞ = 6.566, ∆ε = 16.86 and σ s = 0.3231. The relaxation time τ is fixed as 0.14288 ns for all the materials.
Starting from this head model, we have studied various different imaging scenarios, which will be discussed in next section. Moreover, we have used a simplified experimental head phantom [54] shown in Figure 3 to validate our algorithm with more realistic data. The phantom is made of a 3D-printed plastic mould and is filled with an "average brain tissue" material with ε r ≈ 41.6 − 5.9j at 1.0 GHz. We should note that the permittivity of the brain tissue differs at different positions, i.e., the permittivity near the different surface, in the inner parts or near the bottom has different values [54]. The brain phantom is immersed in the 90-10% glycerol-water mixture inside an imaging tank, which is made of acrylic and is surrounded by absorbing material ECCOSORB MCS covered by a metallic shield. A printed monopole triangular patch antenna is used, as shown in Figure 3c, which operates well in the frequency range 0.5-2 GHz [55]. An eight-antenna array inside the tank captures MWT data with the help of an eight-port VNA system embedded in the Keysight M9019A PXIe Chassis connected to a desktop, shown in Figure 3d,e. Further details of the experimental system and the process of making the phantom can be found in [54].
The number of required antennas for a 2-D microwave imaging problem has been defined in [56] using the approximation, where α is the radius of the investigated region, while for 3-D problems the required number is slightly larger. This requirement may be difficult to satisfy in practice due to experimental limitations; for example, our experimental setup uses eight antenna hosted by an eight-port VNA system, and the same setup is used for each antenna ring in our CST scenarios. We note that the antennas are modelled as point sources in our algorithms, and the received signals by the antennas are normalized with respect to the source. The calibrated scattered fields are a standard total field calibration procedure [57] as, where S b and S t represent the S-parameters for the background case and for the unknown case to be imaged, as measured experimentally or modeled in CST. In (36), E b represents the background incident signals calculated by the FDTD forward model, and (i, j) refers to the transmitter-receiver pair.

Validation of the Proposed DBIM-FISTA Algorithm
This subsection presents results from simple imaging scenarios where only a strokelike target is unknown, with the aim to validate our 3-D DBIM-FISTA implementation. We consider both CST-calculated and experimental data and compare our results with 2-D scenarios where the target is placed at the same height as the antenna ring used by the 2-D algorithm. The 2-D DBIM-FISTA can perform well for these scenarios, as it deals with a detectable scattered signal and an inverse problem of much fewer unknowns. Although the number of iterations required for convergence will be different for each 2-D and 3-D imaging scenario, we have selected a fixed number of 20 iterations for this comparison.
To study the impact of keeping the number of iterations fixed, we have compared results with a much greater number of iterations for "Case I" presented in Section 3.1.1. We note that all of our reconstructions estimate 3D volume distributions, but we only show the 2D axial slice results, as we have limited array coverage along the z-axis, and we also want to compare directly 2D and 3D simulations along the axial planes.
The true dielectric constant values (ε r ) of the numerical and experimental phantoms are shown in Figure 4.
The root mean square error (RMSE) of the xy-slice reconstruction contrast is used to compare the 2-D and 3-D reconstructions, defined as where N r is the number of voxels inside the reconstruction region, p is the index of the voxel, and e ε is the difference between the real part of relative permittivity of reconstruction ε r (p) and true values of ε gt (p), i.e., The max error of the xy-slice reconstruction contrast is defined as and the relative error is defined as the relative norm between the residual errors at the K-th and first DBIM iteration,

Reconstructions with CST Data
"Case I" considers the numerical model in Figure 2 immersed inside "infinite" 90-10% glycerol-water mixture. A cylindrical blood target centered at O tg = (20 mm, 20 mm) with ε r = 61.1 − 28.4j, radius ρ = 15 mm and height h = 30 mm is inserted into the Zubal head phantom. Eight antennas are placed in an elliptical array configuration with semi-major and semi-minor axes equal to 100 mm and 85 mm, respectively. Two CST simulations are performed to obtain the scattered field data, with and without the target, to which we refer as WT and NT, respectively. The initial guess of FDTD model for reconstruction is chosen to be similar to the NT case CST model. The reconstruction area inside the brain is a cubical volume with height along the Z-axis between [10, 70] mm for the 3-D model. The 3-D reconstruction results using the simplified Green's function are shown in Figure 5, where the top shows the reconstructed relative permittivity (ε r ) and the bottom shows the contrast (δ ε ) due to the target.
We have also performed reconstructions using the vectorial Green's function (13), and the reconstruction results of (δ ε ) are shown in Figure 6 with errors RMSE = 3.65, e max = 20.42, RE = 0.27, which are close to the reconstructed values using the simplified Green's function shown in the first row of Table 3.  This comparison suggests that the z-only approximation for the Green's function does not affect the accuracy of the results. We thus argue that this formulation is sufficient to achieve similar accuracy with the more complete vector formulation of (13). To re-enforce this argument, we have performed another comparison of the two formulations for a more challenging imaging scenario in Section 3.2.
The target is also detected in 2-D as shown in Figure 7, but the reconstructed values are lower. Inaccuracies in these 2-D and 3-D results can be attributed to the complex brain structure which leads to a highly non-linear scattering problem, and to the mismatch between the CST model producing the data and the FDTD forward model of our DBIM-FISTA algorithm. For example, our FDTD solver models the antennas as point sources to avoid the additional computational complexity and much finer resolution required for modeling the full antenna structure, which leads to a model mismatch from the full antenna CST model. Moreover, the models in FDTD and CST have discrepancies in the Debye material properties due to their different computational environments. This includes errors in the head boundary between the CST (or experimental model) and the FDTD solver, due to discretization errors and the coarser resolution used in FDTD.
To assess whether a fixed number of 20 iterations for both 2-D and 3-D algorithms leads to fair comparison, we have run Case I with a large number of iterations which can ensure that the residual error RE converges to almost a fixed value (see plots (a) and (b) in Figure 8). The numbers of iterations and errors for these cases are given in Table 3, and the resulting images are shown in Figure 8.
Relative to the reconstruction results with 20 DBIM iterations, reconstructions after a much greater number of iterations (120 in 3-D, 200 in 2-D) lead to rather modest reductions in the errors in Table 3, despite a significant decrease in the data residual RE. Importantly, the ratio between the 3-D and 2-D RMSE errors for the increased iterations is similar to the 20 iterations case, thereby allowing us to draw conclusions for 2-D vs. 3-D results using this limited fixed number. These observations and the results in Figure 8 suggest that using a fixed number of 20 iterations is sufficient for our comparison, although it does not represent complete convergence of the algorithms.

Reconstructions with Experimental Data
We conducted an experiment with the antenna array surrounding the lower half of the phantom shown in Figure 3. We inserted a cylindrical blood target with radius ρ = 15 mm and permittivity ε r ≈ 67.3 − 9.3j at 1.0 GHz into the phantom, centered at O tg = (−30 mm, 30 mm) along the horizontal axes. Reconstructions in 3-D and for the 2-D slice defined by the antenna centers are shown in Figure 9, where the Debye parameters of the brain material ε ∞ = 20, ∆ε = 20 and σ s = 0.147 are used. The 2-D reconstruction has less artifacts than the 3-D image along the same x-y slice, as it solves an inverse problem with a lot fewer unknowns. The y-z and x-z slices from the 3-D reconstructions show that the bottom part of the cylindrical target is detected more clearly than the upper half, as the eight-antenna ring is placed in the lower half. These images also suggest that reflections from the plastic container have created artifacts in the 3-D reconstruction images. Overall, both 2-D and 3-D algorithms have detected the target at the right position, albeit with artifacts that are more pronounced in 3-D for the axial slice. We have also performed reconstructions when the antenna ring is placed at different heights, and the results become worse and even fail to detect the target when the antenna ring is near the upper part of the head phantom due to that the curvature of the surface reflects signals into the air which are not captured by the antenna ring.

Imaging Performance with Limited Data
This subsection focuses on reconstructions of more challenging imaging scenarios. These include imaging with limited prior information, as well as detecting a target positioned in-between the antenna rings of a 3-D array. We compare 2-D and 3-D reconstructions for these cases to examine possible benefits of using the 3-D algorithm. Similar to the previous results, we have performed twenty DBIM iterations in all our reconstructions below.

Reconstructions with Limited Prior Information
To investigate a scenario of limited prior information where only the boundary of the head is known, the model of Figure 2 was filled with white matter only for the NT case. Moreover, taking into account that the dielectric properties of gray and white matter are not very different (see Table 2), gray matter is replaced with white matter in Figure 2 to reduce the model complexity [41]. The resulting head model is then used in two WT cases (WT1 and WT2), "Case II.1" and "Case II.2", which differ by the presence or absence of the CSF layer with ε r ≈ 68.4 − 44.9j. Cross-sectional views of the 3-D models for these cases are shown in Figure 10 while the true values of (ε r ) is shown in Figure 11.
3-D and 2-D reconstruction results are shown in Figures 12 and 13, respectively, with errors shown in Table 4.

Reconstructions of Small Target at an Offset Height Using a Headband
We consider another scenario in which a small target is placed in-between a two-ring array surrounding the model of Figure 2. The head phantom comprises white matter and includes a target at the same x-y position as previously, but with a smaller radius ρ = 10 mm and height h = 20 mm for the WT case. As the immersion liquid cannot extend infinitely in a practical scenario, we consider a headband of finite dimensions filled with the immersion liquid and surrounded by air. The setup is shown in Figure 14. The two rings are placed with an offset in the x-y plane, to obtain information from more angles and reduce coupling. The target covers an area between the x-y planes of the two rings only, as shown in Figure 14a. We have considered different ways of applying our 2-D imaging algorithm to the data by the two-ring array by using: (1) the bottom ring data, (2) the top ring data, and (3) combined data from both rings as if they were on the same plane, effectively creating a sixteen-antenna array for the slice reconstructed in 2-D, which was selected as the slice of the target center.
The head model for the cases considered in this subsection (Cases III.1, 2 and 3) is shown in Figure 15, while the headband used in each case is different.
We note that the "homogeneous white-matter" head model of Figure 15 is unrealistic, but it was selected to focus on comparing 3-D with 2-D results of more realistic arrays and smaller targets. To compare the errors of the target domain more clearly, we have also defined an RMSE in the selected target slice as RMSE-T. Despite including the interface in the 3-D model, the reconstructed images suffer from strong artifacts near the interface, which suggest that the mismatch between the CST and FDTD model is significant. The target is reconstructed to some extent, but it is difficult to detect it with certainty. The 2-D reconstructions fail to detect the target completely regardless of whether we use data from the top, bottom, or both rings. This is not surprising given that the target is not aligned with any of these rings, and that the interface between air and the glycerol-water immersion liquid cannot be fully modeled in 2-D.
To examine the impact of this interface and improve detection performance, "Case III.2" considers the same headband as Case III.1, which is now surrounded by an additional layer of absorbing material. The headband layers are as follows: glycerol-water mixture, plastic, absorbing material (ECCOSORB MCS), and a metallic shield. Finally, "Case III.3" uses a a larger size (ρ = 140 mm and h = 120 mm) headband with the same materials as Case III.1. As with the absorbers in Case III.2, increasing the distance between the antennas and the interface with air can reduce the resulting reflections.
Results Similarly, 3-D and 2-D reconstruction results for Case III.3 are presented in Figure 19 and Figure 20. Detection is improved from Case III.1, as the distance from the interface has been increased, but is less accurate than Case III.2 where the absorbing material was added. Similar to the other two cases, the 2-D imaging results in Figure 20 suggest that the algorithm fails to detect the target in all cases.  Figure 14 (with an increased distance between the array and the interface with air) is used to obtain the data.
The RMSE, RMSE-T, e max and RE of Case III.1, 2 and 3 are shown in Table 5, where the 3-D reconstructions have a lower RMSE-T for all the three cases.  Figure 14 (with an increased distance between the array and the interface with air) is used to obtain the data. Reconstructions of other cases with the vectorial Green's function (13) have also been performed and Case III.2 is shown in Figure 21. The relative errors are RMSE = 2.37, RMSE − T = 20.89, e max = 21.08, RE = 0.62, which again shows similar reconstruction quality with the simplified Green's function even when the antenna array is placed at an offset height. This comparison confirms that the simplified Green's function does not affect the 3-D reconstruction quality significantly, leading to equally accurate results for the scenarios considered in this paper. Finally, we note that we have also performed reconstructions for additional scenarios related to Cases III.1-3, i.e., headbands with radii ρ = 135 and 140 mm, with and without the absorbing layers. As expected, when the interface distance increases, detection accuracy and reconstruction quality improves. Moreover, the absorbing materials improve reconstructions further. Importantly, the 3-D imaging outperforms its 2-D counterpart in all these cases.

Discussion
We have developed and validated a computationally efficient 3-D DBIM algorithm for microwave head imaging. Our DBIM implementation uses the FISTA solver for the linear inverse problem, which has shown advantages over traditional CGLS solvers in our previous work [58]. The 3-D DBIM-FISTA algorithm relies on an in-house 3-D FDTD forward solver implemented on GPU, which is equally accurate but runs considerably faster than previous implementations with commercial software Acceleware. Our implementation combined this 3-D GPU-based FDTD solver with the inversion code on MATLAB via MEX functions.
Inspired by microwave head imaging applications, we used CST-calculated as well as experimental data to validate the 3-D DBIM-FISTA algorithm. We implemented a simplified 3-D inverse algorithm that considers a scalar Green's function for linearly polarized antennas, such as the printed monopoles used by our system. Importantly, we showed that this approximation does not lead to worse performance than a vectorial Green's function implementation for the considered head imaging scenarios. Our numerical studies employed different numerical head models in CST based on the Zubal head phantom. Our experimental study examined a simpler, homogeneous head model, with the purpose to validate the 3-D implementation rather than examining more complex head models.
Our results showed that reconstructions are of similar quality to those produced by a 2-D version of the algorithm for cases where the data processed by the 2-D imaging algorithm is of sufficient quality. To show that this may not always be possible, we considered cases where the target was not centered at the same transverse plane as the antenna ring and reflections from the interface of a headband with air were significant. Comparison of 3-D and 2-D reconstructions for these imaging scenarios showcased the advantages of imaging in 3-D. Moreover, these cases also demonstrated that terminating the imaging headband with absorbing material can improve drastically the 3-D array's imaging performance.
Future work will focus further on reconstructions with experimental data inspired by the outcomes of this study. In particular, using prior information to improve reconstruction accuracy in realistic clinical settings requires further investigation. This prior information is not only required for dealing with the complexity of the brain's heterogeneous tissues, but also to estimate accurately the head boundary, which constitutes the reconstruction domain of our MWT approach. Moreover, the issue of patient movement should also be investigated in clinical settings. Motion artefacts are an important issue in imaging techniques which require long data acquisition times such as MRI, and a MWT system must be designed to minimise their impact by ensuring very short data acquisition times.

Conflicts of Interest:
The authors declare no conflict of interest. The funders had no role in the design of the study; in the collection, analyses, or interpretation of data; in the writing of the manuscript, or in the decision to publish the results.