# Compressive Imaging with Iterative Forward Models

###### Abstract

We propose a new compressive imaging method for reconstructing 2D or 3D objects from their scattered wave-field measurements. Our method relies on a novel, nonlinear measurement model that can account for the multiple scattering phenomenon, which makes the method preferable in applications where linear measurement models are inaccurate. We construct the measurement model by expanding the scattered wave-field with an accelerated-gradient method, which is guaranteed to converge and is suitable for large-scale problems. We provide explicit formulas for computing the gradient of our measurement model with respect to the unknown image, which enables image formation with a sparsity-driven numerical optimization algorithm. We validate the method both analytically and with numerical simulations.

## 1 Introduction

Some of the most difficult, yet important, problems in computational sensing involve imaging objects that are hidden behind an opaque medium. For example, identifying a tumor inside a human body in medical diagnosis, detecting defects within a structure in industrial testing, or visualizing the shape of a multicellular organism in biology are all instances of this fundamental problem of subsurface imaging. The most commonly used approach in such applications is based on probing the medium with a controlled incident wave of a specific frequency range that can penetrate the medium, and then to rely on the physics of wave scattering to infer or visualize the spatial distribution of the refractive index within the medium.

This problem of inferring the refractive index distribution from the scattered wave-field is known as inverse scattering. It is often formulated as a large-scale optimization that relies on models for describing both the physical (forward or measurement model) and signal-related (regularization, prior constraints) aspects of the problem. Traditional approaches to inverse scattering are based on linear measurement models that can be obtained by assuming a straight-ray propagation of waves [1], or by adopting more refined scattering models based on the first Born [2] or Rytov approximations [3]. Recent works have demonstrated impressive imaging capability of the optimization-based approaches that also incorporate prior constraints on the solution [4, 5, 6]. In particular, dramatic improvements were obtained by relying on sparsity-promoting regularization [7, 8, 9], which is an essential component of compressive sensing [10, 11]. The basic motivation is that many natural objects are inherently sparse in some transform domain and can be reconstructed with high accuracy even with a small amount of measured data.

Linear measurement models, though simple and efficient, are only accurate for weakly scattering objects. This limits their applicability for imaging larger objects and/or those with large refractive index contrasts. Recent experimental results also indicate that the resolution and quality of the reconstructed image is improved when nonlinear measurement models are used [12, 13, 14, 15, 16, 17, 18, 19]. In particular, nonlinear models can account for multiple scattering and thus provide a more accurate interpretation of the measured data.

We propose a new method for reconstructing the refractive index from measurements of the scattered wave-field. This method combines our nonlinear forward model with an edge-preserving total variation (TV) regularizer [20] and forms images by solving a large-scale optimization problem. Our measurement model—called series expansion with accelerated gradient descent on Lippmann-Schwinger equation (SEAGLE)—is based on formulating wave-scattering as a smooth optimization subproblem and using Nesterov’s fast gradient method [21] to iteratively approximate the scattered wave. The key advantage of SEAGLE is its guaranteed convergence, even for objects with large refractive index contrasts. We provide explicit formulas for computing the gradient of our measurement model with respect to the refractive index, which enables large-scale 2D and 3D imaging using fast iterative shrinkage/thresholding algorithm (FISTA) [22]. We validate our forward model and reconstruction method analytically and with numerical simulations.

## 2 Main Contribution

Our approach expands the scattered wave field with the iterates of Nesterov’s accelerated gradient descent and efficiently computes the derivative of the field with respect to the object. The inverse problem is formulated as a TV-regularized data fitting with complex scattered-wave measurements.

### 2.1 Problem formulation

Consider a setup illustrated in Fig. 1 where the unknown object resides inside the image domain , where . We want to recover the refractive index of the unknown object given wave measurements at point sensors (on the vertical solid lines) and controllable sources (solid circles). Monochromatic light scattering by nonuniform refractive index can be described by the Lippmann-Schwinger equation

(1) |

