MCMC and variational approaches for Bayesian inversion in diffraction imaging

International audienceThe term “diffraction imaging” is meant, herein, in the sense of an “inverse scattering problem” where the goal is to build up an image of an unknown object from measurements of the scattered field that results from its interaction with a known probing wave. This type of problem occurs in many imaging and non-destructive testing applications. It corresponds to the situation where looking for a good trade-off between the image resolution and the penetration of the incident wave in the probed medium, leads to choosing the frequency of the latter in such a way that its wavelength lies in the “resonance” domain, in the sense that it is approximately of the same order of magnitude as the characteristic dimensions of the inhomogeneities of the inspected object. In this situation the wave-object interaction gives rise to important diffraction phenomena. This is the case for the two applications considered herein, where the interrogating waves are electromagnetic waves with wavelengths in the microwave and optical domains, whereas the characteristic dimensions of the sought object are 1 cm and 1 μm, respectively.The solution of an inverse problem obviously requires previous construction of a forward model that expresses the scattered field as a function of the parameters of the sought object. In this model, diffraction phenomena are taken into account by means of domain integral representations of the electric fields. The forward model is then described by two coupled integral equations, whose discrete versions are obtained using a method of moments and whose inversion leads to a non-linear problem.Concerning inversion, at the beginning of the 1980s, accounting for the diffraction phenomena has been the subject of much attention in the field of acoustic imaging for applications in geophysics, non-destructive testing or biomedical imaging. It led to techniques such as diffraction tomography, a term that denotes “applications that employs diffracting wavefields in the tomographic reconstruction process” , but which generally implies reconstruction processes based on the generalized projection-slice theorem, an extension to the diffraction case of the projection-slice theorem of the classical computed tomography whose forward model is given by a Radon transform . This theorem is based upon first- order linearizing assumptions such as the Born’s or Rytov’s approximations. So, the term diffraction tomography was paradoxically used to describe reconstruction techniques adapted to weakly scattering environments that do not provide quantitative information on highly contrasted dielectric objects such as those encountered in the applications considered herein, where multiple diffraction cannot be ignored.Furthermore, the resolution of these techniques is limited because evanescent waves are not taken into consideration. These limitations have led researchers to develop inversion algorithms able to deal with non-linear problems, at the beginning of the 1990s for microwave imaging and more recently for optical imaging. Many studies have focused on the development of deterministic methods, such as the Newton-Kantorovich algorithm, the modified gradient method (MGM) or the contrast-source inversion technique (CSI), where the solution is sought for by means of an iterative minimization by a gradient method of a cost functional that expresses the difference between the scattered field and the estimated model output. But, in addition to be non-linear, inverse scattering problems are also known to be ill-posed, which means that their resolution requires a regularization which generally consists in introducing prior information on the sought object. In the present case, for example, we look for man-made objects that are composed of homogeneous and compact regions made of a finite number of different materials, and with the aforementioned deterministic methods, it is not easy to take into account such prior information because it must be introduced into the cost functional to be minimized.On the contrary, the probabilistic framework of Bayesian estimation, basis of the model presented herein, is especially well suited for this situation. Prior information is appropriately introduced via a probabilistic Gauss-Markov-Potts model. The marginal contrast distribution is modeled as a mixture of Gaussians, where each Gaussian distribution represents a class of materials and the compactness of the regions is taken into account using a hidden Markov model. Estimation of the unknowns and parameters introduced into the prior model is performed via an unsupervised joint approach.Two iterative algorithms are proposed. The first one, denoted as the MCMC algorithm (Monte-Carlo Markov Chain), is rather classic ; it consists in expressing all the joint posterior or conditional distributions of all the unknowns and, then, using a Gibbs sampling algorithm for estimating the posterior mean of the unknowns. This algorithm yields good results, however, it is computationally intensive mainly because Gibbs sampling requires a significant number of samples.The second algorithm is based upon the variational Bayesian approximation (VBA). The latter was first introduced in the field of Bayesian inference for applications to neural networks, learning graphic models and model parameter estimation. Its appearance in the field of inverse problems is relatively recent, starting with source separation and image restoration. It consists in approximating the joint posterior distribution of all the unknowns by a free-form separable distribution that minimizes, with respect to the posterior law, the Kullback-Leibler divergence which has interesting properties for optimization and leads to an implicit parametric optimization scheme. Once the approximate distribution is built up, the estimator can be easily obtained.A solution to this functional optimization problem can be found in terms of exponential distributions whose shape parameters are estimated iteratively. It can be noted that, at each iteration, the updating expression for these parameters is similar to the one that could be obtained if a gradient method was used to solve the optimization problem. Moreover, the gradient and the step size have an interpretation in terms of statistical moments (means, variances, etc.).Both algorithms introduced herein are applied to two quite different configurations. The one related to microwave imaging is quasi-optimal: data are quasi-complete and frequency diverse. This means that the scattered fields are measured all around the object for several directions of illumination and several frequencies. The configuration used in optical imaging is less favorable since only aspect-limited data are available at a single frequency. This means that illuminations and measurements can only be performed in a limited angular sector. This limited aspect reinforces the ill-posedness of the inverse problem and makes essential the introduction of prior information. However, it will be shown that, in both cases, satisfactory results are obtained