On Algorithms Based on Joint Estimation of Currents and Contrast in Microwave Tomography

This paper deals with improvements to the contrast source inversion method which is widely used in microwave tomography. First, the method is reviewed and weaknesses of both the criterion form and the optimization strategy are underlined. Then, two new algorithms are proposed. Both of them are based on the same criterion, similar but more robust than the one used in contrast source inversion. The first technique keeps the main characteristics of the contrast source inversion optimization scheme but is based on a better exploitation of the conjugate gradient algorithm. The second technique is based on a preconditioned conjugate gradient algorithm and performs simultaneous updates of sets of unknowns that are normally processed sequentially. Both techniques are shown to be more efficient than original contrast source inversion.Comment: 12 pages, 12 figures, 5 table

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

On the 3D electromagnetic quantitative inverse scattering problem: algorithms and regularization

In this thesis, 3D quantitative microwave imaging algorithms are developed with emphasis on efficiency of the algorithms and quality of the reconstruction. First, a fast simulation tool has been implemented which makes use of a volume integral equation (VIE) to solve the forward scattering problem. The solution of the resulting linear system is done iteratively. To do this efficiently, two strategies are combined. First, the matrix-vector multiplications needed in every step of the iterative solution are accelerated using a combination of the Fast Fourier Transform (FFT) method and the Multilevel Fast Multipole Algorithm (MLFMA). It is shown that this hybridMLFMA-FFT method is most suited for large, sparse scattering problems. Secondly, the number of iterations is reduced by using an extrapolation technique to determine suitable initial guesses, which are already close to the solution. This technique combines a marching-on-in-source-position scheme with a linear extrapolation over the permittivity under the form of a Born approximation. It is shown that this forward simulator indeed exhibits a better efficiency. The fast forward simulator is incorporated in an optimization technique which minimizes the discrepancy between measured data and simulated data by adjusting the permittivity profile. A Gauss-Newton optimization method with line search is employed in this dissertation to minimize a least squares data fit cost function with additional regularization. Two different regularization methods were developed in this research. The first regularization method penalizes strong fluctuations in the permittivity by imposing a smoothing constraint, which is a widely used approach in inverse scattering. However, in this thesis, this constraint is incorporated in a multiplicative way instead of in the usual additive way, i.e. its weight in the cost function is reduced with an improving data fit. The second regularization method is Value Picking regularization, which is a new method proposed in this dissertation. This regularization is designed to reconstruct piecewise homogeneous permittivity profiles. Such profiles are hard to reconstruct since sharp interfaces between different permittivity regions have to be preserved, while other strong fluctuations need to be suppressed. Instead of operating on the spatial distribution of the permittivity, as certain existing methods for edge preservation do, it imposes the restriction that only a few different permittivity values should appear in the reconstruction. The permittivity values just mentioned do not have to be known in advance, however, and their number is also updated in a stepwise relaxed VP (SRVP) regularization scheme. Both regularization techniques have been incorporated in the Gauss-Newton optimization framework and yield significantly improved reconstruction quality. The efficiency of the minimization algorithm can also be improved. In every step of the iterative optimization, a linear Gauss-Newton update system has to be solved. This typically is a large system and therefore is solved iteratively. However, these systems are ill-conditioned as a result of the ill-posedness of the inverse scattering problem. Fortunately, the aforementioned regularization techniques allow for the use of a subspace preconditioned LSQR method to solve these systems efficiently, as is shown in this thesis. Finally, the incorporation of constraints on the permittivity through a modified line search path, helps to keep the forward problem well-posed and thus the number of forward iterations low. Another contribution of this thesis is the proposal of a new Consistency Inversion (CI) algorithm. It is based on the same principles as another well known reconstruction algorithm, the Contrast Source Inversion (CSI) method, which considers the contrast currents – equivalent currents that generate a field identical to the scattered field – as fundamental unknowns together with the permittivity. In the CI method, however, the permittivity variables are eliminated from the optimization and are only reconstructed in a final step. This avoids alternating updates of permittivity and contrast currents, which may result in a faster convergence. The CI method has also been supplemented with VP regularization, yielding the VPCI method. The quantitative electromagnetic imaging methods developed in this work have been validated on both synthetic and measured data, for both homogeneous and inhomogeneous objects and yield a high reconstruction quality in all these cases. The successful, completely blind reconstruction of an unknown target from measured data, provided by the Institut Fresnel in Marseille, France, demonstrates at once the validity of the forward scattering code, the performance of the reconstruction algorithm and the quality of the measurements. The reconstruction of a numerical MRI based breast phantom is encouraging for the further development of biomedical microwave imaging and of microwave breast cancer screening in particular

Embedding approach to modeling electromagnetic fields in a complex two-dimensional environment