where is the complex total electric field, is the complex incident electric field, is the free-space Green’s function, is the scattering potential, which is assumed to be real and contains the map of refractive index of the object , and is the wavenumber of the background medium. This integral is only over domain as is zero outside of . The free-space Green’s function in (1) is given by

(2) |

where and is the Hankel function of first kind. The Green’s function is obtained under the outgoing wave boundary condition, Helmholtz equation , and the time-dependence convention where the physical electric field equals to .

Eq. (1) provides a nonlinear relationship between the wave-field and the scattering potential , whereas first Born and Rytov approximations provide simplified linearized versions of this relationship. The inverse problem is to find an estimation of given the measurements of at point sensors.

### 2.2 Algorithmic Expansion of the Wave Model

For points in the domain , the discretized version of (1) can be expressed as

(3) |

where the operator forms a diagonal matrix from its argument, , and are the discretized versions of , and , respectively, and is the matrix representing the convolution with the Green’s function within . Here, denotes the number of sample points. Note that the field outside can be evaluated by using a different matrix , which corresponds to evaluating (1) at sensor points.

As a forward model, we propose to solve (3) for by applying Nesterov’s fast gradient method to the following minimization problem

(4) |

where . The full procedure is summarized in Algorithm 1, which we call SEAGLE. Note the dependence of the solution on the scattering potential and the fact that it can be interpreted as an expansion of the wave-field with the iterates of the accelerated gradient descent method.

### 2.3 Inverse problem

We formulate the inverse problem as the minimization

(5) |

where

(6a) | |||

(6b) |

Here, is the quadratic data-fidelity term that measures the discrepancy between the measured data and the output of SEAGLE forward computation . The functional is -dimensional isotropic TV regularizer, where is the discrete gradient operators along the axis . The regularization parameter controls the strength of regularization, while the set is used for enforcing additional physical constraints on such as, for example, non-negativity.

The two key steps of FISTA for solving the optimization problem (5) are computing the and evaluating the proximal operator

(7) |

for some and [22]. The proximal operator corresponding to TV can be efficiently computed [23, 24].

The iterative structure of the SEAGLE forward computation allows for an efficient computation of the gradient, which can be expressed as

(8) |

with . By differentiating lines 11, 10 and 7 in Algorithm 1 with respect to and keeping constant, we obtain an efficient error-back propagation rule for computing (8) summarized in Algorithm 2. Note that the matrices in Algorithm 2 are implemented as operators on without the need for an explicit storage in memory. The remarkable aspect of Algorithm 2 is that it explicitly provides the gradient that can be used for FISTA-based minimization of (5). While the optimization problem (5) is generally non-convex, we did not observe any practical convergence issues in our simulations (see Section (3)).

### 2.4 Network interpretation

Figure 2 graphically illustrates Algorithms 1 and 2 as feedforward networks. Each square module represents the operation in each iteration of SEAGLE, and the edges represent the vectors as message carriers. A nice feature of SEAGLE is that it can easily incorporate additional modules for modeling other physical phenomena. For example, we can prepend a module representing an initialization with Rytov approximation in front of the module Iter 1. Then, in the back propagation, and are fed into the Rytov module in which the Rytov inverse step is applied. The flexibility of SEAGLE will be explored in future works for speeding up the computation or for dealing with other measurement scenarios where only intensity of the wave-field is preserved.

## 3 Numerical Results

### 3.1 Analytic validation of the forward model

In order to validate the forward model, we first consider two simple scattering experiments where it is possible to derive analytic forms of the scattered wave-fields: a 2D point source scattered by a cylinder, and a 3D point source scattered by a sphere (see Sections 3.8-3.11 in [25]). In both cases the scatterers have diameters of wavelengths (Fig. 3(a) and 3(c)). The wavelength is mm, the grid size is mm ( mm) and there are points ( points) along each axis in 2D (3D). We define the contrast of an object as . In order to evaluate the performance quantitatively, we plot in Fig. 3(b) and 3(c) the normalized error, , where is the analytic solution. We additionally provide errors achieved when using the first Born and Rytov approximations at contrast. In Fig. 3(e) and 3(f), we provide visual comparison between the analytic solution and the result of our model.

