Directional sinogram inpainting for limited angle tomography

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.


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.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 common model in X-ray tomography is that we assume the map from physical samples to observed data, its sinogram, is the Xray transform which will be formally defined in (2.1).This link allows us easily to adapt tools which were developed for one application and translate them to other problems which also inherit structure from this sinogram data (see Table 1).Figure 1 gives a graphical representation of the X-ray transform in relation to the main application of this paper, limited angle tomography.This scenario is very common in imaging techniques such as TEM where unavoidable hardware limitations force strong sampling constraints on the system.In particular, we shall say that the ideal fully sampled inverse problem that we face is: Given data b find u such that Ru = b where R is the X-ray transform as will be defined in (2.1).Taking into account sampling constraints and the presence of noise we shall relax this to: Given data b find optimal pair (u, v) such that Ru ≈ v, Sv ≈ b where S represents the sub-sampling pattern which gives us the flexibility to address each situation described in Table 1.Presenting the problem in this format makes it apparent that we have two (potentially) ill-posed problems, a reconstruction problem and an inpainting problem.Trivially, there are three possible approaches to solving this problem.Constraining v = Ru, solving the inpainting and reconstruction independently, and solving the inpainting and reconstruction problem jointly.The core principle that we follow here is that the flexibility of the joint approach can allow for noticeable improvements when S is particularly 'bad'.Our main theoretical contribution shall be to present a broad joint framework for combining prior knowledge on two representations of the same sample.i.e. our original sample is an infinite dimensional object with a discretised density representation, u, and a discretised sinogram, v, with the consistency condition Ru ≈ v.In the application which we shall focus on, the limited angle problem, the inpainting problem is very ill-posed as there is a large contiguous range of data which cannot be sampled [7,8].Figures 1 and 2 respectively demonstrate what this problem looks like in the data space and why this is common in practical situations.We aim to demonstrate that by leveraging structure jointly in u and v we can minimise the impact of a severe limited angle scenario.

High Noise
In applications such as PET and Single Particle Analysis there is a full range of data available but there is a large amount of noise.[

9, 10] Few Angles
There are many reasons to restrict to a few angles, either for faster acquisition or to protect the sample from long exposure to radiation.
In such cases there are few measurements available but they have a low level of noise and no constraints on which view points are observed.Reconstruction artefacts are dominated by the low number of measurements.[7]

Occlusions
In some applications there can be multiple scales of signal intensity within the same image (possibly infinite intensity), such an example is Heavy Metal artefacts in medical X-ray CT scans.Untreated, this can cause severe artefacts in the smaller scales after reconstruction.The typical treatment is to cut all measurements of the 'heavy' object from the sinogram and inpaint before reconstruction.[11,12] Limited Angle If there are physical constraints on the acquisition process then there will often be a limited range of viewpoints which can be sampled.Reconstruction artefacts are characteristically dominated by the width of the contiguous missing data and so inpainting is required during the reconstruction process.(Figures 3, 4) [8,13] Table 1.Overview of signature reconstruction problems in X-ray Tomography.

Sample partially hidden
Figure 2. Demonstration of the physical constraints which lead to limited angle imaging.The key issue is that some materials block radiation almost completely, sending signal to noise ratio to 0. Typically, this is not a problem at small angles but will become unavoidable past a critical angle range.Left: Beam strength is constantly attenuated going through the sample.At high angles this attenuation becomes more extreme until the beam cannot penetrate the sample.This is particularly the case in Cryo-EM where samples must be contained in a frozen block [14].Right: The sample holders are typically impenetrable to the beam, for instance in Electron Microscopy.At high angles these will shield parts of the sample which can no longer be sampled.

