Variationally Bayesian Medical Image
Reconstruction Applied to X-ray Tomography?
Jonas Latz1[0000−0002−4600−0247] , Sergio Bacallado1[0000−0002−7193−6450] , Claire
Delplancke2[0000−0001−7483−0419] , Matthias J. Ehrhardt2[0000−0001−8523−353X] ,and Carola-Bibiane Schönlieb1[0000−0003−0099−6306]
1

Centre for Mathematical Sciences, University of Cambrige, Wilberforce Road,Cambridge CB3 0WA, United Kingdom
2
Department of Mathematical Sciences, University of Bath, North Rd, Claverton
Down, Bath BA2 7AY, United Kingdom{jl2160, sb2116, cbs31}@cam.ac.uk, {cd902, me549}@bath.ac.uk

Abstract. We present a Bayesian statistical approach to medical image reconstruction. Here, we obtain not a point estimate for the imageof interest, but a posterior distribution that also quantifies the reconstruction uncertainty. We approximate the posterior variationally. We
illustrate the effectiveness of our approach in numerical experiments inwhich we study variationally Bayesian X-ray tomography.
Keywords: Medical Image Reconstruction · Statistical Learning

1

Statistical Image Reconstruction

Medical image reconstruction problems can often be formulated in the followingway: Identify an image θ† ∈ RNim given a data set y † ∈ RNdt which satisfiesy † = Aθ† + η † ,where A ∈ RNdt ×Nim is a matrix and η † is noise. In, e.g., X-ray tomography A
represents the Radon transform and y † is the sinogram acquired by the scanner.
This image reconstruction problem is usually approached by a (regularised) maximum likelihood method, which gives an estimate θ∗ ∈ RNim that approximates
the true θ† .
Due to noise and/or ill-conditioning of A, the estimate θ∗ may contain aconsiderable amount of remaining uncertainty concerning the true image θ† .
Unfortunately, this uncertainty cannot be quantified through θ∗ itself. Instead,a Bayesian statistical approach can be employed; see [3] for details. Here, theimage θ† is represented by a random variable θ that has a certain probability
distribution µprior := P(θ ∈ ·). This so-called prior represents knowledge and?

The authors acknowledge support from the EPSRC grant EP/S026045/1 “PET++:
Improving Localisation, Diagnosis and Quantification in Clinical and Medical PET
Imaging with Randomised Optimisation”.

2J. Latz et al.

uncertainty in the image before the data y † is observed. Then, Bayes’ formula isemployed to determine the posterior µpost := P(θ ∈ ·|Aθ + η † = y † ), which gives
an accurate description of the knowledge gained through learning the data andof the remaining uncertainty in the reconstruction. In practice, this posterior
usually needs to be approximated computationally.

2

Variational Bayes

To approximate the posterior in our medical image reconstruction problem weemploy the so-called variational Bayes approach; see, e.g., Chapter 10 of [1]: wedefine a set M of probability distributions on RNim and findµvar := argminµ∈M DKL (µkµpost ),(1)where DKL denotes the Kullback–Leibler Divergence.
For the computational approximation of µvar , it is vital that DKL (µkµpost )can be evaluated without first approximating µpost . This is the case, if, e.g.,µprior and µ ∈ M describe the distribution of some f (ϕ), where f is a givenfunction and ϕ is a non-degenerate Gaussian random vector.

3

Eulerian and Lagrangian approaches

Inspired by [2], we consider two different models for µprior and M : a Eulerian
and a Lagragian model. In both cases, we aim to represent the uncertain image
Nim ×Nbabaθ through a system (ψk )N, i.e.
k=1 ∈ R
θ=NbaXξk ψkk=1

for suitable random prefactors ξ1 , . . . , ξNba . The difference of Eulerian and Laba
grangian approach appears in the choice of the system (ψk )N
k=1 : In the Eulerian
Nbaapproach, we choose (ψk )k=1 to be a high-dimensional fixed frame or basis, e.g.,baa system of wavelets or curvelets, and the (ξk )N
k=1 to be sparsity-inducing. In theba
Lagragian approach, we assume to have a low-dimensional system (ψk )N
k=1 thatNba
can be adapted. Hence, not only the prefactors (ξk )k=1 are uncertain, but alsoba
the vectors in (ψk )N
k=1 .
3.1

