---
abstract: 'In this work, we consider the inverse problem of reconstructing the internal structure of an object from limited x-ray projections. We use a Gaussian process prior to model the target function and estimate its (hyper)parameters from measured data. In contrast to other established methods, this comes with the advantage of not requiring any manual parameter tuning, which usually arises in classical regularization strategies. Our method uses a basis function expansion technique for the Gaussian process which significantly reduces the computational complexity and avoids the need for numerical integration. The approach also allows for reformulation of come classical regularization methods as Laplacian and Tikhonov regularization as Gaussian process regression, and hence provides an efficient algorithm and principled means for their parameter tuning. Results from simulated and real data indicate that this approach is less sensitive to streak artifacts as compared to the commonly used method of filtered backprojection.'
address: |
    $^1$Department of Electrical Engineering and Automation, Aalto University, Finland\
    $^2$Department of Information Technology, Uppsala University, Sweden\
author:
- 'Zenith Purisha$^1$, Carl Jidling$^2$, Niklas Wahlstr[ö]{}m$^2$, Thomas B. Sch[ö]{}n$^2$, Simo S[ä]{}rkk[ä]{}$^1$'
bibliography:
- 'sample.bib'
title: 'Probabilistic approach to limited-data computed tomography reconstruction'
---

2019

[*Keywords*]{}: computed tomography, limited data; probabilistic method; Gaussian process; Markov chain Monte Carlo

Introduction
============

X-ray computed tomography (CT) imaging is a non-invasive method to recover the internal structure of an object by collecting projection data from multiple angles. The projection data is recorded by a detector array and it represents the attenuation of the x-rays which are transmitted through the object. Since the 1960s, CT has been used to a deluge of applications in medicine [@cormack1963representation; @cormack1964representation; @herman1979image; @kuchment2013radon; @national1996mathematics; @shepp1978computerized] and industry [@akin2003computed; @cartz1995nondestructive; @de2014industrial].

Currently, the so-called filtered back projection (FBP) is the reconstruction algorithm of choice because it is very fast [@avinash1988principles; @buzug2008computed]. This method requires dense sampling of the projection data to obtain a satisfying image reconstruction. However, for some decades, the limited-data x-ray tomography problem has been a major concern in, for instance, the medical imaging community. The limited data case—also referred to as [*sparse projections*]{}—calls for a good solution for several important reasons, including:

-   the needs to examine a patient by using low radiation doses to reduce the risk of malignancy or to [*in vivo*]{} samples to avoid the modification of the properties of living tissues,

-   geometric restrictions in the measurement setting make it difficult to acquire the complete data [@riis2018limited], such as in [*mammography*]{} [@niklason1997digital; @rantala2006wavelet; @wu2003tomographic; @zhang2006comparative] and electron imaging [@fanelli2008electron],

-   the high demand to obtain the data using short acquisition times and to avoid massive memory storage, and

-   the needs to avoid—or at least minimize the impact of—the moving artifacts during the acquisition.

Classical algorithms—such as FBP—fail to generate good image reconstruction when dense sampling is not possible and we only have access to limited data. The under-sampling of the projection data makes the image reconstruction (in classical terms) an [*ill-posed*]{} problem [@natterer1986mathematics]. In other words, the inverse problem is sensitive to measurement noise and modeling errors. Hence, alternative and more powerful methods are required. Statistical estimation methods play an important role in handling the ill-posedness of the problem by restating the inverse problem as a [*well-posed extension*]{} in a larger space of probability distributions [@kaipio2006statistical]. Over the years there have been a lot of work on tomographic reconstruction from limited data using statistical methods (see, e.g., [@rantala2006wavelet; @bouman1996unified; @haario2017shape; @kolehmainen2003statistical; @siltanen2003statistical; @sauer1994bayesian]). In the statistical approach, incorporation of [*a priori*]{} knowledge is a crucial part in improving the quality of the image reconstructed from limited projection data. That can be viewed as an equivalent of the regularization parameter in classical regularization methods. However, statistical methods, unlike classical regularization methods, also provide a principled means to estimate the parameters of the prior (i.e., the hyperparameters) which corresponds to automatic tuning of regularization parameters.

In our work we build the statistical model by using a Gaussian process model [@Rasmussen2006] with a hierarchical prior in which the (hyper)parameters in the prior become part of the inference problem. As this kind of hierarchical prior can be seen as an instance of a Gaussian process (GP) regression model, the computational methods developed for GP regression in machine learning context [@Rasmussen2006] become applicable. It is worth noting that some works on employing GP methods for tomographic problems have also appeared before. An iterative algorithm to compute a maximum likelihood point in which the prior information is represented by GP is introduced in [@tarantola2005inverse]. In [@hendriks2018implementation; @jidling2018probabilistic], tomographic reconstruction using GPs to model the strain field from neutron Bragg-edge measurements has been studied. Tomographic inversion using GP for plasma fusion and soft x-ray tomography have been done in [@li2013bayesian; @svensson2011non]. Nevertheless, the proposed approach is different from the existing work.

Our aim is to employ a hierarchical Gaussian process regression model to reconstruct the x-ray tomographic image from limited projection data. Due to the measurement model involving line integral computations, the direct GP approach does not allow for closed form expressions. The first contribution of this article is to overcome this issue by employing the basis function expansion method proposed in [@SolinSarkka2015], which makes the line integral computations tractable as it detaches the integrals from the model parameters. This approach can be directly used for common GP regression covariance functions such as Matérn or squared exponential. The second contribution of this article is to point out that the we can also reformulate classical regularization, in particular Laplacian and Tikhonov regularization, as Gaussian process regression where only the spectral density of the process (although not the covariance function itself) is well defined. As the basis function expansion only requires the availability of the spectral density, we can build a hierarchical model off a classical regularization model as well and have a principles means to tune the regularization parameters. Finally, the third contribution is to present methods for hyperparameter estimation that arise from the machine learning literature and apply the methodology to the tomographic reconstruction problem. In particular, the proposed methods are applied to simulated 2D chest phantom data available in <span style="font-variant:small-caps;">Matlab</span> and real carved cheese data measured with $\mu$CT system. The results show that the reconstruction images created using the proposed GP method outperforms the FBP reconstructions in terms of image quality measured as relative error and as peak signal to noise ratio.

Constructing the model
======================

The tomographic measurement data
--------------------------------

Consider a physical domain $\Omega \subset {{\mathbb R}}^2$ and an attenuation function $f:\Omega\rightarrow{{\mathbb R}}$. The x-rays travel through $\Omega$ along straight lines and we assume that the initial intensity (photons) of the x-ray is $I_0$ and the exiting x-ray intensity is $I_d$. If we denote a ray through the object as function $s \mapsto (x_1(s),x_2(s))$ Then the formula for the intensity loss of the x-ray within a small distance $ds$ is given as:

$$\label{calibration1}
\frac{dI(s)}{I(s)}= -f(x_1(s),x_2(s)) ds,$$

and by integrating both sides of , the following relationship is obtained $$\label{calibration2}
\int_{-R}^{R} f(x_1(s),x_2(s)) ds = \log\frac{I_0}{I_d},$$ where $R$ is the radius of the object or area being examined.

In x-ray tomographic imaging, the aim is to reconstruct $f$ using measurement data collected from the intensities $I_d$ of x-rays for all lines through the object taken from different angles of view. The problem can be expressed using the Radon transform, which can be expressed as $$\label{Measurement Model}
\mathcal{R} f(r,\theta) = \int f(x_1,x_2) d{\mathbf x}_L,
$$ where $d{\mathbf x}_L$ denotes the $1$-dimensional Lebesgue measure along the line defined by $L=\{(x_1,x_2) \in {{\mathbb R}}^2 : x_1\cos \theta + x_2\sin\theta = r \}$, where $\theta\in[0,\pi)$ is the angle and $r\in{{\mathbb R}}$ is the distance of $L$ from the origin as shown in Figure \[Radon transform\].

(100,80) (0,-65)[![An illustration of the Radon transform. It maps the object $f$ on the $(x_1,x_2)$-domain into $f$ on the $(r,\theta)$ domain. The measurement data is collected from the intensities $I_d$ of x-rays for all lines L through the object $f(x_1,x_2)$ and from different angles of view. []{data-label="Radon transform"}](figure1 "fig:"){width="10cm"}]{}

