1. INTRODUCTION

In inverse scattering, one emits wave fields that are incident on a medium containing scattering objects and seeks to reconstruct the location, shape, material properties, etc., of those scattering objects from measurements of scattered fields. Imaging applications for inverse scattering usually have essential limitations on these measurements. These limitations often lead to the inverse scattering problem being mathematically underdetermined and ill posed. In optics, inverse scattering problems are inherently complicated because measurements are often limited to intensities of the scattered field. As a result, those measurements do not possess the phase information needed for commonly used imaging methods. Phase retrieval methods attempt to recover this missing phase information from intensity measurements, e.g., see [1] and references contained therein. Alternatively, one may take interferometric measurements using a known reference field to obtain relative phase information that can be used to solve the inverse problem. For example, Choi et al. [2] use a Mach–Zehnder heterodyne interferometer to produce quantitative phase images to image cells and multicellular organisms. Their imaging method is based on a filtered back-projection algorithm and requires a phase unwrapping method to remove a $ 2\pi $-phase ambiguity.

There are several noteworthy works for intensity-only inverse scattering. Wolf [3] gives a method for extracting the amplitude and phase of the scattered fields from holograms. Devaney [4] and Maleki et al. [5] use the Rytov approximation for scattering to derive a formula for propagation of the complex phase perturbation incurred by the weakly scattering obstacle. A filtered back-projection method is used on that formula to recover the complex phase perturbation from far-field measurements of the scattered field. Maleki et al. [6] compare this method with a phase-retrieval method and show that the phase-retrieval method performs better, but requires a priori information about the scattering object. Carney et al. [7] make use of the optical theorem to develop an intensity-only imaging method using extinguished power measurements. Gbur and Wolf [8] describe an intensity-only diffraction tomography method that requires two or more measurement planes separated by a known distance. Marengo et al. [9] modify the multiple signal classification (MUSIC) method for intensity-only measurements by linearizing the problem with respect to a proxy source function that maintains the location and shape of the scattering object. The challenge in using this method is that one needs either to exhaustively search in a six-dimensional space, which is computationally challenging, or consider a reduced, three-dimensional space at the expense of resolution and sensitivity to noise. D’Urso et al. [10] describe an imaging method using intensity measurements of the total field based on the contrast source-extended Born approximation. That approximation suggests a cost function to minimize for which they give an effective optimization method. By taking advantage of the sparsity of the image, Brady et al. [11] and Zhang et al. [12] obtain accurate images in in-line digital holography. Tahir et al. [13] has extended those compressive holography results to multiple scattering regimes. Chai et al. [14] and Pham et al. [15] formulate the intensity-only imaging problem as a non-convex optimization problem that they then solve effectively using sparsity-promoting methods. Tian and Waller [16] and Horstmeyer et al. [17] use principles from ptychography and diffraction tomography to reconstruct images of the complex index of refraction based on the Born approximation. In contrast to those works, Ling et al. [18] have developed a slice-based framework for three-dimensional microscopy that reduces the overall number of measurements required for reconstruction.

We present a method to solve the inverse scattering problem that consists of determining the location and shape of scattering objects using intensity-only measurements taken over a single measurement plane. This study is limited to scalar waves and does not take into account the polarization state of light. This imaging method is direct and therefore does not require phase retrieval or any other iterative procedure. It is not restricted in the validity region of the Born or Rytov approximations. As a result, this method can be used on high-contrast scattering objects in thick samples. In the examples considered here, we have accurately reconstructed the location and shape of scattering objects whose sizes and separation distances from one another are on the order of the wavelength and have a relative refractive index of 1.4. The resulting image is reconstructed over the entire imaging region at once without the need for scanning or slicing of the imaging region. Our results indicate that the resolution is on the order of the wavelength and is the same along axes perpendicular and parallel to the measurement plane.