Particular modelling

Let S(x) = (1 + exp(−x))−1 be a non-linearity acting element-wise on vectors.
Then, we choose the following:
Eulerian approach. We set ξk := ζk S(ζk0 ), k = 1, . . . , Nba . In the prior, we choose0ζ1 , . . . , ζNba , ζ10 , . . . , ζN
∼ N(0, 12 ) independent. In M we consider the samebaba
structure but assume that the means of the (ζk , ζk0 )N
k=1 and the variances of theVariationally Bayesian Medical Image Reconstruction

3Fig. 1. Eulerian variational posterior: Point-wise mean (left column) and variance(right column).

Nbaba(ζk )N
k=1 are degrees of freedom. The resulting prior of the prefactors (ξk )k=1 isa smooth spike-and-slab prior that enforces sparsity.
Lagrangian approach. We set ξk := ζk , ψk = S(Σζ 0 − mk ), k = 1, . . . , Nba , whereΣ ∈ RNim ×Nin is a fixed matrix and m1 < . . . < mNba are fixed scalars. In0∼ N(0, 12 ) independent. In M wethe prior, we choose ζ1 , . . . , ζNba , ζ10 , . . . , ζN
inconsider the same structure but assume that mean and covariance matrix of
0 Ninba(ζk )N
k=1 , as well as the mean and variances of the (ζn )n=1 are degrees of freedom.
The basis functions ψk are obtained through a thresholding of the random object
Σζ 0 , which represents a discretisation of a (smooth) Gaussian process.
Importantly, both models represent non-linear transformations of Gaussian
random variables. Thus, they are suitable for a computationally efficient variational approximation, as discussed in §2.

4

Examples of variationally Bayesian reconstructions

We now give computational examples in which we reconstruct a 50 × 50 Shepp–
Logan phantom based on full-angle X-ray data. We assume to have perturbed
the data with scaled Gaussian white noise and quantify the uncertainty in theimage as if this noise had been applied. The results shown, however, are based onnoisefree data. In each case, we employ suitable stochastic optimisation methods
to solve the variational optimisation problems (1).
ba
In the Eulerian approach, we choose (ψk )N
k=1 to be a frame of 5184 coiflets.
We show the results in Figure 1.
In the Lagrangian approach, we consider two settings. In both settings, wechoose Σ to be a matrix the columns of which are Nin = 2500 orthogonalsinusoidal basis vectors. The basis vectors are weighted with prefactors which go
to zero as the frequency of oscillation of the basis vector increases. In the firstsetting (‘long correlation length’) the prefactors converge quickly, in the secondcase (‘short correlation length’) the prefactors converge slower. In both cases, wechoose Nba = 11 adapted basis vectors. We show the results of our estimationsin Figure 2.

4J. Latz et al.

Fig. 2. Lagrangian variational posterior: Point-wise mean (left column) and variance(right column). Long correlation length (top row) and short correlation length (bottomrow).

In the Eulerian approach, we see an overall rougher appearance of the piecewise constant phantom and a very small posterior variance. The Lagrangianposterior means are a much more accurate reconstruction of the phantom andthe variances appear to highlight areas in which we expect a larger reconstruction
uncertainty. Using the long correlation length appears to give us a smoother image and cleaner edges, but loses small details. The short correlation length blurs
the edges slightly, but retains all fine details.

References
1. Bishop, C.: Pattern recognition and machine learning. Springer-Verlag, New York(2006).
2. Candès, E.J., Donoho, D.L.: Curvelets: A surprisingly effective nonadaptive representation for objects with edges. In: Cohen, A., Rabut, C., Schumaker, L.L. (eds.)
Curve and Surface Fitting: Saint-Malo 1999, Vanderbilt Univ. Press, Nashville, TN(2000).
3. Stuart, A.M.: Inverse Problems: A Bayesian perspective. Acta Numerica 19, 451–559 (2010).