Context and Proposed Model
The traditional methods for X-ray Tomography reconstruction attempt to solve SRu ≈ f using prior knowledge of only u or v.There are three main methods which fit into this category: • Filtered back projection (FBP) is a linear inversion method which takes advantage of the fact that R * R is a known convolution operator.Thus, FBP combines this deconvolution operator with smoothing on the sinogram data to account for noise.[15][16][17] • The Simultaneous Iterative Reconstruction Technique (SIRT) can be thought of as a preconditioned gradient descent on the function SRu − b 2 2 .Regularisation is then typically implemented by an early-stopping technique.[3,[18][19][20] • 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 method in this field is Total Variational (TV) regularisation where u is encouraged to have a piecewise constant intensity.This will be introduced formally in Section 2. [2,13,21] 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 3 gives examples where TV can vastly outperform FBP and SIRT when either the noise level is large or the angular range small.However,   Recently, this classical group of methods have received a revival through machine learning methods for reconstruction [22,23].In both of these examples the inpainting problem remains implicitly treated and the main reduction of artefacts is performed post-inversion.
The most common method that has been used to reconstruct pairs (u, v) is to solve each inverse problem sequentially.In general we can express the pipeline for such methods as: This has seen a lot of success in heavy metal artefact reduction where a regularisation functional for the inpainting problem may be constructed from dictionary learning [24], fractional order TV [12] and Euler's Elastica [25].Euler's Elastica has also been used in the limited angle problem [26] along with more customised interpolation methods [7].This approach allows 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 prior knowledge and consistency terms into a single energy functional such that our optimal pair of reconstructions will minimise this joint functional which we shall write as: 1) where α i ≥ 0 are weighting parameters, d i are some 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 and r 1 shall be chosen such that our reconstruction becomes a generalisation of the TV reconstruction.Recent work [9] has shown that such a joint approach can be advantageous in similar situations but closest to our approach is that of [27] where r 1 and r 2 were chosen to encode wavelet sparsity in both u and v.Where we hope to improve is in our choice for r 2 .It has previously been suggested that enhancing lines in a sinogram can also enhance edges in the density reconstruction, u [28].This aligns with our observations in Figure 4.In both examples shown, the reconstructions show characteristic blurring at the top and bottom edge which is seen in the sinogram domain as a loss of structure outside of the given angular range.The idea that lines starting either side of the inpainting region should continue through inpainting region is exactly the assertion that samples should have some rotational symmetries or consistencies.The archetypal example of this is that of concentric rings; it is perfectly radially symmetric and thus its sinogram is a vertical stack of constant horizontal stripes.The Shepp-Logan phantom is more complex but all of the edges are constantly sharp around each boundary level set rather than becoming more or less blurred at particular angles.This type of symmetry can still be seen in the sinogram, for instance the sharp yellow curve at the bottom the sinogram which traces the outer yellow ellipse in the phantom.The theoretical observations of [28] with these visual findings has led us to choose r 2 to encourage sharp line structures in v.The exact form of this will be formalised in Section 2. We shall demonstrate that our joint model, (1.1), has two main advantages over the other methods mentioned here.The first is that by including these weights α i we can enforce a more adaptive consistency between each u, v and b.For instance, if we let α 2 , α 4 → ∞ then we recover the TV reconstruction model.Alternatively, if we let α 3 , α 5 → ∞ then we recover a method which performs a single inpainting step and then calculates u from v and b as in [12,25,26].Another advantage of our proposed model is that it can easily be adapted to any of the problems described in Table 1.In limited angle tomography we have a particular sampling pattern, S, and particular types of samples which we encode with r i .In other applications one could simply substitute these components to customise the model and the numerical scheme that we present shall be valid for a broad class of non-convex and non-smooth functionals.

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, even when these problems become non-smooth and non-convex.
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. [29] 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 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, Ω ⊂ R n , by where In most applications we have n = 3 although in many cases (e.g.EM, X-ray CT) the scanners proceed slice-wise and so we have Consequently, we shall provide our numerical examples in 2D and this can be extended to 3D by a trivial extension of our proposed model.

Total Variation Regularisation
Total Variation (TV) regularisation is extremely common across many fields of image processing [30][31][32].The definition of the (isotropic) TV semi-norm on domains Ω ⊂ R n is stated as The intuition behind this is that when g ∈ W 1,1 (Ω, R) then we have exactly From this point forward we shall write the • 2,1 representation where |∇g| is to be understood as a measure if necessary.We shall denote the space of Bounded Variation as The most important fact that we shall use about the space of Bounded Variation is BV ⊂ L 2 ⊂ L 1 which holds whenever Ω is compact by the Poincaré inequality: 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,21].
The parameter λ is a regularisation parameter, which allows to enforce more or less regularity, depending on the quality of the data b.

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 [33] 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 [34][35][36].The heuristic for our regularisation on the sinogram was described in Figure 4 and we shall encode it in an anisotropic TV penalty which shall be written as 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 Examples of this are presented in Figure 7.Note that the choice c 1 = c 2 recovers the traditional TV regulariser but for |c 1 | << c 2 we allow for much larger (sparse) gradients in the direction of e 1 .This allows for large jumps in the direction of e 1 whilst maintaining flatness in the direction of e 2 .The question which remains is how to choose c i and e i for any given situation.In the case of denoising we can extract this information directly from the image.Following the construction of Weickert [33], given a noisy image, d, we construct It will be convenient to view this from a more natural parametrisation of λ 1 , λ 2 .
From this point on we shall drop the explicit x dependence and reparameterise, with a slight abuse of notation, This choice is natural as Σ will turn out to be a useful indicator for detecting regions with edges whereas ∆ will indicate edges with clear direction, see Figure 6.Typical examples of c 1 , c 2 include The key idea here is that c 1 << c 2 near lines but 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 → A d , now characterised uniquely by our choice of c i .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.
• Property (i) is necessary for A d to be well defined and continuous for all fixed d