### 3.2 Inverse scattering experiment

We next use the proposed technique for reconstructing the Shepp-Logan phantom in the ill-posed and strongly scattering regime. Specifically, we consider the setup in Fig. 1 where the scattered wave measurements are generated by a high-fidelity finite-difference time-domain (FDTD) [26] simulator. The object is of size cm cm and has a contrast of . We put two linear detectors on both sides of the phantom at a distance of cm from the center of the object, and each detector has sensors placed with in-between spacing of cm. The transmitters are put on a line cm left to the left detector. They are spaced uniformly in azimuth with respect to the center of the phantom (every within ). We setup a cm cm square area for reconstructing the object, with pixel size cm. The wavelength of the illuminating light is cm.

Fig. 4 summarizes the performances of the proposed method, as well as two baseline methods based on the first Born and Rytov approximations. All three methods are implemented iteratively with isotropic TV regularizer (). In SEAGLE, we set but the forward algorithm may stop earlier when the objective function (4) is below . It is shown that SEAGLE outperforms first Born and Rytov methods. It can be seen that due to the ill-posed nature of the measurements, the reconstructed images suffer from the missing frequency artifacts [27]. However, our method is still able to accurately capture most features of the object while the linear methods cannot.

## 4 Conclusion

Our method is suitable for compressive imaging in the presence of multiple scattering. It can handle both transmission and reflection data, unlike alternative methods based on beam propagation. It is additionally more stable compared to the methods based on iterative Born approximations.

## References