An approach is presented to combine the response of a two-dimensionally inhomogeneous dielectric object in a homogeneous environment with that of an empty inhomogeneous environment. This allows an efficient computation of the scattering behavior of the dielectric cylinder with the aid of the CGFFT method and a dedicated extrapolation procedure. Since a circular observation contour is adopted, an angular spectral representation can be employed for the embedding. Implementation details are discussed for the case of a closed 434 MHz microwave scanner, and the accuracy and efficiency of all steps in the numerical procedure are investigated. Guidelines are proposed for choosing computational parameters such as truncation limits and tolerances. We show that the embedding approach does not increase the CPU time with respect to the forward problem solution in a homogeneous environment, if only the fields on the observation contour are computed, and that it leads to a relatively small increase when the fields on the mesh are computed as well

On the Adjoint Operator in Photoacoustic Tomography

Photoacoustic Tomography (PAT) is an emerging biomedical "imaging from coupled physics" technique, in which the image contrast is due to optical absorption, but the information is carried to the surface of the tissue as ultrasound pulses. Many algorithms and formulae for PAT image reconstruction have been proposed for the case when a complete data set is available. In many practical imaging scenarios, however, it is not possible to obtain the full data, or the data may be sub-sampled for faster data acquisition. In such cases, image reconstruction algorithms that can incorporate prior knowledge to ameliorate the loss of data are required. Hence, recently there has been an increased interest in using variational image reconstruction. A crucial ingredient for the application of these techniques is the adjoint of the PAT forward operator, which is described in this article from physical, theoretical and numerical perspectives. First, a simple mathematical derivation of the adjoint of the PAT forward operator in the continuous framework is presented. Then, an efficient numerical implementation of the adjoint using a k-space time domain wave propagation model is described and illustrated in the context of variational PAT image reconstruction, on both 2D and 3D examples including inhomogeneous sound speed. The principal advantage of this analytical adjoint over an algebraic adjoint (obtained by taking the direct adjoint of the particular numerical forward scheme used) is that it can be implemented using currently available fast wave propagation solvers.Comment: submitted to "Inverse Problems

Expansion of the nodal-adjoint method for simple and efficient computation of the 2d tomographic imaging jacobian matrix

This paper focuses on the construction of the Jacobian matrix required in tomographic reconstruction algorithms. In microwave tomography, computing the forward solutions during the iterative reconstruction process impacts the accuracy and computational efficiency. Towards this end, we have applied the discrete dipole approximation for the forward solutions with significant time savings. However, while we have discovered that the imaging problem configuration can dramatically impact the computation time required for the forward solver, it can be equally beneficial in constructing the Jacobian matrix calculated in iterative image reconstruction algorithms. Key to this implementation, we propose to use the same simulation grid for both the forward and imaging domain discretizations for the discrete dipole approximation solutions and report in detail the theoretical aspects for this localization. In this way, the computational cost of the nodal adjoint method decreases by several orders of magnitude. Our investigations show that this expansion is a significant enhancement compared to previous implementations and results in a rapid calculation of the Jacobian matrix with a high level of accuracy. The discrete dipole approximation and the newly efficient Jacobian matrices are effectively implemented to produce quantitative images of the simplified breast phantom from the microwave imaging system

Microwave Sensing and Imaging

In recent years, microwave sensing and imaging have acquired an ever-growing importance in several applicative fields, such as non-destructive evaluations in industry and civil engineering, subsurface prospection, security, and biomedical imaging. Indeed, microwave techniques allow, in principle, for information to be obtained directly regarding the physical parameters of the inspected targets (dielectric properties, shape, etc.) by using safe electromagnetic radiations and cost-effective systems. Consequently, a great deal of research activity has recently been devoted to the development of efficient/reliable measurement systems, which are effective data processing algorithms that can be used to solve the underlying electromagnetic inverse scattering problem, and efficient forward solvers to model electromagnetic interactions. Within this framework, this Special Issue aims to provide some insights into recent microwave sensing and imaging systems and techniques

A subspace preconditioned LSQR Gauss-Newton method with a constrained line search path applied to 3D biomedical microwave imaging

Three contributions that can improve the performance of a Newton-type iterative quantitative microwave imaging algorithm in a biomedical context are proposed. (i) To speed up the iterative forward problem solution, we extrapolate the initial guess of the field from a few field solutions corresponding to previous source positions for the same complex permittivity (i.e., “marching on in source position”) as well as from a Born-type approximation that is computed from a field solution corresponding to one previous complex permittivity profile for the same source position. (ii) The regularized Gauss-Newton update system can be ill-conditioned; hence we propose to employ a two-level preconditioned iterative solution method. We apply the subspace preconditioned LSQR algorithm from Jacobsen et al. (2003) and we employ a 3D cosine basis. (iii) We propose a new constrained line search path in the Gauss-Newton optimization, which incorporates in a smooth manner lower and upper bounds on the object permittivity, such that these bounds never can be violated along the search path. Single-frequency reconstructions from bipolarized synthetic data are shown for various three-dimensional numerical biological phantoms, including a realistic breast phantom from the University of Wisconsin-Madison (UWCEM) online repository