The parametrization of the straight line $L$ with respect to the arc length $s$ can be written as: $$\begin{split}
  x_1(s,\theta,r) &= r \, \cos(\theta) - s \, \sin(\theta), \\
  x_2(s,\theta,r) &= r \, \sin(\theta) + s \, \cos(\theta). \\
\end{split}
\label{eq:xystr}$$ In this work, the object is placed inside a circular disk with radius $R$. Then, as a function of $r$ and $\theta$ the line integral in (\[Measurement Model\]) can be written as $$\label{Radon}
\begin{split}
     \mathcal{R} f(r,\theta) &= \int_{-R}^{R} 
    f(x_1(s,\theta,r),x_2(s,\theta,r)) \, ds \\
   &= \int_{-R}^R f({\mathbf x}^0+s\hat{{\mathbf u}}) ds,
\end{split}$$ where $$\begin{aligned}
    {\mathbf x}^0 = \begin{bmatrix} r\cos(\theta) & r\sin(\theta) \end{bmatrix}^{\mathsf{T}}, \qquad
    \hat{{\mathbf u}} = \begin{bmatrix} -\sin(\theta) & \cos(\theta) \end{bmatrix}^{\mathsf{T}}.\end{aligned}$$

In a real x-ray tomography application, the measurement is corrupted by at least two noise types: photons statistics and electronic noise. In x-ray imaging, a massive number of photons are usually recorded at each detector pixel. In such case, a Gaussian approximation for the attenuation data in $\eqref{calibration2}$ can be used [@bouman1992generalized; @sachs19993d]. Recall that a logarithm of the intensity is involved in $\eqref{Radon}$, and so additive noise is a reasonable model for the electronic noise.

We collect a set of measurements as $$\label{noisy measurement}
    y_i =  \int_{-R}^R f({\mathbf x}_i^0+s\hat{{\mathbf u}}_i) ds + \varepsilon_i,$$ where $i$ corresponds to the data point index. The corresponding inverse problem is given the noisy measurement data $\{y_i\}_{i=1}^n$ in to reconstruct the object $f$.

