Highlights of 2015

Recent advances in technology have enabled the combination of positron emission tomography (PET) with magnetic resonance imaging (MRI). These PET-MRI scanners simultaneously acquire functional PET and anatomical or functional MRI data. As function and anatomy are not independent of one another the images to be reconstructed are likely to have shared structures. We aim to exploit this inherent structural similarity by reconstructing from both modalities in a joint reconstruction framework. The structural similarity between two modalities can be modelled in two different ways: edges are more likely to be at similar positions and/or to have similar orientations. We analyse the diffusion process generated by minimizing priors that encapsulate these different models. It turns out that the class of parallel level set priors always corresponds to anisotropic diffusion which is sometimes forward and sometimes backward diffusion. We perform numerical experiments where we jointly reconstruct from blurred Radon data with Poisson noise (PET) and under-sampled Fourier data with Gaussian noise (MRI). Our results show that both modalities benefit from each other in areas of shared edge information. The joint reconstructions have less artefacts and sharper edges compared to separate reconstructions and the ℓ²-error can be reduced in all of the considered cases of under-sampling.

The goal of quantitative elastography is to identify biomechanical parameters from interior displacement data, which are provided by other modalities, such as ultrasound or magnetic resonance imaging. In this paper, we analyze the stability of several linearized problems in quantitative elastography. Our method is based on the theory of redundant systems of linear partial differential equations. We analyze the ellipticity properties of the corresponding PDE systems augmented with the interior displacement data; we explicitly characterize the kernel of the forward operators and show injectivity for particular linearizations. Stability criteria can then be deduced. While joint reconstruction of all parameters suffers from non-ellipticity even for more measurements, our results show stability of the separate reconstruction of shear modulus, pressure and density; they indicate that singular strain fields should be avoided, and show how additional measurements can help in ensuring stability of particular linearized problems.

A limited-view scheme is proposed for multi-source quantitative photoacoustic tomography (MS-QPAT), in which the acoustic measurements following each optical illumination are acquired on the partial boundary near the optical source instead of the entire boundary, namely the limited-view MS-QPAT. The proposed limited-view scheme has an improved signal-to-noise ratio when the data are measured near the optical source, and reduces the acquisition time of the imaging system with a single or limited-view acoustic detector. A limited-view MS-QPAT example is to acquire $4{}^\circ$ acoustic data following each of 90 optical illuminations, in contrast to $360{}^\circ$ acoustic data for each of 90 optical illuminations under the conventional MS-QPAT setting. However, due to the incomplete data, the initial acoustic pressure can no longer be stably reconstructed that serves as an intermediate step in the conventional two-step reconstruction that first reconstructs the initial acoustic pressure and then the optical coefficients. Therefore the direct reconstruction of optical coefficients is considered using the coupled opto-acoustic forward model. The reconstruction algorithm is based on the quasi-Newton method, i.e. limited-memory BFGS with efficient adjoint computations of objective function gradients, and the sparsity-regularized formulation is also considered with tensor framelet sparsity transform and solved by the alternating direction method of multipliers.

We develop a fast algorithm for a Kalman filter applied to the random walk forecast model. The key idea is an efficient representation of the estimate covariance matrix at each time step as a weighted sum of two contributions—the process noise covariance matrix and a low-rank term computed from a generalized eigenvalue problem, which combines information from the noise covariance matrix and the data. We describe an efficient algorithm to update the weights of the preceding terms and the computation of eigenmodes of the generalized eigenvalue problem. The resulting algorithm for the Kalman filter with a random walk forecast model scales as $\mathcal{O}(N)$ in memory and $\mathcal{O}(N{\rm log} N)$ in computational cost, where N is the number of grid points. We show how to efficiently compute measures of uncertainty and conditional realizations from the state distribution at each time step. An extension to the case with nonlinear measurement operators is also discussed. Numerical experiments demonstrate the performance of our algorithms, which are applied to a synthetic example from monitoring CO₂ in the subsurface using travel-time tomography.

In this paper, we propose a new approach for the atmospheric tomography based on the method of the approximate inverse. The image quality of Earth-bound telescopes is severely degraded by turbulences of the atmosphere of the Earth. In order to receive sharp images, the incoming light is corrected for these distortions using deformable mirrors. In atmospheric tomography, the turbulence profile is reconstructed so that the shape of the mirrors can be adjusted in an optimal way. We perform this reconstruction step by applying the method of the approximate inverse, a non-iterative regularization method which requires the solution of an adjoint problem, to Multi-Conjugate Adaptive Optics. We show that the approximate inverse leads to efficient algorithms and give numerical examples. The adjoint problem, which can be precomputed, is solved with an iterative Kaczmarz-type algorithm. In this first study of the new algorithm, tip/tilt effects and spot elongation are not included.