holds trivially
The proof of this theorem is contained in Appendix A.

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, R, is the X-ray transform from (2.1) • The desired reconstructed sample, u ∈ BV(Ω, R) on some domain Ω We extend such that b| Ω c = 0 for notational convenience.
We also define S = S Ω to be the indicator function on Ω .The joint model that we thus propose is where α i , β i > 0 are weighting parameters.α 1 is embedded in the norm as it is a spatially varying weight, taking different values inside and outside of Ω .We note that the norms involving b are determined by our noise model of the acquisition process, in this case Gaussian noise.The final metric pairing Ru 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]).Isotropic TV penalty is chosen for a spatial prior as it represents the current state of the art method.Heuristically, it encourages piece-wise constant intensities which is also a natural assumption in many situations.The motivation for an anisotropic regularisation functional on the sinogram can be seen in Figure 4 where we see that reconstructions fail to continue sharp edge structures into the inpainting region of the sinogram.Finally, we must choose the directional component for DTV.If we do not know this beforehand then it must be estimated along with the reconstruction.We have decided to incorporate this into our directional regulariser to be of the form 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.
This theorem justifies numerical minimisation of E because it tells us that all (locally) descending sequences ( ) have a convergent subsequence and any limit point must minimise E over the original sequence.

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 : X × Y → 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: Sub-level sets of E are weakly compact (4.1) Note that if g is a norm then L x can be chosen to be the standard Lipschitz factor of ∇ x J.If J is twice differentiable then these constants must be finite.In our case we have Finiteness of ∇ 2 A 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. Classical alternating proximal descent would take but because of the complexity of A Ru each sub-problem would have the same complexity as the full problem, making it infeasible.The classical solution to this would have been to find a smooth approximation of • 2,1 with a small maximum error, ε.This now decomposes E(x, y n ) into the sum of a differentiable function and a convex function for which there exist many minimisation algorithms, e.g.[30].There are two problems with this approach, the original energy functional has been changed and the new energy is illconditioned with Lipschitz gradient O( |∂ d A| / ε ).By linearising A d inside • 2,1 we overcome both of these problems.Our energy has not changed and the x n+1 sub-problem is now a convex optimisation problem with O(|∂ d A|) conditioning which is much better.This second approach puts us in the framework of the ProxDescent algorithm [29,37] which we shall extend to cover alternating descent.In this setting we have the equivalent Algorithm 1 Input: Any x 0 ∈ X, τ x = L x + τ X for some τ X ≥ 0 and equivalent for y. n ← 0 while not converged do n ← n + 1 end while Output: (x n , y n ) convergence guarantees as for the original ProxDescent algorithm.With this algorithm our statement of convergence is the following theorem.• A subsequence of (x n , y n ) must converge weakly in X × Y • If {(x n , y n ) s. t. n = 1, . ..} 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 Lemma 4.6) E in both the direction of x and y.
This result is the parallel of Lemma 10, Theorem 18 and Corollary 21 in [29] 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:

Sketch Proof
The proof of Theorem 4.1 will be a consequence of two lemmas.
• (Lemma 4.3) For τ x , τ y sufficiently large, the sequence E n,n is monotonically decreasing and sequences x n − x n−1 X , y n − y n−1 Y converge to 0. This follows by a relatively standard sufficient decrease argument as seen in [29,38,39].
• (Lemma 4.6) First we shall 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 1 shows a uniformly convex energy E which shows this to be sharp.The common technique for overcoming this is assuming differentiability in the terms including both x and y [38,40,41].The necessity of our algorithm is that the non-convex term is also nondifferentiable and so it is unclear what would be required to make this statement stronger.
Example 4.2.Define E(x, y) = max(x, y) + x 2 + y 2 for x, y ∈ 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 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 ).

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 Proof.Note by Equations 4.6, 4.7 (definition of our sequence) we have Further, by application of the triangle inequality for g and the mean value theorem we have By equivalent argument, (4.11) Summing Equations 4.8-4.11gives 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:

