## Abstract

We present a method for inverse scattering that relies on intensity-only measurements of the scattered field on a single measurement plane. By collecting measurements from a suite of experiments in which the sample is illuminated using different incident fields, we create sufficient data diversity to overcome the limitations of the intensity-only measurements. We give an explicit procedure that uses an algebraic relation called the polarization identity to convert intensity measurements of scattered fields to interferometric measurements in which one of the scattered fields serves as the reference. By adjusting the multiple signal classification (MUSIC) method for these interferometric data, we effectively recover the location and shapes of multiple objects contained in the imaging region. This method is effective and robust to noise as long as there is sufficiently high data diversity and the fractional volume of the scattering objects is not too high. We present image reconstructions for several three-dimensional examples with simulated data computed using the Method of Fundamental Solutions that demonstrate the effectiveness of this imaging method.

© 2019 Optical Society of America under the terms of the OSA Open Access Publishing Agreement

## 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

and where $ \textrm{i} = \sqrt { - 1} $ denotes the imaginary constant, $ \textrm{Re}[ \cdot ] $ denotes the real part of the expression, $ \textrm{Im}[ \cdot ] $ denotes the imaginary part of the expression, and $ z_1^ * $ denotes the complex conjugate of $ {z_1} $. By using Eqs. (1) and (2) to solve for $ z_1^ * {z_2} = \textrm{Re}[z_1^ * {z_2}] + \textrm{i}\textrm{Im}[z_1^ * {z_2}] $, we find thatLet $ {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.

- 1. Illuminate the medium with $u_m^i $ for $ m = 1, \cdots ,M $ and measure $ |u_m^s(\textbf{r}{)|^2} = |Lu_m^i(\textbf{r}{)|^2} $ for $ m = 1, \cdots ,M $.
- 2. Illuminate the medium with $ u_1^i + u_m^i $ for $ m = 2, \cdots ,M $ and measure $ |u_1^s(\textbf{r}) + u_m^s(\textbf{r}{)|^2} = |Lu_1^i(\textbf{r}) + Lu_m^i(\textbf{r}{)|^2} $ for $ m = 2, \cdots ,M $.
- 3. Illuminate the medium with $ u_1^i + \textrm{i}u_m^i $ for $ m = 2, \cdots ,M $ and measure $ |u_1^s(\textbf{r}) + \textrm{i}u_m^s(\textbf{r}{)|^2} = |Lu_1^i(\textbf{r}) + \textrm{i}Lu_m^i(\textbf{r}{)|^2} $ for $ m = 2, \cdots ,M $.

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

*et al.*[22]. They show that MUSIC provides the exact support of $ x $ for Eq. (10) when the data are noiseless and remains robust with respect to additive noise, provided the diversity in the data is high enough.

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

with $ B $ given in Eq. (6). This linear system has a special structure in which the unknown vector $ x $ is multiplied by the $ N $ diagonal matrices $ {\Lambda _n} $ for $ n = 1, \cdots ,N $. According to the formulation above, $ {\Lambda _n} $ contains the unknown multiplicative constant, $ u_1^{s * }({\textbf{r}_n}) $.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

with $ {I_M} $ denoting the $ M \times M $ identity matrix. The matrix $ P $ is a projection onto the subspace orthogonal to that spanned by the columns of $ \tilde U $. When we apply $ P $ to the $ k $-th column of $ A $ and find that its length is small, then that column lies in the signal subspace. Let $ {\eta _k} \,=\, \parallel\! P{a_k}\parallel $ for $ k = 1, \cdots ,K $, with $ {a_k} $ denoting the $ k $-th column of $ A $ and let $ {\eta _{\min}} = \min {\eta _k} $. We form an image by plotting The peaks of $ {I_k} $ correspond to the scattering objects. Note that the unknown constants, $ u_1^{s * }({\textbf{r}_n}) $ for $ n = 1, \cdots ,N $, do not affect the signal subspace. Consequently, the imaging method does not rely on their knowledge, and that is why we can image with interferometric measurements with respect to an unknown reference field.To summarize the imaging method, we give the following procedure for forming images.

- 1. Compute the singular value decomposition, $ B = U\Sigma {V^H} $, and determine the significant diagonal entries of $ \Sigma $, denoted by $ {\sigma _j} $, satisfying $ {\sigma _j} \gt \delta \mathop {\max }\nolimits_j {\sigma _j} $ for $ j = 1, \cdots ,\min(M,N) $ with $ \delta $ denoting a user-defined tolerance.
- 2. Compute the projection $ P = {I_M} - \tilde U{\tilde U^H} $ with $ \tilde U $ denoting the columns of $ U $ associated with the significant singular values determined in the first step.
- 3. Compute $ {\eta _k} \,=\, \parallel\!P{a_k}\!\parallel $ for $ k = 1, \cdots ,K $ with $ {a_k} $ denoting the $ k $-th column of $ A $, whose entries are given in Eq. (11), and determine $ {\eta _{\min}} = \mathop {\min }\nolimits_k {\eta _k} $.
- 4. Plot the values of $ {I_k} = {\eta _{\min}}/{\eta _k} $ for $ k = 1, \cdots ,K $.

## 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

where $ \ell \gt 0 $ is a parameter whose value is to be specified, and $ {\nu _{qp}} $ is the unit outward normal for $ \partial {\Omega _q} $ at $ {\textbf{r}_{qp}} $. These points lie*outside*of $\, {\bar \Omega _q} $. Similarly, we introduce the points These points lie

*inside*of $ \,{\Omega _q} $. Using these points, we approximate the solution of the boundary-value problem (18) as

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.

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.

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.

#### 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

In particular, the principal semi-axes of the ellipsoids we consider are $ a = 632\, \,\textrm{nm} $, $ b = 948\, \,\textrm{nm} $, and $ c = 474\, \,\textrm{nm} $. An image reconstruction for this ellipsoidal scattering object with relative refractive index $ {n_1}/{n_0} = 1.4 $ is shown in Fig. 4. The actual ellipsoidal scattering object is plotted in blue, and the isosurface of the reconstructed image is plotted in red. This result shows that the imaging method accurately recovers the location and shape of the ellipsoidal scattering object.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.

## 6. PRACTICAL CONSIDERATIONS

We discuss several practical considerations of this imaging method.

- • This method is for scalar waves and does not consider the vectorial nature of electromagnetic waves. To apply MUSIC to problems that require the consideration of the polarization state of scattered light, we must use the methods of Ammari
*et al.*[26]. Nonetheless, these scalar wave results give valuable insight into the problem and strongly indicate that it can be extended to the full electromagnetic inverse scattering problem. - • This imaging method recovers the locations and shapes of the scattering objects. It does not accurately recover quantitative information such as the complex refractive index. However, using this imaging method as a first step followed by solving a nonlinear least-squares problem similar to what Chen has done [29] is one possible strategy to extend this method for quantitative imaging applications.
- • This method uses the intensity of the scattered field and not the total field. Hence, it is restricted to dark-field measurements, which is a current limitation.
- • Sufficiently high diversity introduced at the source is crucial for the effectiveness of this method. Far fewer sources with more noise and limited range will adversely affect the reconstructed images.
- • The measurement procedure given in Section 2 based on the polarization identity requires controlling the phase of the incident field with sufficient accuracy, which may not be practically feasible. Nonetheless, this approach indicates that there may be alternative methods for leveraging gains in source diversity to circumvent the need for phase retrieval.
- • The results shown above are for very sparse distributions of scattering objects. However, the method is not limited to only these very sparse examples. The method is limited only when the volume fraction of scatterers is so large that the number of nonzero entries of the grid function $ {X_k} $ introduced in Section 3 becomes comparable to the total number of grid points in the mesh of the imaging region.

## 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.

## Funding

Air Force Office of Scientific Research (FA9550-17-1-0238, FA9550-18-1-0519).

## Acknowledgment

C. Tsogka and A. D. Kim are supported by the Air Force Office of Scientific Research.

## REFERENCES

**1. **J. R. Fienup, “Phase retrieval
algorithms: a personal tour [Invited],”
Appl. Opt. **52**,
45–56 (2013). [CrossRef]

**2. **W. Choi, C. Fang-Yen, K. Badizadegan, S. Oh, N. Lue, R. R. Dasari, and M. S. Feld, “Tomographic phase
microscopy,” Nat. Methods **4**, 717–719
(2007). [CrossRef]

**3. **E. Wolf, “Determination of the
amplitude and the phase of scattered fields by
holography,” J. Opt. Soc. Am. **60**, 18–20
(1970). [CrossRef]

**4. **A. J. Devaney, “Diffraction tomographic
reconstruction from intensity data,”
IEEE Trans. Image Process. **1**,
221–228
(1992). [CrossRef]

**5. **M. H. Maleki, A. J. Devaney, and A. Schatzberg, “Tomographic
reconstruction from optical scattered
intensities,” J. Opt. Soc. Am.
A **9**,
1356–1363
(1992). [CrossRef]

**6. **M. H. Maleki and A. J. Devaney, “Phase-retrieval and
intensity-only reconstruction algorithms for optical diffraction
tomography,” J. Opt. Soc. Am.
A **10**,
1086–1092
(1993). [CrossRef]

**7. **P. S. Carney, E. Wolf, and G. Agarwal, “Diffraction tomography
using power extinction measurements,”
J. Opt. Soc. Am. A **16**,
2643–2648
(1999). [CrossRef]

**8. **G. Gbur and E. Wolf, “Hybrid diffraction
tomography without phase information,”
J. Opt. Soc. Am. A **19**,
2194–2202
(2002). [CrossRef]

**9. **E. A. Marengo, R. D. Hernandez, and H. Lev-Ari, “Intensity-only
signal-subspace-based imaging,” J. Opt.
Soc. Am. A **24**,
3619–3635
(2007). [CrossRef]

**10. **M. D’Urso, K. Belkebir, L. Crocco, T. Isernia, and A. Litman, “Phaseless imaging with
experimental data: facts and challenges,”
J. Opt. Soc. Am. A **25**,
271–281
(2008). [CrossRef]

**11. **D. J. Brady, K. Choi, D. L. Marks, R. Horisaki, and S. Lim, “Compressive
holography,” Opt. Express **17**, 13040–13049
(2009). [CrossRef]

**12. **W. Zhang, L. Cao, D. J. Brady, H. Zhang, J. Cang, H. Zhang, and G. Jin, “Twin-image-free
holography: a compressive sensing approach,”
Phys. Rev. Lett. **121**,
093902 (2018). [CrossRef]

**13. **W. Tahir, U. S. Kamilov, and L. Tian, “Holographic particle
localization under multiple scattering,”
Adv. Photonics **1**,
036003 (2019). [CrossRef]

**14. **A. Chai, M. Moscoso, and G. Papanicolaou, “Array imaging using
intensity-only measurements,” Inverse
Prob. **27**, 015005
(2010). [CrossRef]

**15. **T.-A. Pham, E. Soubies, A. Goy, J. Lim, F. Soulez, D. Psaltis, and M. Unser, “Versatile
reconstruction framework for diffraction tomography with intensity
measurements and multiple scattering,”
Biomed. Opt. Express **26**,
2749–2763
(2018). [CrossRef]

**16. **L. Tian and L. Waller, “3D intensity and phase
imaging from light field measurements in an LED array
microscopy,” Optica **2**, 104–111
(2015). [CrossRef]

**17. **R. Horstmeyer, J. Chung, X. Ou, G. Zheng, and C. Yang, “Diffraction tomography
with Fourier ptychography,”
Optica **3**,
827–835
(2016). [CrossRef]

**18. **R. Ling, W. Tahir, H.-Y. Lin, H. Lee, and L. Tian, “High-throughput
intensity diffraction tomography with a computational
microscope,” Biomed. Opt.
Express **9**,
2130–2141
(2018). [CrossRef]

**19. **A. Novikov, M. Moscoso, and G. Papanicolaou, “Illumination strategies
for intensity-only imaging,” SIAM J.
Imaging Sci. **8**,
1547–1573
(2015). [CrossRef]

**20. **M. Moscoso, A. Novikov, and G. Papanicolaou, “Coherent imaging
without phases,” SIAM J. Imaging
Sci. **9**,
1689–1707
(2016). [CrossRef]

**21. **M. Moscoso, A. Novikov, G. Papanicolaou, and C. Tsogka, “Multifrequency
interferometric imaging with intensity-only
measurements,” SIAM J. Imaging
Sci. **10**,
1005–1032
(2017). [CrossRef]

**22. **M. Moscoso, A. Novikov, G. Papanicolaou, and C. Tsogka, “Robust multifrequency
imaging with MUSIC,” Inverse
Prob. **35**, 015007
(2018). [CrossRef]

**23. **H. Ammari, E. Iakovleva, and D. Lesselier, “A MUSIC algorithm for
locating small inclusions buried in a half-space from the scattering
amplitude at a fixed frequency,”
Multiscale Model. Simul. **3**,
597–628
(2005). [CrossRef]

**24. **F. K. Gruber, E. A. Marengo, and A. J. Devaney, “Time-reversal imaging
with multiple signal classification considering multiple scattering
between the targets,” J. Acoust. Soc.
Am. **115**,
3042–3047
(2004). [CrossRef]

**25. **A. Chai, M. Moscoso, and G. Papanicolaou, “Imaging strong
localized scatterers with sparcity promoting
optimization,” SIAM J. Imaging
Sci. **7**,
1358–1387
(2014). [CrossRef]

**26. **H. Ammari, E. Iakovleva, D. Lesselier, and G. Perrusson, “MUSIC-type
electromagnetic imaging of a collection of small three-dimensional
bounded inclusions,” SIAM J. Sci.
Comput. **29**,
674–709
(2007). [CrossRef]

**27. **E. Iakovleva, S. Gdourra, D. Lesselier, and G. Perrusson, “Multistatic response
matrix of a 3-D inclusion in half space and MUSIC
imaging,” IEEE Trans. Antennas
Propag. **55**,
2598–2609
(2007). [CrossRef]

**28. **X. Chen, “Signal-subspace method
approach to the intensity-only electromagnetic inverse scattering
problem,” J. Opt. Soc. Am. A **25**, 2018–2024
(2008). [CrossRef]

**29. **X. Chen, “Application of
signal-subspace and optimization methods in reconstructing extended
scatterers,” J. Opt. Soc. Am.
A **26**,
1022–1026
(2009). [CrossRef]

**30. **G. Zheng, C. Kolner, and C. Yang, “Microscopy refocusing
and dark-field imaging by using a simple LED
array,” Opt. Lett. **36**, 3987–3989
(2011). [CrossRef]

**31. **P. Jordan and J. von Neumann, “On inner products in
linear, metric spaces,” Ann.
Math. **36**,
719–723
(1935). [CrossRef]

**32. **G. Beylkin, “Imaging of
discontinuities in the inverse scattering problem by inversion of a
causal generalized radon transform,” J.
Math. Phys. **26**,
99–108
(1985). [CrossRef]

**33. **G. Beylkin, “Reconstructing
discontinuities in multidimensional inverse scattering problems:
smooth errors vs small errors,” Appl.
Opt. **24**,
4086–4088
(1985). [CrossRef]

**34. **R. Mathon and R. L. Johnston, “The approximate
solution of elliptic boundary-value problems by fundamental
solutions,” SIAM J. Numer.
Anal. **14**,
638–650
(1977). [CrossRef]

**35. **T. Wriedt, A. Doicu, and Y. Eremin, *Acoustic and Electromagnetic
Scattering Analysis Using Discrete Sources*
(Academic,
2000).

**36. **A. D. Kim,
“IntensityOnlyImagingwMUSIC,” https://github.com/arnolddkim/IntensityOnlyImagingwMUSIC.