The key to this method lies in creating data diversity through a suite of experiments in which the imaging region is illuminated using different incident fields. An explicit algebraic relation called the polarization identity is used to process intensity measurements of the scattered field. The polarization identity has been used to study several intensity-only imaging problems [19–21]. Here, we use it to convert intensity measurements of the scattered field to interferometric measurements in which one of the scattered fields acts as the reference. Despite the fact that this reference scattered field is unknown, we show that the MUSIC method data can be formulated for these interferometric data without needing the modifications used by Marengo et al. [9].

MUSIC is a direct sampling method that recovers the location and shape of the scattering objects. Moscoso et al. [22] give a rigorous mathematical analysis of the MUSIC algorithm. Ammari et al. [23] have extended MUSIC to study subsurface imaging problems. Gruber et al. [24] demonstrated that MUSIC can be used beyond the Born linearization regime when multiple scattering is considered. Chai et al. [25] give an alternative method for imaging strong point scatterers relying on sparsity promoting optimization. Ammari et al. [26] have extended the theory of MUSIC to the full electromagnetic inverse scattering problem. This work highlights inherent differences in the full electromagnetic scattering case from the scalar wave approximation. Iakovleva et al. [27] apply this full electromagnetic theory to subsurface imaging of an extended object. Chen [28] extends the intensity-only theory by Marengo et al. [9] to the full electromagnetic scattering problem for point-like scatterers. Chen [29] extends this theory to reconstruct extended objects in two dimensions by first applying the MUSIC method and then using that result to solve a nonlinear least-squares problem.

In the examples we show here, we illuminate the imaging region, i.e., a cube of size $ 20\lambda \times 20\lambda \times 20\lambda $, using plane waves at different angles of incidence similarly to how Carney et al. [7] have done. Varying the angle of incidence has been demonstrated to be practically feasible in many imaging systems. Choi et al. [2] have used a galvanometer-mounted tilting mirror to vary the angle of illumination. Zheng et al. [30] introduce the use of an array of LEDs to create source diversity that effectively illuminates the imaging region with plane waves with different angles of incidence. However, the formulation described in Section 2 shows that the approach is general and can incorporate any source diversity available.

The remainder of this paper is as follows. In Section 2, we introduce the polarization identity and use it in an explicit procedure that converts intensity measurements to interferometric measurements in which one of the scattered fields acts as the reference. We describe a modification of the MUSIC method that uses those interferometric data in Section 3. In Section 4, we describe the specific problem used here to test this imaging method. Additionally, we describe the method of fundamental solutions used to simulate measurements. We give several image reconstructions in Section 5 that demonstrate the effectiveness of this imaging method. We discuss several practical considerations of this imaging method in Section 6. In Section 7, we give our conclusions.

2. POLARIZATION IDENTITY

Consider two complex numbers: $ {z_1} $ and $ {z_2} $. It follows that

Let $ {u^s} $ denote the field scattered by a medium due to $ {u^i} $, the field incident on that medium. These fields are related to one another according to $ {u^s}(\textbf{r}) = L{u^i}(\textbf{r}) $, where $ L $ is the scattering operator giving the linear mapping from $ {u^i} $ to $ {u^s} $ evaluated at position $ \textbf{r} $. Suppose that measurements correspond to $ |{u^s}(\textbf{r}{)|^2} $. These measurements are limited since there is no phase information available in them. However, suppose that we are able to introduce source diversity into the problem. For example, suppose we are able to illuminate the medium with linear combinations of $ M $ different, known incident fields, denoted by $ u_m^i $ for $ m = 1, \cdots ,M $. Then we perform the following suite of experiments.

With these $ 3M - 2 $ measurements, we apply Eq. (3) and obtain

For the imaging method we discuss below, we use the data matrix $ B $ as measurements for reconstructing images of a multiple scattering medium. Using the procedure described above, we explicitly obtain $ B $ using $ 3M - 2 $ illuminations. At this point, we have not specified exactly what $ M $ incident fields are to be used. In fact, one can consider introducing any form of spatial, angular, or wavelength diversity at the source plane and apply the procedure given above to obtain a useful matrix $ B $ of interferometric measurements. In the results below, we consider angular diversity introduced by illuminating the imaging region with plane waves propagating with different angles of incidence.

