Ville Kolehmainen

Statistical inversion and Monte Carlo sampling methods in electrical impedance tomography

This paper discusses the electrical impedance tomography (EIT) problem: electric currents are injected into a body with unknown electromagnetic properties through a set of contact electrodes. The corresponding voltages that are needed to maintain these currents are measured. The objective is to estimate the unknown resistivity, or more generally the impedivity distribution of the body based on this information. The most commonly used method to tackle this problem in practice is to use gradient-based local linearizations. We give a proof for the differentiability of the electrode boundary data with respect to the resistivity distribution and the contact impedances. Due to the ill-posedness of the problem, regularization has to be employed. In this paper, we consider the EIT problem in the framework of Bayesian statistics, where the inverse problem is recast into a form of statistical inference. The problem is to estimate the posterior distribution of the unknown parameters conditioned on measurement data. From the posterior density, various estimates for the resistivity distribution can be calculated as well as a posteriori uncertainties. The search of the maximum a posteriori estimate is typically an optimization problem, while the conditional expectation is computed by integrating the variable with respect to the posterior probability distribution. In practice, especially when the dimension of the parameter space is large, this integration must be done by Monte Carlo methods such as the Markov chain Monte Carlo (MCMC) integration. These methods can also be used for calculation of a posteriori uncertainties for the estimators. In this paper, we concentrate on MCMC integration methods. In particular, we demonstrate by numerical examples the statistical approach when the prior densities are non-differentiable, such as the prior penalizing the total variation or the L1 norm of the resistivity.

Inverse problems with structural prior information

In this paper we propose a method for the regularization of inverse problems whose solutions are known to exhibit anisotropic characteristics. The method is based on the generalized Tikhonov regularization and on the spatial prior information on the underlying solution. We allow the prior information to be only of approximate nature. In the proposed method, the prior information is incorporated into the regularization operator with the aid of a properly constructed matrix-valued field. Although the approach is deterministic it also has a clear statistical interpretation that will be discussed from the Bayesian viewpoint. The method is applied to two examples, the first is the inversion of a Fredholm integral equation of the first kind and the second is a case study of electrical impedance tomography (EIT).

Approximation errors and model reduction with an application in optical diffusion tomography

Model reduction is often required in several applications, typically due to limited available time, computer memory or other restrictions. In problems that are related to partial differential equations, this often means that we are bound to use sparse meshes in the model for the forward problem. Conversely, if we are given more and more accurate measurements, we have to employ increasingly accurate forward problem solvers in order to exploit the information in the measurements. Optical diffusion tomography (ODT) is an example in which the typical required accuracy for the forward problem solver leads to computational times that may be unacceptable both in biomedical and industrial end applications. In this paper we review the approximation error theory and investigate the interplay between the mesh density and measurement accuracy in the case of optical diffusion tomography. We show that if the approximation errors are estimated and employed, it is possible to use mesh densities that would be unacceptable with a conventional measurement model.

Statistical inversion for medical x-ray tomography with few radiographs: I. General theory

In x-ray tomography, the structure of a three-dimensional body is reconstructed from a collection of projection images of the body. Medical CT imaging does this using an extensive set of projections from all around the body. However, in many practical imaging situations only a small number of truncated projections are available from a limited angle of view. Three-dimensional imaging using such data is complicated for two reasons: (i) typically, sparse projection data do not contain sufficient information to completely describe the 3D body, and (ii) traditional CT reconstruction algorithms, such as filtered backprojection, do not work well when applied to few irregularly spaced projections. Concerning (i), existing results about the information content of sparse projection data are reviewed and discussed. Concerning (ii), it is shown how Bayesian inversion methods can be used to incorporate a priori information into the reconstruction method, leading to improved image quality over traditional methods. Based on the discussion, a low-dose three-dimensional x-ray imaging modality is described.

Recovery of region boundaries of piecewise constant coefficients of an elliptic PDE from boundary data

In this study we consider the recovery of smooth region boundaries of piecewise constant coefficients of an elliptic PDE, - a+b = f, from data on the exterior boundary . The assumption made is that the values of the coefficients (a,b) are known a priori but the information about the geometry of the smooth region boundaries where a and b are discontinous is missing. For the full characterization of (a,b) it is then sufficient to find the region boundaries separating different values of the coefficients. This leads to a nonlinear ill-posed inverse problem. In this study we propose a numerical algorithm that is based on the finite-element method and subdivision of the discretization elements. We formulate the forward problem as a mapping from a set of coefficients representing boundary shapes to data on , and derive the Jacobian of this forward mapping. Then an iterative algorithm which seeks a boundary configuration minimizing the residual norm between measured and predicted data is implemented. The method is illustrated first for a general elliptic PDE and then applied to optical tomography where the goal is to find the diffusion and absorption coefficients of the object by transilluminating the object with visible or near-infrared light. Numerical test results for this specific application are given with synthetic data.