Convergence to Critical Points
First we follow the work by Drusvyatskiy et.al. [42] we shall define criticality in terms of the, so called, slope of a function.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 However, at a point of non-differentiability with |∂F | = 0 we can still say a little more than local criticality (a priori a minimum/maximum/saddle point).Due to the lim sup in the definition we can in fact rule out some of these situations in a way that biases us towards minima, this is why we refer to this as a first order minimum rather than simply a first order critical point.
Figure 8. Examples of functions where 0 is/isn't a first order minimum in terms of the slope.The top row show characteristic examples of functions which are not smooth at 0. 0 is a first order minimum iff 0 is a minimum with respect to both left and right limits.In any finite dimension this generalises to minimum with respect to every one-sided direction.The bottom row shows the sharpness of statement when the function is locally some power of x on an interval of x > 0, say F (x) = cx 1+ε .If c > 0 then this is always a first order minimum.If c is negative then 0 is a first order minimum iff ε > 0. Note that for c < 0, 0 is never a local minimum of F but for ε > 0 the decrease is slower than linear and so this is a first order minimum.Hence, 0 is not a first order minimum of F which fits with an intuitive understanding of the term.In fact, due to the lim sup, if any directional derivative is negative then we are not at a first order minimum.More examples of this are shown in Figure 8.
Now we shall show that our iterative sequence converges to a first order minimum in this sense.Lemma 4.6.If (x n , y n ) are as given above and X, Y have finite dimension then every limit point of (x n , y n ), e.g.(x * , y * ), is a first order minimum of E in each coordinate direction.i.e. • Our assumption of finite dimensionality accounts for what is referred to as 'Assumption 3' in [29] 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 at all.The difficulty of showing convergence to critical points is common across both convex [43] and non-convex [29,38] optimisation.
Proof.First we recall that, in finite dimensional spaces, convex functions are continuous on the relative interior of their domain [44].Also note that by our definition of τ x , for all x, x , y we have where existence of such ξ is given by the Mean Value theorem.Hence, for all x we have 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: X for all x (4.12) This completes the proof for x * as The proof for y * follows by symmetry.Note that 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 ) (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 * ) ≤ lim inf E(x n+1 , y n ) by weak lower semi-continuity.However, on the right hand side we require In particular, note that we already require x − x n X to be weakly upper semicontinuous.Topologically, this is the statement that weak and norm convergence are equivalent which will not be the case in most practical (infinite dimensional) examples.
. In general we cannot guarantee that limit points are local minima but we can quantify the speed at which the energy may decay in a neighbourhood.In the case of quadratic penalty functions we have an explicit formula for a supporting cone which bounds the energy from below.