This measurement procedure leverages the diversity created by using multiple illuminations at the source plane to overcome the inherent limitations of intensity-only measurements. The data contained in $ B $ correspond to unusual interferometric measurements because the reference field, $ u_1^s({\textbf{r}_n}) $ for $ n = 1, \cdots ,N $, is not known. Nonetheless, we will show that $ B $ contains sufficient information for reconstructing images. Because $ B $ is formed through the explicit algebraic relation given in Eq. (3), this procedure represents a substantial simplification over other methods that require phase retrieval.

3. MULTIPLE SIGNAL CLASSIFICATION METHOD

To reconstruct an image over an imaging region, we introduce a mesh that covers the imaging region with the set of grid points $ {\boldsymbol{\rho }_k} $ for $ k = 1, \cdots ,K $. We then seek to recover a grid function, $ {X_k} $ for $ k = 1, \cdots ,K $, whose entries are nonzero when the corresponding grid point lies inside a scattering object and zero otherwise. In this framework, a scattering object is represented as the collection of grid points that lie inside that object. When there are relatively few objects in the imaging region, we expect that most of the grid function values $ {X_k} $ are zero so that $ X $ is sparse.

A mapping from $ {X_k} $ for $ k = 1, \cdots ,K $ to the data matrix $ B $ defined in Eq. (6) is required for reconstructing an image. Here, we consider a very simple model for this mapping in which each of the grid points independently acts as a secondary point source. The grid function $ {X_k} $ for $ k = 1, \cdots ,K $ is then the collection of point source strengths associated with each of the grid points, and the measurements correspond to the superposition of these $ K $ point sources. We explain this model in detail below.

Suppose we illuminate the imaging region with incident field $ u_m^i $. The model we use for the scattered field at $ {\textbf{r}_n} $ is given by

When we combine all $ N $ measurements, we obtain the matrix system

The model given in Eq. (7) leading to Eq. (13) is a linearized approximation of the direct scattering problem. Therefore, its inversion provides only an approximate solution to the inverse scattering problem. However, Beylkin [32,33] has shown that the solution of the linearized inverse scattering problem preserves discontinuities. It follows that the linearized problem is sufficient for recovering the location and shape of the scattering objects, which is why we use it here.

In Eq. (13), each column of $ B $ is given as a linear combination of the columns of $ A $. According to Eq. (11), the columns of $ A $ are evaluations of the known incident fields on the grid points. It follows from the special structure of the linear system given in Eq. (13) that our measurements correspond to different linear combinations of the same columns of $ A $ since $ x $ is the same for each of the measurements. Furthermore, because the vector $ x $ of grid function values is sparse, there are only a small number of columns of $ A $ that contribute to the data. If we are able to determine the span of those columns of $ A $, we effectively determine those nonzero entries of $ x $ that, in turn, give the location and shape of the scattering objects.

Our implementation of the MUSIC method uses the fact that measurements are linear combinations of the same columns of $ A $. The span of those columns of $ A $ is called the signal subspace. When we compute the singular value decomposition, $ B = U\Sigma {V^H} $, where the superscript $ H $ denotes a conjugate transpose, the columns of $ U $ corresponding to the significant singular values contained in the diagonal entries of $ \Sigma $ give an orthonormal basis for the signal subspace. Let $ \tilde U $ denote the matrix containing the columns of $ U $ corresponding to the significant singular values. We introduce the matrix

To summarize the imaging method, we give the following procedure for forming images.

4. SIMULATING MEASUREMENTS

To test and evaluate the imaging method described in Section 3, we consider a general, three-dimensional imaging problem. We neglect polarization and consider scalar fields. This problem provides a simple setting to evaluate the effectiveness of the imaging method. This general imaging problem can be easily modified or extended to consider specific imaging systems.

A. Imaging System

Let $ z = - 200\lambda $ denote the source plane and $ z = + 200\lambda $ denote the measurement plane. The imaging region is a $ 20\lambda \times 20\lambda \times 20\lambda $ cube centered at the origin that contains one or more scattering objects. The background refractive index is $ {n_0} $, and the refractive index inside the scattering objects is $ {n_1} $.

