# Directional Sinogram Inpainting for Limited Angle Tomography

Robert Tovey1, Martin Benning2, Christoph Brune3, Marinus J. Lagerwerf4, Sean
M. Collins5, Rowan K. Leary5, Paul A. Midgley5, Carola-Bibiane Schönlieb1 1
Centre for Mathematical Sciences, University of Cambridge  2 School of
Mathematical Sciences, Queen Mary University of London  3 University of Twente
4 Centrum Wiskunde & Informatica  5 Department of Materials Science and
Metallurgy, University of Cambridge

()

###### Abstract

In this paper we propose a new joint model for the reconstruction of
tomography data under limited angle sampling regimes. In many applications of
Tomography, e.g. Electron Microscopy and Mammography, physical limitations on
acquisition lead to regions of data which cannot be sampled. Depending on the
severity of the restriction, reconstructions can contain severe,
characteristic, artefacts. Our model aims to address these artefacts by
inpainting the missing data simultaneously with the reconstruction.
Numerically, this problem naturally evolves to require the minimisation of a
non-convex and non-smooth functional so we review recent work in this topic
and extend results to fit an alternating (block) descent framework. We perform
numerical experiments on two synthetic datasets and one Electron Microscopy
dataset. Our results show consistently that the joint inpainting and
reconstruction framework can recover cleaner and more accurate structural
information than the current state of the art methods.

##  1 Introduction

###  1.1 Problem Formulation

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, 3, 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

$\text{Given data }b\text{ find optimal pair }(u,v)\text{ such that
}Sv=b,\mathcal{R}u=v$

where $\mathcal{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.

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

Figure 1: Diagrammatic representation of the acquisition of 2D X-ray transform
data , the sinogram, in both full range and limited angle acquisition. Note
that measurement at $180^{\circ}$ is exactly a reflection of that at
$0^{\circ}$. This symmetry allows us to consider a $180^{\circ}$ range of the
sinogram as a full sample. In the limited angle setting we can only rotate the
sample a small amount clockwise and anti-clockwise which results in missing
data in the middle of the sinogram.

###  1.2 Context and Proposed Model

Traditional methods for X-ray Tomography reconstruction find approximate
solutions to $S\mathcal{R}u=b$, constraining $\mathcal{R}u=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, 13, 14].

  * •

