SEAGLE: Sparsity-Driven Image Reconstruction under Multiple Scattering

TL;DR

SEAGLE introduces a novel nonlinear model and an accelerated gradient method for robust electromagnetic imaging under multiple scattering, enabling accurate reconstruction even in complex scattering scenarios.

Contribution

It presents a new series expansion model with accelerated gradient descent for improved imaging under multiple scattering, incorporating a total variation regularizer.

Findings

Effective in handling multiple scattering scenarios

Validated on simulated and real data

Achieves fast and accurate image reconstruction

Abstract

Multiple scattering of an electromagnetic wave as it passes through an object is a fundamental problem that limits the performance of current imaging systems. In this paper, we describe a new technique-called Series Expansion with Accelerated Gradient Descent on Lippmann-Schwinger Equation (SEAGLE)-for robust imaging under multiple scattering based on a combination of a new nonlinear forward model and a total variation (TV) regularizer. The proposed forward model can account for multiple scattering, which makes it advantageous in applications where linear models are inaccurate. Specifically, it corresponds to a series expansion of the scattered wave with an accelerated-gradient method. This expansion guarantees the convergence even for strongly scattering objects. One of our key insights is that it is possible to obtain an explicit formula for computing the gradient of our nonlinear forward model with respect to the unknown object, thus enabling fast image reconstruction with the state-of-the-art fast iterative shrinkage/thresholding algorithm (FISTA). The proposed method is validated on both simulated and experimentally measured data.

Figures

Equations

$$u(\mathbf{x})=u_{\text{in}}(\mathbf{x})+{\int }_{\Omega }g(\mathbf{x}-{\mathbf{x}}^{\prime })f({\mathbf{x}}^{\prime })u({\mathbf{x}}^{\prime })\mathrm{d}{\mathbf{x}}^{\prime },(\mathbf{x}\in \Omega )$$