We illuminate the imaging region by plane waves of the form

Measurements of the scattered field are collected on a $ 400\lambda \times 400\lambda $ region of the measurement plane, $ z = 200\lambda $. The scattered field is sampled on this measurement region with an equi-spaced mesh with mesh width $ \Delta x = \Delta y = 10\lambda $. Consequently, we have $ N = 1681 $ scattered field measurements for each illumination.

B. Forward Problem

Suppose that $ Q $ scattering objects in the imaging region correspond to the disjoint regions $ {\Omega _q} $ for $ q = 1, \cdots ,Q $ with corresponding boundaries $ \partial {\Omega _q} $ for $ q = 1, \cdots ,Q $. Let $ {\bar \Omega _q} = {\Omega _q} \cup \partial {\Omega _q} $ and let $ E = {{\mathbb R}^3}\backslash \bigcup\nolimits_{q = 1}^Q\! {\bar \Omega _q} $ denote the exterior to all of the scattering objects. The refractive index in $ E $ is $ {n_0} $, and the refractive index in $ {\Omega _q} $ is $ {n_1} $ for $ q = 1, \cdots ,Q $. We define the corresponding wavenumbers: $ {k_0} = 2\pi {n_0}/\lambda $ and $ {k_1} = 2\pi {n_1}/\lambda $. We seek the solution of the following system of boundary-value problems for the scalar, time-harmonic wave equation:

To compute the solution of boundary-value problem Eq. (18), we use the Method of Fundamental Solutions. This method was introduced by Mathon and Johnston [34]. It provides an accurate and efficient computational method for solving the full scattering problem. The text by Wriedt et al. [35] provides an overview of this method and its applications to various problems.

The fundamental solution

Let $ {\textbf{r}_{qp}} $ for $ p = 1, \cdots ,P $ denote a set of points on $ \partial {\Omega _q} $. We introduce the points

Because Eq. (23a) is a sum of fundamental solutions, each satisfying Eq. (20), and because the points $ \textbf{r}_{qp}^s \in {\Omega _q} $ for $ p = 1, \cdots ,P $, this approximation exactly satisfies Eq. (18a) for all $ \textbf{r} \in E $. Similarly, Eq. (23b) exactly satisfies Eq. (18b) for all $ \textbf{r} \in {\Omega _q} $. All that remains is to require that Eqs. (23a) and (23b) satisfy boundary conditions Eqs. (18c) and (18d), and appropriate radiation conditions. Because each fundamental solution, $ {G_0} $, satisfies the outgoing radiation condition, Eq. (23a) automatically satisfies the outgoing radiation condition. Thus, we need to require only that Eqs. (23a) and (23b) satisfy boundary conditions (18c) and (18d) for $ q = 1, \cdots ,Q $. We cannot require that Eqs. (23a) and (23b) satisfy the boundary conditions exactly. Instead, we require that Eqs. (23a) and (23b) satisfy boundary conditions (18c) and (18d) on the $ P $ boundary points, $ {\textbf{r}_{qp}} $ for $ p = 1, \cdots ,P $. Doing so yields the $ 2QP $ equations

The Method of Fundamental Solutions described above approximates the full boundary-value problem for the wave equation with multiple scattering obstacles. It does not make any approximations with regards to the scattering—it includes all orders of multiple scattering needed for the problem. There are no assumptions on the refractive indices used or the shape of the scattering objects besides some reasonable smoothness requirements for the boundaries. This method approximates only the requirement that the fields satisfy boundary conditions, which depends on the number $ P $ of points used to sample each of the boundaries and the distance $ \ell $ used. For the simulations used in the results below, we have chosen $ P = 512 $ and $ \ell \approx 0.1{\sigma _g^{1/2}} $, with $ {\sigma _g} $ denoting the geometrical cross section of the scattering object. With these choices of user-defined parameters, we obtain less than 0.1% relative error made by this method to compute the scattering cross section by a single, spherical scattering object with radius $ a = 632\,\,\textrm{nm} $ over the visible spectrum compared to the analytical solution for this problem.