The Simultaneous Iterative Reconstruction Technique (SIRT) can be thought of
as a preconditioned gradient descent on the function $\left\lVert
S\mathcal{R}u-b\right\rVert_{2}^{2}$ [15, 3, 16, 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 [18, 2, 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$.

Figure 2: Demonstration of TV reconstruction in comparison to FBP and SIRT.
The top row shows reconstructions from noise-less limited angle data and the
bottom shows reconstructions from noisy limited view data (far left images).
Comparing the columns, we immediately see that FBP and SIRT are much more
prone to angular artefacts than TV. In both cases we notice that the TV
reconstructions better show the broad structure of the phantom.

Figure 3: Examples when TV reconstructions cannot recover the global
structures of samples. When there is a large missing wedge
(${}^{2}\\!\\!/\\!_{3}$ of data unseen) and noise on the projections then
reconstructions exhibit characteristic blurring in the vertical direction.
This can also be seen in the extrapolated region of the sinograms as a loss of
structure.

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:

$\displaystyle v$ $\displaystyle=\text{ optimal inpainted sinogram given $b$}$
$\displaystyle u$ $\displaystyle=\text{ optimal reconstruction given $v$}$

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:

$E(u,v)=\alpha_{1}d_{1}(\mathcal{R}u,v)+\alpha_{2}d_{2}(S\mathcal{R}u,b)+\alpha_{3}d_{3}(Sv,b)+\alpha_{4}r_{1}(u)+\alpha_{5}r_{2}(v)$
(1.1)

where $\alpha_{i}\geq 0$ 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 that the flexibility of our joint model, (1.1), can allow
for a better adaptive fitting to the data.

Figure 4: Demonstration of the improvement which can be achieved by using a
model as in (1.1). Left hand images show state of the art reconstructions
using Total Variational regularisation ($\alpha_{1}=\alpha_{3}=\alpha_{5}=0$).
This reconstruction clearly shows the characteristic blurring artefacts at the
top and bottom. In our proposed joint reconstruction ( right hand) we can
minimise these effects.

###  1.3 Overview and Contributions

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.

##  2 Preliminaries

###  2.1 The X-ray Transform

The principal forward model for X-ray tomography is provided by the X-ray
transform which can be defined for any bounded domain,
$\Omega\subset\mathds{R}^{n}$, by

$\mathcal{R}\colon L^{1}(\Omega,\mathds{R})\to
L^{1}(\mathcal{S}^{n-1}\times\mathds{R}^{n},\mathds{R})\text{ such that
}\mathcal{R}u(\theta,y)=\int_{x=y+t\theta,t\in\mathds{R}}u(x)dt$ (2.1)

where
$\mathcal{S}^{n-1}=\\{s\in\mathds{R}^{n}\operatorname{\;s.t.\;}|s|=1\\}$. In
this work, for simplicity, we will only be using $n=2$ although the case $n=3$
is completely analogous.

###  2.2 Total Variation Regularisation

Total Variation (TV) regularisation is extremely common across many fields of
image processing [32, 33, 34]. The definition of the (isotropic) TV semi-norm
on domains $\Omega\subset\mathds{R}^{n}$ is stated as

$g\in
L^{1}(\Omega,\mathds{R})\implies\operatorname{TV}(g)=\sup\left\\{\int\left\langle
g(x),\operatorname{div}(\varphi)(x)\right\rangle dx\right.\\\
\left.\vphantom{\int}\operatorname{\;s.t.\;}\varphi\in
C_{c}^{\infty}(\Omega,\mathds{R}^{n}),|\varphi(x)|_{2}\leq 1\text{ for all
}x\right\\}$ (2.2)

The intuition behind this is that when $g\in W^{1,1}(\Omega,\mathds{R})$ then
we have exactly

$\operatorname{TV}(g)=\left\lVert\nabla g\right\rVert_{2,1}=\int|\nabla
g(x)|_{2}dx$

From this point forward we shall write the
$\left\lVert\cdot\right\rVert_{2,1}$ representation where $|\nabla g|$ is to
be understood as a measure if necessary. We shall denote the space of Bounded
Variation as

$\operatorname{BV}(\Omega,\mathds{R})=\\{g\in
L^{1}\operatorname{\;s.t.\;}\operatorname{TV}(g)<\infty\\}$

One property that we shall use about the space of Bounded Variation is
$\operatorname{BV}\subset L^{2}\subset L^{1}$ which holds whenever $\Omega$ is
compact by the Poincaré inequality: $\left\lVert g-^{\int g}\\!\\!/\\!_{\int
1}\right\rVert_{2}\lesssim\operatorname{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 [18, 2].

$u=\operatorname*{argmin}_{u\geq 0}\frac{1}{2}\left\lVert
S\mathcal{R}u-b\right\rVert_{2}^{2}+\lambda\operatorname{TV}(u)\text{ for some
}\lambda\geq 0$ (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$.

###  2.3 Directional Total Variation Regularisation

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, 37, 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

$\operatorname{DTV}(v)=\int|A(x)\nabla v(x)|dx=\left\lVert A\nabla
v\right\rVert_{2,1}\text{ for some anisotropic
}A\colon\mathds{R}^{2}\to\mathds{R}^{2\times 2}.$

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{split}A(x)=&\;c_{1}(x)\mathbf{e}_{1}(x)\mathbf{e}_{1}(x)^{T}+c_{2}(x)\mathbf{e}_{2}(x)\mathbf{e}_{2}(x)^{T}\\\
&\text{ such that
}\mathbf{e}_{i}\colon\mathds{R}^{2}\to\mathds{R}^{2},|\mathbf{e}_{i}(x)|=1,\left\langle\mathbf{e}_{1}(x),\mathbf{e}_{2}(x)\right\rangle=0\end{split}$
(2.4)

i.e.

$\operatorname{DTV}(v)=\int\sqrt{c_{1}^{2}|\left\langle\mathbf{e}_{1},\nabla
v\right\rangle|^{2}+c_{2}^{2}|\left\langle\mathbf{e}_{2},\nabla
v\right\rangle|^{2}}dx.$

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 $\mathbf{e}_{1}$. This
allows for large jumps in the direction of $\mathbf{e}_{1}$ whilst maintaining
flatness in the direction of $\mathbf{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:

$(\nabla d_{\rho}\nabla
d_{\rho}^{T})_{\sigma}(x)=\lambda_{1}(x)\mathbf{e}_{1}(x)\mathbf{e}_{1}(x)^{T}+\lambda_{2}(x)\mathbf{e}_{2}(x)\mathbf{e}_{2}(x)^{T}\text{
such that }\lambda_{1}(x)\geq\lambda_{2}(x)\geq 0$

where

$d_{\rho}(x)=[d\star\exp(-^{|\cdot|^{2}}\\!\\!/\\!_{2\rho^{2}})](x)=\int\exp(-^{|y-x|^{2}}\\!\\!/\\!_{2\rho^{2}})d(y)\text{
etc.}$

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

$\Delta(x)=\lambda_{1}(x)-\lambda_{2}(x)\text{ is
coherence}\qquad\Sigma(x)=\lambda_{1}(x)+\lambda_{2}(x)\text{ is energy}$

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

$A_{d}(x)\coloneqq
c_{1}(x|\Delta,\Sigma)\mathbf{e}_{1}(x)\mathbf{e}_{1}(x)^{T}+c_{2}(x|\Delta,\Sigma)\mathbf{e}_{2}(x)\mathbf{e}_{2}(x)^{T}$
(2.5)

where the only parameters left to choose are $c_{i}$. Typical examples of
include

$c_{1}=\frac{1}{\sqrt{1+\Sigma^{2}}},\qquad c_{2}=1$ $c_{1}=\varepsilon,\qquad
c_{2}=\varepsilon+\exp\left(-^{1}\\!\\!/\\!_{\Delta^{2}}\right)\text{ for some
}\varepsilon>0$

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}$.

Figure 5: Examples comparing TV with directional TV for inpainting and
denoising. Both examples have the same edge structure and so parameters in
(2.4) are the same in both. DTV uses $c_{2}=1$ and $c_{1}$ as the indicator (0
or 1) shown in the bottom left plot, TV is the case $c_{1}=c_{2}=1$. Left
block: Top left image is inpainting input where the dark square shows the
inpainting region. The structure of $c_{1}$ allows DTV (bottom right) to
connect lines over arbitrary distances, whereas TV inpainting (top right) will
fail to connect the lines if the horizontal distance is greater than the
vertical separation of the edges. Right block: Top left image is denoising
input. DTV has two advantages. Firstly, the structure of $c_{1}$ recovers a
much straighter line than that in the TV reconstruction. Secondly, TV
penalises jumps equally in each direction and so the contrast is reduced, DTV
is able to penalise noise oscillations independently from edge discontinuities
which allows us to maintain much higher contrast. Figure 6: Surface
representing a characteristic image, $d$, to demonstrate the behaviour of
$\Sigma$ and $\Delta$. Away from edges (A) we have $\Sigma\approx\Delta\approx
0$. On simple edges (B) we have $\Sigma\approx\Delta\gg 0$ and, finally, at
corners (C) we have $\Sigma\gg\Delta$.

######  Theorem 2.1.

If

  1. 1.

$c_{i}$ are $2k$ times continuously differentiable in $\Delta$ and $\Sigma$,
$k\geq 1$

  2. 2.

$c_{1}(x|0,\Sigma)=c_{2}(x|0,\Sigma)$ for all $x$ and $\Sigma\geq 0$

  3. 3.

$\partial_{\Delta}^{2j-1}c_{1}(x|0,\Sigma)=\partial_{\Delta}^{2j-1}c_{2}(x|0,\Sigma)=0$
for all $x$ and $\Sigma\geq 0,j=1\ldots,k$

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

######  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 A.

##  3 The Joint Model

Now that all of the notation and concepts have been defined we can formalise
the statement of our particular joint model of the form in (1.1):

  * •

The forward operator, $\mathcal{R}$, is the X-ray transform from (2.1)

  * •

The desired reconstructed sample, $u\in\operatorname{BV}(\Omega,\mathds{R})$
on some domain $\Omega$

  * •

The noisy sub-sampled data, $b\in L^{1}(\Omega^{\prime},\mathds{R})$ on some
$\Omega^{\prime}\Subset\mathcal{S}^{1}\times\mathds{R}_{\geq 0}$. We extend
such that $b|_{\Omega^{\prime c}}=0$ for notational convenience.

  * •

The full reconstructed sinogram, $v\in
L^{1}(\mathcal{S}^{1}\times\mathds{R}_{\geq 0},\mathds{R})$

We also define $S=S_{\Omega^{\prime}}$ to be the indicator function on
$\Omega^{\prime}$. The joint model that we thus propose is

$\begin{split}(u,v)=\operatorname*{argmin}_{u\geq
0}E(u,v)=&\;{\operatorname*{argmin}_{u\geq
0}}\frac{1}{2}\left\lVert\mathcal{R}u-v\right\rVert^{2}_{\alpha_{1}}+\frac{\alpha_{2}}{2}\left\lVert
S\mathcal{R}u-b\right\rVert_{2}^{2}\\\
&\qquad\quad+\frac{\alpha_{3}}{2}\left\lVert
Sv-b\right\rVert_{2}^{2}+\beta_{1}\operatorname{TV}(u)+\beta_{2}\operatorname{DTV}_{\mathcal{R}u}(v)\end{split}$
(3.1)

where

$\operatorname{DTV}_{\mathcal{R}u}(v)=\left\lVert A_{\mathcal{R}u}\nabla
v\right\rVert_{2,1}$

and $\alpha_{i},\beta_{i}>0$ are weighting parameters, $A_{Ru}$ as defined in
(2.5). $\alpha_{1}$ is embedded in the norm as it is a spatially varying
weight, taking different values inside and outside of $\Omega^{\prime}$. We
note that the norms involving $b$ are determined by the noise model of the
acquisition process, in this case Gaussian noise. The final metric pairing
$\mathcal{R}u$ and $v$ was free to be chosen to promote any structural
similarity between the two quantities. We have chosen the squared $L^{2}$ norm
for simplicity although if some structure is known to be important then there
is a wide choice of specialised functions from which to choose (e.g. [1]).

The choice of regularisation functionals reflects prior assumptions on the
expected type of sample; all of the examples shown later will follow these
assumptions. The isotropic TV penalty is chosen as $u$ is expected to be
piece-wise constant. This will reduce oscillations from $u$ and favour ‘stair-
case’ like changes of intensity over smooth ones. The assumptions of our
regularisation on $v$ must also be derived from expected properties of $u$.
What is known from [28], and can be seen in Figure 3, is that discontinuities
of $u$ along curves in the spatial domain, say $\gamma$, generate a so called
‘dual curve’ in the sinogram. $\mathcal{R}u$ will also have an edge, although
possibly not a discontinuity, along this dual curve. Thus, perpendicular to
the dual curve $v$ should have sharp features and parallel to the dual curve
intensity should vary slowly. The assumption of our regularisation is that if
a dual curve is present in the visible component of the data then it should
correspond to some $\gamma$ in the spatial domain. The extrapolation of this
dual curve must behave like the boundary of a level-set of $u$, i.e. preserve
the sharp edge and slow varying intensities in $v$. The main influence of this
regularisation is in the inpainting region and so any artefacts it introduces
should also only effect edges corresponding to these invisible singularities,
including streaking artefacts. Another bias that is present is an assumption
that dual curves are themselves smooth. In the inpainting region, this will
encourage dual curves with low curvature thus invisible singularities are
likely to follow near-circular arcs in the spatial domain. Final parameter
choices, such as $\alpha_{i},\beta_{i}$ and $c_{i}$, are not necessary at this
point and will be chosen in Section 5.1.

The immediate question to ask is whether this model is well posed. For a non-
convex function we typically cannot expect numerically to find global
minimisers but the following result shows we can expect some convergence to
local minima. We present the following result which justifies looking for
minima of this functional.

######  Theorem 3.1.

If

  * •

$c_{i}$ are bounded away from 0

  * •

$\rho>0$

  * •

$A_{d}$ is differentiable in $d$

then sub-level sets of $E$ are weakly compact in
$L^{2}(\Omega,\mathds{R})\times L^{2}(\mathds{R}^{2},\mathds{R})$ and $E$ is
weakly lower semi-continuous. i.e. for all $(u_{n},v_{n})\in
L^{2}(\Omega,\mathds{R})\times L^{2}(\mathds{R}^{2},\mathds{R})$:

$E(u_{n},v_{n})\text{ uniformly bounded implies a subsequence converges
weakly}$ $\liminf E(u_{n},v_{n})\geq E(u,v)\text{ whenever
}u_{n}\rightharpoonup u,v_{n}\rightharpoonup v$

The proof of this theorem is contained in B. This theorem justifies numerical
minimisation of $E$ because it tells us that all descending sequences
($E(u_{n},v_{n})\leq E(u_{n-1},v_{n-1})$) have a convergent subsequence and
any limit point must minimise $E$ over the original sequence.

##  4 Numerical Solution

To address the issue of convergence we shall first generalise our functional
and prove the result in the general setting. Functionals will be of the form
$E\colon X\times Y\to\mathds{R}$ where $X$ and $Y$ are Banach spaces and $E$
accepts the decomposition

$E(x,y)=f(x,y)+g(J(x,y))$

such that:

$\displaystyle\text{Sub-level sets of }E\text{ are weakly compact}$ (4.1)
$\displaystyle f\colon X\times Y\to\mathds{R}\text{ is jointly convex in
}(x,y)\text{ and bounded below}$ (4.2) $\displaystyle g\colon
Z\to\mathds{R}\text{ is a semi-norm on Banach space }Z\text{, i.e. for all
}t\in\mathds{R},z,z_{1},z_{2}\in Z$ $\displaystyle\qquad g(z)\geq
0,g(tz)=|t|g(z)\text{ and }g(z_{1}+z_{2})\leq g(z_{1})+g(z_{2})$ (4.3)
$\displaystyle J\colon X\times Y\to Z\text{ is }C^{1}\text{ and for all
}{K\Subset X\times Y},\;\exists L_{x},L_{y}<\infty{\text{ such that
}\forall(x,y)\in K}$ $\displaystyle\qquad
g(J(x+dx,y)-J(x,y)-\nabla_{x}J(x,y)dx)\leq L_{x}\left\lVert
dx\right\rVert_{X}^{2}$ (4.4) $\displaystyle\qquad
g(J(x,y+dy)-J(x,y)-\nabla_{y}J(x,y)dy)\leq L_{y}\left\lVert
dy\right\rVert_{Y}^{2}$ (4.5)

Note if $g$ is a norm then $L_{x}$ can be chosen to be the standard Lipschitz
factor of $\nabla_{x}J$. If $J$ is twice Frèchet-differentiable then these
constants must be finite. In our case:

$\displaystyle f(x,y)$ $\displaystyle=\frac{1}{2}\left\lVert
Rx-y\right\rVert_{\alpha_{1}}^{2}+\frac{\alpha_{2}}{2}\left\lVert
S\mathcal{R}x-b\right\rVert^{2}$
$\displaystyle\qquad+\frac{\alpha_{3}}{2}\left\lVert
Sy-b\right\rVert^{2}+\beta_{1}\operatorname{TV}(x)+\left\\{\begin{array}[]{cc}0&x\geq
0\\\ \infty&\text{else}\end{array}\right.$ $\displaystyle g(z)$
$\displaystyle=\beta_{2}\left\lVert z\right\rVert_{2,1}$ $\displaystyle
J(x,y)$ $\displaystyle=A_{Rx}\nabla
y\implies\tau_{x}\sim\beta_{2}\left\lVert\nabla^{2}A_{\cdot}\right\rVert\left\lVert
R\right\rVert\operatorname{TV}(y),\tau_{y}=0$

Finiteness of $\left\lVert\nabla^{2}A\right\rVert$ and weak compactness of
sub-level sets are given by Theorems 2.1 and 3.1 respectively. The alternating
descent algorithm is stated in Algorithm 1.

Algorithm 1

0: Any $x_{0}\in X$, $\tau_{x},\tau_{y}\geq 0$.

$n\leftarrow 0$

while not converged do

$n\leftarrow n+1$

$x_{n}=\operatorname*{argmin}_{x\in X}f(x,y_{n-1})+\tau_{x}\left\lVert
x-x_{n-1}\right\rVert_{X}^{2}\\\
+g(J(x_{n-1},y_{n-1})+\nabla_{x}J(x_{n-1},y_{n-1})(x-x_{n-1}))$ (4.6)

$y_{n}=\operatorname*{argmin}_{y\in Y}f(x_{n},y)+\tau_{y}\left\lVert
y-y_{n-1}\right\rVert_{Y}^{2}\\\
\qquad\qquad+g(J(x_{n},y_{n-1})+\nabla_{y}J(x_{n},y_{n-1})(y-y_{n-1}))$ (4.7)

end while

$(x_{n},y_{n})$

Classical alternating proximal descent would take
$x_{n+1}=\operatorname*{argmin}E(x,y_{n})+\tau_{x}\left\lVert
x-x_{n}\right\rVert_{2}^{2}$ but because of the complexity of
$A_{\mathcal{R}u}$ each sub-problem would have the same complexity as the full
problem, making it computationally infeasible. By linearising $A_{d}$ we
overcome this problem as both sub-problems are convex, for which there are
many standard solvers such as [39, 40]. This second approach is similar to
that of the ProxDescent algorithm [41, 31]. We extend this algorithm to cover
alternating descent and achieve equivalent convergence guarantees. Using
Algorithm 1, our statement of convergence is the following theorem.

######  Theorem 4.1.

Convergence of alternating minimisation: If $E$ satisfies (4.1)-(4.5) and
$(x_{n},y_{n})$ are a sequence generated by Algorithm 1 then

  * •

$E(x_{n+1},y_{n+1})\leq E(x_{n},y_{n})$ for each $n$.

  * •

A subsequence of $(x_{n},y_{n})$ must converge weakly in $X\times Y$

  * •

If $\\{(x_{n},y_{n})\operatorname{\;s.t.\;}n=1,\ldots\\}$ is contained in a
finite dimensional space then every limit point of $(x_{n},y_{n})$ must be a
critical point (as will be defined in Definition 8) of $E$ in both the
direction of $x$ and $y$.

This result is the parallel of Lemma 10, Theorem 18 and Corollary 21 in [31]
without the alternating or block descent setting. There is some overlap in the
analysis for the two settings although we present an independent proof which
is more direct and we feel gives more intuition for our more restricted class
of functionals. The rest of this section is now dedicated to the proof of this
theorem.

For notational convenience we shall compress notation such that:

$f_{n,m}=f(x_{n},y_{m}),\quad J_{n,m}=J(x_{n},y_{m}),\quad
E_{n,m}=E(x_{n},y_{m})\text{ etc.}$

###  4.1 Sketch Proof

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 $\left\lVert
x_{n}-x_{n-1}\right\rVert_{X},\left\lVert y_{n}-y_{n-1}\right\rVert_{Y}$
converge to 0. This follows by a relatively standard sufficient decrease
argument as seen in [42, 31, 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 $x,y\in\mathds{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

$\displaystyle x_{1}$
$\displaystyle=\operatorname*{argmin}E(x,y_{0})+\tau_{x}(x-x_{0})^{2}=0$
$\displaystyle y_{1}$
$\displaystyle=\operatorname*{argmin}E(x_{1},y)+\tau_{y}(y-y_{0})^{2}=0$

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,
$(-^{1}\\!\\!/\\!_{2},-^{1}\\!\\!/\\!_{2})$.

###  4.2 Sufficient Decrease Lemma

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

######  Lemma 4.3.

If for each $n$

$\tau_{x}\geq L_{x}+\tau_{X},\qquad\tau_{y}\geq L_{y}+\tau_{Y}$

for some $\tau_{X},\tau_{Y}\geq 0$ then

$\sum^{\infty}\tau_{X}\left\lVert
x_{n}-x_{n-1}\right\rVert_{X}^{2}+\tau_{Y}\left\lVert
y_{n}-y_{n-1}\right\rVert_{Y}^{2}\leq E(x_{0},y_{0})$

and

$E(x_{n+1},y_{n+1})\leq E(x_{n},y_{n})\text{ for all }n$

###### Proof.

Note by Equations (4.6), (4.7) (definition of our sequence) we have

$\displaystyle
f_{n+1,n}+g(J_{n,n}+\nabla_{x}J_{n,n}(x_{n+1}-x_{n}))+\tau_{x}\left\lVert
x_{n+1}-x_{n}\right\rVert_{X}^{2}$ $\displaystyle\leq E_{n,n}$ (4.8)
$\displaystyle
f_{n+1,n+1}+g(J_{n+1,n}+\nabla_{y}J_{n+1,n}(y_{n+1}-y_{n}))+\tau_{y}\left\lVert
y_{n+1}-y_{n}\right\rVert_{Y}^{2}$ $\displaystyle\leq E_{n+1,n}$ (4.9)

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

$\displaystyle g(J_{n+1,n})-$ $\displaystyle
g(J_{n,n}+\nabla_{x}J_{n,n}(x_{n+1}-x_{n}))+\tau_{X}\left\lVert
x_{n+1}-x_{n}\right\rVert_{X}^{2}$ $\displaystyle\qquad\leq
g(J_{n+1,n}-J_{n,n}-\nabla_{x}J_{n,n}(x_{n+1}-x_{n}))+\tau_{X}\left\lVert
x_{n+1}-x_{n}\right\rVert_{X}^{2}$
$\displaystyle\qquad=g([\nabla_{x}J(\xi)-\nabla_{x}J_{n,n}](x_{n+1}-x_{n}))+\tau_{X}\left\lVert
x_{n+1}-x_{n}\right\rVert_{X}^{2}$
$\displaystyle\qquad\leq\operatorname{Lip}_{X,g}(\nabla_{x}J(\cdot,y_{n}))\left\lVert
x_{n+1}-x_{n}\right\rVert_{X}^{2}+\tau_{X}\left\lVert
x_{n+1}-x_{n}\right\rVert_{X}^{2}$ $\displaystyle\qquad\leq\tau_{x}\left\lVert
x_{n+1}-x_{n}\right\rVert_{X}^{2}$ (4.10)

By equivalent argument,

$g(J_{n+1,n+1})-g(J_{n,n+1}+\nabla_{y}J_{n+1,n}(y_{n+1}-y_{n}))+\tau_{Y}\left\lVert
y_{n+1}-y_{n}\right\rVert_{Y}^{2}\leq\tau_{y}\left\lVert
y_{n+1}-y_{n}\right\rVert_{Y}^{2}$ (4.11)

Summing Equations (4.8)-(4.11) gives

$E_{n+1,n+1}+\tau_{X}\left\lVert
x_{n+1}-x_{n}\right\rVert_{X}^{2}+\tau_{Y}\left\lVert
y_{n+1}-y_{n}\right\rVert_{Y}^{2}\leq E_{n,n}$

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:

$\sum_{n=1}^{\infty}\tau_{X}\left\lVert
x_{n+1}-x_{n}\right\rVert_{X}^{2}+\tau_{Y}\left\lVert
y_{n+1}-y_{n}\right\rVert_{Y}^{2}\leq E_{0,0}-\lim E_{n,n}\leq E_{0,0}$

∎

###  4.3 Convergence to Critical Points

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 $F\colon X\to\mathds{R}$ if

$|\partial F(x_{*})|=0$

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

$|\partial F(x_{*})|=\limsup_{dx\to
0}\frac{\max(0,F(x_{*})-F(x_{*}+dx))}{\left\lVert dx\right\rVert}$

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

$\displaystyle F\in C^{1}\implies$ $\displaystyle|\partial
F(x_{*})|=\max\left(0,\sup_{\left\lVert dx\right\rVert=1}-\left\langle\nabla
F(x_{*}),dx\right\rangle\right)=\left\lVert\nabla F(x_{*})\right\rVert$ Hence
$\displaystyle|\partial F(x_{*})|=0\iff\left\lVert\nabla
F(x_{*})\right\rVert=0\iff\nabla F(x_{*})=0$ $\displaystyle F\text{ convex,
hence }$ $\displaystyle x^{*}\in\operatorname*{argmin}F\iff F(x_{*})\leq
F(x_{*}+dx)\text{ for all }dx$ Hence $\displaystyle|\partial
F(x_{*})|=0\iff\forall dx,0\geq\frac{F(x_{*})-F(x_{*}+dx)}{\left\lVert
dx\right\rVert}\iff x_{*}\in\operatorname*{argmin}F$

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 8, 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 $F(x)=-\left\lVert x\right\rVert$

$|\partial F(0)|=\limsup_{dx\to 0}\;\max\left(0,\frac{0+\left\lVert
0+dx\right\rVert}{\left\lVert dx\right\rVert}\right)=1$

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 doesn’t affect smooth functions as directional derivatives must vanish
continuously to 0 about a critical point.

Some more examples are shown in Figure 7.

$x>0,\varepsilon>0$$F(x)=-x^{1+\varepsilon}$ (a) Examples of Critical Points

$x>0,\varepsilon\leq 0$$F(x)=-x^{1+\varepsilon}$ (b) Examples of Non-Critical
Points

Figure 7: Examples of 1D functions where 0 is/isn’t a critical point by
Definition 8. If a function is piece-wise linear then 0 is a critical point
iff each directional derivative is non-negative. If the function can be
approximated on an interval of $x>0$ to first order terms by
$F(x)=cx^{1+\varepsilon}$ then criticality can be characterised sharply. If
$c\geq 0$ then 0 is always a critical point. If $c<0$ then 0 is a critical
point iff $\varepsilon>0$, however, 0 is never a local minimum.

Now we shall show that our iterative sequence converges to a critical point in
this sense.

######  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.

$|\partial_{x}E(x_{*},y_{*})|=|\partial_{y}E(x_{*},y_{*})|=0$

######  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)-(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
[42, 31] 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^{\prime},y$ we have

$\displaystyle|g(J(x,y)+\nabla_{x}J(x,y)(x^{\prime}-x))$
$\displaystyle-g(J(x^{\prime},y))|$
$\displaystyle\leq|g([J(x,y)-J(x^{\prime},y)]-\nabla_{x}J(x,y)(x^{\prime}-x))|$
$\displaystyle=|g(\nabla_{x}J(\xi,y)(x^{\prime}-x)-\nabla_{x}J(x,y)(x^{\prime}-x))|$
$\displaystyle\leq\tau_{x}\left\lVert\xi-x\right\rVert_{X}\left\lVert
x^{\prime}-x\right\rVert_{X}$ $\displaystyle\leq\tau_{x}\left\lVert
x^{\prime}-x\right\rVert_{X}^{2}$

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

$\displaystyle E(x_{n+1},y_{n})$ $\displaystyle=f_{n+1,n}+g(J_{n+1,n})$
$\displaystyle\leq
f_{n+1,n}+g(J_{n,n}+\nabla_{x}J_{n,n}(x_{n+1}-x_{n}))+\tau_{x}\left\lVert
x_{n+1}-x_{n}\right\rVert_{X}^{2}$ $\displaystyle\leq
f(x,y_{n})+g(J_{n,n}+\nabla_{x}J_{n,n}(x-x_{n}))+\tau_{x}\left\lVert
x-x_{n}\right\rVert_{X}^{2}$ $\displaystyle\leq
f(x,y_{n})+g(J(x,y_{n}))+2\tau_{x}\left\lVert x-x_{n}\right\rVert_{X}^{2}$
$\displaystyle=E(x,y_{n})+2\tau_{x}\left\lVert x-x_{n}\right\rVert_{X}^{2}$

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:

$\implies E(x_{*},y_{*})\leq E(x,y_{*})+2\tau_{x}\left\lVert
x-x_{*}\right\rVert_{X}^{2}\text{ for all }x$ (4.12)

This completes the proof for $x_{*}$ as

$|\partial_{x}E(x_{*},y_{*})|=\limsup_{x\to
x_{*}}\max\left(0,\frac{E(x_{*},y_{*})-E(x,y_{*})}{\left\lVert
x_{*}-x\right\rVert_{X}}\right)\leq\limsup 2\tau_{x}\left\lVert
x-x_{*}\right\rVert_{X}=0$

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,

$E(x_{*},y_{*})\leq\liminf E(x_{n+1},y_{n})\text{ by weak lower semi-
continuity.}$

However, on the right hand side we require:

$\lim E(x,y_{n})+2\tau_{x}\left\lVert x-x_{n}\right\rVert_{X}^{2}\leq
E(x,y_{*})+2\tau_{x}\left\lVert x-x_{*}\right\rVert_{X}^{2}$

In particular, we already require $\left\lVert x-x_{n}\right\rVert_{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 8 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.

$E(x_{*},y_{*})-2\tau_{x}\left\lVert x-x_{*}\right\rVert^{2}$$x$$E(x,y_{*})$
Figure 8:  Theorem 3.1 tells us that $(x_{*},y_{*})$ is a local critical point
but does not qualify the globality of the limit point. Equation (4.12) further
allows us to quantify the idea that if a lower energy critical point exists
then it must lie far from $(x_{*},y_{*})$. In particular, it must lie outside
of the shaded cone given by the supporting quadratic.

###  4.4 Proof of Theorem 4.1

###### Proof.

Fix arbitrary $(x_{0},y_{0})\in X\times Y$ and $\tau_{X},\tau_{Y}\geq 0$. Let
$(x_{n},y_{n})$ be defined as in Algorithm 1. By Lemma 4.3 we know that
$\\{(x_{n},y_{n})\operatorname{\;s.t.\;}n\in\mathds{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. ∎

##  5 Results

For numerical results we present two synthetic examples and one experimental
dataset from Transmission Electron Tomography. The two synthetic examples are
discretised at a resolution of $200\times 200$ then simulated using the X-ray
transform with a parallel beam geometry sampled at $1^{\circ}$ intervals over
a range of $60^{\circ}$ resulting in a full sinogram of size $287\times 180$
and sub-sampled data at $287\times 60$. Gaussian white noise (standard
deviation 5% maximum signal) is then added to give the data. The experimental
dataset was acquired with an annular dark field (parallel beam) Scanning TEM
modality from which we have 46 projections spaced uniformly in $3^{\circ}$
intervals over a range of $135^{\circ}$. Because of the geometry of the
acquisition, we can treat the original 3D dataset as a stack of 2D and thus
extract one of these slices as our example. This 2D dataset is then sub-
sampled to 29 projections over $87^{\circ}$, reducing the size from $173\times
45$ to $173\times 29$. This results in a reconstruction with $u$ of size
$120\times 120$ and $v$ of size $173\times 180$. A more detailed description
of the acquisition and sample properties of the experimental dataset can be
found in [48]. The code, and data, for all examples is available
111https://github.com/robtovey/2018_Directional_Inpainting_for_Limited_Angle
under the Creative Commons Attribution (CC BY) license.

###  5.1 Numerical Details

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.

$u_{0}=\operatorname*{argmin}_{u\geq 0}\frac{1}{2}\left\lVert
S\mathcal{R}u-b\right\rVert_{2}^{2}+\lambda\operatorname{TV}(u),\qquad
v_{0}=\mathcal{R}u_{0}$

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

$E(u,v)=\frac{1}{2}\left\lVert\mathcal{R}u-v\right\rVert^{2}_{\alpha_{1}}+\frac{\alpha_{2}}{2}\left\lVert
S\mathcal{R}u-b\right\rVert_{2}^{2}+\frac{\alpha_{3}}{2}\left\lVert
Sv-b\right\rVert_{2}^{2}+\beta_{1}\operatorname{TV}(u)+\beta_{2}\left\lVert
A_{Ru}\nabla v\right\rVert_{2,1}$ $\displaystyle A_{d}(x)$
$\displaystyle=c_{1}(x|\lambda_{1}-\lambda_{2},\lambda_{1}+\lambda_{2})\mathbf{e}_{1}(x)\mathbf{e}_{1}(x)^{T}+c_{2}(x|\lambda_{1}-\lambda_{2},\lambda_{1}+\lambda_{2})\mathbf{e}_{2}(x)\mathbf{e}_{2}(x)^{T}$
where $\displaystyle(\nabla d_{\rho}\nabla
d_{\rho}^{T})_{\sigma}=\lambda_{1}\mathbf{e}_{1}\mathbf{e}_{1}^{T}+\lambda_{2}\mathbf{e}_{2}\mathbf{e}_{2}^{T}\text{
is a pointwise eigenvalue decomposition}$ $\displaystyle
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))$

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:

  1. 1.

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
$10^{\pm^{1}\\!\\!/\\!_{2}}$.

  2. 2.

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}$.

  3. 3.

$\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.

  4. 4.

$\operatorname{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:

$\operatorname*{argmin}\frac{1}{2}\left\lVert
v-v_{0}\right\rVert_{2}^{2}+\beta_{2}\left\lVert A_{v_{0}}\nabla
v\right\rVert_{2,1}$

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
$\operatorname{DTV}$ prior is for each dataset.

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.

| $\alpha_{1}$ | $\alpha_{2}$ | $\alpha_{3}$ | $\beta_{1}$ | $\beta_{2}$ | $\beta_{3}$ | $\rho$ | $\sigma$
---|---|---|---|---|---|---|---|---
Concentric Rings Phantom | $\frac{1}{2^{2}}\mathds{1}_{\Omega^{\prime c}}$ | 1 | $1\times 10^{-1}$ | $3\times 10^{-5}$ | $3\times 10^{3}$ | $10^{10}$ | $1$ | $8$
Shepp-Logan Phantom | $\frac{1}{4^{2}}\mathds{1}_{\Omega^{\prime c}}$ | 1 | $3\times 10^{-1}$ | $3\times 10^{-5}$ | $3\times 10^{2}$ | $10^{10}$ | $1$ | $8$
Experimental Dataset (Both sampling ratios) | $\frac{1}{2^{2}}\mathds{1}_{\Omega^{\prime c}}$ | 1 | $3\times 10^{2}$ | $1\times 10^{-5}$ | $3\times 10^{1}$ | $10^{6}$ | $1$ | $0$
Table 1: Parameter choices for numerical experiments. Each algorithm was run
for 300 iterations

###  5.2 Canonical Synthetic Dataset

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.

Figure 9: Canonical synthetic example. Top row shows the reconstructions, $u$,
while the bottom row shows the reconstructed sinogram, $v$.

###  5.3 Non-Trivial Synthetic Dataset

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.

Figure 10: Non-trivial synthetic example of the modified Shepp-Logan phantom.
Top row shows the reconstructions, $u$, while the bottom row shows the
reconstructed sinogram, $v$. We regain the large-scale geometry of the shape
without losing much of the interior features.

###  5.4 Experimental Dataset

Figure 11: Raw data for Transmission Electron microscopy example. Projections
at large angles, e.g. $-68^{\circ}$, show the presence of the sample holder
which violates the X-ray modelling assumption that outside of the region of
interest is vacuum. If the violation is too extreme then this can cause strong
artefacts in reconstructions and so the common action is to discard such data.
The plane surface also violates this model but is relatively weak at low
angles and so will cause weaker artefacts. A source of noise in this
acquisition is that over time the surface becomes coated with carbon. This is
first visible as a thin film at $-2^{\circ}$ and progressively gets thicker
through the remaining projections. At $34^{\circ}$ we see a bump of carbon
appear on the top right edge. After pre-processing, we extract a 2D slice of
all projections to form the full range as shown top right artificially sub-
sample to compare TV with our proposed reconstruction method.

Figure 12: Reconstructions from a slice of the experimental data. We have
chosen the slice half-way down through the projections shown in Figure 11 to
coincide with one of the rounded corners. The arc artefact was an anticipated
consequence due to the out-of-view mass, the pre-processing has simply reduced
the intensity. Proposed reconstructions consistently show better homogeneity
inside the particle and sharper boundaries. The missing angles direction is
the bottom-left to top-right diagonal where we see most error in each
reconstruction, in particular, the blurring of the top right corner of the
particle is a limited angle artefact.

Figure 13: Comparison between each reconstruction after thresholding. The
geometrical properties of interest are that each edge should be linear, the
left hand corner is rounded and the remaining corners are not. The particle of
interest is homogeneous so thresholding the images should emphasise this in a
way which is very unsympathetic to blurred edges. Again, the top right corner
of each particle in the sub-sampled reconstructions coincides with the
exacerbated missing wedge direction and so we expect each reconstruction to
make some error here.

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 ${}^{3}\\!\\!/\\!_{4}$ 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.

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.

##  6 Conclusions and Outlook

In this paper we have presented a novel method for tomographic reconstructions
in a limited angle scenario along with a numerical algorithm with convergence
guarantees. We have also tested our approach on synthetic and experimental
data and shown consistent improvement over alternative reconstruction methods.
Even when the X-ray transform model is noticeably violated, as with our
experimental data, we still better recover boundaries of the reconstructed
sample.

There are three main directions which could be explored in future. Firstly, we
think there is great potential to apply our framework to other applications,
such as in tomographic imaging with occlusions and heavy metal artefacts where
the inpainting region is much smaller [22, 23]. Secondly, we would like to
find an alternative numerical algorithm with either faster practical
convergence or one which is more capable of avoiding local minima in this non-
convex landscape. Finally, we would like to explore the potential for an
alternative regularisation functional on the sinogram which is better able to
treat visible and invisible singularities, denoising and inpainting problems,
independently. At the moment, the TV prior alone can reconstruct visible
singularities well however, introducing a sinogram regulariser currently
improves on the invisible region at the expensive of damaging the visible.
Overall, we feel that this presents the natural progression for the current
work although it remains unclear how to regularise these invisible
singularities.

## Acknowledgements

RT would like to thank Ozan Öktem for valuable discussion on the application
of micro-local analysis in this setting. RT acknowledges support from the
EPSRC grant EP/L016516/1 for the University of Cambridge Centre for Doctoral
Training, the Cambridge Centre for Analysis. MB acknowledges support from the
Leverhulme Trust Early Career Fellowship Learning from mistakes: ‘A supervised
feedback-loop for imaging applications’, and the Isaac Newton Trust. CBS
acknowledges support from the Leverhulme Trust project on ‘Breaking the non-
convexity barrier’, EPSRC grant Nr. EP/M00483X/1, the EPSRC Centre Nr.
EP/N014588/1, the RISE projects CHiPS and NoMADS, and the Alan Turing
Institute. CBS, MB and RT also acknowledge the Cantab Capital Institute for
the Mathematics of Information. CB acknowledges support from the EU RISE
project NoMADS, the PIHC innovation fund of the Technical Medical Centre of UT
and the Dutch 4TU programme Precision Medicine. MJL acknowledges financial
support from the Netherlands Organization for Scientific Research (NWO),
project 639.073.506. SMC acknowledges support from the Henslow Research
Fellowship at Girton College, Cambridge. RKL acknowledges support from a Clare
College Junior Research Fellowship. PAM acknowledges EPSRC grant EP/R008779/1.

## References

  * [1] Matthias J. Ehrhardt, Kris Thielemans, Luis Pizarro, David Atkinson, Sébastien Ourselin, Brian F. Hutton, and Simon R. Arridge.  Joint reconstruction of PET-MRI by exploiting structural similarity.  Inverse Problems, 31(1):015001, 2015.
  * [2] Rowan K. Leary, Zineb Saghi, Paul A. Midgley, and Daniel J. Holland.  Compressed sensing electron tomography.  Ultramicroscopy, 131:70–91, aug 2013.
  * [3] Christian Kübel, Andreas Voigt, Remco Schoenmakers, Max Otten, David Su, Tan-Chen Lee, Anna Carlsson, and John Bradley.  Recent advances in electron tomography: TEM and HAADF-STEM tomography for materials science and semiconductor applications.  Microscopy and Microanalysis, 11(5):378–400, oct 2005.
  * [4] Gongpu Zhao, Juan R. Perilla, Ernest L. Yufenyuy, Xin Meng, Bo Chen, Jiying Ning, Jinwoo Ahn, Angela M. Gronenborn, Klaus Schulten, Christopher Aiken, and Peijun Zhang.  Mature HIV-1 capsid structure by cryo-electron microscopy and all-atom molecular dynamics.  Nature, 497(7451):643–646, may 2013.
  * [5] Willi A. Kalender.  X-ray computed tomography.  Physics in Medicine and Biology, 51(13):29–43, jul 2006.
  * [6] Veerle Cnudde and Matthieu Nicolaas Boone.  High-resolution X-ray computed tomography in geosciences: A review of the current technology and applications.  Earth-Science Reviews, 123:1–17, aug 2013.
  * [7] Noboru Kawase, Mitsuro Kato, Hideo Nishioka, and Hiroshi Jinnai.  Transmission electron microtomography without the “missing wedge” for quantitative structural analysis.  Ultramicroscopy, 107(1):8–15, jan 2007.
  * [8] Yiheng Zhang, Heang Ping Chan, Berkman Sahiner, Jun Wei, Mitchell M. Goodsitt, Lubomir M. Hadjiiski, Jun Ge, and Chuan Zhou.  A comparative study of limited-angle cone-beam reconstruction methods for breast tomosynthesis.  Medical Physics, 33(10):3781–3795, sep 2006.
  * [9] E T Quinto.  Singularities of the X-ray transform and limited data tomography in {R^2} and {R^3}.  SIAM Journal on Mathematical Analysis, 24(5):1215–1225, 1993.
  * [10] Jürgen Frikel and Eric Todd Quinto.  Characterization and reduction of artifacts in limited angle tomography.  Inverse Problems, 29(12):21, 2013.
  * [11] Alexander I. Katsevich.  Local tomography for the limited-angle problem.  Journal of Mathematical Analysis and Applications, 213(1):160–182, 1997.
  * [12] Lee A. Feldkamp, L. C. Davis, and James W. Kress.  Practical cone-beam algorithm.  Journal of the Optical Society of America A, 1(6):612, jun 1984\.
  * [13] Thomas Flohr, Karl Stierstorfer, Herbert Bruder, J. Simon, Arkadiusz Polacin, and Stefan Schaller.  Image reconstruction and image quality evaluation for a 16-slice CT scanner.  Medical Physics, 30(5):832–845, apr 2003.
  * [14] Xiangyang Tang, Jiang Hsieh, Roy A. Nilsen, Sandeep Dutta, Dmitry Samsonov, and Akira Hagiwara.  A three-dimensional-weighted cone beam filtered backprojection (CB-FBP) algorithm for image reconstruction in volumetric CT-helical scanning.  Physics in Medicine and Biology, 51(4):855–874, feb 2006.
  * [15] Peter Gilbert.  Iterative methods for the three-dimensional reconstruction of an object from projections.  Journal of Theoretical Biology, 36(1):105–117, jul 1972.
  * [16] Jose. I. Agulleiro, Eduardo. M. Garzón, Inmaculada García, and J. J. Fernández.  Vectorization with SIMD extensions speeds up reconstruction in electron tomography.  Journal of Structural Biology, 170(3):570–575, 2010.
  * [17] Martin Spitzbarth and Malte Drescher.  Simultaneous iterative reconstruction technique software for spectral-spatial EPR imaging.  Journal of Magnetic Resonance, 257:79–88, aug 2015.
  * [18] Bart Goris, Wouter Van den Broek, Kees Joost Batenburg, Hamed Heidari Mezerji, and Sara Bals.  Electron tomography based on a total variation minimization reconstruction technique.  Ultramicroscopy, 113:120–130, feb 2012.
  * [19] Zhiqiang Chen, Xin Jin, Liang Li, and Ge Wang.  A limited-angle CT reconstruction method based on anisotropic TV minimization.  Physics in medicine and biology, 58(7):2119–41, apr 2013.
  * [20] Jawook Gu and Jong Chul Ye.  Multi-Scale Wavelet Domain Residual Learning for Limited-Angle CT Reconstruction.  arXiv submission, 2017.
  * [21] Kerstin Hammernik, Tobias Würfl, Thomas Pock, and Andreas Maier.  A Deep Learning Architecture for Limited-Angle Computed Tomography Reconstruction.  In Bildverarbeitung für die Medizin, pages 92–97. Springer Vieweg, Berlin, Heidelberg, 2017.
  * [22] Harald Köstler, Michael Prümmer, Ulrich Rüde, and Joachim Hornegger.  Adaptive variational sinogram interpolation of sparsely sampled CT data.  In Proceedings - International Conference on Pattern Recognition, volume 3, pages 778–781. IEEE, 2006.
  * [23] Yi Zhang, Yi Fei Pu, Jin Rong Hu, Yan Liu, and Ji Liu Zhou.  A new CT metal artifacts reduction algorithm based on fractional-order sinogram inpainting.  Journal of X-Ray Science and Technology, 19(3):373–384, 2011.
  * [24] Si Li, Qing Cao, Yang Chen, Yining Hu, Limin Luo, and Christine Toumoulin.  Dictionary learning based sinogram inpainting for CT sparse reconstruction.  Optik, 125(12):2862–2867, 2014.
  * [25] Jianwei Gu, Li Zhang, Guoqiang Yu, Yuxiang Xing, and Zhiqiang Chen.  X-ray CT metal artifacts reduction through curvature based sinogram inpainting.  Journal of X-ray Science and Technology, 14(2):73–82, 2006.
  * [26] Hanming Zhang, Linyuan Wang, Yuping Duan, Lei Li, Guoen Hu, and Bin Yan.  Euler’s Elastica Strategy for Limited-Angle Computed Tomography Image Reconstruction.  IEEE Transactions on Nuclear Science, 64(8):2395–2405, 2017.
  * [27] Martti Kalke and Samuli Siltanen.  Sinogram Interpolation Method for Sparse-Angle Tomography.  Applied Mathematics, 05(03):423–441, 2014.
  * [28] Jean-Philippe Thirion.  A Geometric Alternative to Computed Tomography.  In Engineering in Medicine and Biology Society, volume 13, page 34, 1991.
  * [29] Martin Burger, Jahn Müller, Evangelos Papoutsellis, and Carola-Bibiane Schönlieb.  Total Variation Regularisation in Measurement and Image Space for PET reconstruction.  Inverse Problems, 30(10), oct 2014.
  * [30] Bin Dong, Jia Li, and Zuowei Shen.  X-ray CT image reconstruction via wavelet frame based regularization and Radon domain inpainting.  Journal of Scientific Computing, 54(2-3):333–349, feb 2013.
  * [31] Peter Ochs, Jalal M. Fadili, and Thomas Brox.  Non-smooth Non-convex Bregman Minimization: Unification and new Algorithms.  arXiv submission, 2017.
  * [32] Amir Beck and Marc Teboulle.  Fast gradient-based algorithms for constrained total variation image denoising and deblurring problems.  IEEE Transactions on Image Processing, 18(11):2419–2434, 2009.
  * [33] Tony F. Chan, Selim Esedoglu, and Mila Nikolova.  Algorithms for Finding Global Minimizers of Image Segmentation and Denoising Models.  SIAM Journal on Applied Mathematics, 66(5):1632–1648, 2006.
  * [34] Martin Benning, Michael Möller, Raz Z. Nossek, Martin Burger, Daniel Cremers, Guy Gilboa, and Carola-Bibiane Schönlieb.  Nonlinear spectral image fusion.  In Scale Space and Variational Methods in Computer Vision, volume LNCS 10302, pages 41–53, mar 2017.
  * [35] Joachim Weickert.  Anisotropic diffusion in image processing.  Image Rochester NY, 256(3):170, 1998.
  * [36] Benjamin Berkels, Martin Burger, Marc Droske, Oliver Nemitz, and Martin Rumpf.  Cartoon extraction based on anisotropic image classification.  In Vision, Modeling, and Visualization, page 293. IOS Press, 2006\.
  * [37] Virginia Estellers, Stefano Soatto, and Xavier Bresson.  Adaptive Regularization With the Structure Tensor.  IEEE Transactions on Image Processing, 24(6):1777–1790, jun 2015\.
  * [38] Marcelo Bertalmio, Guillermo Sapiro, Vincent Caselles, and Coloma Ballester.  Image inpainting.  In Proceedings of the 27th annual conference on Computer graphics and interactive techniques - SIGGRAPH ’00, pages 417–424, New York, New York, USA, 2000. ACM Press.
  * [39] Antonin Chambolle and Thomas Pock.  A first-order primal-dual algorithm for convex problems with applications to imaging.  Journal of Mathematical Imaging and Vision, 40(1):120–145, may 2011\.
  * [40] Mosek ApS.  The Mosek optimization software.  Online at http://www.mosek.com, 2010.
  * [41] Dmitriy Drusvyatskiy and Adrian S. Lewis.  Error bounds, quadratic growth, and linear convergence of proximal methods.  arXiv submission, page 23, 2016.
  * [42] Thomas Pock and Shoham Sabach.  Inertial Proximal Alternating Linearized Minimization (iPALM) for Nonconvex and Nonsmooth Problems.  SIAM Journal on Imaging Sciences, 9(4):1756–1787, jan 2016.
  * [43] Jingwei Liang, Jalal M. Fadili, and Gabriel Peyré.  A Multi-step Inertial Forward-Backward Splitting Method for Non-convex Optimization.  Advances in Neural Information Processing Systems, (29):1–9, 2016\.
  * [44] Peter Ochs, Yunjin Chen, Thomas Brox, and Thomas Pock.  iPiano: Inertial Proximal Algorithm for Non-Convex Optimization.  SIAM Journal on Imaging Sciences, 7(2):1388–1419, apr 2014.
  * [45] Jérôme Bolte, Shoham Sabach, and Marc Teboulle.  Proximal alternating linearized minimization for nonconvex and nonsmooth problems.  Mathematical Programming, 146(1-2):459–494, 2014.
  * [46] Dmitriy Drusvyatskiy, Alexander D. Ioffe, and Adrian S. Lewis.  Nonsmooth optimization using Taylor-like models: error bounds, convergence, and termination criteria.  arXiv submission, pages 1–21, 2016.
  * [47] John C. Duchi.  Introductory Lectures on Stochastic Population Systems.  arXiv submission, pages 1–84, may 2017.
  * [48] Sean M. Collins, Emilie Ringe, Martial Duchamp, Zineb Saghi, Rafal E. Dunin-Borkowski, and Paul A. Midgley.  Eigenmode Tomography of Surface Charge Oscillations of Plasmonic Nanoparticles by Electron Energy Loss Spectroscopy.  ACS Photonics, 2(11):1628–1635, nov 2015.
  * [49] Michael Grant, Stephen Boyd, and Yinyu Ye.  CVX: Matlab software for disciplined convex programming, 2008.
  * [50] Duccio Fanelli and Ozan Öktem.  Electron tomography: A short overview with an emphasis on the absorption potential model for the forward problem.  Inverse Problems, 24(1):013001, feb 2008.
  * [51] Antonin Chambolle and Pierre-Louis Lions.  Image recovery via total variation minimization and related problems.  Numerische Mathematik, 76(2):167–188, 1997.

##  Appendix A Theorem 2.1

######  Theorem.

If

  1. 1.

$c_{i}$ are $2k$ times continuously differentiable in $\Delta$ and $\Sigma$,
$k\geq 1$

  2. 2.

$c_{1}(x|0,\Sigma)=c_{2}(x|0,\Sigma)$ for all $x$ and $\Sigma\geq 0$

  3. 3.

$\partial_{\Delta}^{2j-1}c_{1}(x|0,\Sigma)=\partial_{\Delta}^{2j-1}c_{2}(x|0,\Sigma)=0$
for all $x$ and $\Sigma\geq 0,j=1\ldots,k$

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

In this proof we will drop the $x$ argument from $c_{i}$ for conciseness of
notation. Define

$M_{d}=(\nabla d_{\rho}\nabla d_{\rho}^{T})_{\sigma}$

Note that convolutions are bounded linear maps and $\nabla d_{\rho}\in L^{2}$
by Young’s inequality hence $M\colon L^{1}(\mathds{R}^{2},\mathds{R})\to
L^{\infty}(\mathds{R}^{2},\operatorname{Sym}_{2})$ is well defined and
differentiable w.r.t. $d$. Hence, it suffices to show that $A$ is
differentiable w.r.t. $M$ where

$M_{d}=\lambda_{1}e_{1}e_{1}^{T}+\lambda_{2}e_{2}e_{2}^{T},\lambda_{1}\geq\lambda_{2}\implies
A=c_{1}(\Delta,\Sigma)e_{1}e_{1}^{T}+c_{2}(\Delta,\Sigma)e_{2}e_{2}^{T}$

where $\Delta=\lambda_{1}-\lambda_{2},\Sigma=\lambda_{1}+\lambda_{2}$. Note
that this is not a trivial statement, we can envisage very simple cases in
which the (ordered) eigenvalue decomposition is not even continuous. For
instance

$M(t)=\begin{pmatrix}1-t&0\\\
0&t\end{pmatrix}\implies\Sigma(t)=1,\Delta(t)=|1-2t|,e_{1}=\left\\{\begin{array}[]{cc}(1,0)^{T}&t<^{1}\\!\\!/\\!_{2}\\\
(0,1)^{T}&t>^{1}\\!\\!/\\!_{2}\end{array}\right.$

Hence we can see that despite having $M\in C^{\infty}$ we cannot even
guarantee that the decomposition is continuous and so cannot employ any chain
rule to say that $A$ is smooth.

The structure of this proof breaks into 4 parts:

  1. 1.

If $c_{1}(0,\Sigma)=c_{2}(0,\Sigma)$ then $A$ is well defined

  2. 2.

If $c_{i}\in C^{2}$ for some open neighbourhood of point $x$ such that
$\lambda_{1}(x)>\lambda_{2}(x)$ then $A$ is differentiable w.r.t. $M$ on an
open neighbourhood of $x$

  3. 3.

If $\partial_{\Delta}c_{1}(0,\Sigma)=\partial_{\Delta}c_{2}(0,\Sigma)=0$ at a
point, $x$, where $\lambda_{1}(x)=\lambda_{2}(x)$ then $A$ is differentiable
on an open neighbourhood of $x$

  4. 4.

If
$\partial_{\Delta}^{2j-1}c_{1}(0,\Sigma)=\partial_{\Delta}^{2j-1}c_{2}(0,\Sigma)=0$
at a point $x$ where $\lambda_{1}(x)=\lambda_{2}(x)$ and for all $j=1\ldots,k$
then $A$ is $C^{2k-1}$ on an open neighbourhood of $x$

###### Proof.

Proof of part i: Note that when $\lambda_{1}=\lambda_{2}$ the choice of
$e_{i}$ is non-unique subject to
$e_{1}e_{1}^{T}+e_{2}e_{2}^{T}=\operatorname{id}$ and so we get

$A=c_{1}(0,\Sigma)\operatorname{id}+(c_{2}(0,\Sigma)-c_{1}(0,\Sigma))e_{2}e_{2}^{T}$

Hence $A$ is well defined if and only if $c_{1}(0,\Sigma)=c_{2}(0,\Sigma)$ for
all $\Sigma\geq 0$.

As we are decomposing $2\times 2$ matrices it will simplify the proof to write
explicit forms for the values under consideration.

$M=\begin{pmatrix}M_{11}&M_{12}\\\
M_{12}&M_{22}\end{pmatrix}\implies\lambda_{i}=\frac{M_{11}+M_{22}\pm\sqrt{(M_{11}-M_{22})^{2}+4M_{12}^{2}}}{2}$
$\Sigma=M_{11}+M_{22},\quad\Delta=\sqrt{(M_{11}-M_{22})^{2}+4M_{12}^{2}}$
$\Delta\neq 0\implies e_{1}=\frac{(2M_{12},\Delta-
M_{11}+M_{22})^{T}}{\sqrt{(\Delta-
M_{11}+M_{22})^{2}+4M_{12}^{2}}}=\frac{(\Delta+M_{11}-M_{22},2M_{12})^{T}}{\sqrt{(\Delta+M_{11}-M_{22})^{2}+4M_{12}^{2}}},$
$\quad e_{2}=\begin{pmatrix}0&-1\\\ 1&0\end{pmatrix}e_{1}$

We give two equations for $e_{1}$ to account for the case when we get
$\frac{(0,0)^{T}}{0}$.

Proof of part ii: Note that $\Sigma$ is always smooth and $\Delta$ is smooth
on the set $\\{\Delta>0\\}$
Case $M_{12}(x)\neq 0$: Now both equations of $e_{1}$ are valid (and equal)
and the denominators non-zero on a neighbourhood. Hence, we can apply the
standard chain rule and product rule to conclude.
Case $M_{12}(x)=0$: In this case $M(x)$ is diagonal but as
$\Delta=|M_{11}-M_{22}|>0$ we know that one of our formulae for $e_{1}$ must
be valid with non-vanishing denominator in a neighbourhood and so we can
conclude as in the first case.

Proof of part iii: Given the Neumann condition on $c_{i}$ we shall express
$c_{i}$ locally by Taylor’s theorem.

$c_{i}(\Delta,\Sigma)=c_{i}(0,\Sigma)+O(\Delta^{2})=c_{1}(0,\Sigma)+O(\Delta^{2})$

Now we consider a perturbation:

$M=\begin{pmatrix}m&0\\\
0&m\end{pmatrix},\quad\varepsilon=\begin{pmatrix}\varepsilon_{11}&\varepsilon_{12}\\\
\varepsilon_{12}&\varepsilon_{22}\end{pmatrix}$ $\implies
A(M+\varepsilon)-A(M)=(c_{1}(0,2m+\varepsilon_{11}+\varepsilon_{22})-c_{1}(0,2m))\operatorname{id}+O(\Delta^{2})$
$\Delta^{2}=(\varepsilon_{11}-\varepsilon_{22})^{2}+4\varepsilon_{12}^{2}=O(\left\lVert\varepsilon\right\rVert^{2})\implies
O(\Delta^{2})\leq O(\left\lVert\varepsilon\right\rVert^{2})$
$\implies\frac{A(M+\varepsilon)-A(M)}{\left\lVert\varepsilon\right\rVert}=\frac{\partial_{\Sigma}c_{1}(0,2m)\operatorname{tr}(\varepsilon)}{\left\lVert\varepsilon\right\rVert}+O(\left\lVert\varepsilon\right\rVert)$

In particular, $A$ is differentiable w.r.t. $M$ at $x$.

Proof of part iv: The proof of this follows exactly as the previous part,

$c_{i}(\Delta,\Sigma)=\sum_{0}^{k-1}\frac{\Delta^{2j}}{j!}\partial_{\Delta}^{2j}c_{i}(0,\Sigma)+O(\Delta^{2k})$

where the remainder term is sufficiently smooth. Hence $c_{i}$ is at least
$2k-1$ times differentiable w.r.t. $M$ ∎

##  Appendix B Theorem 3.1

######  Theorem.

If

  * •

$c_{i}$ are bounded away from 0

  * •

$\rho>0$

  * •

$A_{d}$ is differentiable in $d$

then sub-level sets of $E$ are weakly compact in
$L^{2}(\Omega,\mathds{R})\times L^{2}(\mathds{R}^{2},\mathds{R})$ and $E$ is
weakly lower semi-continuous. i.e. for all $(u_{n},v_{n})\in
L^{2}(\Omega,\mathds{R})\times L^{2}(\mathds{R}^{2},\mathds{R})$:

$E(u_{n},v_{n})\text{ uniformly bounded implies a subsequence converges
weakly}$ $\liminf E(u_{n},v_{n})\geq E(u,v)\text{ whenever
}u_{n}\rightharpoonup u,v_{n}\rightharpoonup v$

###### Proof.

If $c_{i}$ are bounded away from 0 then in particular we have
$A_{\mathcal{R}u_{n}}\gtrsim 1$ so
$\operatorname{DTV}_{u_{n}}(v_{n})=\left\lVert A_{\mathcal{R}u_{n}}\nabla
v_{n}\right\rVert\gtrsim\left\lVert\nabla
v_{n}\right\rVert=\operatorname{TV}(v_{n})$. Hence,

$\begin{split}E(u_{n},v_{n})\text{ uniformly bounded }&\implies\\\
&\left\lVert S_{\Omega^{\prime
c}}(\mathcal{R}u_{n}-v_{n})\right\rVert_{2}^{2}+\left\lVert
S_{\Omega^{\prime}}\mathcal{R}u_{n}-b\right\rVert_{2}^{2}+\left\lVert
S_{\Omega^{\prime}}v_{n}-b\right\rVert_{2}^{2}\\\
&\qquad+\operatorname{TV}(u_{n})+\operatorname{TV}(v_{n})\text{ uniformly
bounded}\\\ &\implies\left\lVert
A(u,v)^{T}-b\right\rVert_{2}^{2}+\operatorname{TV}\left((u,v)\right)\text{
uniformly bounded}\end{split}$

for some linear $A$ and constant $b$. Thus we are in a very classical setting
where weak compactness can be shown in both the
$\left\lVert(u,v)\right\rVert_{2}$ norm and
$\left\lVert(u,v)\right\rVert_{1}+\operatorname{TV}((u,v))$ [51].

We now proceed to the second claim, that $E$ is weakly lower semi-continuous.
Note that all of the convex terms in our energy are lower semi-continuous by
classical arguments so it remains to show that the non-convex term is too.
i.e.

$(u_{n},v_{n})\rightharpoonup(u,v)\stackrel{{\scriptstyle?}}{{\implies}}\liminf\left\lVert
A_{\mathcal{R}u_{n}}\nabla v_{n}\right\rVert_{2,1}\geq\left\lVert
A_{\mathcal{R}u}\nabla v\right\rVert_{2,1}$

We shall first present a lemma from distribution theory.

######  Lemma B.1.

If $\Omega$ is compact, $\varphi\in C^{\infty}(\Omega,\mathds{R})$ and
$w_{n}\stackrel{{\scriptstyle L^{p}}}{{\rightharpoonup}}w$ then

$w_{n}\star\varphi\to w\star\varphi\text{ convergence in
}C^{k}(\Omega,\mathds{R})\text{ for all }k<\infty$

###### Proof.

Recall that $w_{n}\rightharpoonup w\implies\left\lVert
w_{n}\right\rVert_{p}\leq W$ for some $W<\infty$. By Hölder’s inequality:

$|w_{n}\star\varphi(x)-w\star\varphi(y)|\leq\int|w_{n}(z)||\varphi(x-z)-\varphi(y-z)|\lesssim_{p,\Omega}|x-y|W\left\lVert\varphi^{\prime}\right\rVert_{\infty}$
$|w_{n}\star\varphi(x)|\lesssim_{p,\Omega}W\left\lVert\varphi\right\rVert_{\infty}$

i.e. $\\{w_{n}\operatorname{\;s.t.\;}n\in\mathds{N}\\}$ is uniformly bounded
and uniformly (Lipschitz) continuous.

$w_{n}\rightharpoonup w\implies
w_{n}\star\varphi(x)-w\star\varphi(x)=\left\langle
w_{n}-w,\varphi(x-\cdot)\right\rangle\to 0\text{ for every }x$

Hence, we also know $w_{n}\star\varphi$ converges point-wise to a unique
continuous function. Suppose $\left\lVert
w_{n_{k}}\star\varphi-w\star\varphi\right\rVert_{\infty}\geq\varepsilon>0$ for
some $\varepsilon$ and sub-sequence $n_{k}\to\infty$. By the Arzela-Ascoli
theorem some further subsequence must converge uniformly to the point-wise
limit, $w\star\varphi$, which gives us the required contradiction. Hence,
$w_{n}\star\varphi\to w\star\varphi$ in $L^{\infty}=C^{0}$. The general
theorem follows by induction. ∎

This lemma gives us two direct results.

$\rho>0\implies(Ru_{n})_{\rho}\to(Ru)_{\rho}\text{ in }L^{\infty}$
$\\{(Ru_{n})_{\rho}\\}\operatorname{\cup}\\{(Ru)_{\rho}\\}\text{ compact,
}A_{d}\in C^{1}(d)\implies A_{Ru_{n}}\to A_{Ru}\text{ in
}\left\lVert\cdot\right\rVert_{2,\infty}$

Hence, whenever $w\in W^{1,1}$ we have

$\displaystyle\left\lVert A_{\mathcal{R}u_{n}}\nabla w\right\rVert$
$\displaystyle\geq\left\lVert A_{\mathcal{R}u}\nabla
w\right\rVert-\left\lVert(A_{\mathcal{R}u_{n}}-A_{\mathcal{R}u})\nabla
w\right\rVert$ $\displaystyle\geq\left\lVert A_{\mathcal{R}u}\nabla
w\right\rVert-\left\lVert
A_{\mathcal{R}u_{n}}-A_{\mathcal{R}u}\right\rVert_{2,\infty}\operatorname{TV}(w)$

By density of $W^{1,1}$ in the space of Bounded Variation, we know the same
holds for $w=v_{n}$. We also know $\operatorname{TV}(v_{n})$ is uniformly
bounded thus

$\liminf\left\lVert A_{\mathcal{R}u_{n}}\nabla
v_{n}\right\rVert=\liminf\left\lVert A_{\mathcal{R}u}\nabla v_{n}\right\rVert$

Hence, for all $\left\lVert\varphi\right\rVert_{2,\infty}\leq 1$ we have

$\left\langle
v,\nabla\cdot(A_{\mathcal{R}u}\varphi)\right\rangle=\liminf_{n}\left\langle
v_{n},\nabla\cdot(A_{\mathcal{R}u}\varphi)\right\rangle\leq\liminf\left\lVert
A_{\mathcal{R}u}\nabla v_{n}\right\rVert\leq\liminf\left\lVert
A_{\mathcal{R}u_{n}}\nabla v_{n}\right\rVert$

as required. ∎