- [1] A. C. Kak and M. Slaney, Principles of Computerized Tomographic Imaging, IEEE, 1988.
- [2] M. Born and E. Wolf, Principles of Optics, Cambridge Univ. Press, 7th (expanded) edition edition, 1999.
- [3] A. J. Devaney, “Inverse-scattering theory within the Rytov approximation,” Opt. Lett., vol. 6, no. 8, pp. 374–376, August 1981.
- [4] W. Choi, C. Fang-Yen, K. Badizadegan, S. Oh, N. Lue, R. R. Dasari, and M. S. Feld, “Tomographic phase microscopy,” Nat. Methods, vol. 4, no. 9, pp. 717–719, September 2007.
- [5] Y. Sung, W. Choi, C. Fang-Yen, K. Badizadegan, R. R. Dasari, and M. S. Feld, “Optical diffraction tomography for high resolution live cell imaging,” Opt. Express, vol. 17, no. 1, pp. 266–277, December 2009.
- [6] T. Kim, R. Zhou, M. Mir, S. D. Babacan, P. S. Carney, L. L. Goddard, and G. Popescu, “White-light diffraction tomography of unlabelled live cells,” Nat. Photonics, vol. 8, pp. 256–263, March 2014.
- [7] M. M. Bronstein, A. M. Bronstein, M. Zibulevsky, and H. Azhari, “Reconstruction in diffraction ultrasound tomography using nonuniform FFT,” IEEE Trans. Med. Imag., vol. 21, no. 11, pp. 1395–1401, November 2002.
- [8] Y. Sung and R. R. Dasari, “Deterministic regularization of three-dimensional optical diffraction tomography,” J. Opt. Soc. Am. A, vol. 28, no. 8, pp. 1554–1561, August 2011.
- [9] J. W. Lim, K. R. Lee, K. H. Jin, S. Shin, S. E. Lee, Y. K. Park, and J. C. Ye, “Comparative study of iterative reconstruction algorithms for missing cone problems in optical diffraction tomography,” Opt. Express, vol. 23, no. 13, pp. 16933–16948, June 2015.
- [10] E. J. Candès, J. Romberg, and T. Tao, “Robust uncertainty principles: Exact signal reconstruction from highly incomplete frequency information,” IEEE Trans. Inf. Theory, vol. 52, no. 2, pp. 489–509, February 2006.
- [11] D. L. Donoho, “Compressed sensing,” IEEE Trans. Inf. Theory, vol. 52, no. 4, pp. 1289–1306, April 2006.
- [12] F. Simonetti, “Multiple scattering: The key to unravel the subwavelength world from the far-field pattern of scattered wave,” Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., vol. 73, no. 3, pp. 036619, March 2006.
- [13] F. Simonetti, M. Fleming, and E. A. Marengo, “Illustration of the role of multiple scattering in subwavelength imaging from far-field measurements,” J. Opt. Soc. Am. A, vol. 25, no. 2, pp. 292–303, February 2008.
- [14] G. Maire, F. Drsek, J. Girard, H. Giovaninni, A. Talneau, D. Konan, K. Belkebir, P. C. Chaumet, and A. Sentenac, “Experimental demonstration of quantitative imaging beyond Abbe’s limit with optical diffraction tomography,” Phys. Rev. Lett., vol. 102, pp. 213905, May 2009.
- [15] L. Tian and L. Waller, “3D intensity and phase imaging from light field measurements in an LED array microscope,” Optica., vol. 2, pp. 104–111, February 2015.
- [16] U. S. Kamilov, I. N. Papadopoulos, M. H. Shoreh, A. Goy, C. Vonesch, M. Unser, and D. Psaltis, “Learning approach to optical tomography,” Optica, vol. 2, no. 6, pp. 517–522, June 2015.
- [17] U. S. Kamilov, I. N. Papadopoulos, M. H. Shoreh, A. Goy, C. Vonesch, M. Unser, and D. Psaltis, “Optical tomographic image reconstruction based on beam propagation and sparse regularization,” IEEE Trans. Comp. Imag., vol. 2, no. 1, pp. 59–70,, March 2016.
- [18] U. S. Kamilov, D. Liu, H. Mansour, and P. T. Boufounos, “A recursive Born approach to nonlinear inverse scattering,” IEEE Signal Process. Lett., 2016, arXiv:1603.03768 [cs.LG].
- [19] T. Zhang, C. Godavarthi, P. C. Chaumet, G. Maire, H. Giovannini, A. Talneau, M. Allain, K. Belkebir, and A. Sentenac, “Far-field diffraction microscopy at resolution,” Optica, vol. 3, no. 6, pp. 609–612, June 2016.
- [20] L. I. Rudin, S. Osher, and E. Fatemi, “Nonlinear total variation based noise removal algorithms,” Physica D, vol. 60, no. 1–4, pp. 259–268, November 1992.
- [21] Y. E. Nesterov, “A method for solving the convex programming problem with convergence rate ,” Dokl. Akad. Nauk SSSR, vol. 269, pp. 543–547, 1983, (in Russian).
- [22] A. Beck and M. Teboulle, “A fast iterative shrinkage-thresholding algorithm for linear inverse problems,” SIAM J. Imaging Sciences, vol. 2, no. 1, pp. 183–202, 2009.
- [23] A. Beck and M. Teboulle, “Fast gradient-based algorithm for constrained total variation image denoising and deblurring problems,” IEEE Trans. Image Process., vol. 18, no. 11, pp. 2419–2434, November 2009.
- [24] U. S. Kamilov, “Parallel proximal methods for total variation minimization,” in IEEE Int. Conf. Acoustics, Speech and Signal Process. (ICASSP 2016), Shanghai, China, March 19-25, 2016, pp. 4697–4701.
- [25] J. D. Jackson, Classical electrodynamics, Wiley, 3 edition, 1999.
- [26] A. Taflove and S. C. Hagness, Computational Electrodynamics: The Finite-Difference Time-Domain Method, Artech House Publishers, 3 edition, 2005.
- [27] C. J. R. Sheppard and C.J. Cogswell, “Three-dimensional image formation in confocal microscopy,” Journal of microscopy, vol. 159, no. 2, pp. 179–194, 1990.