For a particular configuration of scattering objects, we compute the scattered field $ u_m^s $ for each of the $ M = 625 $ incident plane waves given in Eq. (16) and evaluate those results at the points $ {\textbf{r}_n} $ for $ n = 1, \cdots ,N $ on the measurement plane. We then use those results to form the matrix $ B $ defined in Eq. (6). This $ B $ matrix is then used in the MUSIC imaging method described in Section 3 to form reconstructions of those scattering objects.

5. RESULTS

We present results of images constructed using the MUSIC method described in Section 3 with simulated measurements computed using the Method of Fundamental Solutions described in Section 4. For all the results that follow, the wavelength is $ \lambda = 632\,\,\textrm{nm} $, the refractive index of the background is $ {n_0} = 1 $, and the refractive index of the scattering objects is $ {n_1} = 1.4 $. We illuminated the $ 20\lambda \times 20\lambda \times 20\lambda $ imaging region containing the scattering objects using $ M = 625 $ plane waves with different angles of incidence on the imaging region. We collect $ N = 1681 $ measurements over the $ 200\lambda \times 200\lambda $ region of the measurement plane $ z = 200\lambda $, sampling at $ \Delta x = \Delta y = 10\lambda $.

In what follows, we show an isosurface of $ {I_k} $ for $ k = 1, \cdots ,K $ given in Eq. (15) plotted over the surface of the actual scattering object. To compute $ {I_k} $ using the procedure given at the end of Section 3, we set the user-defined threshold in step 1 to be $ \delta { = 10^{ - 8}} $. Note that $ 0 \le {I_k} \le 1 $. We find that the reconstructed images are quite accurate. In the examples that follow, we plot the isosurface corresponding to $ {I_k} = 0.1 $.

The MATLAB codes used to generate all of the examples that follow are available in a GitHub repository [36].

A. Reconstructing Spherical Objects

We first consider spherical scattering objects each with radius $ a = 632\, \,\textrm{nm} $ and relative refractive index of $ {n_1}/{n_0} = 1.4 $. In Fig. 1, we show a single, spherical scattering object in blue and an isosurface corresponding to $ {I_k} = 0.1 $ of the reconstructed image in red. This result shows that the reconstruction accurately recovers the location and shape of the scattering object.

 

Fig. 1. Spherical scattering object plotted in blue with radius $ a = 632\, \,\textrm{nm} $ and relative refractive index $ {n_1}/{n_0} = 1.4 $ and an isosurface corresponding to $ {I_k} = 0.1 $ of the reconstructed image plotted in red.

Download Full Size | PDF

To test the ability of the imaging method to recover multiple scattering objects, we plot in Fig. 2 an image reconstruction for three spheres, each with radius $ a = 632\, \,\textrm{nm} $ and relative refractive index $ {n_1}/{n_0} = 1.4 $. The distance between each of the spheres is approximately $ 1.5\lambda $. Thus, this problem corresponds to high-contrast scattering objects on the order of the wavelength separated by a distance that is on the order of the wavelength. It is a particularly challenging case since the multiple scattering between these scattering objects is significant and complex. Nonetheless, we find that the reconstructed image accurately recovers the locations and shapes of the three, distinct scattering objects. We have found that when these same spheres are closer than this distance, the reconstructed image will deteriorate and does not distinguish the three distinct scattering objects from one another. However, when the relative refractive index is decreased, we find that the reconstructed images do distinguish three distinct scattering objects, and accurately recover their locations and shapes.

 

Fig. 2. Three spherical scattering objects plotted in blue, each with radius $ a = 632\,\, \textrm{nm} $ and relative refractive index $ {n_1}/{n_0} = 1.4 $ and an isosurface corresponding to $ {I_k} = 0.1 $ of the reconstructed image plotted in red.

Download Full Size | PDF