$$g(\mathbf{x})\triangleq \{\begin{array}{ll}\frac{\mathrm{j}}{4}{H}_0^{(1)}(k_b\Vert \mathbf{x}{\Vert }_{{\ell }_2})& \text{in 2D}\\ \frac{{\mathrm{e}}^{\mathrm{j}{k}_b\Vert \mathbf{x}{\Vert }_{{\ell }_2}}}{4\pi \Vert \mathbf{x}{\Vert }_{{\ell }_2}}& \text{in 3D},\end{array}$$

$$\left({\nabla }^2+k_b^2\text{I}\right)g(\mathbf{x})=-\delta (\mathbf{x}),$$

$$u_{\text{sc}}(\mathbf{x})={\int }_{\Omega }g(\mathbf{x}-{\mathbf{x}}^{\prime })f({\mathbf{x}}^{\prime })u({\mathbf{x}}^{\prime })\mathrm{d}{\mathbf{x}}^{\prime }.(\mathbf{x}\in \Gamma )$$

$$\mathbf{y}=\mathbf{H}(\mathbf{u}\bullet \mathbf{f})+\mathbf{e}$$

$$\mathbf{u}={\mathbf{u}}_{\text{in}}+\mathbf{G}(\mathbf{u}\bullet \mathbf{f}),$$

$$\hat{\mathbf{u}}(\mathbf{f})\triangleq {\mathrm{argmin}}_{\mathbf{u}\in {\mathbb{C}}^N}\left\{\mathcal{S}(\mathbf{u})\right\}$$

$$\text{with}\mathcal{S}(\mathbf{u})\triangleq \frac{1}{2}\Vert \mathbf{A}\mathbf{u}-{\mathbf{u}}_{\text{in}}{\Vert }_{{\ell }_2}^2,$$

$$\nabla \mathcal{S}(\mathbf{u})={\mathbf{A}}^{\text{H}}(\mathbf{A}\mathbf{u}-{\mathbf{u}}_{\text{in}}).$$

$$t_k\leftarrow \frac{1}{2}\left(1+\sqrt{1+t_{k-1}^2}\right)$$

$${\mathbf{s}}^k\leftarrow {\mathbf{u}}^{k-1}+((t_{k-1}-1)/t_k)({\mathbf{u}}^{k-1}-{\mathbf{u}}^{k-2})$$

$${\mathbf{u}}^k\leftarrow {\mathbf{s}}^k-\nu {\mathbf{A}}^{\text{H}}(\mathbf{A}{\mathbf{s}}^k-{\mathbf{u}}_{\text{in}}),$$

$$\hat{\mathbf{f}}={\mathrm{argmin}}_{\mathbf{f}\in {\mathbb{R}}^N}\left\{\mathcal{D}(\mathbf{f})+\mathcal{R}(\mathbf{f})\right\},$$

$$\mathcal{D}(\mathbf{f})\triangleq \frac{1}{2}\Vert \mathbf{y}-\mathbf{z}(\mathbf{f}){\Vert }_{{\ell }_2}^2\text{and}$$

$$\mathcal{R}(\mathbf{f})\triangleq \tau \sum _{n=1}^N\Vert [\mathbf{D}\mathbf{f}{]}_n{\Vert }_{{\ell }_2}=\tau \sum _{n=1}^N\sqrt{\sum _{d=1}^D|[{\mathbf{D}}_d\mathbf{f}{]}_n{|}^2}.$$

$${\mathbf{f}}^t\leftarrow {\mathrm{prox}}_{\gamma \mathcal{R}}\left({\mathbf{f}}^{t-1}-\gamma \nabla \mathcal{D}({\mathbf{f}}^{t-1})\right),$$

$$\nabla \mathcal{D}(\mathbf{f})=\mathrm{Re}\left\{{\left[\frac{\partial \mathbf{z}(\mathbf{f})}{\partial \mathbf{f}}\right]}^{\text{H}}\left(\mathbf{z}(\mathbf{f})-\mathbf{y}\right)\right\},$$

$$\text{normalized error}\triangleq \frac{\Vert \hat{\mathbf{u}}-{\mathbf{u}}_{\text{true}}{\Vert }_{{\ell }_2}^2}{\Vert {\mathbf{u}}_{\text{true}}{\Vert }_{{\ell }_2}^2},$$

$$\text{normalized data fit}\triangleq \frac{\mathcal{D}(\hat{\mathbf{f}})}{\mathcal{D}(0)}=\frac{\Vert \mathbf{z}(\hat{\mathbf{f}})-\mathbf{y}{\Vert }_{{\ell }_2}^2}{\Vert \mathbf{y}{\Vert }_{{\ell }_2}^2}$$

$$\text{normalized reconst. error}\triangleq \frac{\Vert \hat{\mathbf{f}}-\mathbf{f}{\Vert }_2^{{\ell }_2}}{\Vert \mathbf{f}{\Vert }_2^{{\ell }_2}},$$

$$\text{SNR (dB) }\triangleq 10{\log }_{10}\left(\frac{\Vert {\mathbf{f}}_{\text{ref}}{\Vert }_{{\ell }_2}^2}{\Vert \hat{\mathbf{f}}-{\mathbf{f}}_{\text{ref}}{\Vert }_{{\ell }_2}^2}\right),$$

$$\mathcal{R}(\mathbf{f})\triangleq \sum _{n=1}^N\Vert [\mathbf{D}\mathbf{f}{]}_n{\Vert }_{{\ell }_1}=\sum _{n=1}^N\sum _{d=1}^D|[{\mathbf{D}}_d\mathbf{f}{]}_n|$$

$$\mathcal{R}(\mathbf{f})\triangleq \sum _{n=1}^N\Vert [\mathbf{D}\mathbf{f}{]}_n{\Vert }_{{\ell }_2}=\sum _{n=1}^N\sqrt{\sum _{d=1}^D|[{\mathbf{D}}_d\mathbf{f}{]}_n{|}^2}$$

$${\mathrm{prox}}_{\tau \mathcal{R}}(\mathbf{z})\triangleq {\mathrm{argmin}}_{\mathbf{f}\in \mathcal{F}}\left\{\frac{1}{2}\Vert \mathbf{f}-\mathbf{z}{\Vert }_{{\ell }_2}^2+\tau \mathcal{R}(\mathbf{f})\right\}.$$

$$\hat{\mathbf{g}}={\mathrm{argmin}}_{\mathbf{g}\in \mathcal{G}}\left\{\mathcal{Q}(\mathbf{g})\right\},$$

$$\mathcal{Q}(\mathbf{g})$$

$$+\frac{1}{2}\Vert \mathbf{z}-\tau {\mathbf{D}}^{\text{T}}\mathbf{g}{\Vert }_{{\ell }_2}^2.$$

$${\mathbf{f}}^t={\mathrm{proj}}_{\mathcal{F}}(\mathbf{z}-\tau {\mathbf{D}}^{\text{T}}{\mathbf{g}}^t).$$

$$\mathcal{F}\triangleq \left\{\mathbf{f}\in {\mathbb{R}}^N:a\le f_n\le b,\forall n\in [1,\dots ,N]\right\},$$

$$[{\mathrm{proj}}_{\mathcal{F}}(\mathbf{f}){]}_n=\{\begin{array}{ll}a& \text{if }{f}_n<a\\ f_n& \text{if }a\le f_n\le b\\ b& \text{if }{f}_n>b,\end{array}$$

$$\mathcal{G}\triangleq \{\mathbf{g}\in {\mathbb{R}}^{N\times D}:\Vert [\mathbf{g}{]}_n{\Vert }_{{\ell }_{\infty }}\le 1,\forall n\in [1,\dots ,N]\}$$

$$[{\mathrm{proj}}_{\mathcal{G}}(\mathbf{g}){]}_n=\left(\begin{array}{c}\frac{[{\mathbf{g}}_x{]}_n}{\max (1,|[{\mathbf{g}}_x{]}_n|)}\\ \frac{[{\mathbf{g}}_y{]}_n}{\max (1,|[{\mathbf{g}}_y{]}_n|)}\end{array}\right),$$

$$\mathcal{G}\triangleq \{\mathbf{g}\in {\mathbb{R}}^{N\times D}:\Vert [\mathbf{g}{]}_n{\Vert }_{{\ell }_2}\le 1,\forall n\in [1,\dots ,N]\}$$

$$[{\mathrm{proj}}_{\mathcal{G}}(\mathbf{g}){]}_n=\frac{[\mathbf{g}{]}_n}{\max \left(1,\Vert [\mathbf{g}{]}_n{\Vert }_{{\ell }_2}\right)},$$

Full text

Introduction

Reconstruction of the spatial permittivity distribution of an unknown object from the measurements of the scattered wave is common in numerous applications. Traditional formulations of the problem are based on linearizing the relationship between the permittivity and the measured wave. For example, if one assumes a straight-ray propagation of waves, the phase of the transmitted wave can be interpreted as a line integral of the permittivity along the propagation direction. This approximation leads to an efficient reconstruction with the filtered back-projection algorithm [1]. Diffraction tomography is a more refined linear scattering model based on the first Born or Rytov approximations [2, 3, 4]. It establishes a Fourier transform-based relationship between the measured wave and the permittivity, and thus enables the reconstruction of the latter with a single numerical application of the inverse Fourier transform.

Recent research in compressive sensing and sparse signal processing has established that sparse regularization can dramatically improve the quality of reconstructed images, even when the amount of measured data is severely limited [5, 6]. This has popularized optimization-based inverse scattering approaches that combine linear forward models with regularizers that mitigate ill-posedness by promoting solutions that are sparse in a suitable transform domain. One class of such regularizers is total variation (TV) [7], which was demonstrated to improve imaging quality by substantially reducing undesired artifacts due to missing data [8, 9, 10].

The main advantage of imaging with linear forward models is that the reconstruction can be reduced to a convex optimization problem that is relatively simple and efficient to solve [11, 12, 13, 14]. However, such models are only accurate when scattering is weak, which is true when the object is relatively small or has a low permittivity contrast compared to the background [15]. When scattering is strong, multiple scattering of waves as they pass through the object limit the validity and accuracy of linear forward models. Multiple scattering is a fundamental problem in tomographic imaging and its complete resolution would enable imaging through strongly scattering objects such as human tissue [16]. As multiple scattering leads to nonlinear forward models, the challenge is in finding computationally tractable methods that can account for the nonlinearity while also accommodating sparsity-driven priors for compressive imaging. To that end, we propose a new method that efficiently combines a nonlinear forward model with the TV regularizer, thus enabling high-quality imaging from a limited number of measurements. Figure 1 provides an example, illustrating the quality of reconstruction for a strongly scattering object with limited data. In particular, the $320\times 320$ image in Figure 1(d) was obtained from only $16$ experimentally collected measurements at $3$ GHz from the public dataset [17].

Contributions

Our work builds upon prior work on inverse scattering that has been applied to a variety of practical problems in optical, microwave, and radar imaging. The proposed method—called Series Expansion with Accelerated Gradient Descent on Lippmann-Schwinger Equation (SEAGLE)—further extends this work with a new forward model that enables efficient sparsity-driven inversion. The performance of SEAGLE is robust to large permittivity contrasts, data reduction, and measurement noise.

The key contributions of this paper are summarized as follows:

• •

A new forward model that iteratively approximates the scattered wave with Nesterov’s accelerated-gradient method [18]. The key benefit is the guaranteed convergence even for objects with large permittivity contrasts.

• •

An efficient optimization strategy based on explicit evaluation of the gradient of our forward model with respect to the permittivity. The key benefit of this strategy is that it enables efficient sparsity-driven reconstruction for the nonlinear forward model.

• •

An extensive validation of our approach on analytical, simulated, and experimental data. The experimental data used in our evaluations comes from a public dataset [17], which enables easy comparisons with several other related approaches.

Related Work

The vast majority of current imaging systems—such as optical projection tomography (OPT), diffraction tomography, optical coherence tomography (OCT), digital holography, and subsurface radar—still rely on the linearization of object-wave interaction [19, 20, 8, 21, 22, 10, 23, 9, 24, 25, 26, 27, 28, 29, 30, 31]. However, early work in microwave imaging have shown the promise of accounting for the nonlinear nature of scattering [32, 33, 34, 35]. These have been extended by a large number of techniques incorporating the nonlinear nature of scattering. Several recent publications have reviewed these methods [36, 37, 38], which include conjugate gradient method (CGM) [39, 40], contrast source inversion method (CSIM) [41], hybrid method (HM) [37], beam propagation method (BPM) [42, 43, 44], recursive Born method (RBM) [45], and bounded inversion method (BIM) [46]. Some recent works have explored the idea of statistical modeling of scattering for imaging through diffusive or turbid media [47, 48, 49].

The key innovation in this paper is in the new forward model, which enables tractable optimization of a cost function that incorporates nonlinear scattering and a non-differentiable regularizer such as TV. As corroborated by experiments, this formulation enables fast, stable, and reliable convergence, even when working with limited amount of data. This paper extends our preliminary work [50] by including key mathematical derivations, as well as more extensive validation on experimentally collected data.

The experimental data used in this paper comes from a public dataset [17]. Several other methods have been tested on this dataset [51, 52, 53, 41, 54]. This enables qualitative evaluation of the performance of the proposed technique against some other algorithms based on nonlinear forward models. While many of the methods tested on the data use multi-frequency measurements, the results in this paper rely on a single frequency. However, the method uses the latest techniques in large-scale optimization with sparse regularization, which enables subsampling and leads to improvements in imaging performance. We expect the performance of the proposed method to improve further if multi-frequency measurements are incorporated.

Outline

Section 2 describes the problem formulation and our forward model. Our algorithmic strategy for reconstructing the image is described in Section 3. Experimental results validating our approach are presented in Section 4. Finally, Section 5 discusses our findings and concludes. We additionally provide three appendices that disclose technical aspects of the proposed method.

Forward model

The forward problem computes the scattered field given a distribution of inhomogeneous permittivity, while the inverse problem reconstructs this distribution. The model we propose here can be interpreted as a series expansion based on the iterates of the gradient method. This expansion can be made arbitrarily accurate and is stable for high permittivity objects. Additionally, it enables efficient computation of the gradient of the cost function, which is essential for fast image reconstruction. In this section we focus on introducing the forward model, leaving the gradient evaluation and inversion to the next section. Our derivations are for the scenario of a single illumination, but the generalization to an arbitrary number of illuminations is straightforward.

Problem formulation

Consider the scattering problem in Figure 2, where an object of the permittivity distribution $\epsilon({\bm{x}})$ in the bounded domain $\Omega\subseteq\mathbb{R}^{D}$, with $D\in\{2,3\}$, is immersed into a background medium of permittivity $\epsilon_{b}$, and illuminated with the incident electric field ${u_{\text{\tiny in}}}({\bm{x}})$. We assume that the incident field is monochromatic and coherent, and it is known inside $\Omega$ and at the locations of the sensors $\Gamma$. The result of object-wave interaction is measured at the location of the sensors as a scattered field ${u_{\text{\tiny sc}}}({\bm{x}})$. The scattering of light can be accurately described by the Lippmann-Schwinger equation inside the image domain [55]

[TABLE]

where $u({\bm{x}})={u_{\text{\tiny in}}}({\bm{x}})+{u_{\text{\tiny sc}}}({\bm{x}})$ is the total electric field, ${f({\bm{x}})\,\triangleq\,k^{2}(\epsilon({\bm{x}})-\epsilon_{b})}$ is the scattering potential, which is assumed to be real, and $k=2\pi/\lambda$ is the wavenumber in vacuum. The function $g({\bm{x}})$ is the Green’s function defined as

[TABLE]

where $k_{b}\,\triangleq\,k\sqrt{\epsilon_{b}}$ is the wavenumber of the background medium and $H_{0}^{(1)}$ is the zero-order Hankel function of the first kind. Note that the Green’s function satisfies the Helmholtz equation

[TABLE]

as well as the outgoing-wave boundary condition, and a time-dependence convention under which the physical electric field equals to $\mathrm{Re}\left\{u({\bm{x}}){\mathrm{e}}^{-{\mathrm{j}}\omega t}\right\}$. The knowledge of the total-field $u$ inside the image domain $\Omega$ enables the prediction of the scattered field at the sensor area

[TABLE]

The computation of the scattered wave is equivalent to solving (1) for $u({\bm{x}})$ inside the image and evaluating (3) for ${u_{\text{\tiny sc}}}({\bm{x}})$ at the sensor locations. Note that $u$ is present on both sides of (1) and that the relation between $u$ and the scattering potential $f$ is nonlinear. The first Born and the Rytov approximations [2, 3, 4] are linear models that replace $u$ in (3) with a suitable approximation that decouples the nonlinear dependence of $u$ on $f$. However, such linearization imposes a strong assumption that the object is weakly scattering, which makes the corresponding reconstruction methods not applicable to a large variety of imaging problems [15]. Our forward model described next is a fast nonlinear method for solving (1) that overcomes the weakly scattering object assumptions.

Algorithmic Expansion of the Scattered Waves

We separate the computation of the electric field into two parts: the total field $u({\bm{x}})$ in the image domain and the scattered field ${u_{\text{\tiny sc}}}({\bm{x}})$ at the sensors. The discretization and combination of (1) and (3) leads to the following matrix-vector description of the forward problem

[TABLE]

where ${\mathbf{f}}\in\mathbb{R}^{N}$ is the discretized scattering potential $f$, ${\mathbf{y}}\in\mathbb{C}^{M}$ is the measured scattered field ${u_{\text{\tiny sc}}}$ at $\Gamma$, ${\mathbf{u}_{\text{\tiny in}}}\in\mathbb{C}^{N}$ is the input field ${u_{\text{\tiny in}}}$ inside $\Omega$, ${\mathbf{H}}\in\mathbb{C}^{M\times N}$ is the discretization of the Green’s function at $\Gamma$, ${\mathbf{G}}\in\mathbb{C}^{N\times N}$ is the discretization of the Green’s function inside $\Omega$, the symbol $\bullet$ denotes a component-wise multiplication between two vectors, and ${\mathbf{e}}\in\mathbb{C}^{M}$ models the random noise at the measurements. Using the shorthand notation ${\mathbf{A}}\,\triangleq\,{\mathbf{I}}-{\mathbf{G}}{\mathrm{diag}}\{{\mathbf{f}}\}$, where ${\mathbf{I}}\in\mathbb{R}^{N\times N}$ is the identity matrix and ${\mathrm{diag}}\{\cdot\}$ is an operator that forms a diagonal matrix from its argument, we can represent the forward scattering in (4b) as a minimization problem

[TABLE]

where the matrix ${\mathbf{A}}$ is a function of ${\mathbf{f}}$. The gradient of $\mathcal{S}$ can be computed as

[TABLE]

Since (5) corresponds to the optimization of a differentiable function, it is possible to compute the total field ${\widehat{\mathbf{u}}}$ iteratively using Nesterov’s accelerated-gradient method [18]

[TABLE]

for $k=1,2,\dots,K$, where ${\mathbf{u}}^{0}={\mathbf{u}}^{-1}={\mathbf{u}_{\text{\tiny in}}}$, $q_{0}=1$, and $\nu>0$ is the step-size. At any moment, the predicted scattered field can be set to ${\mathbf{z}}^{k}={\mathbf{H}}({\mathbf{u}}^{k}\bullet{\mathbf{f}})$ with ${\mathbf{u}}^{k}$ given by (7c). Note that the resulting set of fields $\{{\mathbf{u}}^{k}\}_{k\in[1\dots K]}$ and $\{{\mathbf{z}}^{k}\}_{k\in[1\dots K]}$ can be interpreted as a $K$-term series expansion of the wave-fields inside the object and at the sensor locations, respectively. The full procedure for forward computation with a convenient adaptive step-size is summarized in Algorithm 1.

Inverse Problem

We now present the overall image reconstruction algorithm, based on the state-of-the-art fast iterative shrinkage/thresholding algorithm (FISTA) [56]. The application of FISTA to nonlinear inverse scattering is, however, nontrivial due to the requirement of the gradient of the scattered field with respect to the object. We solve this by providing an explicit formula, based on the error backpropagation method [57], for computing this gradient.

Image Reconstruction

We formulate image reconstruction as the following optimization problem

[TABLE]

The data-fidelity term $\mathcal{D}$ measures the discrepancy between the actual measurements ${\mathbf{y}}$ and the ones predicted by our scattering model ${\mathbf{z}}$. The function $\mathcal{R}$ is the isotropic TV regularizer with the parameter $\tau>0$ controlling the strength of the regularization, where ${\mathbf{D}}:\mathbb{R}^{N}\rightarrow\mathbb{R}^{N\times D}$ is the discrete gradient operator with matrix ${\mathbf{D}}_{d}$ denoting the finite difference operation along dimension $d$ (see Appendix 6 for details).

The image can then be formed iteratively using a first order method such as ISTA [58, 59, 60]

[TABLE]

for $t=1,2,3,\dots$ or its accelerated variant FISTA [56] summarized in Algorithm 2. Here, $\gamma>0$ is a step-size. The operator $\mathrm{prox}_{\gamma\mathcal{R}}$ denotes the proximity operator, and for isotropic TV it can be efficiently evaluated with an algorithm described in Appendix 6. Finally, an efficient implementation of the imaging algorithm requires the gradient of the data-fidelity term

[TABLE]

which can be evaluated using the algorithm summarized in Algorithm 3. The mathematical derivation of this algorithm is similar to that of the standard error backpropagation used in deep learning [57, 61], and its full derivation for multiple scattering is provided in Appendix 7. As illustrated in Figure 3, the proposed algorithm uses the recursive structure of accelerated-gradient expansion in order to obtain an efficient formula for evaluating the gradient of the data-fidelity term. While the theoretical convergence of FISTA is difficult to analyze for nonconvex functions, it is often used as a faster alternative to the standard gradient-based methods in the context deep learning and broader machine learning [62, 63, 64]. In fact, we observed that our method reliably converges and achieves excellent results on a wide array of problems, as reported in Section 4.

Experimental Evaluation

We now present the results of validating our method on analytically obtained scattering data for simple scenarios, scattering data obtained with a high-fidelity finite-difference time-domain (FDTD) simulator [65], and experimentally collected data from the public dataset [66].

Validation on analytic data

In the first set of experiments, we validated our forward model for two simple objects where analytic expressions of the scattered wave exist: a two-dimensional point source scattered by a cylinder, and a three-dimensional point source scattered by a sphere. The expressions are derived following the mathematical formalism in [67], which we review in Appendix 8 for completeness. As illustrated in Figures 4(a) and (b), in both cases, the objects have diameters equal to $6$ wavelengths. The wavelength is set to $74.9$ mm, the source is placed $1$ m away from the center of the objects, the grid size is set to $4.8$ mm ($6$ mm), and there are $250$ points ($128$ points) along each axis in 2D (3D). The contrast of an object is defined as $\max(|{\mathbf{f}}|)/k_{b}^{2}$. In Figures 4(b) and (d), we quantitatively evaluate the performance of our forward model with the normalized error defined as

[TABLE]

where ${\widehat{\mathbf{u}}}$ is the solution of (5) and ${\mathbf{u}}_{\text{\tiny true}}$ is the analytic expression. For comparison, we additionally provide the performance of the first Born (FB) and Rytov approximations at $20\%$ contrast. In Fig. 4(e) and 4(f), we demonstrate a visual comparison between the analytic expression and the result of our model. Overall, we observed that, by allowing for a large enough value of $K$, our forward model can match the analytically obtained field with arbitrarily high precision. The actual value of $K$ depends on the severity of multiple scattering and must be adapted on the basis of the application of interest. For example, we observed that for objects closely resembling biological samples, one generally requires $10\leq K\leq 30$.

Validation on FDTD simulations

We next validated the proposed technique for reconstructing the Shepp-Logan phantom in an ill-posed, strongly scattering, and compressive regime ($M=25\times 338$ and $N=250\times 250$). Specifically, we consider the setup in Fig. 2 where the scattered wave measurements are generated with an FDTD simulator. The object is of size $84.9$ cm $\times$ $113$ cm and has a contrast of $20\%$. We place two linear detectors on either side of the phantom at a distance of $95.9$ cm from the center of the object. Each detector has $169$ sensors placed with a spacing of $3.84$ cm. The transmitters are positioned on a line $48.0$ cm left to the left detector. They are spaced uniformly in azimuth with respect to the center of the phantom (every $5^{\circ}$ within $\pm 60^{\circ}$). We set up a $120$ cm $\times$ $120$ cm square area for reconstructing the object, with pixel size $0.479$ cm. The wavelength of the illuminating light is $7.49$ cm.

Figure 5 summarizes the performance of the proposed method, as well as that of the first Born and the Rytov approximations—the two baseline methods. All three approaches rely on FISTA with TV regularizer that was empirically set to $\tau=1.5\times 10^{-9}\|{\mathbf{y}}\|_{\ell_{2}}^{2}$ for the best performance. The order of our forward model is set to $K=120$, but Algorithm 1 may terminate earlier when the objective function (5) is below $\delta_{\text{\tiny tol}}=5\times 10^{-7}\|{\mathbf{u}_{\text{\tiny in}}}\|_{\ell_{2}}^{2}$. The figure additionally provides quantitative performance evaluation in terms of two quantities defined as

[TABLE]

where ${\mathbf{f}}$ and ${\widehat{\mathbf{f}}}$ denote the true and estimated object, respectively. Simulation results corroborate the benefit of using the proposed method for strongly scattering objects. It can be seen that, due to the ill-posed nature of the measurements, the reconstructed images suffer from missing frequency artifacts [68]. However, the proposed method is still able to accurately capture most features of the object while the linear methods fail to do so. Note also, that our method was initialized with the background value of the dielectric permittivity, $\epsilon_{\text{\tiny b}}=1$, and that it takes fewer than $50$ FISTA iterations for converging to a stationary point (see convergence plot in Figure 5(f)).

Validation on experimental data

We apply our method to three objects from the public dataset provided by the Fresnel institute [17]: FoamDielExtTM, FoamDielIntTM, and FoamTwinDielTM. These objects are placed in a region of size 15 cm $\times$ 15 cm at the center of a circular rim of radius 1.67 m and measured using 360 detectors and 8 transmitters evenly distributed on the rim. The number of transmitters is increased to 18 for FoamTwinDielTM and are also uniformly spaced. In all cases, only one transmitter is turned on at a time, while 241 detectors are used for each transmitter by excluding 119 detectors that are closest to the transmitter. While the full data contains multiple frequency measurements, we only use the data corresponding to the $3$ GHz. As before, we compare the result of our method with first Born and Rytov approximations that are also regularized with TV. We set the highest order of forward scattering to $K=200$ and the TV regularization to $\tau=0.25\times 10^{-8}\|{\mathbf{y}}\|_{\ell_{2}}^{2}$, and run the image formation algorithm for $40$ FISTA iterations. The reconstruction was initialized with the background value of the dielectric permittivity, which in this case corresponds to $\epsilon_{\text{\tiny b}}=1$.

Figure 6 summarizes the imaging results on the experimental data. The quantitative evaluation is performed using the same metrics defined in (12). The results show that our method successfully captures the shape of the samples, as well as the value of the permittivity. Both first Born and Rytov approximations underestimate the permittivity. One can also see that the data-fit error for both of the linear forward models remain high as iterations progress. On the other hand, the object reconstructed by the proposed method closely agrees with the measured data (see rightmost column in Figure 6).

Figure 7 illustrates the performance of our method when using a limited number of measured data. In particular, we consider the reconstruction of the same three objects, but reduce the number of measurements for each transmission using regular downsampling by factors of $2,4,8,16,32,64,$ and $128$. The full dataset consists of $8$ transmissions with $241$ measurements each; a factor of $128$ downsampling reduces to $8$ transmissions with $2$ measurements each. The size of the reconstructed image is set to $320\times 320$ pixels. The reconstruction performance is quantified as

[TABLE]

where ${\mathbf{f}}_{\text{\tiny ref}}$ is the reconstructed image with all the measured data (see Figure 6). The visual illustration is provided for FoamDielExtTM in Figure 1. This result highlights the stability of the proposed method to subsampling and experimental noise, even at highly nonlinear scattering scenarios.

Note that several other methods have been tested on this dataset [51, 52, 53, 41, 54]. Qualitative comparison of our results in Figure 6 with the results of those methods indicates that our approach achieves comparable performance using only a fraction of data (i.e., a single frequency with possible subsampling). Additionally, we observe a reliably stable and fast convergence starting from the initialization to the background permittivity, which is desirable in strongly scattering regimes.

Conclusion

In conclusion, we have demonstrated a nonconvex optimization technique for solving nonlinear inverse scattering problems. We have applied the technique to simulated and experimentally measured data in microwave frequencies. The scattering was modeled as a series expansion with Nesterov’s accelerated-gradient method. By structuring the expansion as a recursive feedforward network, we derived an error backpropagation formula for evaluating the gradient that can be used for iterative image reconstruction. The algorithm yields images of better quality than methods using linear forward models and is competitive with state-of-the-art inverse scattering approaches, tested on the same dataset. While the optimization problem is not convex, we have observed that the algorithm converges reliably within 100 iterations from a constant initialization of the permittivity. Our approach provides a promising framework for active correction of scattering in various applications and has the potential of significantly increasing the resolution and robustness when imaging strongly scattering objects.

Total variation minimization

In this section, we review ideas and algorithms behind total variation (TV) regularized reconstruction. This section reviews the algorithm that was originally proposed by Beck and Teboulle [14].

Two common variants of TV regularizers are the anisotropic TV regularizer

[TABLE]

and the isotropic TV regularizer

[TABLE]

Here, ${\mathbf{D}}:\mathbb{R}^{N}\rightarrow\mathbb{R}^{N\times D}$ is the discrete gradient operator and $D$ is the number of dimensions in the signal. The matrix ${\mathbf{D}}_{d}$ denotes the finite difference operator along the dimension $d$ with appropriate boundary conditions (periodization, Neumann boundary conditions, etc.).

The anisotropic TV regularizer (14) can be interpreted as a sparsity-promoting $\ell_{1}$-penalty on the image gradient, while its isotropic counterpart (15) as an $\ell_{1}$-penalty on the magnitudes of the image gradient, which can also be viewed as a penalty promoting joint-sparsity of the gradient components. By promoting signals with sparse gradients, TV minimization recovers images that are piecewise-smooth, which means that they consist of smooth regions separated by sharp edges. Isotropic TV regularizer (15) is rotation invariant, which makes it preferable in the context of image reconstruction.

FISTA summarized in Algorithm 2, is a popular algorithm for solving inverse problems. FISTA relies on the efficient evaluation of the gradient $\nabla\mathcal{D}$ and of the proximal operator

[TABLE]

The error backpropagation algorithm makes application of FISTA straightforward for solving inverse scattering with regularizers that admit closed form poximal operators such as the $\ell_{1}$-penalty. However, some regularizers including TV do not have closed form proximals and require an additional iterative algorithm for solving (16).

Our implementation solves (16) with the dual approach that was proposed by Beck and Teboulle [14]. The fast-gradient projection (FGP) approach, summarized in Algorithm 4, is based on iteratively solving the dual optimization problem

[TABLE]

where

[TABLE]

Given the dual iterate ${\mathbf{g}}^{t}$, the corresponding primal iterate can be computed as

[TABLE]

The operator $\mathrm{proj}_{\mathcal{F}}$ represents an orthogonal projection onto the convex set $\mathcal{F}$. For example, a projection onto $N$-dimensional cube

[TABLE]

with bounds $a$, $b>0$, is given by

[TABLE]

for all $n\in[1,\dots,N]$.

The set $\mathcal{G}\subseteq\mathbb{R}^{N\times D}$ in (17) depends on the variant of TV used for regularization. For anisotropic TV (14), the set corresponds to

[TABLE]

with the corresponding projection

[TABLE]

for all $n\in[1,\dots,N]$. Similarly, for isotropic TV (15), the set corresponds to

[TABLE]

with the corresponding projection

[TABLE]

for all $n\in[1,\dots,N]$.

While the theoretical convergence of FISTA requires the full convergence of inner Algorithm 4, in practice, it is sufficient to run about $5$-$10$ iterations with an initializer that corresponds to the dual variable from the previous outer iteration. In our implementation, we thus fix the maximal number of inner iterations to $t_{\text{\tiny in}}=10$ and enforce an additional stopping criterion based on measuring the relative change of the solution in two successive iterations as $\|{\mathbf{g}}^{t}-{\mathbf{g}}^{t-1}\|_{\ell_{2}}/\|{\mathbf{g}}^{t-1}\|_{\ell_{2}}\leq\delta_{\text{in}}$, where $\delta_{\text{in}}=10^{-4}$ in all the experiments.

Derivation of Error Backpropagation

In this appendix, we provide the derivation of error backpropagation applied to our method. The method essentially computes the gradient of the data-fidelity term. This gradient is a key step of updating the scattering potential in solving the inverse problem. We now present the mathematical derivation of the gradient computation and relate it to Algorithm 3.

The inputs of the error-back propagation are the data mismatch and the intermediate variables ($\{{\mathbf{s}}^{k}\}$, $\{\gamma_{k}\}$, $\{\mu_{k}\}$, ${\widehat{\mathbf{u}}}={\mathbf{u}}^{K}$) of the forward model computation, and the output is the gradient. Here we follow the differentiation conventions for vectors: $(\frac{\partial{\mathbf{u}}}{\partial{\mathbf{f}}})_{ij}=\frac{\partial{\mathbf{u}}_{i}}{\partial{\mathbf{f}}_{j}}$ and $(\nabla_{\mathbf{f}}{\mathbf{u}})_{ij}=[(\frac{\partial{\mathbf{u}}}{\partial{\mathbf{f}}})^{\text{H}}]_{ij}=\frac{\partial{\mathbf{u}}^{*}_{j}}{\partial{\mathbf{f}}_{i}}$. All boldface lower-case variables are column vectors.

Let us begin with the gradient of $\mathcal{D}=\frac{1}{2}\|{\mathbf{z}}-{\mathbf{y}}\|^{2}_{2}$.

[TABLE]

This can be evaluated by applying the chain rule to $\nabla_{\mathbf{f}}{\mathbf{z}}$ and all the variables composing ${\mathbf{z}}$. The equations leading from the initialization all the way to ${\mathbf{z}}$ are listed below:

[TABLE]

for $k=1,\ldots,K$, where ${\mathbf{A}}\,\triangleq\,{\mathbf{I}}-{\mathbf{G}}{\mathrm{diag}}\{{\mathbf{f}}\}$, and ${\mathbf{u}}^{-1}={\mathbf{u}}^{0}$. It is worth noting that, while the step-size $\gamma_{k}$ also depends on ${\mathbf{f}}$, we ignore this dependency to simplify the computation. The rationale for this simplification is that the step-size can be replaced by a fixed one. Furthermore, in practice, $\gamma_{k}$ attains a stationary value for large enough $k$, which indicates that this simplification has a negligible effect on backpropagation.

Initialization of backpropagation

The initialization in Algorithm 3 is obtained by differentiating the first of the above equations with respect to ${\mathbf{f}}$. With ${{\mathrm{diag}}\{{\mathbf{f}}\}{\mathbf{u}}={\mathrm{diag}}\{{\mathbf{u}}\}{\mathbf{f}}}$, we have

[TABLE]

The first term gives the remainder that contributes to the final result while the second term gives the vector that multiplies with $\nabla_{\mathbf{f}}{\mathbf{u}}^{K}$. For convenience, we define two sets of vectors:

• •

${\mathbf{q}}^{k}$: the vector that multiplies with $\nabla_{\mathbf{f}}{\mathbf{u}}^{k}$

• •

${\mathbf{r}}^{k}$: the remainder before computing $(\nabla_{\mathbf{f}}{\mathbf{u}}^{k}){\mathbf{q}}^{k}$.

In addition, due to the acceleration step in the forward computation, we expect subsequent ${\mathbf{q}}^{k-1}$ to have a contribution from $(\nabla_{\mathbf{f}}{\mathbf{u}}^{k+1}){\mathbf{q}}^{k+1}$ in addition to the contribution from its direct neighbor $(\nabla_{\mathbf{f}}{\mathbf{u}}^{k}){\mathbf{q}}^{k}$. This leads to the third set of vectors:

• •

${\mathbf{p}}^{k}$: the explicit contribution of $(\nabla_{\mathbf{f}}{\mathbf{u}}^{k+1}){\mathbf{q}}^{k+1}$ to ${\mathbf{q}}^{k-1}$.

Finally, multiplying (27) with $({\mathbf{z}}-{\mathbf{y}})$ we identify

[TABLE]

Since there is no term multiplying with $\nabla_{\mathbf{f}}{\mathbf{u}}^{K-1}$ explicitly (hence nothing to pass to ${\mathbf{q}}^{K-1}$), we have

[TABLE]

Recursive updates for ${\mathbf{u}}^{k}$

The computation of $(\nabla_{\mathbf{f}}{\mathbf{u}}^{k}){\mathbf{q}}^{k}$ is the key step in error-back propagation. We evaluate the gradient $\nabla_{\mathbf{f}}{\mathbf{u}}^{k}$ by taking the Hermitian of the derivative, and the multiplication with ${\mathbf{q}}^{k}$ follows. The result should be passed onto another gradient with smaller $k$. Before we start, let us write out the gradient of ${\mathbf{s}}^{k}$ which is straightforward from its definition,

[TABLE]

The derivative of ${\mathbf{u}}^{k}$ is

[TABLE]

The first term becomes

[TABLE]

and the second term becomes

[TABLE]

By taking Hermitian transpose of (32), we have

[TABLE]

where

[TABLE]

By multiplying (33) with ${\mathbf{q}}^{k}$ and substituting with (31), we obtain the expression for $(\nabla_{\mathbf{f}}{\mathbf{u}}^{k}){\mathbf{q}}^{k}$,

[TABLE]

Note that because we set ${\mathbf{u}}^{-1}={\mathbf{u}}^{0}$,

[TABLE]

Error backpropagation equations

From equations (27) to (30), we have

[TABLE]

Substituting (36) into (38), we have the following expressions

[TABLE]

and the recursion relations for $k=2\ldots K$

[TABLE]

For the case $k=1$, or namely ${\mathbf{r}}_{0}$ and ${\mathbf{q}}_{0}$ (note that ${\mathbf{p}}_{0}$ does not exist due to ${\mathbf{u}}^{-1}={\mathbf{u}}^{0}$), we plug (37) into (39),

[TABLE]

Hence we have

[TABLE]

In the initialization of our forward model, ${\mathbf{u}}^{0}$ is the incident field and does not depend on ${\mathbf{f}}$. Therefore $\nabla_{\mathbf{f}}{\mathbf{u}}^{0}=\mathbf{0}$ and the gradient of data-fidelity is

[TABLE]

We summarize these recursion relations of error backpropagation in Algorithm 3 of the main text.

Analytic Solutions to Scattering

In this section, our aim is to present the analytic expressions for scalar electric fields resulting from a point source outside a dielectric sphere in 2D and 3D (strictly speaking, the 2D case should be understood as an infinitely long line source illuminating a cylinder parallel to it and looking at the cross-section). A sketch of the derivation is provided after the actual expressions. A more complete description can be found in a number of standard textbooks such as [67].

Expressions

Consider a sphere of a radius $r_{\text{\tiny sph}}$ and a refractive index $n=\sqrt{\epsilon}$. The source is located $r_{\text{\tiny s}}$ distance away from the center of the sphere and the wavenumber of the source outside the sphere is $k_{b}$.

2D case:

We consider the polar coordinates

[TABLE]

and, without loss of generality, assume that the source is at $\theta_{\text{\tiny s}}=0$. The field can be expressed as

[TABLE]

where $\rho=k_{b}r$, $\rho_{\text{\tiny sph}}=k_{b}r_{\text{\tiny sph}}$ and $\rho_{\text{\tiny s}}=k_{b}r_{\text{\tiny s}}$,

[TABLE]

and $J_{m}$ and $Y_{m}$ are the $m$’th order Bessel functions of the first kind and the second kind, and $H^{(1)}_{m}=J_{m}+\mathrm{j}Y_{m}$ is the $m$’th order Hankel’s function of the first kind.

3D case:

We consider the spherical coordinates

[TABLE]

and, without loss of generality, assume that the source has zenith angle $\theta_{\text{\tiny s}}=0$ and azimuthal angle $\phi_{\text{\tiny s}}=0$. The field then reads

[TABLE]

where, with $\rho=k_{b}r$, $\rho_{\text{\tiny sph}}=k_{b}r_{\text{\tiny sph}}$ and $\rho_{\text{\tiny s}}=k_{b}r_{\text{\tiny s}}$,

[TABLE]

and $j_{l}$ and $n_{l}$ are the $l$’th order spherical Bessel function of the first kind and the second kind, $h^{(1)}_{l}=j_{l}+\mathrm{j}n_{l}$ is the $l$’th order of spherical Hankel function of the first kind, and $P_{l}(x)$ is the Legendre polynomial defined as

[TABLE]

Helmholtz equation

The Helmholtz equation for a point source is

[TABLE]

where $E$ is the complex electric field at position ${\mathbf{r}}$ when the point source is at ${\mathbf{r}}_{\text{\tiny s}}$ and $k^{2}({\mathbf{r}})$ is defined as

[TABLE]

Note that both 2D and 3D cases follow the same form. Their difference is that ${\mathbf{r}}$ and ${\mathbf{r}}_{\text{\tiny s}}$ have 2 or 3 coordinates.

Derivation for 2D case

We consider the polar coordinates, assume $\theta_{\text{\tiny s}}=0$ without loss of generality, and use an ansatz for the electric field in (46). With the Laplacian in the polar coordinate and the following expansion of a 2D delta function [67]

[TABLE]

eq. (71) becomes a sequence of equations on $R_{m}(r,r_{\text{\tiny s}})$

[TABLE]

for each $m$. Each equations is a Bessel differential equation so the solution can be composed of Bessel functions of order $m$. The boundary conditions for $R_{m}$ are as follows

finite value at $r=0$ 2. 2.

only outgoing component at $r=\infty$ 3. 3.

continuous and first-derivative-continuous at $r=r_{\text{\tiny sph}}$ 4. 4.

continuous at $r=r_{\text{\tiny s}}$ 5. 5.

$\left.\frac{\partial R_{m}}{\partial r}\right|_{r_{\text{\tiny s}}^{+}}-\left.\frac{\partial R_{m}}{\partial r}\right|_{r_{\text{\tiny s}}^{-}}=-\frac{1}{r_{\text{\tiny s}}}$ (integrate (LABEL:eq:analytical:2Dderive3) around $r_{\text{\tiny s}}$)

The above condition and equations lead to Eqs. (46)-(57).

Derivation for 3D case

We consider the spherical coordinates and assume that the source lies on the zenith axis. The ansatz for the electric field is (58), the expansion of a 3D delta function is

[TABLE]

and eq. (71) becomes

[TABLE]

for each $l$. These equations are spherical Bessel equations and there are corresponding spherical Bessel functions to compose the solution. The boundary conditions for the solution are the same as listed above except the last one becoming

[TABLE]