Consider time-harmonic acoustic scattering problems governed by the Helmholtz equation in two and three dimensions. We prove that bounded penetrable obstacles with corners or edges scatter every incident wave nontrivially, provided the function of refractive index is real-analytic. Moreover, if such a penetrable obstacle is a convex polyhedron or polygon, then its shape can be uniquely determined by the far-field pattern over all observation directions incited by a single incident wave. Our arguments are elementary and rely on the expansion of solutions to the Helmholtz equation.

A demanding challenge in Bayesian inversion is to efficiently characterize the posterior distribution. This task is problematic especially in high-dimensional non-Gaussian problems, where the structure of the posterior can be very chaotic and difficult to analyse. Current inverse problem literature often approaches the problem by considering suitable point estimators for the task. Typically the choice is made between the maximum a posteriori (MAP) or the conditional mean (CM) estimate. The benefits of either choice are not well-understood from the perspective of infinite-dimensional theory. Most importantly, there exists no general scheme regarding how to connect the topological description of a MAP estimate to a variational problem. The recent results by Dashti and others (Dashti et al 2013 Inverse Problems 29 095017) resolve this issue for nonlinear inverse problems in Gaussian framework. In this work we improve the current understanding by introducing a novel concept called the weak MAP (wMAP) estimate. We show that any MAP estimate in the sense of Dashti et al (2013 Inverse Problems 29 095017) is a wMAP estimate and, moreover, how the wMAP estimate connects to a variational formulation in general infinite-dimensional non-Gaussian problems. The variational formulation enables to study many properties of the infinite-dimensional MAP estimate that were earlier impossible to study. In a recent work by the authors (Burger and Lucka 2014 Maximum a posteriori estimates in linear inverse problems with logconcave priors are proper bayes estimators preprint) the MAP estimator was studied in the context of the Bayes cost method. Using Bregman distances, proper convex Bayes cost functions were introduced for which the MAP estimator is the Bayes estimator. Here, we generalize these results to the infinite-dimensional setting. Moreover, we discuss the implications of our results for some examples of prior models such as the Besov prior and hierarchical prior.

The inverse problem of MEG aims at estimating electromagnetic cerebral activity from measurements of the magnetic fields outside the head. After formulating the problem within the Bayesian framework, a hierarchical conditionally Gaussian prior model is introduced, including a physiologically inspired prior model that takes into account the preferred directions of the source currents. The hyperparameter vector consists of prior variances of the dipole moments, assumed to follow a non-conjugate gamma distribution with variable scaling and shape parameters. A point estimate of both dipole moments and their variances can be computed using an iterative alternating sequential updating algorithm, which is shown to be globally convergent. The numerical solution is based on computing an approximation of the dipole moments using a Krylov subspace iterative linear solver equipped with statistically inspired preconditioning and a suitable termination rule. The shape parameters of the model are shown to control the focality, and furthermore, using an empirical Bayes argument, it is shown that the scaling parameters can be naturally adjusted to provide a statistically well justified depth sensitivity scaling. The validity of this interpretation is verified through computed numerical examples. Also, a computed example showing the applicability of the algorithm to analyze realistic time series data is presented.

The Mumford–Shah model is a very powerful variational approach for edge preserving regularization of image reconstruction processes. However, it is algorithmically challenging because one has to deal with a non-smooth and non-convex functional. In this paper, we propose a new efficient algorithmic framework for Mumford–Shah regularization of inverse problems in imaging. It is based on a splitting into specific subproblems that can be solved exactly. We derive fast solvers for the subproblems which are key for an efficient overall algorithm. Our method neither requires a priori knowledge of the gray or color levels nor of the shape of the discontinuity set. We demonstrate the wide applicability of the method for different modalities. In particular, we consider the reconstruction from Radon data, inpainting, and deconvolution. Our method can be easily adapted to many further imaging setups. The relevant condition is that the proximal mapping of the data fidelity can be evaluated a within reasonable time. In other words, it can be used whenever classical Tikhonov regularization is possible.