Gaussian processes as functional priors {#functional priors}
---------------------------------------

A Gaussian process (GP) [@Rasmussen2006] can be viewed as a distribution over functions, where the function value in each point is treated as a Gaussian random variable. To denote that the function $f$ is modeled as a GP, we formally write $$\begin{aligned}
\label{eq:fGP}
    f({\mathbf x}) \sim {\GP \left(m({\mathbf x}),\,\, k({\mathbf x},{\mathbf x}') \right)}.\end{aligned}$$ The GP is uniquely specified by the *mean function* $m({\mathbf x})=\mathbb{E}[f({\mathbf x})]$ and the *covariance function* $k({\mathbf x},{\mathbf x}')=\mathbb{E}[(f({\mathbf x})-m({\mathbf x}))(f({\mathbf x}')-m({\mathbf x}'))]$. The mean function encodes our prior belief of the value of $f$ in any point. In lack of better knowledge it is common to pick $m({\mathbf x})=0$, a choice that we will stick to also in this paper.

The covariance function on the other hand describes the covariance between two different function values $f({\mathbf x})$ and $f({\mathbf x}')$. The choice of covariance function is the most important part in the GP model, as it stipulates the properties assigned to $f$. A few different options are discussed in Section \[sec:cov\_func\].

As data is collected our belief about $f$ is updated. The aim of regression is to predict the function value $f({\mathbf x}_*)$ at an unseen test point ${\mathbf x}_*$ by conditioning on the seen data. Consider direct function measurements on the form $$y_i=f({\mathbf x}_i)+\varepsilon_i,$$ where $\varepsilon_i$ is independent and identically distributed (iid) Gaussian noise with variance $\sigma^2$, that is, $\varepsilon_i\sim{\mathcal{N}}(0,\sigma^2)$. Let the measurements be stored in the vector ${\mathbf y}$. Then the mean value and the variance of the predictive distribution $p(f({\mathbf x}_*) \mid {\mathbf y})$ are given by [@Rasmussen2006]

\[eq:GPreg\] $$\begin{aligned}
    \mathbb{E}[f({\mathbf x}_*) \mid {\mathbf y}] &=
    {\mathbf k}_*^{\mathsf{T}}(K +\sigma^2I)^{-1}{\mathbf y}, \\
    \mathbb{V}[f({\mathbf x}_*) \mid {\mathbf y}] &=
    k({\mathbf x}_*,{\mathbf x}_*)-
    {\mathbf k}_*^{\mathsf{T}}(K +\sigma^2I)^{-1}{\mathbf k}_*.
    \end{aligned}$$

Here the vector ${\mathbf k}_*$ contains the covariances between $f({\mathbf x}_*)$ and each measurement while the matrix $K$ contains the covariance between all measurements, such that

$$\begin{aligned}
    ({\mathbf k}_*)_i &= k({\mathbf x}_i,{\mathbf x}_*), \\
    K_{ij} &= k({\mathbf x}_i,{\mathbf x}_j).\end{aligned}$$

An example of GP regression for a two-dimensional input is given in Figure \[fig:se\_ill\_post\]. The red stars indicate the measurements, while the shaded surface is the GP prediction. The blue line highlights a slice of the plot that is shown explicitly to the right, including the $95\%$ credibility region.

![Left: GP prediction (shaded surface) obtained from the measurements (red stars, also indicated by their deviation from the prediction). Right: slice plot of the blue line in the left figure, including the $95\%$ credibility region.[]{data-label="fig:se_ill_post"}](gp_ill.eps){width="\textwidth"}

The Gaussian process for x-ray tomography {#GP Xray}
-----------------------------------------

In this section, we show how to apply the functional priors presented in Section \[functional priors\] to x-ray tomography application. Since the x-ray measurements are line integrals of the unknown function $f({\mathbf x})$, they are linear functionals of the Gaussian process. Hence, we can define a linear functional $\mathcal{H}_{{\mathbf x},i}$ as follows: $$\begin{aligned}
    \mathcal{H}_{{\mathbf x},i} f({\mathbf x}) = \int_{-R}^R f({\mathbf x}^0_i+s\hat{{\mathbf u}}_i) ds.\end{aligned}$$ and thus the GP regression problem becomes

\[eq:invprob\] $$\begin{aligned}
        f({\mathbf x}) &\sim {\GP \left(m({\mathbf x}),\,\, k({\mathbf x},{\mathbf x}') \right)}, \\
        y_i &= \mathcal{H}_{{\mathbf x},i} f({\mathbf x})+\varepsilon_i.
    \end{aligned}$$

As discussed, for example, in [@Sarkka:2011; @SolinSarkka2015] the GP regression equations can be extended to this kind of models, which in this case leads to the following:

\[eq:GPreg\] $$\begin{aligned}
    \mathbb{E}[f({\mathbf x}_*)|{\mathbf y}] &=
    {\mathbf q}_*^{\mathsf{T}}(K +\sigma^2I)^{-1}{\mathbf y}, \\
    \mathbb{V}[f({\mathbf x}_*)|{\mathbf y}] &=
    k({\mathbf x}_*,{\mathbf x}_*)-
    {\mathbf q}_*^{\mathsf{T}}(Q +\sigma^2I)^{-1}{\mathbf q}_*,
    \end{aligned}$$

where $\mathbf{y} = \begin{bmatrix} y_1 & \cdots & y_n \end{bmatrix}^{\mathsf{T}}$ and

$$\begin{aligned}
    \label{eq:crosscov}
    ({\mathbf q}_*)_i &= \int_{-R}^R k({\mathbf x}_i^0+s\hat{{\mathbf u}}_i,{\mathbf x}_*) ds, \\
    \label{eq:Gram}
    Q_{ij} &= \int_{-R}^R\int_{-R}^R k({\mathbf x}_i^0+s\hat{{\mathbf u}}_i,{\mathbf x}_j^0+s'\hat{{\mathbf u}}_j) ds ds'.\end{aligned}$$

In general we can not expect closed form solutions to – and numerical computations are then required. However, even with efficient numerical methods, the process of selecting the hyperparameters is tedious since the hyperparameters are in general not decoupled from the integrand and the integrals need to be computed repeatedly in several iterations. In this paper, we avoid this by using the basis function expansion that will be described in Section \[sec:approx\].

Squared exponential and Matérn covariance functions {#sec:cov_func}
---------------------------------------------------

An important modeling parameter in Gaussian process regression is the covariance function $k({\mathbf x},{\mathbf x}')$ which can be selected in various ways. Because the basis function expansion described in Section \[sec:approx\] requires the covariance function to be *stationary*, we here limit our discussion to covariance functions of this form. *Stationarity* means that $k({\mathbf x},{\mathbf x}')=k({\boldsymbol{\mathrm{r}}})$ where ${\boldsymbol{\mathrm{r}}}={\mathbf x}-{\mathbf x}'$, so the covariance only depends on the distance between the input points. In that case we can also work with the spectral density, which is the Fourier transform of the stationary covariance function $$\label{eq:fourier}
S({\boldsymbol{\omega}}) =
\mathcal{F}[k]
= \int k({\boldsymbol{\mathrm{r}}})e^{-\text{i}{\boldsymbol{\omega}}^{\mathsf{T}}{\boldsymbol{\mathrm{r}}}}d{\boldsymbol{\mathrm{r}}},$$ where again ${\boldsymbol{\mathrm{r}}}={\mathbf x}-{\mathbf x}'$.

The perhaps most commonly used covariance function within the machine learning context [@Rasmussen2006] is the *squared exponential* (SE) covariance function $$\begin{aligned}
\label{eq:SE}
    k_{\textrm{SE}}({\boldsymbol{\mathrm{r}}}) &= \sigma_f^2\exp\left[ -\frac{1}{2l^2}\|{\boldsymbol{\mathrm{r}}}\|_2^2 \right],\end{aligned}$$ which has the following spectral density $$\begin{aligned}
\label{eq:SES}
    S_\textrm{SE}({\boldsymbol{\omega}})&=\sigma_f^2 (2\pi)^{d/2} l^d
    \exp\left[
    -\frac{l^2 \| {\boldsymbol{\omega}} \|_2^2}{2}
    \right],\end{aligned}$$ where $d$ is the dimensionality of ${\mathbf x}$ (in our case $d=2$). The SE covariance function is characterized by the magnitude parameter $\sigma_f$ and the *length scale* $l$. The squared exponential covariance function is popular due to its simplicity and ease of implementation. It corresponds to a process whose sample paths are infinitely many times differentiable and thus the functions modeled by it are very smooth.

Another common family of covariance functions is given by the Matérn class

$$\begin{aligned}
k_{\textrm{Matern}}({\boldsymbol{\mathrm{r}}})&=
\sigma_f^2
\frac{2^{1-\nu}}{\Gamma(\nu)}
\left(
\frac{\sqrt{2\nu}\|{\boldsymbol{\mathrm{r}}}\|_2}{l}
\right)^\nu
K_\nu\left(
\frac{\sqrt{2\nu}\|{\boldsymbol{\mathrm{r}}}\|_2}{l}
\right),
\\
S_{\textrm{Matern}}({\boldsymbol{\omega}})&=
\sigma_f^2\frac{2^d \pi^{d/2}\Gamma(\nu+d/2)(2\nu)^\nu}{\Gamma(\nu)l^{2\nu}}
\left(
\frac{2\nu}{l^2} + \| {\boldsymbol{\omega}} \|_2^2
\right)^{-(\nu+d/2)},\end{aligned}$$

where $K_\nu$ is a modified Bessel function [@Rasmussen2006]. The smoothness of the process is increased with the parameter $\nu$: in the limit $\nu\rightarrow\infty$ we recover the squared exponential covariance function.

Gaussian processes are also closely connected to classical spline smoothing [@kimeldorf1970] as well as other classical regularization methods [@kaipio2006statistical; @mueller2012linear] for inverse problems. Although the construction of the corresponding covariance function is hard (or impossible), it is still possible to construct the corresponding spectral density in many cases. With these spectral densities and the basis function method of Section \[sec:approx\], we can construct probabilistic versions of the classical regularization methods as discussed in the next section.

Covariance functions arising from classical regularization
----------------------------------------------------------

Let us recall that a classical way to seek for solutions to inverse problems is via optimization of a functional of the form $$\mathcal{J}[f] = \frac{1}{2 \sigma^2} \sum_i (y_i - \mathcal{H}_{{\mathbf x},i} f({\mathbf x}))^2
  + \frac{1}{2 \sigma_f^2} \int | \mathcal{L} f({\mathbf x}) |^2 \, d{\mathbf x},$$ where $\mathcal{L}$ is a linear operator. This is equivalent to a Gaussian process regression problem, where the covariance operator is formally chosen to be $\mathcal{K} = [\mathcal{L}^* \mathcal{L}]^{-1}$. In (classical) Tikhonov regularization we have $\mathcal{L} = \mathcal{I}$ (identity operator) which corresponds to penalizing the norm of the solution. Another option is to penalize the Laplacian which gives $\mathcal{L} = \nabla^2$.

Although the kernel of this covariance operator is ill-defined with the classical choices of $\mathcal{L}$ and thus it is not possible to form the corresponding covariance function, we can still compute the corresponding spectral density function by computing the Fourier transform of $\mathcal{L}^* \mathcal{L}$ and then inverting it to form the spectral density: $$S({\boldsymbol{\omega}}) = \frac{\sigma_f^2}{\mathcal{F}[\mathcal{L}^* \mathcal{L}]}.$$ In particular, the minimum norm or (classical) Tikhonov regularization can be recovered by using a white noise prior which is given by the constant spectral density $$\label{eq:STikhonov}
  S_{\textrm{Tikhonov}}({\boldsymbol{\omega}}) = \sigma_f^2,$$ where $\sigma_f$ is a scaling parameter. Another interesting case is the Laplacian operator based regularization which corresponds to $$\label{eq:SLaplacian}
  S_{\textrm{Laplacian}}({\boldsymbol{\omega}}) = \frac{\sigma_f^2}{\| {\boldsymbol{\omega}} \|^4_2}.$$ It is useful to note that the latter spectral density corresponds to a $l \to \infty$ limit of the Matérn covariance function with $\nu + d/2 = 2$ and the white noise to $l \to 0$ in either the SE or the Matérn covariance functions. The covariance functions corresponding to the above spectral densities would be degenerate, but this does not prevent us from using the spectral densities in the basis function expansion method described in Section \[sec:approx\] as the method only requires the availability of the spectral density.

Basis function expansion {#sec:approx}
------------------------

To overcome the computational hazard described in Section \[GP Xray\], we consider the approximation method proposed in [@SolinSarkka2015], which relies on the following truncated basis function expansion $$\begin{aligned}
    \label{eq:BFE}
    k({\mathbf x},{\mathbf x}')\approx \sum_{i=1}^{m}S({\sqrt{\lambda}_i})\phi_i({\mathbf x})\phi_i({\mathbf x}'),  \end{aligned}$$ where $S$ denotes the spectral density of the covariance function, and $m$ is the truncation number. The basis functions $\phi_i({\mathbf x})$ and eigenvalues ${\lambda}_i$ are obtained from the solution to the Laplace eigenvalue problem on the domain $\Omega$ $$\begin{aligned}
\label{eq:lapl_eig}
\begin{cases}
\hspace{-4mm}
\begin{split}
-\Delta\phi_i({\mathbf x})&  ={\lambda}_i\phi_i({\mathbf x}),  \\
\phi_i({\mathbf x})&         =0,         
\end{split}
\end{cases}
\quad
\begin{split}
{\mathbf x}&\in\Omega, \\
{\mathbf x}&\in\partial\Omega. 
\end{split}\end{aligned}$$ In two dimensions with $\Omega=[-L_1,L_1]\times[-L_2,L_2]$ we introduce the positive integers $i_1\le m_1$ and $i_2\le m_2$. The number of basis functions is then $m=m_1m_2$ and the solution to is given by

$$\begin{aligned}
\phi_i({\mathbf x})  &= \frac{1}{\sqrt{L_1L_2}}\sin\big(\varphi_{i_1}(x_1+L_1)\big)\sin\big(\varphi_{i_2}(x_2+L_2)\big) ,
\\
\lambda_i &= \varphi_{i_1}^2+\varphi_{i_2}^2, \quad \varphi_{i_1}=\frac{\pi i_1}{2L_1}, \quad \varphi_{i_2}=\frac{\pi i_2}{2L_2}, \end{aligned}$$

where $i=i_1+m_1(i_2-1)$. Let us now build the vector ${\boldsymbol\phi}_*\in{{\mathbb R}}^{m\times1}$, the matrix $\Phi\in{{\mathbb R}}^{m\times M}$ and the diagonal matrix $\Lambda\in{{\mathbb R}}^{m\times m}$ as

$$\begin{aligned}
    ({\boldsymbol\phi}_*)_i&=\phi_i({\mathbf x}_*), \\
    \label{eq:Phi_entry}
    \Phi_{ij} &= \int_{-R}^R \phi_i({\mathbf x}^0_j+s\hat{{\mathbf u}}_j)ds, \\
    \Lambda_{ii} &= S(\sqrt{\lambda_i}). \end{aligned}$$

The entries $\Phi_{ij}$ can be computed in closed form with details given in \[app:compdet\]. Now we substitute $Q\approx\Phi^{\mathsf{T}}\Lambda\Phi$ and ${\mathbf q}_*\approx\Phi^{\mathsf{T}}\Lambda{\boldsymbol\phi}_*$ to obtain

\[eq:pred\_appr2\] $$\begin{aligned}
\mathbb{E}[f({\mathbf x}_*) \mid {\mathbf y}]
&\approx{\boldsymbol\phi}_{*}^{\mathsf{T}}\Lambda \Phi (\Phi^{\mathsf{T}}\Lambda \Phi + \sigma^2 I)^{-1}{\mathbf y},
\\
\mathbb{V}[f({\mathbf x}_*) \mid {\mathbf y}]
&\approx {\boldsymbol\phi}_{*}^{\mathsf{T}}\Lambda {\boldsymbol\phi}_{*} -  {\boldsymbol\phi}_{*}^{\mathsf{T}}\Lambda \Phi (\Phi^{\mathsf{T}}\Lambda \Phi + \sigma^2I)^{-1} \Phi^{\mathsf{T}}\Lambda  {\boldsymbol\phi}_{*}.\end{aligned}$$

When using the spectral densities corresponding to the classical regularization methods in and , the mean equation reduces to the classical solution (on the given basis). However, also for the classical regularization methods we can compute the variance function which gives uncertainty estimate for the solution which in the classical formulation is not available. Furthermore, the hyperparameter estimation methods outlined in the next section provide principled means to estimate the parameters also in the classical regularization methods.

Hyperparameter estimation
=========================

In this section, we will consider some methods for estimating the *hyperparameters*. The free parameters of the covariance function, for example, the parameters $\sigma_f$ and $l$ in the squared exponential covariance function, are together with the noise parameter $\sigma$ referred to as the hyperparameters of the model. In this work, we employ a Bayesian approach to estimate the hyperparameters, and comparisons with standard parameter estimation methods such as L-curve and cross-validation methods are given as well.

Posterior distribution of hyperparameters
-----------------------------------------

The marginal likelihood function corresponding to the model is given as $$\label{likelihood}
 p(\mathbf{y} \mid \sigma_f,l,\sigma)
 = \mathcal{N}(\mathbf{y} \mid \mathbf{0}, Q(\sigma_f,l) + \sigma^2 \, I),$$ where $Q(\sigma_f,l)$ is defined by . The posterior distribution of parameters can now be written as follows:

$$\label{posterior}
    p(\sigma_f,l,\sigma \mid \mathbf{y}) \propto
    p(\mathbf{y} \mid \sigma_f,l,\sigma) 
    p(\sigma_f) p(l) p(\sigma),$$

where non-informative priors are used: $p(\sigma_f)\propto \frac{1}{\sigma_f}$, $p(l)\propto \frac{1}{l}$ and $p(\sigma)\propto \frac{1}{\sigma}$. The logarithm of  can be written as $$\label{logposterior}
    \log\,p(\sigma_f,l,\sigma \mid \mathbf{y}) =
    \text{const.} - \frac{1}{2} \, \log \det (Q + \sigma^2 I) 
    - \frac{1}{2} \mathbf{y}^{\mathsf{T}}(Q + \sigma^2 I)^{-1} \mathbf{y} 
    - \log \frac{1}{\sigma_f}
    - \log \frac{1}{l} - \log \frac{1}{\sigma}.$$ Given the posterior distribution we have a wide selection of methods from statistics to estimate the parameters. One approach is to compute the maximum a posteriori (MAP) estimate of the parameters by using, for example, gradient-based optimization methods [@Rasmussen2006]. However, using this kind of point estimate loses the uncertainty information of the hyperparameters and therefore in this article we use Markov chain Monte Carlo (MCMC) methods [@brooks2011handbook] which retain the information about the uncertainty in the final result.

Metropolis–Hastings sampling of hyperparameters
-----------------------------------------------

As discussed in the previous section, the statistical formulation of the inverse problem gives a posterior distribution of the hyperparameters ${\boldsymbol{\varphi}} = (\sigma_f,l,\sigma)$ as the solution rather than single estimates. The MCMC methods are capable of generating samples from the distribution. The Monte Carlo samples can then be used for computing the mean, the variance, or some other statistics of the posterior distribution [@Gelman_et_al:2013]. In this work, we employ the Metropolis–Hastings algorithm to sample from the posterior distribution.

The L-curve method
------------------

One of the classical methods to obtain information about the optimum value for $\sigma$ is the L-curve method [@hansen1992analysis], which operats by plotting the norm of the solution $\lVert f_{\sigma}({\mathbf x}) \rVert _2$ versus the residual norm $\lVert \mathcal{H}_{{\mathbf x},i} f_{\sigma}({\mathbf x}) - y_i \rVert _2$. The associated L-curve is defined as the continous curve consisting of all the points $(\lVert \mathcal{H}_{{\mathbf x},i} f_{\sigma}({\mathbf x}) - y_i \rVert _2,\lVert f_{\sigma}({\mathbf x}) \rVert _2)$ for $\sigma \in [0,\infty)$.

Cross-validation
----------------

As a comparison, we also consider to use methods of cross-validation (CV) for model selection. In $k$-fold CV, the data are partitioned into $k$ disjoint sets $\mathbf{y}_{j}$, and at each round $j$ of CV, the predictive likelihood of the set $\mathbf{y}_{j}$ is computed given the rest of the data $\mathbf{y}_{-j}$. These likelihoods are used to monitor the predictive performance of the model. This performance is used to estimate the generalization error, and it can be used to carry out model selection [@kohavi1995study; @Rasmussen2006; @vehtari2017practical].

The Bayesian CV estimate of the predictive fit with given parameters ${\boldsymbol{\varphi}}$ is $$\mbox{CV} = \sum_{j=1}^n
    \log p(\mathbf{y}_j \mid \mathbf{y}_{-j},{\boldsymbol{\varphi}}),$$ where $p(\mathbf{y}_j \mid \mathbf{y}_{-j},{\boldsymbol{\varphi}})$ is the predictive likelihood of the data $\mathbf{y}_j$ given the rest of the data. The best parameter values with respect to CV can be computed by enumerating the possible parameter values and selecting the one which gives the best fit in terms of CV.

Experimental results {#sec:expResults}
====================

In this section, we present numerical results using the GP model for limited x-ray tomography problems. All the computations were implemented in <span style="font-variant:small-caps;">Matlab</span> 9.4 (R2018a) and performed on an Intel Core i5 at 2.3 GHz and CPU 8GB 2133MHz LPDDR3 memory.

For both simulated data (see Section \[SimulatedData\]) and real data (see Section \[RealData\]) we use $m=10^4$ basis functions in . The measurements are obtained from the line integral of each x-ray over the attenuation coefficient of the measured objects. The measurements are taken for each direction (angle of view), and later they will be referred to as projections. The same number of rays in each direction is used. The computation of the hyperparameters is carried out using the Metropolis–Hastings algorithms with $5\,000$ samples, and the first $1\,000$ samples are thrown away ([*burn-in*]{} period). The reconstruction is computed by taking the conditional mean of the object estimate.

Simulated data: 2D Chest phantom {#SimulatedData}
--------------------------------

As for the simulated data, we use one slice of <span style="font-variant:small-caps;">Matlab</span>’s 3D Chest dataset [@matlab_chest] as a ground truth, $f_{\text{true}}$, which is shown in Figure \[fig:ChestPhantomRec\](a). The size of the phantom is $N \times N$, with $N = 128$. The black region indicates zero values and lighter regions indicate higher attenuation function values. The measurements (i.e. sinogram) of the chest phantom are computed using the [radon]{} command in <span style="font-variant:small-caps;">Matlab</span> and corrupted by additive white Gaussian noise with zero mean and $0.1$ variance ($\sigma_{\text{true}} = 0.32$).

Several reconstructions of the chest phantom using different covariance functions, namely squared exponential (SE), Matérn, Laplacian, and Tikhonov, are presented. For the SE, Matérn, and Laplacian covariance functions, the paramameters $\sigma_f$, $l$, and $\sigma$ are estimated using the proposed method. We use $\nu = 1$ for the Matérn covariance. As for the Tikhonov covariance, it is not characterized by the length scale $l$, and hence only $\sigma_f$ and $\sigma$ are estimated. All the estimated parameters are reported in Table \[GP parameters\]. Figure \[fig:Histogram parameters Chest\] presents the histograms of the 1-d marginal posterior distribution of each parameters using different covariance functions. The histograms show the distribution of the parameters samples in the Metropolis–Hastings samples. The results show that the $\sigma_f$ estimate for SE and Matérn covariances is $0.12$, while for Laplacian and Tikhonov, the estimates are $0.05$ and $0.64$. For Matérn, Laplacian, and Tikhonov covariance functions, the $\sigma$ estimates are concentrated around the same values $0.34-0.39$ with standard deviation (SD) between $0.02 - 0.03$. These noise estimates are well-estimated the ground-truth noise, $\sigma_{\text{true}} = 0.32$, with the absolute error is between $0.02 - 0.07$. The estimate of the SE kernel appears to overestimate the noise, $\sigma = 0.60$. It is reported that the length-scale parameter, $l$, for Laplacian and SE covariance functions are concentrated in the same values, while for Matérn yields higher estimate, $l = 10.14$.

Figure \[fig:ChestPhantomRec\](c)-(f) shows GP reconstructions of the 2D chest phantom using different covariance functions from 9 projections (uniformly spaced) out of 180$^\circ$ angle of view and $185$ number of rays for each projection. The computation times for all numerical tests are reported in Table \[Computation time\]. The Metropolis–Hastings reconstruction shows longer computational time due to the need for generation of a large number of samples from the posterior distribution. However, the benefit of this algorithm is that it is easy to implement and it is reliable for sampling from high dimensional distributions.

      Target       FBP    SE     Matérn   Laplacian   Tikhonov
  --------------- ----- ------- -------- ----------- ----------
   Chest phantom   0.5   11210    9676      9615        9615

  : Computation times of chest phantom (in seconds)[]{data-label="Computation time"}

The numerical test of the simulated data reconstructions is compared against figures of merit, namely:

-   the relative error (RE) $$\begin{aligned}
            \frac{\|f_{\text{true}} - f_{\text{rec}} \|_2}{\| f_{\text{true}}\|_2},
        \end{aligned}$$ where $f_{\text{rec}}$ is the image reconstruction, and

-   the peak-signal-to-noise ratio (PSNR) $$\begin{aligned}
      10\log_{10}\left(\frac{\mathrm{peakval}^2}{\mathrm{MSE}}\right),
        \end{aligned}$$ where $\mathrm{peakval}$ is the maximum possible value of the image and $\mathrm{MSE}$ is the mean square error between $f_{\text{true}}$ and $f_{\text{rec}}$,

as shown in Table \[Figures of merit\].

In practice, image quality in CT depends on other parameters as well, such as image contrast, spatial resolution, and image noise [@goldman2007principles]. These parameters can be evaluated when the CT device is equipped with CT numbers for various materials, high-resolution image is available, and statistical fluctuations of image noise which require several times of measurement to record random variations in detected x-ray intensity are acquired. However, in this work, the collected datasets are not supported by the aforementioned factors and they fall outside the scope of this paper. The results presented here are focusing on the implementation of a new algorithm to limited-data CT reconstruction and are reported as a preliminary study.

Reconstruction using a conventional method is computed as well with the built-in <span style="font-variant:small-caps;">Matlab</span> function [iradon]{}, which uses the FBP to invert the Radon transform. It reconstructs a two-dimensional slice of the sample from the corresponding projections. The angles for which the projections are available are given as an argument to the function. Linear interpolation is applied during the backprojection and a Ram–Lak or ramp filter is used. The FBP reconstruction of the chest phantom is shown in Figure \[fig:ChestPhantomRec\](b). For comparison, FBP reconstructions computed using some other filters are seen in Figure \[fig:FBPs\].

  ------------ ----------------- -------------- ---------------
   Covariance   $\sigma_f$ (SD)     $l$ (SD)     $\sigma$ (SD)
    function                                    
       SE         0.12 (0.04)     5.03 (0.03)     0.60 (0.02)
     Matérn       0.12 (0.07)     10.14 (0.08)    0.34 (0.03)
   Laplacian      0.05 (0.10)     4.49 (0.02)     0.39 (0.03)
    Tikhonov      0.64 (0.02)          -          0.35 (0.03)
  ------------ ----------------- -------------- ---------------

  : The GP parameter estimates for the chest phantom. The estimates are calculated using the conditional mean, and the standard deviation (SD) values are also reported in parentheses.[]{data-label="GP parameters"}

(100,600) (52,415)[![Histogram of the 1-d marginal distribution of the GP parameters. Left, middle and right columns are the marginal distribution for parameter $\sigma_f$, $l$ and $\sigma$ with corresponding covariance functions indicated in the vertical text in the left of the figure. The estimate of the parameter $l$ is not available for Tikhonov covariance.[]{data-label="fig:Histogram parameters Chest"}](samples_1_ChestSE "fig:"){width="4.25cm"}]{} (190,415)[![Histogram of the 1-d marginal distribution of the GP parameters. Left, middle and right columns are the marginal distribution for parameter $\sigma_f$, $l$ and $\sigma$ with corresponding covariance functions indicated in the vertical text in the left of the figure. The estimate of the parameter $l$ is not available for Tikhonov covariance.[]{data-label="fig:Histogram parameters Chest"}](samples_2_ChestSE "fig:"){width="4.25cm"}]{} (326,415)[![Histogram of the 1-d marginal distribution of the GP parameters. Left, middle and right columns are the marginal distribution for parameter $\sigma_f$, $l$ and $\sigma$ with corresponding covariance functions indicated in the vertical text in the left of the figure. The estimate of the parameter $l$ is not available for Tikhonov covariance.[]{data-label="fig:Histogram parameters Chest"}](samples_3_ChestSE "fig:"){width="4.25cm"}]{} (52,265)[![Histogram of the 1-d marginal distribution of the GP parameters. Left, middle and right columns are the marginal distribution for parameter $\sigma_f$, $l$ and $\sigma$ with corresponding covariance functions indicated in the vertical text in the left of the figure. The estimate of the parameter $l$ is not available for Tikhonov covariance.[]{data-label="fig:Histogram parameters Chest"}](samples_1_ChestMatern "fig:"){width="4.25cm"}]{} (190,265)[![Histogram of the 1-d marginal distribution of the GP parameters. Left, middle and right columns are the marginal distribution for parameter $\sigma_f$, $l$ and $\sigma$ with corresponding covariance functions indicated in the vertical text in the left of the figure. The estimate of the parameter $l$ is not available for Tikhonov covariance.[]{data-label="fig:Histogram parameters Chest"}](samples_2_ChestMatern "fig:"){width="4.25cm"}]{} (326,265)[![Histogram of the 1-d marginal distribution of the GP parameters. Left, middle and right columns are the marginal distribution for parameter $\sigma_f$, $l$ and $\sigma$ with corresponding covariance functions indicated in the vertical text in the left of the figure. The estimate of the parameter $l$ is not available for Tikhonov covariance.[]{data-label="fig:Histogram parameters Chest"}](samples_3_ChestMatern "fig:"){width="4.25cm"}]{} (52,125)[![Histogram of the 1-d marginal distribution of the GP parameters. Left, middle and right columns are the marginal distribution for parameter $\sigma_f$, $l$ and $\sigma$ with corresponding covariance functions indicated in the vertical text in the left of the figure. The estimate of the parameter $l$ is not available for Tikhonov covariance.[]{data-label="fig:Histogram parameters Chest"}](samples_1_ChestLaplacian "fig:"){width="4.25cm"}]{} (185,125)[![Histogram of the 1-d marginal distribution of the GP parameters. Left, middle and right columns are the marginal distribution for parameter $\sigma_f$, $l$ and $\sigma$ with corresponding covariance functions indicated in the vertical text in the left of the figure. The estimate of the parameter $l$ is not available for Tikhonov covariance.[]{data-label="fig:Histogram parameters Chest"}](samples_2_ChestLaplacian "fig:"){width="4.25cm"}]{} (326,125)[![Histogram of the 1-d marginal distribution of the GP parameters. Left, middle and right columns are the marginal distribution for parameter $\sigma_f$, $l$ and $\sigma$ with corresponding covariance functions indicated in the vertical text in the left of the figure. The estimate of the parameter $l$ is not available for Tikhonov covariance.[]{data-label="fig:Histogram parameters Chest"}](samples_3_ChestLaplacian "fig:"){width="4.25cm"}]{} (52,-10)[![Histogram of the 1-d marginal distribution of the GP parameters. Left, middle and right columns are the marginal distribution for parameter $\sigma_f$, $l$ and $\sigma$ with corresponding covariance functions indicated in the vertical text in the left of the figure. The estimate of the parameter $l$ is not available for Tikhonov covariance.[]{data-label="fig:Histogram parameters Chest"}](samples_1_ChestTikhonov "fig:"){width="4.25cm"}]{} (326,-10)[![Histogram of the 1-d marginal distribution of the GP parameters. Left, middle and right columns are the marginal distribution for parameter $\sigma_f$, $l$ and $\sigma$ with corresponding covariance functions indicated in the vertical text in the left of the figure. The estimate of the parameter $l$ is not available for Tikhonov covariance.[]{data-label="fig:Histogram parameters Chest"}](samples_3_ChestTikhonov "fig:"){width="4.25cm"}]{} (30,480)

[90]{} SE

(30,330)

[90]{} Matérn

(30,180)

[90]{} Laplacian

(30,30)

[90]{} Tikhonov

(125,560)[$\sigma_f$]{} (250,560)[$l$]{} (393,560)[$\sigma$]{}

(100,260) (5,138)[![(a) A ground truth of 2D chest phantom. (b) Filtered backprojection reconstruction (Ram–Lak filter) from 9 projections. (c) GP reconstruction using SE covariance, (d) GP reconstruction using Matérn covariance, (e) GP reconstruction using Laplacian covariance, (f) GP reconstruction using Tikhonov covariance. The GP reconstructions are using 9 projections.[]{data-label="fig:ChestPhantomRec"}](ChestImage "fig:"){width="6.3cm"}]{} (158,138)[![(a) A ground truth of 2D chest phantom. (b) Filtered backprojection reconstruction (Ram–Lak filter) from 9 projections. (c) GP reconstruction using SE covariance, (d) GP reconstruction using Matérn covariance, (e) GP reconstruction using Laplacian covariance, (f) GP reconstruction using Tikhonov covariance. The GP reconstructions are using 9 projections.[]{data-label="fig:ChestPhantomRec"}](SliceChest9AngFBP "fig:"){width="6.3cm"}]{} (306,138)[![(a) A ground truth of 2D chest phantom. (b) Filtered backprojection reconstruction (Ram–Lak filter) from 9 projections. (c) GP reconstruction using SE covariance, (d) GP reconstruction using Matérn covariance, (e) GP reconstruction using Laplacian covariance, (f) GP reconstruction using Tikhonov covariance. The GP reconstructions are using 9 projections.[]{data-label="fig:ChestPhantomRec"}](SliceChest9AngSE "fig:"){width="6.3cm"}]{} (5,3)[![(a) A ground truth of 2D chest phantom. (b) Filtered backprojection reconstruction (Ram–Lak filter) from 9 projections. (c) GP reconstruction using SE covariance, (d) GP reconstruction using Matérn covariance, (e) GP reconstruction using Laplacian covariance, (f) GP reconstruction using Tikhonov covariance. The GP reconstructions are using 9 projections.[]{data-label="fig:ChestPhantomRec"}](SliceChest9AngMatern "fig:"){width="6.3cm"}]{} (158,3)[![(a) A ground truth of 2D chest phantom. (b) Filtered backprojection reconstruction (Ram–Lak filter) from 9 projections. (c) GP reconstruction using SE covariance, (d) GP reconstruction using Matérn covariance, (e) GP reconstruction using Laplacian covariance, (f) GP reconstruction using Tikhonov covariance. The GP reconstructions are using 9 projections.[]{data-label="fig:ChestPhantomRec"}](SliceChest9AngLaplacian "fig:"){width="6.3cm"}]{} (306,3)[![(a) A ground truth of 2D chest phantom. (b) Filtered backprojection reconstruction (Ram–Lak filter) from 9 projections. (c) GP reconstruction using SE covariance, (d) GP reconstruction using Matérn covariance, (e) GP reconstruction using Laplacian covariance, (f) GP reconstruction using Tikhonov covariance. The GP reconstructions are using 9 projections.[]{data-label="fig:ChestPhantomRec"}](SliceChest9AngTikhonov "fig:"){width="6.3cm"}]{} (87,137)[(a)]{} (240,137)[(b)]{} (393,137)[(c)]{} (87,0)[(d)]{} (240,0)[(e)]{} (393,0)[(f)]{}

(270,260) (5,138)[![Filtered backprojection reconstructions using (a) Shepp–Logan filter, (b) Cosine filter (c) Hamming filter, (d) Hann filter. Values of relative error (RE) are between $23.6 - 25.2$ and PSNR values are between $18.1 - 19.9\%$. []{data-label="fig:FBPs"}](SheppLogan "fig:"){width="6.3cm"}]{} (158,138)[![Filtered backprojection reconstructions using (a) Shepp–Logan filter, (b) Cosine filter (c) Hamming filter, (d) Hann filter. Values of relative error (RE) are between $23.6 - 25.2$ and PSNR values are between $18.1 - 19.9\%$. []{data-label="fig:FBPs"}](Cosine "fig:"){width="6.3cm"}]{} (5,3)[![Filtered backprojection reconstructions using (a) Shepp–Logan filter, (b) Cosine filter (c) Hamming filter, (d) Hann filter. Values of relative error (RE) are between $23.6 - 25.2$ and PSNR values are between $18.1 - 19.9\%$. []{data-label="fig:FBPs"}](Hamming "fig:"){width="6.3cm"}]{} (158,3)[![Filtered backprojection reconstructions using (a) Shepp–Logan filter, (b) Cosine filter (c) Hamming filter, (d) Hann filter. Values of relative error (RE) are between $23.6 - 25.2$ and PSNR values are between $18.1 - 19.9\%$. []{data-label="fig:FBPs"}](Hann "fig:"){width="6.3cm"}]{} (87,137)[(a)]{} (240,137)[(b)]{} (87,0)[(c)]{} (240,0)[(d)]{}

  Method                  RE (%)   PSNR   
  ---------------------- -------- ------- --
  FBP (Ram–Lak filter)    25.86    18.44  
  GP-SE                   29.41    21.76  
  GP-Matérn               23.26    22.76  
  GP-Laplacian            29.18    21.79  
  GP-Tikhonov             23.39    22.73  
  Lcurve-Laplacian        23.38    22.62  
  Lcurve-Tikhonov         23.26    22.63  
  CV-Laplacian            25.18    22.31  
  CV-Tikhonov             23.47    22.75  

  : Figures of merit for chest phantom reconstructions.[]{data-label="Figures of merit"}

We also compared the results to the L-curve method and the CV:

-   The L-curve method is applied to the Laplacian and the Tikhonov covariances and the L-curve plots from different values of parameter $10^{-1}\leq\sigma\leq10$ for both covariances are shown in Figure \[L-curve chest data\]. Both plots show that the corner of the L-curve is located in between $0.2\leq \sigma \leq 1$.

    (100,250) (-30,50)[![The L-curve for (a) Tikhonov and (b) Laplacian covariance from the chest phantom reconstruction.[]{data-label="L-curve chest data"}](LcurveChestTikhonov "fig:"){width="9cm"}]{} (210,50)[![The L-curve for (a) Tikhonov and (b) Laplacian covariance from the chest phantom reconstruction.[]{data-label="L-curve chest data"}](LcurveChestLaplacian "fig:"){width="9cm"}]{} (65,45)[$\lVert \mathcal{H}_{{\mathbf x},i} f_{\sigma}({\mathbf x}) - y_i \rVert _2$]{} (-13,130)

    [90]{} $\lVert f_{\sigma}({\mathbf x}) \rVert _2$

    (225,130)

    [90]{} $\lVert f_{\sigma}({\mathbf x}) \rVert _2$

    (305,45)[$\lVert \mathcal{H}_{{\mathbf x},i} f_{\sigma}({\mathbf x}) - y_i \rVert _2$]{} (100,15)[(a)]{} (345,15)[(b)]{} (10,220) (10,150) (50,127) (175,80) (267,225) (269,94) (290,95) (410,95)

-   The CV is tested for the Laplacian and Tikhonov covariances using point-wise evaluation of $10^{-2} \leq \sigma \leq 1$ and $10^{-2}\leq \sigma_f \leq 1$. For the Laplacian covariance, several points of length scale $1\leq \ell \leq 100$ are tested as well. The minimum prediction error was obtained for $\sigma_f = 0.8$, $\sigma = 0.8$ and $\ell = 10$. For the Tikhonov covariance, the minimum prediction error was obtained for $\sigma = 0.5$ and $\sigma_f = 0.5$. The estimates of $\sigma_f$ and $\sigma$ for Laplacian are $0.8$ and $0.5$, respectively, and they give the same estimates for the Tikhonov covariance function. The estimates of $\sigma$ for both kernels appear to overestimate the $\sigma_{\text{true}}$. The absolute error is in between $0.18 - 0.48$. The length-scale estimate from Laplacian covariance, $l = 10$, appears to close to the estimate in Matérn covariance.

Image reconstructions for both L-curve and CV methods are shown in Figure \[fig:Parameter Choice Methods\].

(100,270) (7,138)[![(a) A ground truth of 2D chest phantom. (b) & (c) reconstructions using L-curve parameter choice method with Laplacian (using $\sigma = 1$) and Tikhonov (using $\sigma = 0.2$) covariance functions, respectively. (d) & (e) reconstructions using CV with Laplacian and Tikhonov covariance functions, respectively []{data-label="fig:Parameter Choice Methods"}](ChestImage "fig:"){width="6.3cm"}]{} (155,138)[![(a) A ground truth of 2D chest phantom. (b) & (c) reconstructions using L-curve parameter choice method with Laplacian (using $\sigma = 1$) and Tikhonov (using $\sigma = 0.2$) covariance functions, respectively. (d) & (e) reconstructions using CV with Laplacian and Tikhonov covariance functions, respectively []{data-label="fig:Parameter Choice Methods"}](LcurveChestLaplacian_sigman1 "fig:"){width="6.3cm"}]{} (306,138)[![(a) A ground truth of 2D chest phantom. (b) & (c) reconstructions using L-curve parameter choice method with Laplacian (using $\sigma = 1$) and Tikhonov (using $\sigma = 0.2$) covariance functions, respectively. (d) & (e) reconstructions using CV with Laplacian and Tikhonov covariance functions, respectively []{data-label="fig:Parameter Choice Methods"}](LcurveChestTikhonov_sigman02 "fig:"){width="6.3cm"}]{} (8,3)[![(a) A ground truth of 2D chest phantom. (b) & (c) reconstructions using L-curve parameter choice method with Laplacian (using $\sigma = 1$) and Tikhonov (using $\sigma = 0.2$) covariance functions, respectively. (d) & (e) reconstructions using CV with Laplacian and Tikhonov covariance functions, respectively []{data-label="fig:Parameter Choice Methods"}](CVChestLaplacian "fig:"){width="6.3cm"}]{} (155,3)[![(a) A ground truth of 2D chest phantom. (b) & (c) reconstructions using L-curve parameter choice method with Laplacian (using $\sigma = 1$) and Tikhonov (using $\sigma = 0.2$) covariance functions, respectively. (d) & (e) reconstructions using CV with Laplacian and Tikhonov covariance functions, respectively []{data-label="fig:Parameter Choice Methods"}](CVChestTikhonov "fig:"){width="6.3cm"}]{} (87,137)[(a)]{} (240,137)[(b)]{} (393,137)[(c)]{} (87,0)[(d)]{} (240,0)[(e)]{}

Real data: Carved cheese {#RealData}
------------------------

We now consider a real-world example using the tomographic x-ray data of a carved cheese slice measured with a custom-built CT device available at the University of Helsinki, Finland. The dataset is available online [@fips_dataset]. For a detailed documentation of the acquisition setup—including the specifications of the x-ray systems—see [@bubba2017tomographic]. We use the downsampled sinogram with $140$ rays and 15 projections from $360^\circ$ angle of view. In the computations, the size of the target is set to $120 \times 120$.

Figure \[CheeseRec\](c) shows the GP reconstruction (Matérn covariance function) of the cross section of the carved cheese slice using 15 projections (uniformly spaced) out of 360$^\circ$ angle of view. For comparison, the FBP reconstruction is shown in Figure \[CheeseRec\](b).

  ------------ ----------------- -------------- ---------------
   Covariance   $\sigma_f$ (SD)     $l$ (SD)     $\sigma$ (SD)
    function                                    
     Matérn      0.012 (0.07)     11.00 (0.08)    0.02 (0.04)
  ------------ ----------------- -------------- ---------------

  : Estimated GP parameters for the carved cheese using Matérn covariance function. The estimates are calculated using the conditional mean, and the standard deviation (SD) values are also reported in parentheses.[]{data-label="GP parameter cheese"}

The computation times for the carved cheese are reported in Table \[Computation time cheese\].

(100,200) (150,38)[![(a) FBP reconstruction (Ram–Lak filter) of the carved cheese using dense $360$ projections. (b) Filtered backprojection reconstruction from 15 projections. (c) GP reconstruction using Matérn covariance from 15 projections.[]{data-label="CheeseRec"}](Cheese15AngFBP "fig:"){width="6.21cm"}]{} (0,-31)[![(a) FBP reconstruction (Ram–Lak filter) of the carved cheese using dense $360$ projections. (b) Filtered backprojection reconstruction from 15 projections. (c) GP reconstruction using Matérn covariance from 15 projections.[]{data-label="CheeseRec"}](Cheese_slice_1120b.pdf "fig:"){width="6.8cm"}]{} (295,38)[![(a) FBP reconstruction (Ram–Lak filter) of the carved cheese using dense $360$ projections. (b) Filtered backprojection reconstruction from 15 projections. (c) GP reconstruction using Matérn covariance from 15 projections.[]{data-label="CheeseRec"}](Cheese15AngMatern "fig:"){width="6.21cm"}]{} (90,30)[(a)]{} (232,30)[(b)]{} (395,30)[(c)]{}

      Target       FBP   Matérn        
  --------------- ----- -------- -- -- --
   Carved cheese   0.1   12604         

  : Computation times of the carved cheese (in seconds)[]{data-label="Computation time cheese"}

Discussion
----------

We have presented x-ray tomography reconstructions from both simulated and real data for limited projections (i.e. sparse sampling) using an approach based on the Gaussian process. However, other limited-data problems such as limited angle tomography could be explored as well. The quality of GP reconstructions using different covariance functions looks rather the same qualitatively. However, quantitatively, the reconstruction using Matérn covariance is the best one: it has the lowest RE $23.26\%$ and the highest PSNR $22.76$. PSNR describes the similarity of the original target with the reconstructioned image (the higher value, the better of the reconstruction). Figures of merit estimates are not available for the real cheese data since there is no comparable ground truth. Nevertheless, the quality of the reconstruction can be observed qualitatively by comparing with the FBP reconstruction obtained with dense $360$ projections from $360$ degrees shown in Figure \[CheeseRec\](a). The corresponding parameter estimates for the chest phantom and the cheese are reported in Table \[GP parameters\] and \[GP parameter cheese\]. For the chest phantom case, the estimate of parameter $\sigma$ using Matérn, Laplacian and Tikhonov kernels tend to be close to the true value $\sigma_{true}$. As for the SE covariance, the standard deviation of noise is overestimated.

The reconstructions produced by the FBP benchmark algorithm using sparse projections are overwhelmed by streak artefacts due to the nature of backprojection reconstruction, as shown in Figure \[fig:ChestPhantomRec\](b) for the chest phantom and Figure \[CheeseRec\](b) the for cheese target. The edges of the target are badly reconstructed. Due to the artefacts, especially for the chest phantom, it is difficult to distinguish the lighter region (which is assumed to be tissue) and the black region (air). The FBP reconstruction has the worst quality and it is confirmed in Table \[Figures of merit\] that it has a high RE value ($25.86\%$) and the lowest PSNR ($18.44$). FBP reconstructions computed with different filters are shown in Figure \[fig:FBPs\]. However, there is no significant improvement in the images as it is clarified by the RE and PSNR values in the caption as well as by qualitative investigation. On the other hand, the GP reconstructions outperform the FBP algorithm in terms of image quality as reported in the figures of merit. The PSNR values of the GP-based reconstructions are higher than that of the FBP reconstruction. Nevertheless, in GP reconstructions, sharp boundaries are difficult to achieve due to the smoothness assumptions embedded in the model.

The GP prior clearly suppresses the artifacts in the reconstructions as shown in Figure \[fig:ChestPhantomRec\](c) and \[CheeseRec\](c). In Figure \[fig:ChestPhantomRec\](c), the air and tissue region are recovered much better, since the prominent artefacts are much less. In Figure \[CheeseRec\](c), the air region (outside the cheese and the C and T letters) are much sharper than in the FBP reconstruction. Overall, the results indicate that the image quality can be improved significantly by employing the GP method.

In Figure \[fig:Parameter Choice Methods\] the image reconstructions using L-curve and CV methods are presented. The quality of the reconstructions is reported in Table \[Figures of merit\] as well. In these methods, finer point-wise evaluations might help to improve the quality of the reconstructions.

We emphasize that in the proposed GP-approach, some parameters in the prior is a part of the inference problem (see Equation ). Henceforth, we can avoid the difficulty in choosing the prior parameters. This problem corresponds to the classical regularization methods, in which selecting the regularization parameters is a very crucial step to produce a good reconstruction.

Conclusions
===========

We have employed the Gaussian process with a hierarchical prior to computed tomography using limited projection data. The method was implemented to estimate the x-ray attenuation function from the measured data produced by the Radon transform. The performance has been tested on simulated and real data, with promising results shown. Unlike algorithms commonly used for the limited x-ray tomography problem that require manual tuning of prior parameters, the proposed GP method offers an easier set up as it takes into account the prior parameters as a part of the estimation. Henceforth, it constitutes a promising and user-friendly strategy.

The most important part of the GP model is the selection of the covariance function, since it stipulates the properties of the unknown function. As such, it also leaves most room for improvement. Considering the examples in Section \[sec:expResults\], a common feature of the target functions is that they consist of a number of well-defined, separate regions. The function values are similar and thus highly correlated within the regions, while the correlation is low at the edges where rapid changes occur. This kind of behavior is hard to capture with a stationary covariance function that models the correlation as completely dependent on the distance between the input locations. A non-stationary alternative is provided by, for example, the neural network covariance function, which is known for its ability to model functions with non-smooth features [@Rasmussen2006]. The basis function approximation method employed in this work is only applicable to stationary covariance functions, but other approximations can of course be considered.

Despite its success, the computational burden of the proposed algorithm is relatively high. To solve this problem, speed-up strategies are available, such as implementing parallelized GPU code, optimizing the covariances of the sampling strategy, or by changing the MCMC algorithm to another one. Investigating finer resolution images and statistical records would also be interesting future research to evaluate other image quality parameters. Moreover, the proposed method can be applied to multidetector CT imaging [@mookiah2018multidetector; @flohr2005multi] as well as 3D CT problems using sparse data [@sidky2008image; @purisha2018automatic].

Details on the computation of $\Phi$ {#app:compdet}
====================================

Here we derive the closed-form expression of the entries $\Phi_{ij}$ stated in . We get that $$\begin{split}
\Phi_{ij} &= \int_{-R}^R \phi_i({\mathbf x}^0_j+s\hat{{\mathbf u}}_j)ds \\
&=\frac{1}{\sqrt{L_1L_2}}\int_{-R}^{R} \sin(\varphi_{i_1} r_j\cos\theta_j - \varphi_{i_1} s \sin\theta_j +\varphi_{i_1}L_1)\sin(\varphi_{i_2} r_j\sin\theta_j + \varphi_{i_2} s \cos\theta_j+\varphi_{i_2} L_2) ds 
\\
& = \frac{1}{\sqrt{L_1L_2}}\int_{-R}^{R} \sin(\alpha_{ij} s + \beta_{ij})\sin(\gamma_{ij} s + \delta_{ij}) ds \\ 
& = \frac{1}{2\sqrt{L_1L_1}} \int_{-R}^R \cos((\alpha_{ij}-\gamma_{ij})s+\beta_{ij}-\delta_{ij}) - \cos((\alpha_{ij}+\gamma_{ij})s + \beta_{ij} + \delta_{ij})ds  \\
& = \frac{1}{2\sqrt{L_1L_1}} \Big[ \frac{1}{\alpha_{ij}-\gamma_{ij}}\sin((\alpha_{ij}-\gamma_{ij})s+\beta_{ij}-\delta_{ij}) - \frac{1}{\alpha_{ij}+\gamma_{ij}}\sin((\alpha_{ij}+\gamma_{ij})s + \beta_{ij} + \delta_{ij})\Big]_{-R}^R 
\\
& = \frac{1}{2\sqrt{L_1L_1}} \Big( \frac{1}{\alpha_{ij}-\gamma_{ij}}\sin((\alpha_{ij}-\gamma_{ij})R+\beta_{ij}-\delta_{ij}) - \frac{1}{\alpha_{ij}+\gamma_{ij}}\sin((\alpha_{ij}+\gamma_{ij})R + \beta_{ij} + \delta_{ij}) 
\\
& \qquad - \frac{1}{\alpha_{ij}-\gamma_{ij}}\sin(-(\alpha_{ij}-\gamma_{ij})R+\beta_{ij}-\delta_{ij}) + \frac{1}{\alpha_{ij}+\gamma_{ij}}\sin(-(\alpha+\gamma_{ij})R + \beta_{ij} + \delta_{ij})\Big),
\end{split}$$ where

$$\begin{aligned}
    \alpha_{ij} &= \varphi_{i_1}\sin\theta_j, \\
    \beta_{ij}  &= \varphi_{i_1}r_j\cos\theta_j+\varphi_{i_1}L_1, \\
    \gamma_{ij} &= \varphi_{i_2}\cos\theta_j, \\
    \delta_{ij} &= \varphi_{i_2}r_j\sin\theta_j+\varphi_{i_2}L_2.\end{aligned}$$

References {#references .unnumbered}
==========
