Directional sinogram inpainting for limited angle tomography

Many applications in materials science and medical imaging rely on the x-ray transform as a mathematical model for performing 3D volume reconstructions of a sample from 2D data. We shall refer to any modality using the x-ray transform forward model as x-ray tomography. This encompasses a huge range of applications including positron emission tomography (PET) in front-line medical imaging [1], transmission electron microscopy (TEM) in materials or biological research [2–4], and x-ray computed tomography (CT) which enjoys success across many fields [5, 6].

The limited angle problem is common in x-ray tomography, for instance in TEM [7] and mammography [8], and is caused by a particular limited data scenario. Algebraically, we search for approximate solutions to the inverse problem

$$\begin{align*} \newcommand{\e}{{\rm e}} \newcommand{\C}[1]{\mathcal{#1}} \displaystyle {{\rm Given~data~}} b{{\rm ~find~optimal~pair~}} (u,v) {{\rm ~such~that~}}Sv = b, \C Ru = v\nonumber \end{align*}$$

where $\newcommand{\C}[1]{\mathcal{#1}} \C R$ is the x-ray transform to be defined in (2.1), $S$ represents the limited angle sub-sampling pattern described in figure 1. Typically, limited angle problems can occur due to having a large sample or because equipment does not allow the sample to be fully rotated. Mathematically, microlocal analysis can be used to categorise the limited angle problem and characterise artefacts that occur. Viewed through the Fourier slice theorem, it becomes clear that the Fourier coefficients of $u$ are partitioned into those 'visible' in $b$ and those contained in a 'missing wedge' [9]. These coefficients are referred to respectively as the visible and invisible singularities of $u$. The limited angle problem then is both a denoising and inpainting inverse problem, on the visible and invisible singularities respectively. The artefacts caused by the missing wedge can be explicitly characterised [10, 11] and examples of such streak artefacts and blurred boundaries can be seen in figures 2 and 3.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Whilst the techniques developed here can apply to any limited angle tomography problem, we focus on the application of TEM for specific examples.

Traditional methods for x-ray tomography reconstruction find approximate solutions to $\newcommand{\C}[1]{\mathcal{#1}} S\C Ru = b$, constraining $\newcommand{\C}[1]{\mathcal{#1}} \C Ru=v$ and only using prior knowledge of $u$ or the sinogram, $v$. There are three main methods which fit into this category:

• Filtered back projection (FBP) is a linear inversion method with smoothing on the sinogram, $v$, to account for noise [12–14].
• The simultaneous iterative reconstruction technique (SIRT) can be thought of as a preconditioned gradient descent on the function $\newcommand{\C}[1]{\mathcal{#1}} \newcommand{\norm}[1]{\left\Vert #1 \right\Vert} \norm{S\C Ru-b}_2^2$ [3, 15–17]. Regularisation is then typically implemented by an early-stopping technique.
• Variational methods where prior knowledge is encoded in regularisation functionals have now been applied in this field for nearly a decade. In particular, the current state of the art in electron tomography is total variation (TV) regularisation [2, 18, 19] where $u$ is encouraged to have a piecewise constant intensity. This will be introduced formally in section 2.

FBP and SIRT are commonly used for their speed although variational methods like TV have quickly gained popularity as they enable prior physical knowledge to be explicitly incorporated in reconstructions. This added prior knowledge tends to stabilise the reconstruction process and figure 2 gives examples where TV can vastly outperform FBP and SIRT when either the noise level is large or the angular range small. However, figure 3 further shows the limitations of TV when high noise and limited angles are combined. The only difference between the Shepp–Logan phantom data shown in figures 2 and 3 is that the former is clean data, in the image of the forward operator, whilst the latter has Gaussian white noise added. We see that as soon as there is a combined denoising/inpainting reconstruction problem, the TV prior on $u$ becomes insufficient to recover the structure of the sample.

Recently, these traditional methods have received a revival through machine learning methods, see for instance [20, 21]. In both of these examples the main artefact reduction is a learned denoising step which only enforces prior knowledge on $u$.

The most common method that has been used to reconstruct pairs $(u,v)$ is to solve each inverse problem sequentially. Typically, we can express the pipeline for such methods as:

$$\begin{align*} \newcommand{\e}{{\rm e}} \displaystyle v&= ~{\rm optimal~ inpainted~ sinogram~ given}~ b \nonumber \\ u&=~ {\rm optimal~ reconstruction~ given}~ v. \nonumber \end{align*}$$

This has seen much success in heavy metal artefact reduction [22, 23] where a regularisation functional for the inpainting problem may be constructed from dictionary learning [24], fractional order TV [23], and Euler's Elastica [25]. Euler's Elastica has also been used in the limited angle problem [26] along with more customised interpolation methods [27]. These approaches allow us to use prior knowledge on the sinogram to calculate $v$ and then spatial prior knowledge to calculate $u$ from $v$; at no point is the choice of $v$ influenced by our spatial prior. A full joint approach allows us to go beyond this and use all of our prior knowledge to inform the choice of both $u$ and $v$. If our prior knowledge is consistent with the true data then this extra utilisation of our prior must have the potential to improve the recovery of both $u$ and $v$. To build a model for this framework we shall encode our In this paper, therefore, we propose a full joint approach which allows us to use all of our prior knowledge at once. To realise this idea we encode prior knowledge and consistency terms into a single energy functional such that an optimal pair of reconstructions will minimise this joint functional, which we shall write as:

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\C}[1]{\mathcal{#1}} \displaystyle \label{basic model} E(u,v) = \alpha_1d_1(\C Ru,v) + \alpha_2d_2(S\C Ru,b) + \alpha_3d_3(Sv,b) + \alpha_4r_1(u) +\alpha_5r_2(v) \nonumber \end{align} \tag{ 1.1 }$$

where $\alpha_i\geqslant0$ are weighting parameters, $d_i$ are appropriate distance functionals and $r_i$ are regularisation functionals which encode our prior knowledge. Note that choice of $d_2$ and $d_3$ are dictated by the data noise model. In what follows, $r_1$ is chosen to be the total variation.

Our choice for $r_2$, the sinogram regularisation, is based on theoretical and heuristic observations. Thirion [28] has shown that discontinuities in $u$ correspond to sharp edges in the sinogram. In figure 3 we also see that blurred reconstructions correspond to loss of structure in the sinogram. Therefore, $r_2$ will be chosen to detect sharp features in the given data and preserve these through the inpainting region. The exact form of $r_2$ will be formalised in section 3.

A typical advantage of joint models is that they generalise previous ones. For instance, if we let $\alpha_2,\alpha_4\to\infty$ then we recover the TV reconstruction model. Alternatively, if we let $\alpha_3,\alpha_5\to\infty$ then we recover a method which performs the inpainting and then the reconstruction sequentially, as in [23, 25, 26]. Recent work [29] has shown that such a joint approach can be advantageous in similar situations but closest to our approach is that of [30] where $r_1$ and $r_2$ were chosen to encode wavelet sparsity in both $u$ and $v$. We shall demonstrate, as seen in figure 4, that the flexibility of our joint model, (1.1), can allow for a better adaptive fitting to the data.

Zoom In Zoom Out Reset image size

The main contribution of this work is to provide a framework for building models of the form described in (1.1) and provide new proofs for a numerical scheme for minimising these functionals. This numerical scheme is valid for a very large class of non-smooth and non-convex functionals $r_i$ and thus could be used in many other applications.

Section 2 first outlines the necessary concepts and notation needed to formalise the statement of our specific joint model in section 3. It will become apparent that the main numerical requirements of this work will require minimising a functional which is neither convex or smooth. Section 4 will start by reviewing recent work from Ochs et al [31] and we then provide alternative concise and self-contained proofs. Our main contribution here is to extend the existing results to an alternating (block) descent scenario. Finally, we present numerical results including two synthetic phantoms and experimental electron microscopy data where the limited angle situation occurs naturally.

The principal forward model for x-ray tomography is provided by the x-ray transform which can be defined for any bounded domain, $\newcommand{\F}[1]{\mathbb{#1}} \Omega\subset\F R^n$, by

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\C}[1]{\mathcal{#1}} \newcommand{\F}[1]{\mathbb{#1}} \displaystyle \label{radon} \C R\colon L^1(\Omega,\F R)\to L^1(\C S^{n-1}\times \F R^n,\F R) ~{\rm such~ that }~ \C Ru(\theta,y) = \int_{x = y + t\theta, t\in \F R} u(x){\rm d}t \nonumber \end{align} \tag{ 2.1 }$$

where $\newcommand{\F}[1]{\mathbb{#1}} \newcommand{\C}[1]{\mathcal{#1}} \newcommand{\st}{{\rm \;s.t.\;}} \C S^{n-1}=\{s\in\F R^n\st |s| = 1\}$. In this work, for simplicity, we will only be using $n=2$ although the case $n=3$ is completely analogous.

Total variation (TV) regularisation is extremely common across many fields of image processing [32–34]. The definition of the (isotropic) TV semi-norm on domains $\newcommand{\F}[1]{\mathbb{#1}} \Omega\subset\F R^n$ is stated as

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\st}{{\rm \;s.t.\;}} \newcommand{\IP}[2]{\left\langle #1,#2 \right\rangle} \newcommand{\re}{{\rm Re}} \newcommand{\TV}{{\rm TV}} \newcommand{\op}[1]{{\rm #1}} \newcommand{\F}[1]{\mathbb{#1}} \displaystyle \label{TV} &g\in L^1(\Omega,\F R) \Longrightarrow \TV(g) = \sup\left\{\int \IP{g(x)}{\op{div}(\varphi)(x)}{\rm d}x\right.\nonumber \\ &\qquad \qquad \qquad \qquad \left.\vphantom{\int}\st \varphi\in C_c^\infty(\Omega,\F R^n), |\varphi(x)|_2\leqslant 1{{\rm ~for~all~}} x\right\}. \nonumber \end{align} \tag{ 2.2 }$$

The intuition behind this is that when $\newcommand{\F}[1]{\mathbb{#1}} g\in W^{1,1}(\Omega,\F R)$ then we have exactly

$$\begin{align*} \newcommand{\e}{{\rm e}} \newcommand{\TV}{{\rm TV}} \newcommand{\norm}[1]{\left\Vert #1 \right\Vert} \displaystyle \TV(g) = \norm{\nabla g}_{2,1} = \int |\nabla g(x)|_2{\rm d}x.\nonumber \end{align*}$$

From this point forward we shall write the $\newcommand{\norm}[1]{\left\Vert #1 \right\Vert} \norm{\cdot}_{2,1}$ representation where $|\nabla g|$ is to be understood as a measure if necessary. We shall denote the space of bounded variation as

$$\begin{align*} \newcommand{\e}{{\rm e}} \newcommand{\st}{{\rm \;s.t.\;}} \newcommand{\TV}{{\rm TV}} \newcommand{\op}[1]{{\rm #1}} \newcommand{\F}[1]{\mathbb{#1}} \displaystyle \op{BV}(\Omega, \F R) = \{g\in L^1\st \TV(g) < \infty\} .\nonumber \end{align*}$$

One property that we shall use about the space of bounded variation is $\newcommand{\op}[1]{{\rm #1}} \op{BV}\subset L^2\subset L^1$ which holds whenever $\Omega$ is compact by the Poincaré inequality: $\newcommand{\norm}[1]{\left\Vert #1 \right\Vert} \newcommand{\op}[1]{{\rm #1}} \newcommand{\sfrac}[2]{^{#1}\!\!/\!_{#2}} \norm{g-\sfrac{\int g}{\int 1}}_2 \lesssim \op{TV}(g)$.

The most common way to enforce this prior in reconstruction or inpainting problems is generalised Tikhonov regularisation, which gives us the basic reconstruction method [2, 18].

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\argmin}{{\rm argmin}} \newcommand{\TV}{{\rm TV}} \newcommand{\norm}[1]{\left\Vert #1 \right\Vert} \newcommand{\C}[1]{\mathcal{#1}} \displaystyle u = \mathop{\argmin}\limits_{u\geqslant 0} \frac12\norm{S\C Ru-b}_2^2 + \lambda\TV(u) ~{\rm for~ some }~ \lambda\geqslant0. \nonumber \end{align} \tag{ 2.3 }$$

The parameter $\lambda$ is a regularisation parameter, which allows to enforce more or less regularity, depending on the quality of the data $b$.

For our sinogram regularisation functional we shall use a directionally weighted TV penalty, motivated by the TV diffusion model developed by Joachim Weickert [35] for various imaging techniques including denoising, inpainting and compression. Such an approach has already shown great ability for enhancing edges in noisy or blurred images, and preserves line structures across inpainting regions [36–38]. The heuristic for our regularisation on the sinogram was described in figure 3 and we shall encode it in an anisotropic TV penalty which shall be written as

$$\begin{align*} \newcommand{\e}{{\rm e}} \newcommand{\DTV}{{\rm DTV}} \newcommand{\norm}[1]{\left\Vert #1 \right\Vert} \newcommand{\F}[1]{\mathbb{#1}} \displaystyle \DTV(v) = \int |A(x)\nabla v(x)| {\rm d}x = \norm{A\nabla v}_{2,1} ~{\rm for ~some~ anisotropic }~ A\colon\F R^2\to \F R^{2\times 2}.\nonumber \end{align*}$$

The power of such a weighted extension of TV is that once a line is detected, either known beforehand or detected adaptively, we can embed this in $A$ and enhance or sharpen that line structure in the final result. The general form that we choose for $A$ is

$$\begin{align} \newcommand{\e}{{\rm e}} \renewcommand{\vec }{\mathbf} \newcommand{\re}{{\rm Re}} \newcommand{\IP}[2]{\left\langle #1,#2 \right\rangle} \newcommand{\F}[1]{\mathbb{#1}} \displaystyle \label{matrix form} A(x) = &\;c_1(x)\vec {e}_1(x)\vec {e}_1(x)^T + c_2(x)\vec {e}_2(x)\vec {e}_2(x)^T \nonumber \\&{\rm such~ that }~ \vec {e}_i\colon \F R^2\to \F R^2, |\vec {e}_i(x)| = 1, \IP{\vec {e}_1(x)}{\vec {e}_2(x)} = 0 \nonumber \end{align} \tag{ 2.4 }$$

i.e.

$$\begin{align*} \newcommand{\e}{{\rm e}} \renewcommand{\vec }{\mathbf} \newcommand{\re}{{\rm Re}} \newcommand{\IP}[2]{\left\langle #1,#2 \right\rangle} \newcommand{\DTV}{{\rm DTV}} \displaystyle \DTV(v) = \int \sqrt{c_1^2|\IP{\vec {e}_1}{\nabla v}|^2 + c_2^2|\IP{\vec {e}_2}{\nabla v}|^2} {\rm d}x.\nonumber \end{align*}$$

Examples of this are presented in figure 5. Note that the choice $c_1=c_2$ recovers the traditional TV regulariser but for $|c_1|\ll c_2$ we allow for much larger (sparse) gradients in the direction of $\newcommand{\re}{{\rm Re}} \renewcommand{\vec }{\mathbf} \vec {e}_1$. This allows for large jumps in the direction of $\newcommand{\re}{{\rm Re}} \renewcommand{\vec }{\mathbf} \vec {e}_1$ whilst maintaining flatness in the direction of $\newcommand{\re}{{\rm Re}} \renewcommand{\vec }{\mathbf} \vec {e}_2$. In order to generate these parameters we follow the construction of Weickert [35]. Given a noisy image, $d$, we can construct the structure tensor:

$$\begin{align*} \newcommand{\e}{{\rm e}} \renewcommand{\vec }{\mathbf} \newcommand{\re}{{\rm Re}} \displaystyle (\nabla d_\rho \nabla d_\rho^T)_\sigma(x) = \lambda_1(x)\vec {e}_1(x)\vec {e}_1(x)^T + \lambda_2(x)\vec {e}_2(x)\vec {e}_2(x)^T {{\rm ~such~that~}} \lambda_1(x)\geqslant \lambda_2(x)\geqslant0\nonumber \end{align*}$$

where

$$\begin{align*} \newcommand{\e}{{\rm e}} \newcommand{\st}{{\rm \;s.t.\;}} \newcommand{\sfrac}[2]{^{#1}\!\!/\!_{#2}} \displaystyle d_\rho(x) = [d\star \exp(-\sfrac{|\cdot|^2}{2\rho^2})] (x)= \int \exp(-\sfrac{|y-x|^2}{2\rho^2})d(y) {{\rm ~etc.}}\nonumber \end{align*}$$

denotes convolution with the heat kernel. This eigenvalue decomposition is typically very informative in constructing $A$. If we define

$$\begin{align*} \newcommand{\e}{{\rm e}} \displaystyle \Delta(x) = \lambda_1(x)-\lambda_2(x) ~{\rm is~ coherence}~ \qquad \Sigma(x) = \lambda_1(x)+\lambda_2(x) ~{\rm is ~energy}\nonumber \end{align*}$$

then the eigenvectors give the alignment of edges and $\Delta, \Sigma$ characterise the local behaviour, as in figure 6. In particular, we simplify the model to

$$\begin{align} \newcommand{\e}{{\rm e}} \renewcommand{\vec }{\mathbf} \newcommand{\re}{{\rm Re}} \displaystyle \label{direction tensor} A_d(x) := c_1(x|\Delta,\Sigma)\vec {e}_1(x)\vec {e}_1(x)^T + c_2(x|\Delta,\Sigma)\vec {e}_2(x)\vec {e}_2(x)^T \nonumber \end{align} \tag{ 2.5 }$$

where the only parameters left to choose are $c_i$. Typical examples of include

$$\begin{align*} \newcommand{\e}{{\rm e}} \displaystyle c_1 = \frac{1}{\sqrt{1+\Sigma^2}},\qquad c_2 = 1\nonumber \end{align*}$$

$$\begin{align*} \newcommand{\e}{{\rm e}} \renewcommand{\epsilon}{\varepsilon} \newcommand{\re}{{\rm Re}} \newcommand{\sfrac}[2]{^{#1}\!\!/\!_{#2}} \displaystyle c_1 = \epsilon,\qquad c_2 = \epsilon+\exp\left(-\sfrac{1}{\Delta^2}\right) {{\rm ~for~some~}} \epsilon>0.\nonumber \end{align*}$$

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The key idea here is that $c_1\ll c_2$ near edges and $c_1=c_2$ on flat regions. In practice $d$ will also be an optimisation parameter and so we shall require a regularity result on our choice of $d\mapsto A_d$, now characterised uniquely by our choice of $c_i$.

Theorem 2.1. If

• (i)  
  $c_i$ are $2k$ times continuously differentiable in $\Delta$ and $\Sigma$, $k\geqslant1$
• (ii)  
  $c_1(x|0,\Sigma) = c_2(x|0,\Sigma)$ for all $x$ and $\Sigma\geqslant0$
• (iii)  
  $\partial_\Delta^{2j-1} c_1(x|0,\Sigma) = \partial_\Delta^{2j-1} c_2(x|0,\Sigma) = 0$ for all $x$ and $\Sigma\geqslant0, j=1\ldots,k$

Then $A_d$ is $C^{2k-1}$ with respect to $d$ for all $\rho>0,\sigma\geqslant0$.

Remark 2.2. 

• Property (ii) is necessary for $A_d$ to be well defined and continuous for all fixed $d$
• If we can write $c_i = c_i(\Delta^2,\Sigma)$ then property (iii) holds trivially

The proof of this theorem is contained in appendix A.

The proof of theorem 4.1 will be a consequence of two lemmas.

• In lemma 4.3 we show for $\tau_x,\tau_y$ (algorithm 1) sufficiently large, the sequence $E_{n,n}$ is monotonically decreasing and sequences $\newcommand{\norm}[1]{\left\Vert #1 \right\Vert} \norm{x_n-x_{n-1}}_X, \norm{y_n-y_{n-1}}_Y$ converge to 0. This follows by a relatively standard sufficient decrease argument as seen in [31, 42, 43].
• In lemma 4.6 we first define a critical point for functions which are non-convex and non-differentiable. If a subsequence of our iterates converges in norm then the limit must be a critical point in each of the two axes. Note that it is very difficult to expect more than this in such a general setting, for instance example 4.2 shows a uniformly convex energy which shows this to be sharp. The common technique for overcoming this is assuming differentiability in the terms including both $x$ and $y$ [42, 44, 45]. These previous results and algorithms are not available to us as we allow non-convex terms which are also non-differentiable.

Example 4.2. Define $E(x,y) = \max(x,y)+x^2+y^2$ for $\newcommand{\F}[1]{\mathbb{#1}} x,y\in\F R$. This is clearly jointly convex in $(x,y)$ and thus a simple case of functions considered in theorem 4.1. Suppose $(x_0,y_0) = (0,0)$ then

$$\begin{align*} \newcommand{\e}{{\rm e}} \newcommand{\argmin}{{\rm argmin}} \displaystyle x_1 &= \argmin~ E(x,y_0) + \tau_x(x-x_0)^2 = 0 \nonumber \\y_1 &= \argmin ~E(x_1,y) + \tau_y(y-y_0)^2 = 0. \nonumber \end{align*}$$

Hence $(0,0)$ is a limit point of the algorithm but it is not a critical point, $E$ is uniformly convex and so it has only one critical point, $\newcommand{\sfrac}[2]{^{#1}\!\!/\!_{#2}} (-\sfrac12,-\sfrac12)$.

In the following we prove the monotone decrease property of our energy functional between iterations.

Lemma 4.3. If for each $n$

$$\begin{align*} \newcommand{\e}{{\rm e}} \displaystyle \tau_x \geqslant L_x+\tau_X, \qquad \tau_y \geqslant L_y+\tau_Y\nonumber \end{align*}$$

for some $\tau_X,\tau_Y\geqslant0$ then

$$\begin{align*} \newcommand{\e}{{\rm e}} \newcommand{\norm}[1]{\left\Vert #1 \right\Vert} \displaystyle \sum^\infty\tau_X\norm{x_n-x_{n-1}}_{X}^2 + \tau_Y\norm{y_n-y_{n-1}}_{Y}^2 \leqslant E(x_0,y_0)\nonumber \end{align*}$$

and

$$\begin{align*} \newcommand{\e}{{\rm e}} \displaystyle E(x_{n+1},y_{n+1})\leqslant E(x_n,y_n) ~{\rm for~ all }~ n.\nonumber \end{align*}$$

Proof. Note by equations (4.6) and (4.7) (definition of our sequence) we have

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\norm}[1]{\left\Vert #1 \right\Vert} \displaystyle f_{n+1,n} + g(J_{n,n}+\nabla_xJ_{n,n}(x_{n+1}-x_n)) + \tau_x\norm{x_{n+1}-x_n}_{X}^2 \leqslant E_{n,n} \label{gc1} \nonumber \end{align} \tag{ 4.8 }$$

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\norm}[1]{\left\Vert #1 \right\Vert} \displaystyle f_{n+1,n+1} + g(J_{n+1,n}+\nabla_yJ_{n+1,n}(y_{n+1}-y_{n})) + \tau_y\norm{y_{n+1}-y_{n}}_{Y}^2 \leqslant E_{n+1,n} .\label{gc3} \nonumber \end{align} \tag{ 4.9 }$$

Further, by application of the triangle inequality for $g$ and the mean value theorem we have

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\op}[1]{{\rm #1}} \newcommand{\norm}[1]{\left\Vert #1 \right\Vert} \displaystyle g(J_{n+1,n}) - &g(J_{n,n}+\nabla_xJ_{n,n}(x_{n+1}-x_{n}))+ \tau_X\norm{x_{n+1}-x_{n}}_{X}^2\nonumber \\&\qquad \leqslant g(J_{n+1,n} -J_{n,n}-\nabla_xJ_{n,n}(x_{n+1}-x_{n}))+ \tau_X\norm{x_{n+1}-x_{n}}_{X}^2\nonumber \\&\qquad = g([\nabla_x J(\xi)-\nabla_xJ_{n,n}](x_{n+1}-x_{n}))+ \tau_X\norm{x_{n+1}-x_{n}}_{X}^2\nonumber \\&\qquad \leqslant \op{Lip}_{X,g}(\nabla_x J(\cdot,y_n))\norm{x_{n+1}-x_{n}}_{X}^2+ \tau_X\norm{x_{n+1}-x_{n}}_{X}^2 \nonumber \\&\qquad \leqslant \tau_x\norm{x_{n+1}-x_{n}}_{X}^2 .\label{gc2} \nonumber \end{align} \tag{ 4.10 }$$

By equivalent argument,

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\norm}[1]{\left\Vert #1 \right\Vert} \displaystyle g(J_{n+1,n+1}) - g(J_{n,n+1}+\nabla_yJ_{n+1,n}(y_{n+1}-y_{n}))+ \tau_Y\norm{y_{n+1}-y_{n}}_{Y}^2 \leqslant \tau_y\norm{y_{n+1}-y_n}_{Y}^2. \label{gc4} \nonumber \end{align} \tag{ 4.11 }$$

Summing equations (4.8)–(4.11) gives

$$\begin{align*} \newcommand{\e}{{\rm e}} \newcommand{\norm}[1]{\left\Vert #1 \right\Vert} \displaystyle E_{n+1,n+1} + \tau_X\norm{x_{n+1}-x_n}_{X}^2+\tau_Y\norm{y_{n+1}-y_n}_{Y}^2\leqslant E_{n,n}.\nonumber \end{align*}$$

This immediately gives the monotone decrease property of $E_{n,n}$. If we also sum this over $n$ then we achieve the final statement of the theorem:

$$\begin{align*} \newcommand{\e}{{\rm e}} \newcommand{\norm}[1]{\left\Vert #1 \right\Vert} \displaystyle \sum_{n=1}^\infty\tau_X\norm{x_{n+1}-x_n}_{X}^2+\tau_Y\norm{y_{n+1}-y_n}_{Y}^2\leqslant E_{0,0} - \lim E_{n,n}\leqslant E_{0,0}\nonumber \end{align*}$$

□

First we follow the work by Drusvyatskiy et al [46] we shall define criticality in terms of the slope of a function.

Definition 4.4. We shall say that $x_*$ is a critical point of $\newcommand{\F}[1]{\mathbb{#1}} F\colon X\to \F R$ if

$$\begin{align*} \newcommand{\e}{{\rm e}} \displaystyle |\partial F(x_*)|= 0\nonumber \end{align*}$$

where we define the slope of $F$ at $x_*$ to be

$$\begin{align*} \newcommand{\e}{{\rm e}} \newcommand{\norm}[1]{\left\Vert #1 \right\Vert} \displaystyle |\partial F(x_*)| = \limsup_{{\rm d}x\to 0} \frac{\max(0,F(x_*)-F(x_*+{\rm d}x))}{\norm{{\rm d}x}}.\nonumber \end{align*}$$

The first point to note is that this definition generalises the concept of a first order critical point for both smooth functions and convex functions (in terms of the convex sub-differential). In particular

$$\begin{align*} \newcommand{\e}{{\rm e}} \newcommand{\IP}[2]{\left\langle #1,#2 \right\rangle} \newcommand{\re}{{\rm Re}} \newcommand{\norm}[1]{\left\Vert #1 \right\Vert} \displaystyle F\in C^1 \Longrightarrow &|\partial F(x_*)| = \max\left(0,\sup_{\norm{{\rm d}x}= 1}-\IP{\nabla F(x_*)}{{\rm d}x}\right) = \norm{\nabla F(x_*)} \nonumber \\{{\rm Hence~}} &|\partial F(x_*)| = 0 \iff \norm{\nabla F(x_*)} = 0\iff \nabla F(x_*) = 0 \nonumber \end{align*}$$

$$\begin{align*} \newcommand{\e}{{\rm e}} \newcommand{\argmin}{{\rm argmin}} \newcommand{\norm}[1]{\left\Vert #1 \right\Vert} \displaystyle F{{\rm ~convex,~hence~}} & x^*\in \argmin ~F \iff F(x_*) \leqslant F(x_*+{\rm d}x) {{\rm ~for~all~}} {\rm d}x \nonumber \\{{\rm Hence~}} &\hspace{-1pt}|\partial F(x_*)| = 0 \iff \forall {\rm d}x, 0\geqslant \frac{F(x_*)-F(x_*+{\rm d}x)}{\norm{{\rm d}x}} \iff x_*\in\argmin~ F. \nonumber \end{align*}$$

For a differentiable function we cannot tell whether a critical point is a local minimum, maximum or saddle point. In general, this is also true for definition 4.4, however, at points of non-differentiability there is a bias towards local minima. This can be seen in the following example.

Example 4.5. Consider $\newcommand{\norm}[1]{\left\Vert #1 \right\Vert} F(x) = -\norm{x}$

$$\begin{align*} \newcommand{\e}{{\rm e}} \newcommand{\norm}[1]{\left\Vert #1 \right\Vert} \displaystyle |\partial F(0)| = \limsup_{{\rm d}x\to 0}\;\max\left(0,\frac{0+\norm{0+{\rm d}x}}{\norm{{\rm d}x}}\right) = 1\nonumber \end{align*}$$

Hence, 0 is not a critical point of $F$. This bias is due to the $\limsup$ in the definition which detects the strictly negative directional derivatives. This does not affect smooth functions as directional derivatives must vanish continuously to 0 about a critical point.

Some more examples are shown in figure 7. Now we shall show that our iterative sequence converges to a critical point in this sense.

Zoom In Zoom Out Reset image size

Lemma 4.6. If $(x_n,y_n)$ are as given by algorithm 1 and $X, Y$ are finite dimensional spaces then every limit point of $(x_n,y_n)$, e.g. $(x_*,y_*)$, is a critical point of $E$ in each coordinate direction. i.e.

$$\begin{align*} \newcommand{\e}{{\rm e}} \displaystyle |\partial_xE(x_*,y_*)| = |\partial_yE(x_*,y_*)| = 0.\nonumber \end{align*}$$

Remark 4.7. 

• Finite dimensionality of $X$ and $Y$ accounts for what is referred to as 'assumption 3' in [31] and is some minimal condition which ensures that the limit is also a stationary point of our iteration (equations (4.6) and (4.7)).
• The condition for finite dimensionality, as we shall see, does not directly relate to the non-convexity. The difficulty of showing convergence to critical points in infinite dimensions is common across both convex [39] and non-convex [31, 42] optimisation.

Proof. First we recall that, in finite dimensional spaces, convex functions are continuous on the relative interior of their domain [47]. Also note that by our choice of $\tau_x$ in lemma 4.6, for all $x,x',y$ we have

$$\begin{align*} \newcommand{\e}{{\rm e}} \newcommand{\norm}[1]{\left\Vert #1 \right\Vert} \displaystyle |g(J(x,y)+\nabla_x J(x,y)(x'-x)) &- g(J(x',y))| \nonumber \\&\leqslant |g([J(x,y)-J(x',y)]-\nabla_x J(x,y)(x'-x))| \nonumber \\&= |g(\nabla_x J(\xi,y)(x'-x)-\nabla_x J(x,y)(x'-x))| \nonumber \\&\leqslant \tau_x\norm{\xi-x}_X\norm{x'-x}_X \nonumber \\&\leqslant \tau_x\norm{x'-x}_X^2 \nonumber \end{align*}$$

where existence of such $\xi$ is given by the mean value theorem. Hence, for all $x$ we have

$$\begin{align*} \newcommand{\e}{{\rm e}} \newcommand{\norm}[1]{\left\Vert #1 \right\Vert} \displaystyle E(x_{n+1},y_n) & = f_{n+1,n} + g(J_{n+1,n}) \nonumber \\&\leqslant f_{n+1,n} + g(J_{n,n}+\nabla_xJ_{n,n}(x_{n+1}-x_n)) + \tau_x\norm{x_{n+1}-x_n}_X^2 \nonumber \\&\leqslant f(x,y_n)+g(J_{n,n}+\nabla_xJ_{n,n}(x-x_n))+\tau_x\norm{x-x_n}_X^2 \nonumber \\&\leqslant f(x,y_n)+g(J(x,y_n))+2\tau_x\norm{x-x_n}_X^2 \nonumber \\&= E(x,y_n) + 2\tau_x\norm{x-x_n}_X^2 \nonumber \end{align*}$$

where the first and third inequality are due to the condition shown above and the second is due to the definition of $x_{n+1}$ in (4.6). Finally, by continuity of $f$, $J$ and $g$ we can take limits on both sides of this inequality:

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\norm}[1]{\left\Vert #1 \right\Vert} \displaystyle \label{conv to stationary} \Longrightarrow E(x_*,y_*) \leqslant E(x,y_*) +2\tau_x\norm{x-x_*}_X^2 {{\rm ~for~all~}} x. \nonumber \end{align} \tag{ 4.12 }$$

This completes the proof for $x_*$ as

$$\begin{align*} \newcommand{\e}{{\rm e}} \newcommand{\norm}[1]{\left\Vert #1 \right\Vert} \displaystyle |\partial_xE(x_*,y_*)| = \limsup_{x\to x_*} \max\left(0,\frac{E(x_*,y_*)-E(x,y_*)}{\norm{x_*-x}_X}\right) \leqslant \limsup 2\tau_x\norm{x-x_*}_X = 0.\nonumber \end{align*}$$

The proof for $y_*$ follows by symmetry. □

Remark 4.8. 

• The important line in this proof, and where we need finite dimensionality, is being able to pass to the limit for (4.12). In the general case we can only expect to have $(x_n,y_n)\rightharpoonup (x_*,y_*)$, typically guaranteed by choice of regularisation functionals as in our theorem 3.1. In this reduced setting the left hand limit of (4.12) still remains valid,
  $$\begin{align*} \newcommand{\e}{{\rm e}} \displaystyle E(x_*,y_*) \leqslant \liminf E(x_{n+1},y_n) {{\rm ~by~weak~lower~semi\mbox{-}continuity.}}\nonumber \end{align*}$$
  However, on the right hand side we require:
  $$\begin{align*} \newcommand{\e}{{\rm e}} \newcommand{\norm}[1]{\left\Vert #1 \right\Vert} \displaystyle \lim E(x,y_n) + 2\tau_x\norm{x-x_n}_X^2 \leqslant E(x,y_*) + 2\tau_x\norm{x-x_*}_X^2\nonumber \end{align*}$$
  In particular, we already require $\newcommand{\norm}[1]{\left\Vert #1 \right\Vert} \norm{x-x_n}_X$ to be weakly upper semi-continuous. Topologically, this is the statement that weak and norm convergence are equivalent which will not be the case in most practical (infinite dimensional) examples.
• The properties we derive for $(x_*,y_*)$ are actually slightly stronger than that of definition 4.4 which only depends on an infinitesimal ball about $(x_*,y_*)$. However, (4.12) gives us a quantification for the more global optimality of this point. This is seen in figure 8.

Zoom In Zoom Out Reset image size

Proof. Fix arbitrary $(x_0,y_0)\in X\times Y$ and $\tau_X,\tau_Y\geqslant0$. Let $(x_n,y_n)$ be defined as in algorithm 1. By lemma 4.3 we know that $\newcommand{\F}[1]{\mathbb{#1}} \newcommand{\st}{{\rm \;s.t.\;}} \{(x_n,y_n)\st n\in\F N\}$ is contained in a sub-level set of $E$ which in turn must be weakly compact by (4.1). The assumption of lemma 4.6 is that we are in a finite dimensional space and so weak compactness is equivalent to norm compactness, i.e. some subsequence of $(x_n,y_n)$ converges in norm. Also by lemma 4.6 we know that the limit point of this sequence must be an appropriate critical point. □

All numerics are implemented in MATLAB 2016b. The sub-problem for $u$ is solved with a PDHG algorithm [39] while the sub-problem for $v$ is solved using the MOSEK solver via CVX [40, 49], the step size $\tau$ is adaptively calculated. The initial point of our algorithm is always chosen to be a good TV reconstruction, i.e.

$$\begin{align*} \newcommand{\e}{{\rm e}} \newcommand{\argmin}{{\rm argmin}} \newcommand{\TV}{{\rm TV}} \newcommand{\norm}[1]{\left\Vert #1 \right\Vert} \newcommand{\C}[1]{\mathcal{#1}} \displaystyle u_0 = \mathop{\argmin}\limits_{u\geqslant 0}\frac12\norm{S\C Ru-b}_2^2 + \lambda\TV(u), \qquad v_0 = \C R u_0.\nonumber \end{align*}$$

For clarity, we shall restate our full model with all of the parameters it includes. We seek to minimise the functional (3.1):

$$\begin{align*} \newcommand{\e}{{\rm e}} \newcommand{\op}[1]{{\rm #1}} \newcommand{\norm}[1]{\left\Vert #1 \right\Vert} \newcommand{\C}[1]{\mathcal{#1}} \displaystyle E(u,v) = \frac{1}{2}\norm{\C Ru-v}^2_{\alpha_1} + \frac{\alpha_2}{2}\norm{S\C Ru-b}_2^2 + \frac{\alpha_3}{2}\norm{Sv-b}_2^2 + \beta_1\op{TV}(u) + \beta_2\norm{A_{Ru}\nabla v}_{2,1} \nonumber \end{align*}$$

$$\begin{align*} \newcommand{\e}{{\rm e}} \renewcommand{\vec }{\mathbf} \newcommand{\re}{{\rm Re}} \displaystyle A_d(x) &= c_1(x|\lambda_1-\lambda_2, \lambda_1+\lambda_2)\vec {e}_1(x)\vec {e}_1(x)^T + c_2(x|\lambda_1-\lambda_2, \lambda_1+\lambda_2)\vec {e}_2(x)\vec {e}_2(x)^T \nonumber \\{{\rm where~}} &(\nabla d_\rho\nabla d_\rho^T)_\sigma = \lambda_1 \vec {e}_1\vec {e}_1^T + \lambda_2 \vec {e}_2\vec {e}_2^T {{\rm ~is~a~pointwise~eigenvalue~decomposition}} \nonumber \\& c_1(x|\Delta,\Sigma) = 10^{-6}+\frac{\tanh(\Sigma(x))}{1+\beta_3\Delta(x)^2},\quad c_2(x|\Delta,\Sigma) = 10^{-6}+\tanh(\Sigma(x)). \nonumber \end{align*}$$

We chose these particular $c_i$ according to two simple heuristics. If $\Sigma$ is large (steep gradients) then it is likely a region with edges and so the regularisation should be largest but still bounded above. If $\Delta=0^+$ then there is a small or blurred 'edge' present and so we want to encourage it to become a sharp jump, i.e. $\partial_\Delta c_1<0$. Theorem 2.1 tells us that choosing $c_i$ as functions of $\Delta^2$ will guarantee accordance with our later convergence results; this leads to our natural choice above. The number of iterations for algorithm 1 was chosen to be 200 and 100 for the synthetic and experimental datasets, respectively. To simplify the process of choosing values for the remaining hyper-parameters we made several observations:

• (i)  
  The choice of $\alpha_i$ and $\beta_i$ appeared to be quite insensitive about the optimum. It is clear within 2–3 iterations whether values are of the correct order of magnitude. After this, values were only tuned coarsely. For example, $\alpha_3$ and $\beta_i$ are optimal within a factor of $\newcommand{\sfrac}[2]{^{#1}\!\!/\!_{#2}} 10^{\pm\sfrac12}$.
• (ii)  
  We can chose $\alpha_2=1$ without any loss of generality. In which case, in general, $\beta_1$ should the same order of magnitude as when performing the TV reconstruction to get $u_0, v_0$.
• (iii)  
  $\alpha_2$ pairs $u$ to the given data and $\alpha_1$ pairs $u$ to the inpainted data, $v$. As such, $\alpha_1$ is spatially varying but should be something like a distance to the non-inpainting region. We chose the binary metric so that $u$ is paired to $v$ uniformly on the inpainting region and not at all outside.
• (iv)  
  $\newcommand{\DTV}{{\rm DTV}} \DTV$ specific parameters ($\beta_2, \beta_3, \rho, \sigma$) can be chosen outside of the main reconstruction. These were chosen by solving the toy problem:

  $$\begin{align*} \newcommand{\e}{{\rm e}} \newcommand{\argmin}{{\rm argmin}} \newcommand{\norm}[1]{\left\Vert #1 \right\Vert} \displaystyle \argmin ~\frac12\norm{v-v_0}_2^2 + \beta_2\norm{A_{v_0}\nabla v}_{2,1}\nonumber \end{align*}$$

  which is a lot faster to solve. $\rho>0$ is required for the analysis and so this was fixed at 1. $\sigma$ is a length-scale parameter which indicates the separation between distinct edges. $\beta_3$ relies on the normalisation of the data. As can be seen in table 1, for the two synthetic examples, with same discretisation and scaling, these values are also consistent. The only value which changes is $\beta_2$, as expected, which weights how valid the $\newcommand{\DTV}{{\rm DTV}} \DTV$ prior is for each dataset.

Table 1. Parameter choices for numerical experiments. Each algorithm was run for 300 iterations.

	$\alpha_1$	$\alpha_2$	$\alpha_3$	$\beta_1$	$\beta_2$	$\beta_3$	$\rho$	$\sigma$
Concentric rings phantom	$\newcommand{\F}[1]{\mathbb{#1}} \newcommand{\1}[0]{\F{1}} \frac1{2^2}\1_{\Omega'^c}$	1	$1\times10^{-1}$	$3\times10^{-5}$	$3\times10^{3}$	$10^{10}$	$1$	$8$
Shepp–Logan phantom	$\newcommand{\F}[1]{\mathbb{#1}} \newcommand{\1}[0]{\F{1}} \frac1{4^2}\1_{\Omega'^c}$	1	$3\times10^{-1}$	$3\times10^{-5}$	$3\times10^{2}$	$10^{10}$	$1$	$8$
Experimental dataset (both sampling ratios)	$\newcommand{\F}[1]{\mathbb{#1}} \newcommand{\1}[0]{\F{1}} \frac1{2^2}\1_{\Omega'^c}$	1	$3\times10^{2}$	$1\times10^{-5}$	$3\times10^{1}$	$10^6$	$1$	$0$

It is unclear whether a gridsearch may provide better results although, due to the number of parameters involved, this would definitely take a lot longer and mask some interpretability of the parameters. A further comparison of different choices of the main parameters can be found in the supplementary material.

This example shows two concentric rings. This is the canonical example for our model because the exact sinogram is perfectly radially symmetric which should trivialise the directional inpainting procedure, even with noise present. As can clearly be seen in figure 9, the TV reconstruction is poor in the missing wedge direction which can be seen as a blurring out of the sinogram. By enforcing better structure in the sinogram, our proposed joint model is capable of extrapolating these local structures from the given data domain to recover the global structure and gives an accurate reconstruction.

Zoom In Zoom Out Reset image size

This example shows the modified Shepp–Logan phantom which is built up as a sum of ellipses. This example has a much more complex geometry although the sinogram still has a clear geometry. In figure 10 we see that the largest scale feature, the shape of the largest ellipse, is recovered in our proposed reconstruction with minimal loss of contrast in the interior. One artefact we have not been able to remove is the two rays extending from the top of the reconstructed sample. Looking more closely we found that it was due to a small misalignment of the edge at the bottom of the sinogram as it crosses between the data to the inpainting region. Numerically, this happens because of the convolutions which take place inside the directional TV regularisation functional. Having a non-zero blurring is essential for regularity of the regularisation (theorem 2.1) but the effect of this is that it does not heavily penalise misalignment on such a small scale. This means that at the interface between the fixed data-term there is a slight kink, the line is continuous but not $C^1$. The effect of this on the reconstruction is the two lines which extend from the sample at this point. Looking at quantitative measures, the PSNR value rises from 17.33 to 17.36 whereas the SSIM decreases from 0.76 to 0.62, from TV to the proposed reconstruction, respectively. These measures are inconclusive and the authors feel that they fail to balance the improvement to global geometry verses more local artefacts in the reconstructions.

Zoom In Zoom Out Reset image size

The sample is a silver bipyramidal crystal placed on a planar surface, and the challenges of this dataset are shown in figure 11. We immediately see that the wide angle projections have large artefacts which produces a very low signal to noise ratio. Another issue present is that there is mass seen in some of the projections which cannot be represented within the reconstruction volume. Both of these issues violate the simple x-ray model that is used. Exact modelling would require estimation of parameters which are not available a priori and so the preferred acquisition is one which automatically minimises these modelling errors. Another artefact is that over time each surface becomes coated with carbon. This is a necessary consequence of the sample preparation and this coating is known to occur during the microscopy. The result of modelling errors and time dependent noise is to prefer an acquisition with limited angular range and earliest acquired projections. Because of this, in numerical experiments we compare both TV and our proposed reconstruction using only $\newcommand{\sfrac}[2]{^{#1}\!\!/\!_{#2}} \sfrac34$ of the available data, 29 projections over an $87^\circ$ interval, with a bias towards earlier projections. The artefacts due to the out-of-view mass are unavoidable but we can perform some further pre-processing to minimise the effect. In particular, if we shrink the field of view of the detector then the 'heaviest' part of the data will be the particle of interest and the model violations will be relatively small, increasing the signal to noise ratio. This can be seen as the sharp horizontal cut-off in the pre-processed sinograms seen on the right of figure 11. The effect of this on the reconstruction is going to be that there is a thin ring of mass placed at the edge of the (shrunken) detector view which will be clearly identifiable in the reconstruction. As a ground truth approximation we shall use a TV reconstruction on the full data for the location of the particle alongside prior knowledge of this sample for more precise geometrical features. We also note that the particle should be very homogeneous so this is another example where we expect the TV reconstruction to be very good.

Zoom In Zoom Out Reset image size

The sample is a single crystal of silver and so we know it must have very uniform intensity and we are interested in locating the sharp facets which bound the crystal [48]. In figure 12 we immediately see that the combination of homogeneity and sharp edges is better reconstructed in our proposed reconstruction. Because we expect the reconstruction to be constant on the background and the particle, thresholding the reconstruction allows us to easily locate the boundaries and estimate interior angles of the particle. Figure 13 shows such images where the threshold is chosen to be the approximate midpoint of these two levels. We see that the proposed reconstruction consistently has less jagged edges and the left hand corner is better curved, as is consistent with our knowledge of the sample. Looking back at the full colour images we see that this is a result of lack of sharp decay at the boundary and homogeneity inside the sample. Looking for location error we see the biggest error in both TV and joint reconstruction is on the bottom-left edge where both reconstructions pull the line inwards. However, looking particularly at points $(40,80)$ and $(20,60)$, we see that this was less severe in the proposed method. The other missing wedge artefact is in the top right corner which has been extended in both reconstructions although it is thinner in the proposed reconstruction. This indicates that it was better able to continue the straight edges either side of the corner and the blurring in the missing wedge direction is more localised than in the TV reconstruction. Overall, we see see that the proposed reconstruction method is much more robust to an decrease in angular sampling range.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size