Assessment of errors in static electrical impedance tomography with adjacent and trigonometric current patterns

In electrical impedance tomography (EIT), difference imaging is often preferred over static imaging. This is because of the many unknowns in the forward modelling which make it difficult to obtain reliable absolute resistivity estimates. However, static imaging and absolute resistivity values are needed in some potential applications of EIT. In this paper we demonstrate by simulation the effects of different error components that are included in the reconstruction of static EIT images. All simulations are carried out in two dimensions with the so-called complete electrode model. Errors that are considered are the modelling error in the boundary shape of an object, errors in the electrode sizes and localizations and errors in the contact impedances under the electrodes. Results using both adjacent and trigonometric current patterns are given.

Time-series estimation of biological factors in optical diffusion tomography

We apply state space estimation techniques to the time-varying reconstruction problem in optical tomography. We develop a stochastic model for describing the evolution of quasi-sinusoidal medical signals such as the heartbeat, assuming these are represented as a known frequency with randomly varying amplitude and phase. We use the extended Kalman filter in combination with spatial regularization techniques to reconstruct images from highly under-determined time-series data. This system also naturally segments activity belonging to different biological processes. We present reconstructions of simulated data and of real data recorded from the human motor cortex (Franceschini et al 2000 Optics Express 6 49-57). It is argued that the application of these time-series techniques improves both the fidelity and temporal resolution of reconstruction in optical tomography.

Dynamic physiological modeling for functional diffuse optical tomography

Diffuse optical tomography (DOT) is a noninvasive imaging technology that is sensitive to local concentration changes in oxy- and deoxyhemoglobin. When applied to functional neuroimaging, DOT measures hemodynamics in the scalp and brain that reflect competing metabolic demands and cardiovascular dynamics. The diffuse nature of near-infrared photon migration in tissue and the multitude of physiological systems that affect hemodynamics motivate the use of anatomical and physiological models to improve estimates of the functional hemodynamic response. In this paper, we present a linear state-space model for DOT analysis that models the physiological fluctuations present in the data with either static or dynamic estimation. We demonstrate the approach by using auxiliary measurements of blood pressure variability and heart rate variability as inputs to model the background physiology in DOT data. We evaluate the improvements accorded by modeling this physiology on ten human subjects with simulated functional hemodynamic responses added to the baseline physiology. Adding physiological modeling with a static estimator significantly improved estimates of the simulated functional response, and further significant improvements were achieved with a dynamic Kalman filter estimator (paired t tests, n=10, P<0.05). These results suggest that physiological modeling can improve DOT analysis. The further improvement with the Kalman filter encourages continued research into dynamic linear modeling of the physiology present in DOT. Cardiovascular dynamics also affect the blood-oxygen-dependent (BOLD) signal in functional magnetic resonance imaging (fMRI). This state-space approach to DOT analysis could be extended to BOLD fMRI analysis, multimodal studies and real-time analysis.

Statistical inversion for medical x-ray tomography with few radiographs: II. Application to dental radiology

Diagnostic and operational tasks in dental radiology often require three-dimensional information that is difficult or impossible to see in a projection image. A CT-scan provides the dentist with comprehensive three-dimensional data. However, often CT-scan is impractical and, instead, only a few projection radiographs with sparsely distributed projection directions are available. Statistical (Bayesian) inversion is well-suited approach for reconstruction from such incomplete data. In statistical inversion, a priori information is used to compensate for the incomplete information of the data. The inverse problem is recast in the form of statistical inference from the posterior probability distribution that is based on statistical models of the projection data and the a priori information of the tissue. In this paper, a statistical model for three-dimensional imaging of dentomaxillofacial structures is proposed. Optimization and MCMC algorithms are implemented for the computation of posterior statistics. Results are given with in vitro projection data that were taken with a commercial intraoral x-ray sensor. Examples include limited-angle tomography and full-angle tomography with sparse projection data. Reconstructions with traditional tomographic reconstruction methods are given as reference for the assessment of the estimates that are based on the statistical model.

Three-dimensional reconstruction of shape and piecewise constant region values for optical tomography using spherical harmonic parametrization and a boundary element method

We consider the recovery of smooth 3D region boundaries with piecewise constant coefficients in optical tomography. The method is based on a parametrization of the closed boundaries of the regions by spherical harmonic coefficients, and a Newton type optimization process. A boundary integral formulation is used for the forward modelling. The calculation of the Jacobian is based on an adjoint scheme for calculating the corresponding shape derivatives. We show reconstructions for 3D situations. In addition we show the extension of the method for cases where the constant optical coefficients are also unknown. An advantage of the proposed method is the implicit regularization effect arising from the reduced dimensionality of the inverse problem.