4D cone-beam computed tomography (4DCBCT) reconstructs a temporal sequence of CBCT images for the purpose of motion management or 4D treatment in radiotherapy. However the image reconstruction often involves the binning of projection data to each temporal phase, and therefore suffers from deteriorated image quality due to inaccurate or uneven binning in phase, e.g., under the non-periodic breathing. A 5D model has been developed as an accurate model of (periodic and non-periodic) respiratory motion. That is, given the measurements of breathing amplitude and its time derivative, the 5D model parametrizes the respiratory motion by three time-independent variables, i.e., one reference image and two vector fields. In this work we aim to develop a new 4DCBCT reconstruction method based on 5D model. Instead of reconstructing a temporal sequence of images after the projection binning, the new method reconstructs time-independent reference image and vector fields with no requirement of binning. The image reconstruction is formulated as a optimization problem with total-variation regularization on both reference image and vector fields, and the problem is solved by the proximal alternating minimization algorithm, during which the split Bregman method is used to reconstruct the reference image, and the Chambolle's duality-based algorithm is used to reconstruct the vector fields. The convergence analysis of the proposed algorithm is provided for this nonconvex problem. Validated by the simulation studies, the new method has significantly improved image reconstruction accuracy due to no binning and reduced number of unknowns via the use of the 5D model.

This paper concerns thermoacoustic tomography and photoacoustic tomography, two couple-physics imaging modalities that attempt to combine the high resolution of ultrasound and the high contrast capabilities of electromagnetic waves. We give sufficient conditions to recover both the sound speed of the medium being probed and the source.

Quantitative photoacoustic tomography (QPAT) is a recent hybrid imaging modality that couples optical tomography with ultrasound imaging to achieve high resolution imaging of optical properties of scattering media. Image reconstruction in QPAT is usually a two-step process. In the first step, the initial pressure field inside the medium, generated by the photoacoustic effect, is reconstructed using measured acoustic data. In the second step, this initial ultrasound pressure field datum is used to reconstruct optical properties of the medium. We propose in this work a one-step inversion algorithm for image reconstruction in QPAT that reconstructs the optical absorption coefficient directly from measured acoustic data. The algorithm can be used to recover simultaneously the absorption coefficient and the ultrasound speed of the medium from multiple acoustic data sets, with appropriate a priori bounds on the unknowns. We demonstrate, through numerical simulations based on synthetic data, the feasibility of the proposed reconstruction method.

We propose a new algorithmic approach to the non-smooth and non-convex Potts problem (also called piecewise-constant Mumford–Shah problem) for inverse imaging problems. We derive a suitable splitting into specific subproblems that can all be solved efficiently. Our method does not require a priori knowledge on the gray levels nor on the number of segments of the reconstruction. Further, it avoids anisotropic artifacts such as geometric staircasing. We demonstrate the suitability of our method for joint image reconstruction and segmentation. We focus on Radon data, where we in particular consider limited data situations. For instance, our method is able to recover all segments of the Shepp–Logan phantom from seven angular views only. We illustrate the practical applicability on a real positron emission tomography dataset. As further applications, we consider spherical Radon data as well as blurred data.

Recent years have seen growing interest in exploiting dual- and multi-energy measurements in computed tomography (CT) in order to characterize material properties as well as object shape. Materials characterization is performed by decomposing the scene into constitutive basis functions, such as Compton scatter and photoelectric absorption functions. While well motivated physically, the joint recovery of the spatial distribution of photoelectric and Compton properties is severely complicated by the fact that the data are several orders of magnitude more sensitive to Compton scatter coefficients than to photoelectric absorption, so small errors in Compton estimates can create large artifacts in the photoelectric estimate. To address these issues, we propose a model-based iterative approach which uses patch-based regularization terms to stabilize inversion of photoelectric coefficients, and solve the resulting problem through use of computationally attractive alternating direction method of multipliers (ADMM) solution techniques. Using simulations and experimental data acquired on a commercial scanner, we demonstrate that the proposed processing can lead to more stable material property estimates which should aid materials characterization in future dual- and multi-energy CT systems.

The ill-posedness of the attenuated Radon transform is a challenging issue in practice due to the Poisson noise and the high level of attenuation. The investigation of the smoothing properties of the underlying operator is essential for developing a stable inversion. In this paper, we consider the framework of Sobolev spaces and derive analytically a reconstruction algorithm based on the method of the approximate inverse. The derived method inherits the efficiency and stability of the approximate inverse and supplies a method of extraction of contours. These algorithms appear to be efficient for an attenuation of human body type. However, for higher attenuations the ill-posedness increases exponentially what deteriorates accordingly the quality of reconstructions. Nevertheless, a high attenuation map affects less the contour extraction of a high contrast function and so can be neglected. This leads to simplifying the proposed method and circumvents in this case the artifacts due to the attenuation as attested by simulation results.