Proof of Theorem 4.1
Proof.Fix arbitrary (x 0 , y 0 ) ∈ X × Y and τ X , τ Y ≥ 0. Let (x n , y n ) be defined as in Algorithm 1.By Lemma 4.3 we know that {(x n , y n ) s. t. n ∈ 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.

Results
For numerical results we present two synthetic examples and one experimental dataset from Transmission Electron Tomography.The two synthetic examples are simulated by the X-ray transform with additive Gaussian white noise (standard deviation 5% maximum signal).Angles are sampled at 1 • intervals over a range of 60 • , i.e.only 1 / 3 of the full data is sampled.The experimental dataset was acquired with an annular dark field Scanning TEM modality from which we have 46 projections spaced uniformly in 3 • intervals over a range of 135 • , these are sub-sampled to 30 projections over 87 • .A more detailed description of the acquisition and sample properties can be found in [45].The code for the following examples is available ‡ under the Creative Commons Attribution (CC BY) license.

Numerical Details
All numerics are implemented in MATLAB 2016b.The sub-problem for u is solved with a PDHG algorithm [43] while the sub-problem for v is solved using the MOSEK solver via CVX [46,47], the step size τ is adaptively calculated.The initial point of our algorithm is always chosen to be a good TV reconstruction, i.e.
Our construction for the directional regularisation was chosen to be as the property ∂ ∆ c 1 is maximally negative for Σ, ∆ ≈ 0 which encourages ∆ to grow on regions which are near flat.In particular, this encourages edge continuation in the sinogram when they would otherwise be blurred out.Our synthetic experiments suggested that the parameter choices were relatively intuitive at a fixed resolution and scaling (both synthetic images satisfied u ∞ = 1) although switching to the synthetic dataset (different resolution, scaling, angular distribution) the optimal parameter values shifted significantly.Despite the large number of tuning parameters, we observed that exact tuning was unnecessary and so most were tuned within a factor of 10 ± 1 / 2 .The final parameter choices are found in Table 2 and a further comparison of different parameter choices can be seen in the supplementary material.

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 10, 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.

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

Experimental Dataset
The sample is a silver bipyramidal crystal placed on a planar surface, and the challenges of this dataset are shown in Figure 12.We immediately see that the wide angle projections have large artefacts from the planar surface which masks the sample at these high angles.This suggests that a more limited angle acquisition is desirable.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 these two factors is that we find a pressure to use both a limited angular range and the 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, 30 projections over an 87 • interval, with a bias towards earlier projections.The final artefact in the data is that there is mass seen in some of the projections which cannot be represented within the reconstruction volume.This breaks the X-ray model assumption that the field of view is surrounded by a vacuum and thus will also lead to reconstruction artefacts.To minimise this problem we shall perform preliminary background subtraction which will minimise the total mass contributing to this violation.Effectively, we pretend our detector was narrower such that we have data only in a small neighbourhood of the region of interest, this can be seen as the sharp horizontal cut-off in the pre-processed sinograms seen on the right of Figure 12.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 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 , contain views of the supporting grid as was described in Figure 2 suggesting that a more limited angle acquisition is desirable.The plane surface will also cause artefacts in the reconstruction as it continues outside the field of view and so more of it is visible in some projections than others.This violates the X-ray modelling assumption that outside the field of view is vacuum.The more subtle artefact is that over time each surface becomes coated with carbon.This first becomes visible as a thin film around −2 • and progressively gets thicker through the remaining projections.At 34 • we see a bump of carbon appear on the top right edge.Because of these factors, there is both a limited angle pressure and pressure for only few projections.After pre-processing, we have the full dataset as shown top right and artificially sub-sample to compare TV with our proposed reconstruction method.
intensity and we are interested in locating the sharp facets which bound the crystal [45].In Figure 13 we immediately see that the combination of homogeneity and sharp edges is better reconstructed in our proposed reconstruction.In Figure 14 we show the thresholded images to visually emphasise the geometrical properties of each reconstruction.This is particularly helpful when locating boundaries or estimating interior angles of the particle.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.Reconstructions from a slice of the experimental data.We have chosen the slice half-way down through the projections shown in Figure 12 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.

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 guarrantees.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 noticabley violated, as with our experimental data, we still better recover boundaries of the reconsructed 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.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 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 potential for an alternative regularisation functional on the sinogram which is better able to encode structure of the X-ray transform.At the moment we treat sinograms as simple functions on R 2 ignoring the fact that it must also lie in the image of R. The discrepancy between the two is seen in Figure 11 as the reconstruction with our proposed method introduces an artefact due to the model mis-match with the X-ray transform.The question of forming regularisation functions consistent with the global structure of a sinogram naturally leads us to the field of micro-local analysis.From this viewpoint the inpainting problem is separated into what are known as visible and invisible singularities; formalising the idea that edges at particular orientations are 'visible' but noisy in the limited angle data while the rest must be inpainted.Ideally, our two regularisers would operate indpendently on just the visible and invisible domains respectively whereas, at the moment, this separation cannot be made.Overall, we feel that this presents the natural progression for the current work although it remains unclear how to regularise these invisible singularities.
As we are decomposing 2 × 2 matrices it will simplify the proof to write explicit forms for the values under consideration.
We give two equations for e 1 to account for the case when we get (0,0) T 0 .
Proof of part ii: Note that Σ is always smooth and ∆ is smooth on the set {∆ > 0} Case M 12 (x) = 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 ∆ = |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.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, where the remainder term is sufficiently smooth.Hence c i is at least 2k − 1 times differentiable w.r.t.M Hence, we also know w n ϕ converges point-wise to a unique continuous function.Suppose w n k ϕ − w ϕ ∞ ≥ ε > 0 for some ε and sub-sequence n k → ∞.By the Arzela-Ascoli theorem some further subsequence must converge uniformly to the pointwise limit, w ϕ, which gives us the required contradiction.Hence, w n ϕ → w ϕ in L ∞ = C 0 .The general theorem follows by induction.
This lemma gives us two direct results.
Hence, whenever w ∈ W 1,1 we have As has been noted in the main text, there are many hyper-parameters to tune for the best reconstruction.We commonly found that reconstructions were qualitatively insensitive near the optimal parameter choice but we include here some illustrations of the typical effect of each parameter.The full model is To remove a degree of equivalence we have always normalised such that α 2 = 1.To construct the directional TV functional we need 2 smoothing parameters, ρ and σ A d (x) := c 1 (λ 1 (x), λ 2 (x))e 1 (x)e 1 (x) T + c 2 (λ 1 (x), λ 2 (x))e 2 (x)e 2 (x) T such that (∇d ρ ∇d T ρ ) σ (x) = λ 1 (x)e 1 (x)e 1 (x) T + λ 2 (x)e 2 (x)e 2 (x) T λ 1 (x) ≥ λ 2 (x) ≥ 0 Again, we kept ρ = 1 fixed and only show reconstructions for different values of σ.The optimal parameters for the Shepp-Logan phantom referred to below were      is a factor of 0.5 lower than 'optimal' and 'high' a factor of 2 higher.

Figure 1 .
Figure 1.Diagrammatic representation of the acquisition of X-ray transform data in Electron Microscopy in both full range and limited angle acquisition.Note that measurement at 180 • is exactly a reflection of that at 0 • .This symmetry allows us to consider a 180 • 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.

Figure 3 .
Figure 3. 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 4 .
Figure 4. 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.

Figure 4
Figure 4 further shows the limitations of TV when high noise and limited angles are combined.Recently, this classical group of methods have received a revival through machine learning methods for reconstruction[22,23].In both of these examples the inpainting problem remains implicitly treated and the main reduction of artefacts is performed post-inversion.The most common method that has been used to reconstruct pairs (u, v) is to

Figure 5 .
Figure 5. Demonstration of the improvement which can be achieved by using a model as in (1.1).Top images show state of the art reconstructions using Total Variational regularisation (α 1 = α 3 = α 5 = 0).This reconstruction clearly shows the characteristic blurring artefacts at the top and bottom.In our proposed joint reconstruction (bottom line) we can minimise these effects.

Figure 7 .
Figure 7. 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 10 .
Figure 10.Canonical synthetic example.Top row shows the reconstructions, u, while the bottom row shows the reconstructed sinogram, v.

Figure 11 .
Figure 11.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.

Figure 12 .
Figure 12.Raw data for Transmission Electron microscopy example.The wide angle projections, e.g.68 • , contain views of the supporting grid as was described in Figure2suggesting that a more limited angle acquisition is desirable.The plane surface will also cause artefacts in the reconstruction as it continues outside the field of view and so more of it is visible in some projections than others.This violates the X-ray modelling assumption that outside the field of view is vacuum.The more subtle artefact is that over time each surface becomes coated with carbon.This first becomes visible as a thin film around −2 • and progressively gets thicker through the remaining projections.At 34 • we see a bump of carbon appear on the top right edge.Because of these factors, there is both a limited angle pressure and pressure for only few projections.After pre-processing, we have the full dataset as shown top right and artificially sub-sample to compare TV with our proposed reconstruction method.

Figure 13 .
Figure 13.Reconstructions from a slice of the experimental data.We have chosen the slice half-way down through the projections shown in Figure12to 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 14 .
Figure 14.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.
Proof of part iii: Given the Neumann condition on c i we shall express c i locally by Taylor's theorem.

Figure 1 .
Figure 1.Varying reconstruction for low (first column), optimal (middle column) and high (right column) values of β 1 (TV regularisation parameter).'low' is a factor of 0.1 lower than 'optimal' and 'high' a factor of 10 higher.

Figure 2 .
Figure 2. Varying reconstruction for low (first column), optimal (middle column) and high (right column) values of β 2 (DTV regularisation parameter).'low' is a factor of 0.1 lower than 'optimal' and 'high' a factor of 10 higher.

Figure 3 .
Figure3.Varying reconstruction for low (first column), optimal (middle column) and high (right column) values of α 1 (pairing term between u and v).'low' is a factor of 0.1 lower than 'optimal' and 'high' a factor of 10 higher.

Figure 4 .
Figure 4. Varying reconstruction for low (first column), optimal (middle column) and high (right column) values of α 3 (sinogram noise parameter).'low' is a factor of 0.1 lower than 'optimal' and 'high' a factor of 10 higher.

Figure 5 .
Figure 5. Varying reconstruction for low (first column), optimal (middle column) and high (right column) values of σ (smoothing parameter inside DTV functional).'low' is a factor of 0.5 lower than 'optimal' and 'high' a factor of 2 higher.

Table 2 .
Parameter choices for numerical experiments.