In this paper we demonstrate a computational method to solve the inverse scattering problem for a star-shaped, smooth, penetrable obstacle in 2D. Broadly speaking, inverse problems in physical systems are statistical inference problems restricted by an initial and/or boundary value problem for a partial differential equation. Consequently, inverse problems are difficult since the quantity to be inferred is defined in a function space, while data is finite and noisy. Stuart studied conditions for the well-posedness of the Bayesian formulation of inverse problems. Recently, Christen  showed that that the Bayesian formulation of inverse problems is well posed under more general conditions. In either case, Bayesian inverse problems defined on function spaces have two main difficulties: the computationally efficient representation of the function to be estimated and the sampling of a posterior measure defined on an infinite-dimensional space. Although there are many theoretical approaches to overcome this issues [3, 4, 1], the development of competitive computational methods for inverse problems remains a very active research field.
We rely on classical computational geometry methods to deal with the issues of function representation and sampling. First, instead of posing the problem in an infinite dimensional setting, we approximate the support of each scatterer using a point cloud . The rationale of our approach is that the direct problem has a finite rank property 
, i.e., the Frechet derivative of the mapping from scattering obstacle to near field is a compact operator, and its singular values have zero as a limit point. Secondly, we use the Bayesian paradigm to model the joint conditional probability of the non-convex hull of the point cloud and the scatterer constant refractive index given measurements of the near field (in the remainder of this paper, we shall refer to this conditional probability as the posterior distribution). Of note, we use the non-convex hull of the point cloud as spline control points to evaluate, on a finer mesh, the volume potential arising in the integral equation formulation of the direct problem. Finally, in order to sample the arising posterior distribution, we construct a probability transition kernel based on moves that commute with affine transformations of space. This property, first introduced by Christen and Fox 
, implies that the performance of the Markov Chain Monte Carlo (MCMC) method is independent of the aspect ratio of anisotropic distributions.
The outline of the paper is as follows. In Section 2, we introduce the direct scattering problem and propose a numerical solution. Section 3 is devoted to the representation of the shape of an obstacle based on a point cloud. This representation uses the
-shape algorithm and allows us to parameterize the obstacle shape with a parameter and a cloud of points. In Section 4, we present the Bayesian formulation of the inverse problem. The solution of the inverse problem is the posterior probability distribution. Thereby, an affine invariant MCMC method to sample the posterior distribution is designed. Numerical results and discussion are presented in Section 5. Finally, Section 6 summarizes our findings, pointing out limitations and possible generalizations of our inference method.
2 The direct scattering problem
In this Section we describe the direct scattering problem and propose a method for its numerical solution. We refer the reader to  for a more detailed discussion about the mathematical formulation of the scattering problem for a penetrable obstacle.
Let us assume that the scatterer is bounded and penetrable. If the incident wave is a plane wave, the acoustic scattering problem can be cast as the following problem for the Helmholtz equation (see ),
where is the total field, is the incident field, is the scattered field, , and is the wavenumber. The refractive index is assumed to be constant on ,
The Helmholtz equation is equivalent to the Lippmann-Schwinger equation
and is the Hankel function of the first class of order zero.
2.1 Numerical solution of the direct scattering problem
In order to obtain an efficient solution of the forward model, we approximate the volume potential in equation (2) using a corrected trapezoidal rule. The correction coefficients are calculated in . Taking the first correction coefficient , we obtain a fourth order approximation
and is Euler’s
The discrete approximation of equation (2) can be written in matrix form
Note that is a sparse matrix that depends on the obstacle and the refractive index. Of note, we can readily index the entries where the matrix is different from zero using the definition of . Thus, the matrix structure allows us to solve the system efficiently (5). In Figure 1, we compare computational times to solve the Lippmann-Schwinger equation numerically using equation (5) with the numerical solution of Vainikko . It is apparent that the method introduced above is considerably faster than the method of Vainikko.
3 Point cloud approximation of a scatterer
Given a point cloud in two or three dimensions, many methods of computational geometry are available to represent its boundary [11, 12, 13, 14, 15]. The structural simplicity of these representations allows efficient handling for geometric deformations and topology changes. The results of this paper are based on these ideas. In order to construct this representation, we use the -shape algorithm that computes the non-convex hull of a set of points in the plane. This algorithm produces a polygonal hull (-shape) that depends on the value of a real parameter (called
). To obtain a smooth boundary, we interpolate on the points that define the-shape with a third order B-spline. B-splines are used due to its local stability properties on the control points, see Figure 3 (a), (b). Palafox et al.  introduced this procedure to estimate the boundary of a smooth sound obstacle from far field data. A detailed description of the algorithm is as follows:
Of note, the -shape algorithm uses the Delaunay triangulation with vertices Q to compute (see Figure 2 (a)).
The algorithm defined above accounts for a well defined mapping from and parameter to a boundary . Indeed, the polygon , and thus the smooth boundary , is uniquely determined by . The positive parameter determines how non-convex the smooth boundary is (see Figure 2 (b)).
4 Bayesian Formulation of the Inverse Problem
In the inverse scattering problem, the task is to estimate the shape of the scatterer and the refractive index given some available observations, e.g., near field data observed at some points of the domain . The scattered field is obtained from several incoming plane waves with different directions . We shall formulate the inverse problem in a finite dimensional space using the parametric representation of the obstacle described above. Thereby, the parameter space to be estimated is .
We shall pose the inverse scattering problem using the Bayesian approach.
This formulation allows us restating the inverse problem as “…a well-posed
extension in a larger space of probability distributions…" .
The solution of the inverse problem is the posterior probability distribution.
Let us consider the following additive noise pointwise observation model:
where , is the set of points at which the scattered field is observed, is additive noise, and noisy observations. Concatenating all the observation, one can rewrite (6) as
with denoting the mapping from the shape of (defined by Q and ) and the refractive index to the noise free observables,
being normally distributed as, and .
The Bayesian solution to our inverse problem is 
where the likelihood is given by the noise distribution (7),
assuming independent measurements. This function shows how measurements
affect knowledge about .
The prior distribution expresses the information about independent of the measurements . In our case, it is modeled as follows:
The refractive index , the point cloud and the shape parameter are taken as independent parameters; therefore
The resulting posterior distribution is explored using MCMC method with Metropolis Hastings (MH) by generating a Markov chain with equilibrium distribution (8). In the next Section we propose a the probability transition kernel to implement the MH algorithm.
4.1 An affine invariant probability transition kernel
By construction, the performance of the Metropolis-Hastings (MH) algorithm depends on a proposal probability distribution used to construct the probability transition kernel. Often, proposal probability distributions are gradient informed [4, 3]. However, there is a class of derivative-free proposals, first introduced by Christen and Fox , for arbitrary continuous distributions and requires no tuning, that gives rise to “…an algorithm that is invariant to scale, and approximately invariant to affine transformations of the state space…”. This property implies that the performance of the MCMC method is independent of the aspect ratio of anisotropic target distributions.
In Daza et al. , a MCMC method with the affine invariant property was used to approach the inverse scattering problem in an inhomogeneous medium with a convex obstacle. In this paper we build on previous results introducing an affine invariant proposal also for the parameter of the non-convex hull of a point cloud. This proposal is defined from the circumradius of the Delaunay triangulation of the point cloud at the iteration t. The circumradius is used for two reasons: it contains information about the scale of the space and the -shape is obtained with decision criteria that involve the radius of the Delaunay triangulation (see details in ).
The probability transition kernel for the parameter space , where is given by
and is the weight to choose the proposal . Thus, the new proposal is defined as follows:
Point move: Move randomly a single point
where , and is the mean of the pairwise distances
Translate move: Move randomly the point cloud
where , and is the mean of the pairwise distances
b move: We move the refractive index using its prior distribution (10)
move: The move for the -shape parameter is defined as follows:
where , are the minimum and maximum circumradius
of the Delaunay triangulation, respectively. The value is
taken in , because for values of
the -shape of is the convex hull and for values of
the -shape is the set of points .
The instrumental proposals for the obstacle boundary (1-2) are symmetric because when computing the pairwise distances without the point in the point move, we produce a symmetric proposal on the MCMC method. Similarly, when the translation move is chosen, the entire cloud is moved with the same direction , and therefore the proposal is also symmetric . A detailed description of the algorithm is as follows:
Data: . uniformly distributed on .
Number of samples
Output: A prescribed number of samples of the posterior distribution (8) obtained
with the affine invariant sampler.
Step 1. Take . Compute the energy
Step 2. Update randomly choosing one of the moves:
Move a point ,
Move every point in ,
Step 3. Compute
and accept the proposal with probability .
Step 5. Go back to Step 2 if , else stop.
5 Numerical results and discussion
In the examples presented below we recover a star-shaper scatterer defined by
Of note, scatterer (12) was introduced by Kirsch  and is often used in the inverse scattering literature as a benchmark. Recovering scatterer (12) is challenging for two reasons: first because it is non-convex and secondly, because it possesses detailed structures (the wings) that are small compared with the scale of the whole figure. In the setting for our numerical examples we take fixed refractive index . Furthermore, we use the area of the scatterer as a quantity of interest to explore the robustness of our inference method.
Data is measured at grid points of the numerical domain. In order to avoid the inverse crime we solve the direct scattering problem with the method introduced by Vainikko 
to create synthetic data. We add Gaussian noise with mean zero and standard deviation, i.e., the signal to noise ratio (SNR) is 100.
We propose a design on the distribution of the incident field over the unit circle. This design involves incident fields in two wavenumbers (, ). We take eight incident wave directions uniformly distributed
In directions , , the incident field is applied with wavenumber and in the other directions with wavenumber . Data is assumed independent, then the likelihood is given by
where and denote the scattered field for the wavenumbers and respectively.
In Example 5 we take for the configuration of the incident field directions (13). The mesh of the numerical domain is and . In Example 5 we explore the robustness of our method by rotating the direction of the incident fields with respect to Example 5. Finally, in the third example we verify that by refining the mesh in the solution of the direct problem, the results improve accordingly. All numerical experiments were performed in python. For the sake of reproducibility, all code is available in Github.
We computed 2,000,000 samples of the posterior distribution. Figure 4 shows a trace plot of samples of the logarithm of the posterior distribution at equilibrium, indicating convergence in measure of the MCMC method. On the other hand, histograms of the prior and posterior distributions of both, the refractive index and the area of the obstacle are shown in Figure 5. True values are shown for comparison.
In order to show that the results obtained in Example 5 are independent of the directions of the incident field, in the next example we repeat the numerical experiment except that none of the incident directions coincide with an axis of symmetry of the obstacle.
We consider the same setting of Example 1, except that the directions in the incident field are rotated letting . Figure 6 compares histograms of both, the refractive index and the area of the obstacle obtained in Examples 5 and 5.
Finally, Example 5
illustrates the effect of using a finer mesh in the numerical solution of the direct problem. As expected, the variance of both, the refractive indexand the quantity of interest decreases if the mesh is reduced.
In this example, the configuration is taken as in Example 5, and we refined the mesh taking , . We computed 1,000,000 samples of the posterior distribution. Histograms of the posterior of both, the refractive index and the area of the obstacle are shown in Figure 7 together with the respective results of Example 5. As before, the true values are shown for comparison. Figure 8 summarizes the retrieval of the scatterer in all examples above. Numerical evidence suggest that the posterior distribution is insensitive with respect to a rotation of the direction of the incident fields. On the other hand, there is a trade off between precision and computational cost as expected. Table 1 compares the true value, maximum a posteriori and conditional mean of our quantity of interest. All estimators are correct with roughly two digits. Of note, the area of the true scatterer was computed with the same resolution of the point cloud interpolation.
We may discuss the results of this manuscript as follows: we have introduced a method, based on standard ideas of computational geometry, capable of simultaneously recovering both, the support and constant refractive index of a star-shaped scatterer. In Figures 4 through 8 and Table 1, we offer numerical evidence of the robustness of our method using a star-shaped scatterer commonly used as benchmark in the literature of inverse scattering problems under a particular data design.
On the other hand, although our method does accomplish the task of reliably solving the inverse problem, it is computationally expensive. We remark that our method is derivative free. We believe that further research on affine invariant proposals is justified i.e., a natural extension is to explore the performance of derivative informed proposals of the probability transition kernel that fulfill the affine invariant property. In general, we believe that it might be possible to further extend the idea of axiomatically imposing geometric properties on move proposals for the probability transition kernel, such as moves that commute with affine transformations of space, and then introducing examples of moves that fulfill those conditions.
|Area of the scatterer|
|True value||Conditional Mean||Maximum a posteriori|
In this paper we have demonstrated a reliable method to approximate the solution of the inverse scattering problem in the case of penetrable obstacles. Of note, we can recover simultaneously the support and constant refractive index of the scatterer. Our results are based on standard ideas of computational geometry, which is a very active area of computer science with numerous applications. We consider that our findings are amenable to numerous extensions.
We believe that future work should be directed towards improving the efficiency of the method in terms of computational cost. Three possible research avenues, not necessarily in order of importance, are as follows. First, explore the performance of derivative informed proposals of the probability transition kernel that fulfill the affine invariant property. Second, defining a delayed acceptance method on the frecuency (lower wavenumber and higher wavenumber) spaces with two stages; the lower-model with a coarse mesh and the higher-model with a fine mesh. Finally, another approach might be, changing the number of points in the cloud in each level (delayed acceptance), for instance adding points in the cloud, taking more points for the interpolation of the point cloud.
The authors would like to thank Colin Fox, Robert Scheichl, Tony Shardlow and Antonio Capella for useful discussions and insight. J. A. I. R. acknowledges the financial support of the Spanish Ministry of Economy and Competitiveness under the project MTM2015-64865-P and the research group MOMAT (Ref.910480) supported by Universidad Complutense de Madrid.
-  A. M. Stuart. Inverse problems: A bayesian perspective. Acta Numerica, 19:451–559, 2010.
-  J Andrés Christen. Posterior distribution existence and error control in banach spaces. arXiv preprint arXiv:1712.03299, 2017.
-  Tan Bui-Thanh, Omar Ghattas, James Martin, and Georg Stadler. A computational framework for infinite-dimensional bayesian inverse problems part i: The linearized case, with application to global seismic inversion. SIAM Journal on Scientific Computing, 35(6):A2494–A2523, 2013.
-  Mark Girolami and Ben Calderhead. Riemann manifold langevin and hamiltonian monte carlo methods. Journal of the Royal Statistical Society: Series B, 73(2):123–214, 2011.
-  Herbert Edelsbrunner, David Kirkpatrick, and Raimund Seidel. On the shape of a set of points in the plane. IEEE Transactions on information theory, 29(4):551–559, 1983.
-  M Sh Birman and Mikhail Zakharovich Solomyak. Estimates of singular numbers of integral operators. Russian Mathematical Surveys, 32(1):15, 1977.
-  David Colton and Rainer Kress. Inverse acoustic and electromagnetic scattering theory, volume 93. Springer Science & Business Media, 2012.
-  J Andrés Christen and Colin Fox. A general purpose sampling algorithm for continuous distributions (the t-walk). Bayesian Analysis, 5(2):263–281, 2010.
-  Juan C Aguilar and Yu Chen. High-order corrected trapezoidal quadrature rules for functions with a logarithmic singularity in 2-d. Computers & Mathematics with Applications, 44(8):1031–1039, 2002.
-  Gennadi Vainikko. Fast solvers of the lippmann-schwinger equation. In Direct and inverse problems of mathematical physics, pages 423–440. Springer, 2000.
-  Lars Linsen. Point cloud representation. Univ., Fak. für Informatik, Bibliothek, 2001.
-  Mark Pauly, Richard Keiser, Leif P Kobbelt, and Markus Gross. Shape modeling with point-sampled geometry. In ACM Transactions on Graphics (TOG), pages 641–650. ACM, 2003.
-  Mark Pauly, Niloy J Mitra, and Leonidas J Guibas. Uncertainty and variability in point cloud surface data. SPBG, 4:77–84, 2004.
-  Mark Pauly, Leif P Kobbelt, and Markus Gross. Point-based multiscale surface representation. ACM Transactions on Graphics (TOG), 25(2):177–193, 2006.
-  Matthew Berger, Andrea Tagliasacchi, Lee M Seversky, Pierre Alliez, Gael Guennebaud, Joshua A Levine, Andrei Sharf, and Claudio T Silva. A survey of surface reconstruction from point clouds. In Computer Graphics Forum, pages 301–329. Wiley Online Library, 2017.
-  Abel Palafox, Marcos A Capistrán, and J Andrés Christen. Effective parameter dimension via bayesian model selection in the inverse acoustic scattering problem. Mathematical Problems in Engineering, 2014, 2014.
-  Jari Kaipio and Erkki Somersalo. Statistical and computational inverse problems, volume 160. Springer Science & Business Media, 2006.
-  James F Greenleaf, Mostafa Fatemi, and Michael Insana. Selected methods for imaging elastic properties of biological tissues. Annual review of biomedical engineering, 5(1):57–78, 2003.
-  Maria L Daza, Marcos A Capistrán, J Andrés Christen, and Lilí Guadarrama. Solution of the inverse scattering problem from inhomogeneous media using affine invariant sampling. Mathematical Methods in the Applied Sciences, 2016.
-  Abel Palafox, Marcos Capistrán, and J Andrés Christen. Point cloud-based scatterer approximation and affine invariant sampling in the inverse scattering problem. Mathematical Methods in the Applied Sciences, 40(9):3393–3403, 2017.
-  Andreas Kirsch and R Kress. Two methods for solving the inverse acoustic scattering problem. Inverse problems, 4(3):749, 1988.