A computational geometry method for the inverse scattering problem

07/23/2018 ∙ by Maria L. Daza-Torres, et al. ∙ CIMAT Universidad Complutense de Madrid 0

In this paper we demonstrate a computational method to solve the inverse scattering problem for a star-shaped, smooth, penetrable obstacle in 2D. Our method is based on classical ideas from computational geometry. First, we approximate the support of a scatterer by a point cloud. Secondly, we use the Bayesian paradigm to model the joint conditional probability distribution of the non-convex hull of the point cloud and the constant refractive index of the scatterer given near field data. 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 propose a probability transition kernel that commutes with affine transformations of space. Our findings indicate that our method is reliable to retrieve the support and constant refractive index of the scatterer simultaneously. Indeed, our sampling method is robust to estimate a quantity of interest such as the area of the scatterer. We conclude pointing out a series of generalizations of our method.

READ FULL TEXT VIEW PDF
POST COMMENT

Comments

There are no comments yet.

Authors

This week in AI

Get the week's most popular data science and artificial intelligence research sent straight to your inbox every Saturday.

1 Introduction

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

[1] studied conditions for the well-posedness of the Bayesian formulation of inverse problems. Recently, Christen [2] 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 [5]. The rationale of our approach is that the direct problem has a finite rank property [6]

, 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

[7]. 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 [8]

, 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 [7] 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 [7]),

(1)

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 ,

for .

The Helmholtz equation is equivalent to the Lippmann-Schwinger equation

(2)

where

(3)

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 [9]. Taking the first correction coefficient , we obtain a fourth order approximation

(4)

where ,

, and is Euler’s constant.

The discrete approximation of equation (2) can be written in matrix form

(5)

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 [10]. It is apparent that the method introduced above is considerably faster than the method of Vainikko.

Figure 1: Comparison of computational times between (in seconds) the method of Vainikko (red) and the method propose (blue) for solving the Lippmann-Schwinger equation with the configuration depicted in Example 1 of Section 5.

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. [16] 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:

Data: Let us a set of points within the domain of the problem, , and a positive real number.

Output: A third order B-spline using the point cloud non-convex hull as control points.

Step 1. Compute the -shape of , denoted by , using the -shape algorithm [5].

Step 2. If is a valid polygon (i.e non-empty, non-self-intersected and non-disconnected):
Compute interpolating the points that define the polygon with a third-order B-spline.

Algorithm 1 Spline interpolation of the non-convex hull of a point cloud in 2D

Of note, the -shape algorithm uses the Delaunay triangulation with vertices Q to compute (see Figure 2 (a)).

Figure 2: Part (a) shows the -shape (green), the interpolating spline (blue), , with , and the correspondent Delaunay triangulation (grey). In part (b) shows for different values of .

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

Figure 3: Part (a) shows the interpolating spline , with , when a boundary point is perturbed, while part (b) shows the case when an interior point is perturbed. Local stability of the representation is apparent.

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…" [17]. The solution of the inverse problem is the posterior probability distribution.

Let us consider the following additive noise pointwise observation model:

(6)

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

(7)

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 [17]

(8)

where the likelihood is given by the noise distribution (7),

(9)

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:

  1. The refractive index , the point cloud and the shape parameter are taken as independent parameters; therefore

  2. Uniform distributions for and

    are imposed, as a default selection. Consequently, the prior probabilities

    , are constant on the domain .

  3. The choice of a gamma distribution for

    is motivated by the application of the inverse scattering problem in elasticity imaging methods. The refractive index represents the Young modulus [18]. In this manner, the refractive index is assumed as a positive parameter.

    (10)

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 [8], 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. [19], 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 [5]).

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:

  1. Point move: Move randomly a single point

    where , and is the mean of the pairwise distances

  2. Translate move: Move randomly the point cloud

    where , and is the mean of the pairwise distances

  3. b move: We move the refractive index using its prior distribution (10)

  4. move: The move for the -shape parameter is defined as follows:

(11)

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 [20]. 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:

  1. Move a point ,

  2. Move every point in ,

  3. Draw

  4. Move ,

Step 3. Compute

and accept the proposal with probability .

Step 4.

Step 5. Go back to Step 2 if , else stop.

Algorithm 2 Point Cloud Metropolis-Hastings Random Walk

5 Numerical results and discussion

In the examples presented below we recover a star-shaper scatterer defined by

(12)

Of note, scatterer (12) was introduced by Kirsch [21] 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 [10]

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

(13)

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.

Figure 4: Trace plot of the logarithm of the posterior distribution with a burn-in of 50,000 iterations. It is apparent that the MCMC sampling is at equilibrium.
Figure 5: Parts (a) and (c) show the prior distribution in cyan, the posterior distribution in blue, and the true value in green of both, the area of the scatterer and the refractive index for Example 1. Parts (b) and (d) show a close up along the horizontal axis near the true values.

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.

Figure 6: For examples 5 and 5; posterior distributions of the scatterer area and (b) the refractive index.

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 index

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

Figure 7: Example 5; posterior distributions of the scatterer area and (b) the refractive index.

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
Example 1
Example 2
Example 3
Table 1: 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.
(a) Example 5
(b) Example 5
(c) Example 5
Figure 8: Numerical evidence in Examples 5 and 5 suggests that the posterior distribution is insensitive with respect to a rotation of the direction of the incident fields. On the other hand, Example 5 indicates there is a trade off between precision and computational cost as expected.

6 Conclusions

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.

Acknowledgements

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.

References

  • [1] A. M. Stuart. Inverse problems: A bayesian perspective. Acta Numerica, 19:451–559, 2010.
  • [2] J Andrés Christen. Posterior distribution existence and error control in banach spaces. arXiv preprint arXiv:1712.03299, 2017.
  • [3] 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.
  • [4] 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.
  • [5] 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.
  • [6] M Sh Birman and Mikhail Zakharovich Solomyak. Estimates of singular numbers of integral operators. Russian Mathematical Surveys, 32(1):15, 1977.
  • [7] David Colton and Rainer Kress. Inverse acoustic and electromagnetic scattering theory, volume 93. Springer Science & Business Media, 2012.
  • [8] 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.
  • [9] 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.
  • [10] Gennadi Vainikko. Fast solvers of the lippmann-schwinger equation. In Direct and inverse problems of mathematical physics, pages 423–440. Springer, 2000.
  • [11] Lars Linsen. Point cloud representation. Univ., Fak. für Informatik, Bibliothek, 2001.
  • [12] 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.
  • [13] Mark Pauly, Niloy J Mitra, and Leonidas J Guibas. Uncertainty and variability in point cloud surface data. SPBG, 4:77–84, 2004.
  • [14] Mark Pauly, Leif P Kobbelt, and Markus Gross. Point-based multiscale surface representation. ACM Transactions on Graphics (TOG), 25(2):177–193, 2006.
  • [15] 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.
  • [16] 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.
  • [17] Jari Kaipio and Erkki Somersalo. Statistical and computational inverse problems, volume 160. Springer Science & Business Media, 2006.
  • [18] 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.
  • [19] 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.
  • [20] 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.
  • [21] Andreas Kirsch and R Kress. Two methods for solving the inverse acoustic scattering problem. Inverse problems, 4(3):749, 1988.