## 2 Materials and Methods  
### 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,  
\[\mathbf{\tilde{E}_{s}}(\mathbf{r_{n}},\mathbf{r_{m}})=\mathbf{\tilde{E}_{t}}(\mathbf{r_{n} },\mathbf{r_{m}})-\mathbf{\tilde{E}_{b}}(\mathbf{r_{n}},\mathbf{r_{m}})=\omega^{2}\mu_{0} \varepsilon_{0}\int_{V}\mathbf{\tilde{G}_{b}}(\mathbf{r_{n}},\mathbf{r})\mathbf{\tilde {E}_{t}}(\mathbf{r},\mathbf{r_{m}})\delta_{\varepsilon}(\mathbf{r})d\mathbf{r}, \tag{1}\]  
where \(\mathbf{\tilde{E}_{t}}\), \(\mathbf{\tilde{E}_{s}}\), \(\mathbf{\tilde{E}_{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 \(\mathbf{E_{b}}\) under the Born approximation. The vectors \(\mathbf{r_{m}}\) and \(\mathbf{r_{n}}\) denote the transmitting and receiving antenna locations, \(\omega\) is the angular frequency, \(\mu_{0}\) and \(\varepsilon_{0}\) are the permeability and permittivity of free space, and \(\mathbf{\tilde{G}_{b}}\) is the dyadic Green's function for the background medium. The difference \(\delta_{\varepsilon}\) between the relative complex permittivities of reconstruction \(\varepsilon(\mathbf{r})\) and background \(\varepsilon_{b}(\mathbf{r})\) is defined as  
\[\delta_{\varepsilon}(\mathbf{r})=\varepsilon_{r}(\mathbf{r})-\varepsilon_{b}(\mathbf{r}). \tag{2}\]  
Assuming an equal discretization step \(\Delta t\) for each dimension and constant permittivity inside each voxel, (1) can be written as,  
\[\mathbf{\tilde{E}_{s}}(\mathbf{r_{n}},\mathbf{r_{m}})=k_{0}^{2}\Delta t^{3}\sum_{r\in V }\mathbf{\tilde{G}_{b}}(\mathbf{r_{n}},\mathbf{r})\mathbf{\tilde{E}_{b}}(\mathbf{r},\mathbf{r _{m}})\delta_{\varepsilon}(\mathbf{r}), \tag{3}\]where  
\[\bar{\mathbf{E}}_{\mathbf{s}}(\mathbf{r}_{n},\mathbf{r}_{m})= \begin{bmatrix}E_{s}^{x}(\mathbf{r }_{n},\mathbf{r}_{m})\\ E_{s}^{y}(\mathbf{r}_{n},\mathbf{r}_{m})\\ E_{s}^{z}(\mathbf{r}_{n},\mathbf{r}_{m})\end{bmatrix} ,\ \bar{\mathbf{E}}_{\mathbf{b}}(\mathbf{r},\mathbf{r}_{m})= \begin{bmatrix}E_{x}(\mathbf{r},\mathbf{r}_{m})\\ E_{y}(\mathbf{r},\mathbf{r}_{m})\\ E_{z}(\mathbf{r},\mathbf{r}_{m})\end{bmatrix} , \tag{4}\]  
with \(E_{x}\), \(E_{y}\), \(E_{z}\) representing the \(x\)-, \(y\)-, \(z\)-direction fields of the background field \(\bar{\mathbf{E}}_{\mathbf{b}}\) respectively, and \(k_{0}\) is the free-space wavenumber. The relation between the background field \(\bar{\mathbf{E}}_{\mathbf{b}}\) and the Green's function \(\bar{\mathbf{G}}_{\mathbf{b}}\) due to a source with current density \(\bar{\mathbf{J}}_{0}\) at \(\mathbf{r}_{0}\) holds as [46],  
\[\bar{\mathbf{E}}_{\mathbf{b}}(\mathbf{r},\mathbf{r}_{0}) =-j\omega\mu_{0}\int_{V}\bar{\mathbf{J}}_{0}(\mathbf{r}^{\prime}) \cdot\bar{\mathbf{G}}_{\mathbf{b}}(\mathbf{r}^{\prime},\mathbf{r})d\mathbf{r}^{\prime} \tag{5}\] \[ \begin{bmatrix}E_{x}(\mathbf{r},\mathbf{r}_{0})\\ E_{y}(\mathbf{r},\mathbf{r}_{0})\\ E_{z}(\mathbf{r},\mathbf{r}_{0})\end{bmatrix}  =-j\omega\mu_{0}\int_{V} \begin{bmatrix}J_{x}(\mathbf{r}^{\prime})\\ J_{y}(\mathbf{r}^{\prime})\\ J_{z}(\mathbf{r}^{\prime})\end{bmatrix} \cdot \begin{bmatrix}G_{xx}(\mathbf{r}^{\prime}, \mathbf{r}),G_{xy}(\mathbf{r}^{\prime},\mathbf{r}),G_{xz}(\mathbf{r}^{\prime},\mathbf{r})\\ G_{yx}(\mathbf{r}^{\prime},\mathbf{r}),G_{yy}(\mathbf{r}^{\prime},\mathbf{r}),G_{yz}(\mathbf{r}^{ \prime},\mathbf{r})\\ G_{zx}(\mathbf{r}^{\prime},\mathbf{r}),G_{zy}(\mathbf{r}^{\prime},\mathbf{r}),G_{zz}(\mathbf{r}^{ \prime},\mathbf{r}),G_{zz}(\mathbf{r}^{\prime},\mathbf{r})\end{bmatrix}