To consider even more scattering objects, we show in Fig. 3 a reconstructed image of eight spherical scatterers, each with radius $ a = 632\, \,\textrm{nm} $ and relative refractive index $ {n_1}/{n_0} = 1.4 $. Just as with the single sphere and the three spheres results, we find that the imaging method accurately recovers the locations and shapes of these scattering objects.

 

Fig. 3. Eight spherical scattering objects plotted in blue, each with radius $ a = 632\, \,\textrm{nm} $ and relative refractive index $ {n_1}/{n_0} = 1.4 $ and an isosurface corresponding to $ {I_k} = 0.1 $ of the reconstructed image plotted in red.

Download Full Size | PDF

B. Reconstructing Ellipsoidal Objects

To test the imaging method on scattering objects that have a more complex shape than a sphere, we consider an ellipsoid defined according to

 

Fig. 4. An ellipsoidal scattering object plotted in blue with principal semi-axes: $ a = 632\,\, \textrm{nm} $, $ b = 948\,\, \textrm{nm} $, and $ c = 474\,\, \textrm{nm} $, and relative refractive index $ {n_1}/{n_0} = 1.4 $. The isosurface corresponding to $ {I_k} = 0.1 $ of the reconstructed image plotted in red.

Download Full Size | PDF

In Fig. 5, we plot the reconstructed image of eight ellipsoidal scattering objects. For this result, we have randomly chosen the orientation for each ellipsoid. Here, we take $ a = 500\, \,\textrm{nm} $, $ b = 400\, \,\textrm{nm} $, and $ c = 300\, \,\textrm{nm} $, and relative refractive index $ {n_1}/{n_0} = 1.4 $. The imaging method accurately determines the location and shape of each of these eight scattering objects. There are no imaging artifacts that appear in this reconstruction.

 

Fig. 5. Eight ellipsoidal scattering objects plotted in blue with principal semi-axes: $ a = 500\, \,\textrm{nm} $, $ b = 400\,\, \textrm{nm} $, and $ c = 300\, \,\textrm{nm} $, and relative refractive index $ {n_1}/{n_0} = 1.4 $. The orientations of these ellipsoidal scattering objects have been chosen randomly. The isosurface corresponding to $ {I_k} = 0.1 $ of the reconstructed image plotted in red.

Download Full Size | PDF

6. PRACTICAL CONSIDERATIONS

We discuss several practical considerations of this imaging method.

7. CONCLUSION

We have presented a simple and effective inverse scattering method limited to intensity measurements taken on a single measurement plane. By leveraging diversity created by adequate multiple illuminations of the imaging region and using the explicit algebraic relation given by Eq. (3), we convert intensity measurements to the interferometric measurements contained in the matrix $ B $ defined in Eq. (6). The explicit procedure is given at the end of Section 2. We then apply to this $ B $ matrix a modification of the MUSIC method as described in Section 3. This leads to a simple imaging method that requires only some elementary linear algebra computations. It is efficient because the image reconstruction is the result of a direct computation—there is no iterative procedure to be done. There is no need for phase retrieval. Moreover, this method reconstructs images over the entire imaging region simultaneously. There is no scanning or slicing of the region required. This imaging method has no inherently different resolution on axes perpendicular and parallel to the measurement plane. Provided that the model for the measurements can be modified for a specific imaging system, this method is widely applicable to a broad variety of settings.

By simulating three-dimensional scattering data using the Method of Fundamental Solutions described in Section 4, we have tested this imaging method for a variety of scattering objects. All of the sizes of the objects we have considered are on the order of the wavelength of the illuminating fields. Additionally, when we considered multiple objects, the distances between them were on the order of a wavelength. Finally, their refractive indices were substantially different from the background. This situation is challenging because the multiple scattering by these objects is complex and not subject to simplifying approximations.

In our reconstructions using these simulated data, we have found that the imaging method presented here is remarkably effective at determining the location and shape of the scattering objects. The method is general and easily extends to other problems. Because the imaging method involves an elementary set of direct procedures with explicit formulas given at the end of Section 3, it should be useful for a broad variety of intensity-only inverse scattering problems.