Template-Based Image Reconstruction from Sparse Tomographic Data

We propose a variational regularisation approach for the problem of template-based image reconstruction from indirect, noisy measurements as given, for instance, in X-ray computed tomography. An image is reconstructed from such measurements by deforming a given template image. The image registration is directly incorporated into the variational regularisation approach in the form of a partial differential equation that models the registration as either mass- or intensity-preserving transport from the template to the unknown reconstruction. We provide theoretical results for the proposed variational regularisation for both cases. In particular, we prove existence of a minimiser, stability with respect to the data, and convergence for vanishing noise when either of the abovementioned equations is imposed and more general distance functions are used. Numerically, we solve the problem by extending existing Lagrangian methods and propose a multilevel approach that is applicable whenever a suitable downsampling procedure for the operator and the measured data can be provided. Finally, we demonstrate the performance of our method for template-based image reconstruction from highly undersampled and noisy Radon transform data. We compare results for mass- and intensity-preserving image registration, various regularisation functionals, and different distance functions. Our results show that very reasonable reconstructions can be obtained when only few measurements are available and demonstrate that the use of a normalised cross correlation-based distance is advantageous when the image intensities between the template and the unknown image differ substantially.


Introduction
In medical imaging, an image can typically not be observed directly but only through indirect and potentially noisy measurements, as it is the case, for example, in computed tomography (CT) [41]. Due to the severe ill-posedness of the problem, reconstructing an image from measurements is rendered particularly challenging when only few or partial measurements are available. This is, for instance, the case in limited-angle computed tomography [22,41], where limited-angle data is acquired in order to minimise exposure time of organisms to X-radiation. Therefore, it can be beneficial to impose a priori information on the reconstruction, for instance, in the form of a template image. However, typically neither its exact position nor its exact shape is known.
In image registration, the goal is to find a reasonable deformation of a given template image so that it matches a given target image as closely as possible according to a predefined similarity measure, see [40,39] for an introduction. When the target image is unknown and only given through indirect measurements, it is referred to as indirect image registration and has been explored only recently [12,24,31,45]. As a result, a deformation together with a transformed template can be computed from tomographic data. The prescribed template acts as a prior for the reconstruction and, when chosen reasonably close in a deformation sense, gives outstanding reconstructions in situations where only few measurements are available and competing methods such as filtered backprojection [41] or total variation regularisation [47] fail, see [12,Sec. 10].
In our setting, deformations are maps from the image domain Ω ⊂ R n , n ∈ N, to itself together with an action that specifies exactly how such a map deforms elements in the shape space, which in this work is the space L 2 (Ω, R) of greyscale images supported in the image domain. Natural problems are to characterise admissible deformations and to compute these numerically in an efficient manner.
One possible approach is diffeomorphic image registration, where the set of admissible deformations is restricted to diffeomorphisms in order to preserve the topology of structures within an image [58]. One can, for instance, consider the group of diffeomorphisms together with the composition as group operation. Elements in this group act on greyscale images by means of the group action and thereby allow for a rich set of non-rigid deformations, as required in many applications. For instance, the geometric group action transforms greyscale images in a way such that its intensity values are preserved, whereas the masspreserving group action ensures that, when the image is regarded as a density, the integral over the density is preserved.
A computational challenge in using the above group formalism is that it lacks a natural vector space structure, which is typically desired for the numerical realisation of the scheme. Hence, it is convenient to further restrict the set of admissible deformations. One way to obtain diffeomorphic deformations is to perturb the identity map with a displacement vector field. Provided that the vector field is reasonably small and sufficiently regular, the resulting map is invertible [58,Prop. 8.6]. For indirect image registration this idea was pursued in [45].
The basic idea of the large deformation diffeomorphic mapping (LDDMM) [4,18,36,37,50,53,58] framework is to generate large deformations by considering flows of diffeomorphisms that arise as the solution of an ordinary differential equation (ODE), the so-called flow equation, with velocity fields that stem from a reproducing kernel Hilbert space. In order to ensure that the flow equation admits a unique solution, one typically chooses this vector space so that it can be continuously embedded into C 1 (Ω, R n ), allowing the application of existence and uniqueness results from Cauchy-Lipschitz theory for ODEs, see [15,Chap. 1] for a brief introduction. In [12], the LDDMM framework is adapted for indirect image registration and the authors prove existence, stability, and convergence of solutions for their variational formulation. Numerically, the problem is solved by gradient descent.
The variational problem associated with LDDMM is typically formulated as an ODEconstrained optimisation problem. As the flow equation can be directly related to hyperbolic partial differential equations (PDE) via the method of characteristics [21,Chap. 3.2], the problem can equivalently be rephrased as a PDE-constrained optimisation problem [33]. The resulting PDE is determined by the chosen group action, see [12, Sec. 6.1.1] for a brief discussion. For instance, the geometric group action is associated with the transport (or advection) equation, while the mass-preserving group action is associated with the continuity equation. It is important to highlight that the PDE constraint implements both the flow equation and the chosen diffeomorphic group action.
Such an optimal control approach was also pursued for motion estimation and image interpolation [2,5,6,9,13,29,44]. In the terminology of optimal control, the PDE represents the state equation, the velocity field the control, and the transformed image the resulting state. We refer to the books [7,16,26,32] and to the article [30] for a general introduction to PDE-constrained optimisation and suitable numerical methods. Let us mention that other methods, such as geodesic shooting [3,36,49,56], exist and constitute particularly efficient numerical approaches. In particular, this direction has recently been combined with machine learning methods [57].
A particularly challenging scenario for diffeomorphic image registration occurs when the target image is not contained in the orbit of the template image under abovementioned group action of diffeomorphisms. For instance, this could happen in the case of the geometric group action due to the appearance of new structures in the target image or due to a discrepancy between the image intensities of the template and the target image. A possible solution is provided by the metamorphosis framework [38,46,51,52], which is an extension to LDDMM that allows for modulations of the image intensities along characteristics of the flow. The image intensities change according to an additional flow equation with an unknown source. See [58,Chap. 13] for a general introduction and, for instance, [33] for an application to magnetic resonance imaging. Let us also mention [43], which adopts a discrete geodesic path model for the purpose of image reconstruction, and [34], in which the metamorphosis model is combined with optimal transport.
In [24], the metamorphosis framework is adapted for indirect image registration. The authors prove that their formulation constitutes a well-defined regularisation method by showing existence, stability, and convergence of solutions. However, in the setting where only few measurements-e.g. a few directions in CT-are available, reconstruction of appearing or disappearing structures seems very challenging. Therefore, in order to obtain robustness with respect to differences in the intensities between the transformed template and the sought target image, we follow a different approach. We consider not only the standard sum-of-squared differences (SSD) but also a distance that is based on the normalised cross correlation (NCC) [40,Chap. 7.2], as it is invariant with respect to a scaling of the image intensities.
While image registration itself is already an ill-posed inverse problem that requires regularisation [20], the indirect setting as described above is intrinsically more challenging. It can be phrased as an inverse problem, where measurements (or observations) g ∈ Y are related to an unknown quantity f ∈ X via the operator equation Here, K : X → Y is a (not necessarily linear) operator that models the data acquisition, often by means of a physical process, n δ are measurement errors such as noise, and X and Y are Banach spaces. When f constitutes an image and g are tomographic measurements, solving (1) is often referred to as image reconstruction.
We use a variational scheme [48] to solve the inverse problem of indirect image registration, which can be formulated as a PDE-constrained optimisation problem [12,Sec. 6.1.1]. It is given by min where Here, V is an admissible vector space with norm · V , D : Y × Y → R ≥0 is a data fidelity term that quantifies the misfit of the solution against the measurements, and γ > 0 is a regularisation parameter. Moreover, f v (T, ·) : Ω → R denotes the evaluation at time T > 0 of the (weak) solution of C(v), which is either the Cauchy problem governed by the transport equation, or involving the continuity equation. Here, f 0 ∈ L 2 (Ω, R) denotes an initial condition, which in our case is the template image.
The main goals of this article are the following. First, to study variational and regularising properties of problem (2), and to develop efficient numerical methods for solving it. Second, to investigate alternative choices of distance functions D, such as the abovementioned NCCbased distance. Third, to demonstrate experimentally that excellent reconstructions can be computed from highly undersampled and noisy Radon transform data.
Our numerical approach is based on the Lagrangian methods developed in [35], called LagLDDMM. In contrast to most existing approaches, which are mainly first-order methods (see [35] for a brief classification and discussion), LagLDDMM uses a Gauss-Newton-Krylov method paired with Lagrangian solvers for the hyperbolic PDEs listed above. The characteristics associated with these PDEs are computed with an explicit Runge-Kutta method. One of the main advantages of this approach is that Lagrangian methods are unconditionally stable with regard to the admissible step size. Furthermore, the approach limits numerical diffusion and, in order to evaluate the gradient or the Hessian required for optimisation, does not require the storage of multiple space-time vector fields or images at intermediate time instants. The scheme can also be implemented matrix-free.
In comparison to abovementioned existing methods for indirect image registration, such as [12,24,31,45], our method is conceptually different in several ways. The first difference concerns the discretisation. While [12,24,45] are mainly based on small deformations and use reproducing kernel Hilbert spaces, our method relies on nonparametric registration.
The main advantages are that it directly allows for a multilevel approach and no kernel parameters need to be chosen. Moreover, due to the flexibility of the underlying framework it is straightforward to extend our method to parametric registration. Second, our approach relies on second-order methods for optimisation by using a Gauss-Newton method paired with line search, while the other methods mainly rely on gradient descent. This allows for a fast decrease of the objective within only few iterations. Third, our method allows to easily exchange the underlying PDE solver. Essentially, any solver can be used as long as it can be differentiated efficiently. The used explicit Runge-Kutta method has the advantage that it does not require the storage of multiple images or repeated interpolation of the template, which can potentially lead to a blurred solution. Finally, let us mention that [31] is conceptually different since both a deformation and a template image are computed. Our main focus, however, are applications where only very few and noisy measurements are available and the problem of estimating an additional template seems highly underdetermined in such situations.

Contributions
The contributions of this article are as follows. First, we provide the necessary theoretical background on (weak) solutions of the continuity and the transport equation, and recapitulate existence and uniqueness theory for characteristic curves for the associated ODE. In contrast to the results derived in [12], where the template image is assumed to be contained in the space SBV (Ω, R) ∩ L ∞ (Ω, R) of essentially bounded functions with special bounded variation, our results only require L 2 (Ω, R) regularity. In addition, by using results from [17], we are able to consider the transport equation in the setting with H 1 regularity of vector fields in space as well as in time and with bounded divergence. Moreover, we show the existence of a minimiser of the problem (2), stability with respect to the data, and convergence for vanishing noise.
Second, in order to solve the problem numerically, we follow a discretise-then-optimise approach and extend the LagLDDMM framework [35] to the indirect setting. The library itself is an extension of FAIR [40] and, as a result, our implementation provides great flexibility regarding the selected PDE, and can easily be extended to other distances as well as to other regularisation functionals. The source code of our MATLAB implementation is available online. 1 Finally, we present numerical results for the abovementioned distances and PDEs. To the best of our knowledge, the results obtained for indirect image reconstruction based on the continuity equation are entirely novel. Moreover, we propose to use the NCC-based distance instead of SSD whenever the image intensities of the template and the unknown target are far apart, and show its numerical feasibility.

Theoretical results on the transport and continuity equation
In this section, we review the necessary theoretical background, and state results on the existence and stability of weak solutions of the transport and the continuity equation. Compared to [12], our results are stronger since we do not require space regularity of the template image.

Continuity equation
In what follows, we consider well-posedness of the continuity equation that arises in the LD-DMM framework using the mass-preserving group action via the method of characteristics. The regularity assumptions on v are such that we can apply the theory from [51].
Let Ω ⊂ R n be a bounded, open, convex domain with Lipschitz boundary and let T > 0. In the following, we examine the continuity equation Note that such velocity fields can be continuously extended to the boundary. Clearly, equation (4) has to be understood in a weak sense, i.e. a function f In this notation, the first argument of X is the time dependence, the second the initial time, and the third the initial space coordinate. The following theorem is a reformulation of [51, Thms. 1 and 9] and characterises solutions of (6).
in weak sense (absolutely continuous solutions). The solution operator can be used to characterise solutions of (4) as follows.
, then the unique weak solution of (4), as defined in (5), is given by Proof. The proof is divided in three steps. First, we want to show that f satisfies the regularity conditions of weak solutions. For this purpose, the first step is to show X(0, ·, ·) ∈ C 0 ([0, T ], C 0 (Ω, R n )), i.e. that the flow is continuous in the initial values. Clearly, X(0, t, ·) ∈ C 0 (Ω, R n ) for every t ∈ [0, T ]. For an arbitrary sequence t i → t we get where the first factor is bounded due to [51, Lemma 9]. Next, using the sequence X i (·) = X(0, t i , ·), it follows f 0 (X(0, ·, ·)) ∈ C 0 ([0, T ], L 2 (Ω, R)), where the continuity in time follows from [42,Cor. 3]. Then, by differentiating X(0, t, X(t, 0, x)) = x and rearranging the terms we obtain which follows from (7) since both summands converge to zero.
The second step is to show that (5) is satisfied. Note that X(·, 0, x) is differentiable in t for a.e. t ∈ [0, T ], since it is absolutely continuous by definition. By inserting f into (5) and using the transformation formula, we get For the last equality we used that η(t, X(t, 0, x)) is absolutely continuous. The last step is to prove uniqueness of weak solutions, i.e. that every solution has the given form. Let f 1 , f 2 be two different solutions, then we can find a t such that f 1 (t, ·) − f 2 (t, ·) L 2 (Ω) > 0. By continuity in time we can find an interval I of length δ > 0 that contains t, and a constant c > 0 such that This yields a contradiction.
Additionally, we can state and prove the following stability result for solutions of (4).
Proof. The solution of (6) corresponding to v i is denoted by It is left to show that a subsequence, also denoted by X i , exists such that X i (0, t, ·) → X(0, t, ·) in C 1 (Ω, R n ). This concludes the proof since it also implies the convergence of

Transport equation with H 1 regularity
Here, we prove well-posedness of the transport equation that arises in the LDDMM framework using the geometric group action. Compared to the previous section, the space regularity assumptions on v are weaker and fit the setting in [17].
The transport equation reads as with coefficients for some fixed constant C and initial value f 0 ∈ L 2 (Ω, R). The admissible set A consists of all H 1 functions that are zero on the boundary of the spatial domain and have bounded divergence in the L ∞ norm. Note that the set A is a subset of H 1 ([0, T ] × Ω) n , which is closed and convex so that it is a weakly closed subset of a reflexive Banach space. In the following, we only check that A is closed. Let v i be a convergent sequence in A with limit v. Since the two involved spaces are Banach spaces, we only have to check that v satisfies the constraint. Assume that div x v L ∞ ([0,T ]×Ω) > C, then there exists a set B with positive measure and an > 0 such that for all , which contradicts the convergence in H 1 . Again, equation (9) has to be understood in the weak sense so that f ∈ C 0 ([0, T ], L 2 (Ω, R)) is said to be a solution of (9) if it satisfies Proof. The existence and uniqueness of weak solutions follows from [17, Corrs. II.1 and II.2]. Note that these solutions are also renormalised due to [17,Thm. II.3].
We recast the second part of the theorem such that it has the exact form of [17, Thm. II.5]. First, note that both the velocity fields and the initial condition can be extended to R n by zero outside of Ω due to boundary condition of A. Due to the conditions on v, the weak formulation is equivalent to the one for the extension in the R n setting. The uniform boundedness condition on f i is satisfied since Ω is bounded.
Proof. Combine the previous theorem with the compact embedding of

Regularising properties of template-based image reconstruction
In this section, we prove regularising properties of template-based reconstruction as defined in (2). Recall that the problem reads where C(v) is the Cauchy problem with either the transport or the continuity equation.
The admissible set V is chosen such that the regularity requirements stated in the previous section are satisfied. For the following considerations we require these assumptions on K and D: 1. The operator K is continuous, D(·, g) is lower semi-continuous for each g ∈ Y , and D(g, ·) is continuous for each g ∈ Y .
2. If f n , g n are two convergent sequences with limits f and g, respectively, then D must satisfy lim inf n→∞ D(f n , g) ≤ lim inf n→∞ D(f n , g n ).
Note that the requirements on D are satisfied if D is a metric. The obtained results are along the lines of [12] but are adapted to our setting and to our notation. For simplicity we stick to the notation of the continuity equation, but want to mention that the same derivations hold for the transport equation with coefficients in the set A. First, we prove that a minimiser of the problem exists. Proof. The idea of the proof is to construct a minimising sequence which is weakly convergent and then use that the functional is weakly lower semi-continuous. Let us consider a sequence v n such that J γ,g (v n ) converges to inf v J γ,g (v). By construction of the functional, v n is bounded in L 2 ([0, T ], V) and hence there exists a subsequence, also denoted with v n , such that v n v ∞ . By Prop. 2.1.3, there exists a subsequence, also denoted with v n , such that f vn (T, ·) → f v∞ (T, ·) in L 2 (Ω, R). With this at hand, we are able to prove weak lower semi-continuity of the data term. Indeed, as K is continuous, from f vn (T, ·) → f v∞ (T, ·) we get K(f vn (T, ·)) → K(f v∞ (T, ·)). Since D(·, g) is lower semi-continuous, we obtain that D (K(f v∞ (T, ·)), g) ≤ lim inf n→∞ D (K(f vn (T, ·)), g). This concludes the proof, since the whole functional is (weakly) lower semi-continuous, and hence Next, we state a stability result.
Hence, v n is bounded in L 2 ([0, T ], V) and there exists a subsequence, also denoted with v n , such that v n v. From the weak convergence we obtain γ v 2 V ≤ γ lim inf n→∞ v n 2 V . By passing to a subsequence and by using Prop. 2.1.3, we deduce that f vn (T, ·) → f v (T, ·). Together with the convergence of g n and the convergence property of D this implies Thus, for anyṽ, it holds that because v n minimises J γ,gn . Then, as J γ,gn (ṽ) converges to J γ,g (ṽ) by the assumptions on D, we deduce J γ,g (v) ≤ J γ,g (ṽ) and hence that v minimises J γ,g .
Let v n be a minimiser of J γn,gn , where γ n = γ(δ n ). Then, there exists a subsequence of v n which weakly converges towards an element v such that K(f v (T, ·)) = g.
From the requirements on γ and δ we deduce that v n is bounded in L 2 ([0, T ], V) and then that up to an extraction, v n weakly converges to some v in L 2 ([0, T ], V). Further, it holds D(K(f v (T, ·)), g) ≤ lim inf n→∞ D(K(f vn (T, ·)), g n ) with the same arguments as in the previous proposition. Finally, for every n, it holds that where the two rightmost terms both converge to zero. Thus, K(f v (T, ·)) = g by the assumptions on D.
We conclude with a remark on data discrepancy functionals that satisfy the conditions and will be used in our numerical experiments in Sec. 5.
We will only check the convergence condition. It holds where the last two terms converge to zero since convergent sequences are bounded.

Another function that satisfies the conditions is
which is based on the normalised cross correlation. First, note thatD(·, g) = ·,g 2 g 2 Y and the function · −2 Y are continuous. Thus, we get that D NCC (·, g) is continuous. By symmetry, this also holds for D NCC (g, ·). It remains to check the convergence property: = lim From this we conclude lim inf n→∞ D NCC (f n , g) = lim inf n→∞ D NCC (f n , g n ). Unfortunately, D NCC (f, g) = 0 only implies f = cg, with c ∈ R.

Numerical solution
The focus of this section is to approximately solve problem (2). Our approach is based on the Lagrangian methods developed in [35] and the inexact multilevel Gauss-Newton method used in [40]. Both methods and their necessary modifications are briefly outlined here. As customary in PDE-constrained optimisation [16, Chap. 3], we eliminate the state equation by defining a control-to-state operator, which parametrises the final state f v (T, ·) in terms of the unknown velocities v. With a slight abuse of notation, we define this solution map as Here, f v denotes the unique solution to either the transport or the continuity equation, as defined in Sec. 2. As a result, we obtain the reduced form of (2): Here, R : V → R ≥0 is a regularisation functional that can be written as with B denoting a linear (vectorial) differential operator.
In this work, we consider the operators B = ∇ x and B = ∆ x , which are also used in [35]. We refer to the resulting functionals R as diffusion and curvature regularisation functionals, respectively. Note that B can as well be chosen to incorporate derivatives with respect to time.
Amongst the operators above, we also consider a regularisation functional that resembles the norm of the space V = L 2 ([0, T ], H 3 0 (Ω, R n )). This particular choice is motivated by the fact that, for n = {2, 3}, the space H 3 0 (Ω, R n ) can be continuously embedded in C 1,α 0 (Ω, R n ), for some α > 0, so that the results in Sec. 2 hold. The norm of V is given by Here, |·| H k (Ω,R n ) denotes the usual H k -seminorm including only the highest-order partial derivatives. By the Gagliardo-Nirenberg inequality, (19) is equivalent to the usual norm of L 2 ([0, T ], H 3 0 (Ω, R n )). To simplify numerical optimisation, we omit the requirement that v is compactly supported in Ω and minimise over L 2 ([0, T ], H 3 (Ω, R n )).
In order to solve problem (17), we follow a discretise-then-optimise strategy. Without loss of generality, we assume that the domain is Ω = (0, 1) n . We partition it into a regular grid consisting of m n equally sized cells of edge length h X = 1/m in every coordinate direction.
The template image f 0 ∈ L 2 (Ω, R) is assumed to be sampled at cell-centred locations x c ∈ R m n , giving rise to its discrete version f 0 (x c ) ∈ R m n . The template image is interpolated on the cell-centred grid by means of cubic B-spline interpolation as outlined in [40,Chap. 3.4]. Similarly, the time domain is assumed to be [0, 1] and is partitioned into m t equally sized cells of length h t = 1/m t . We assume that the unknown velocities v : [0, 1] × Ω → R n are sampled at cell-centred locations in space as well as at cell-centred locations in time, leading to a vector of unknowns v ∈ R N , where N = (m t +1)·n·m n is the total number of unknowns of the finite-dimensional minimisation problem.
Lagrangian solver In order to compute the solution map f (v) numerically, i.e. to solve the hyperbolic PDEs (4) and (9), the Lagrangian solver in [35] follows a two-step approach. First, given a vector v ∈ R N of velocities, the ODE (6) is solved approximately using a fourth-order Runge-Kutta (RK4) method with N t equally spaced time steps of size ∆t. For simplicity, we follow the presentation in [35] based on an explicit first-order Euler method and refer to [35,Sec. 3.1] for the full details.
Given initial points x ∈ R m n and velocities v ∈ R N , an approximation X v : [0, 1] 2 ×R m n → R m n of the solution X v is given recursively by for all k = 0, 1, . . . , N t − 1, with initial condition X v (0, 0, x) = x. Here, I(v, t k , X v (0, t k , x)) denotes a componentwise interpolation of v at time t k = k∆t and at the points X v (0, t k , x). Note that, since the characteristic curves for both PDEs coincide, this step is identical regardless of which PDE we impose. The second step computes approximate intensities of the final state f v (1, ·). This step depends on the particular PDE. For the transport equation, in order to compute the intensities at the grid points x c , we follow characteristic curves backwards in time, which is achieved by setting ∆t = −1/N t in (20). The deformed template is then given by where f 0 ∈ R m n is the interpolation of the discrete template image.
For the continuity equation, [35] proposes to use a particle-in-cell (PIC) method, see [14] for details. The density of particles which are initially located at grid points x c is represented by a linear combination of basis functions, which are then shifted by following the characteristics computed in the first step. To determine the final density at grid points, exact integration over the grid cells is performed. By setting ∆t = 1/N t in (20), the transformed template can be computed as where F ∈ R N ×N is the pushforward matrix that computes the integrals over the shifted basis functions. See [35,Sec. 3.1] for its detailed specification using linear, compactly supported basis functions. By design, the method is mass-preserving at the discrete level.
Numerical optimisation Let us denote by K : R N → R M , M ∈ N, a finite-dimensional, Fréchet differentiable approximation of the (not necessarily linear) operator K : L 2 (Ω, R) → Y . With the application to CT in mind, we will outline a discretisation of (17) suitable for the n-dimensional Radon transform, which maps a function on R n into the set of its integrals over the hyperplanes in R n [41,Chap. 2].
An element K(f (v)) ∈ Y is a function on the unit cylinder S n−1 × R of R n+1 , where S n−1 is the (n − 1)-dimensional unit sphere. We discretise this unit cylinder as follows. First, we sample p ∈ N directions from S n−1 . When n = 2, as it is the case in our experiments in Sec. 5, directions are parametrised by angles from the interval [0, 180] degrees. For simplicity, we say (slightly imprecise) that we take one measurement in each direction. Second, similarly to the sampling of Ω, we use an interval (0, 1) instead of R and partition it into q equally sized cells of length h Y = 1/q. Depending on n and the diameter of Ω, the interval length may require adjustment. Each measurement i is then sampled at cell-centred points y c ∈ R q and denoted by g i (y c ) ∈ R q . All measurements are then concatenated into a vector g : The finite-dimensional optimisation problem in abstract form is then given by where D and R are chosen to be discretisations of a distance and of (18), respectively. In further consequence, we approximate integrals using a midpoint quadrature rule. As we are mainly interested in the setting where only few directions are given, we disregard integration over the unit sphere. For vectors x, y ∈ R M , the corresponding approximations of the distance based on sum-of-squared-differences and the normalised cross correlationbased distance are then respectively. See [40, Chaps. 6.2 and 7.2] for details. Note that, due to cancellation, no (spatial) discretisation parameter occurs in the approximation of the NCC above. Moreover, we approximate the regularisation functional in (18) with where B ∈ R N ×N is a finite-difference discretisation of the differential operator in (18), analogous to [39,Chap. 8.5]. In our implementation we use zero Neumann boundary conditions and pad the spatial domain to mitigate boundary effects arising from the discretisation of the operator. In order to apply (inexact) Gauss-Newton optimisation to problem (23), we require firstand (approximate) second-order derivatives of J γ,g (v). By application of the chain rule, we obtain where ∂K/∂f is the Fréchet derivative of K and ∂f (v)/∂v is the derivative of the solution map (16) with respect to the velocities, which is given below. The partial derivatives of the distance functions (24) with respect to its first argument are given by where I N ∈ R N ×N is the identity matrix of size N , and ∂ ∂x D NCC (x, y) = − 2(x y)y x 2 y 2 + 2(x y) 2 x x 4 y 2 respectively. Moreover, the derivatives of (25) are given by In order to obtain an efficient iterative second-order method for solving (23), one requires an approximation of the Hessian H ∈ R N ×N that balances the following tradeoff. Ideally, it is reasonably efficient to compute, consumes limited memory (sparsity is desired), and has sufficient structure so that preconditioning can be used. However, each iteration of the Gauss-Newton method should also provide a suitable descent direction. For these reasons, we approximate the Hessian by where > 0 ensures positive semidefiniteness. For simplicity, the term involving ∂ 2 f (v)/∂v 2 is omitted and, regardless of the chosen distance, we use the second derivative in (26) as an approximation of ∂ 2 D(x, y)/∂x 2 . In our numerical experiments we found that this choice works well for the problem considered in Sec. 5. It remains to discuss the derivative of the solution map. For the transport equation, the application of the chain rule to (21) yields where ∇ x f 0 denotes the gradient of the interpolation of the template image and ∂X v /∂v is the derivative of the endpoints of the characteristic curves with respect to the velocities, see below. Similarly, for the solution map (22) that corresponds to the continuity equation, we obtain ∂ ∂v Here, ∂F/∂X v is the derivative of the pushforward matrix with respect to the endpoints of the characteristics, again see [35,Sec. 3.1].
If explicit time stepping methods are used to solve the ODE (6), the partial derivative ∂X v /∂v can be computed recursively. For example, for the forward Euler approach in (20) it is given by for all k = 0, 1, . . . , N t −1, with ∂I/∂v and ∂I/∂X v being the derivatives of the interpolation schemes with respect to the velocities and with respect to the endpoints of the characteristics, respectively. We refer to [40,Chap. 3.5] for details. The case where characteristics are computed backwards in time can be handled similarly.
In order to solve the finite-dimensional minimisation problem (23), we apply a inexact Gauss-Newton-Krylov method, which proceeds as follows. Given an initial guess v (0) = 0, we update the velocities in each iteration i = 0, 1, . . . by v (i+1) = v (i) + µδv until a termination criterion is satisfied. Here, µ ∈ R denotes a step size that is determined via Armijo line search and δv ∈ R N is the solution to the linear system For details on the stopping criteria and the line search we refer to [40,Chap. 6.3.3]. We solve the system (27) approximately by means of a preconditioned conjugate gradient (PCG) method, which can be implemented matrix-free whenever the derivative of K and its adjoint can be computed matrix-free. See [35,Sec. 3.2] for further details on the preconditioning. Due to the non-convexity of (17) and to speed up computation, we use a multilevel strategy in order to reduce the risk of ending up in a local minimum, see [27]. On each level, we use a subsampled version of the velocities that were computed on the previous, more coarser, discretisation as initial guess. (c) Measured (sinogram) data without noise. Figure 1: Synthetic example based on an artificial brain image [25] that has been deformed manually. We generated six Radon transform measurements that correspond to six equally spaced angles from the interval [0, 60] degrees.
While standard image registration typically uses resampling of the template and the target image [40,Chap. 3.7], the approach described here requires multilevel versions of the operator K together with a suitable method for resampling the measurements g. We stress that, if these are not available, optimisation can as well just be performed on the finest discretisation level.
In the following, we assume that K is a discretisation of the Radon transform [41], which is a linear operator, and outline a suitable procedure for creating multilevel versions of the operator and the measured data. The former is easily achieved with a computational backend such as ASTRA [54,55], which allows to explicitly specify the number of grid cells used to discretise the measurement geometry. For the sake of simplicity, we restrict the presentation here to the case where n = 2, i.e. Ω ⊂ R 2 , and K is linear.
Let us assume that the number of grid cells used to discretise Ω at the finest level is m = 2 , ∈ N. In our experiments, we set the number of grid cells of the one-dimensional measurement domain (0, 1) at the current level k ≤ to q (k) = 1.5 · 2 (k) and set the length of each cell to h (k) Y = 1/q (k) . Then, a multilevel representation of each measurement g i , i ≤ p, at cell-centred grid points y j = (j − 1/2)h (k−1) Y is given by where the denominator arises from averaging over two neighbouring grid points and dividing the edge length of the imaging domain Ω in each coordinate direction in half. The approach can be extended to higher dimensions.

Numerical examples
In our numerical experiments we use the Radon transform [41] as operator. Other choices are possible and, assuming that one has access to a suitable resampling procedure for the measured data, the multilevel strategy can be applied as well. The aim here is to investigate (a) Reconstruction using filtered back-projection.
(c) Reconstruction using R 2 with given template.
(d) Reconstruction using the metamorphosis approach [24]. Figure 2: Comparison of different reconstruction models applied to an artificial brain image [40] that has been deformed manually. We generated six measurements that correspond to six equally spaced angles from the interval [0, 60] degrees.
the reconstruction quality with different regularisation functionals, distances, and noise levels for both PDE constraints. We show synthetic examples for the settings n = 2 and n = 3, and non-synthetic examples for n = 2 using real X-ray tomography data. In the synthetic case, all shown reconstructions were computed from measurements taken from at most 10 directions (i.e. angles) sampled from intervals within [0, 180] degrees. All computations were performed using an Intel Xeon E5-2630 v4 2.2 GHz server equipped with 128 GB RAM and an NVIDIA Quadro P6000 GPU featuring 24 GB of memory. The GPU was only used for computing the Radon transform of 3D volumes.
Before we proceed, we give a brief idea of suitable parameter choices. For the multilevel approach we used in each synthetic example 32×32 pixels at the coarsest level and 128×128 pixels at the finest level, i.e. = 7. The size of the reconstructed images in the nonsynthetic examples was 128 × 128. Again, three levels were used. In the synthetic 3D example the reconstructed volume was 32 × 32 × 32 and the coarsest level was 8 × 8 × 8.
We used time dependent velocity fields with only one time step, i.e. n t = 1, since this keeps the computational cost reasonable and sufficed for our examples. The characteristics were computed using five Runge-Kutta steps, i.e. N t = 5.
The spatial regularisation parameter depends on the chosen regularisation functional and the noise level, and was chosen in the order of 10 −3 , 10 0 , and 10 3 for third-order, curvature, and diffusion regularisation, respectively, in the noisefree case and using the NCC-based distance. The temporal regularisation parameter is less sensitive and was chosen in the order of 10 2 . Furthermore, the parameter corresponding to the norm of L 2 (Ω, R n ) in (19) was set to 10 −6 .
In our first example, we investigate different regularisation functionals with different noise levels together with the the transport equation. The target is 2D Radon transform data based on a digital brain image and the template is a deformed version thereof, see Fig. 1.
Since we want to focus on the behaviour of the regularisation functionals, we do not treat the continuity equation here. The data was generated using parallel beam tomography with only six equally distributed angles from the interval [0, 60] degrees and was corrupted with Gaussian white noise of different levels. Fig. 2 shows results obtained from the generated noisefree measurements using four existing methods. In Fig. 2a filtered backprojection was used. In Figs. 2b and 2c, the following two total variation regularisation-based models, see e.g. [10], with R 1 (u) := TV(u), R 2 (u) := TV(u−f 0 ), and γ > 0 were used. Here, R 2 (u) incorporates template information. Approximate minimisers of both functionals were computed using the primal-dual hybrid gradient method [11]. For the case of filtered backprojection, the standard MATLAB implementation was used. The results in Figs. 2a to 2c highlight why more sophisticated methods, such as the proposed template-based approach, are necessary to obtain satisfying reconstructions in this setting, and illustrate the challenges when dealing with very sparse data.
As outlined in Sec. 1, one possibility is the metamorphosis approach [24]. In Fig. 2d we show a result obtained with this method using the recommended parameters. However, 200 iterations of gradient descent were performed, and the regularisation parameters were set to γ = 10 −5 and τ = 1. Observe the change in image intensities compared to Fig. 1a and the blur in the heavily deformed regions.
In Fig. 3, we show results for the different noise levels and different regularisation functionals computed with our approach. All results were obtained using the NCC-based distance. As expected, the quality of the reconstruction gets worse for higher noise levels and, consequentially, larger regularisation parameters were necessary. Since data is acquired from only six directions, the influence of the noise is very strong. Especially for the diffusive regularisation we needed to choose large regularisation parameters for higher noise levels, see Fig. 3c. Since diffusion corresponds to first-order regularisation, it is much easier to reconstruct the noise with "rough" deformations. Overall, we found that second-and third-order regularisation performed similar when appropriate regularisation parameters were chosen. Even though some theoretical results only hold for higher-order regularity, second-order regularisation seems sufficient for our use case. The computation time for the results in Fig. 3 was between 200 and 700 seconds.
In the second example, see Fig. 4, we compare the behaviour of the SSD and the NCCbased distance. The example consists of two different hands which, in addition, are rotated relative to each other. Here, the deformation is much larger than in the previous example, but still fairly regular. The data was generated similarly to the previous example, but with only five angles from the interval [0, 75] degrees. Note also that the intensities of the template and target image are different (roughly by a factor of two). First, we discuss the transport equation. The intensity difference is a serious issue if we use the SSD distance, as we can see in Fig. 4e. The hand is deformed into a smaller version in order to compensate the differences. If we use the NCC-based distance instead, which can deal with such discrepancies, the result is much better from a visual point of view. The shapes are well-aligned. The resulting SSIM value is still low, which is not surprising since SSIM is not invariant with respect to intensity differences between perfectly aligned images. However, neither of the two approaches is able to remove or create any of the additional (noise) structures in the images. For the combination SSD with continuity equation, no satisfactory results could be obtained. Since no change of intensity is possible by changing the size of the hand, part of it is moved outside of the image. This behaviour could potentially be corrected if    [40] images with different image intensity levels using our method. We generated five measurements that correspond to five equally spaced angles from the interval [0, 75] degrees and added five percent noise.
other boundary conditions are used in the implementation. Therefore, we do not provide an example image for this case. Using the NCC-based distance, the results look similar as for the transport equation with slightly worse SSIM value. These results suggest that the NCC-based distance is a more robust choice that avoids unnatural deformations, which would be required in the case of SSD to compensate intensity differences. In this example, the computation time was between 50 and 325 seconds.
In the next example, see Fig. 5, we compare the continuity equation with the transport equation as constraint together with the NCC-based distance. The continuity equation allows for limited change of mass along the deformation path. Since the intensity change scales with the determinant of the Jacobian, bigger changes are only possible if areas are compressed or extended a lot. In the presented example this occurs only to a mild extent. For this example, the continuity equation and the transport equation yield visually similar results with minor differences in the SSIM value. As in the previous examples, higherorder regularisation is beneficial and artefacts occur for the diffusion regularisation. The computation time amounted to roughly 64 to 360 seconds in this example.
In Fig. 6, we created an artificial pair of images consisting of a disk to show the possibilities of intensity changes when using the continuity equation as a constraint. Both template and unknown image were constructed so that their total mass is equal. The measurements were generated as before using only five angles uniformly distributed in the interval [0, 90] degrees. Furthermore, we used curvature regularisation. For the transport equation we observe that the shape is matched, but the intensity is not correct, see Fig. 6d. If we use the continuity equation instead, intensity changes are possible, which can be observed in Fig. 6e. The computation time for the two results was 90 and 500 seconds.
In order to demonstrate the practicality of our method, we computed results from nonsynthetic X-ray tomography data [8,28], which are available online. 2,3 See Figure 7 for these two examples ('lotus' and 'walnut'). The template was generated by applying filtered backprojection to the full measurements and by subsequently deforming it. Then, this deformed templated was used in our method to compute a reconstruction from only few measurement directions. The computation time amounted to roughly 80 and 600 seconds in these examples. In both nonsynthetic examples the use of the NCC-based distance proved crucial and no satisfactory result could be obtained using SSD.
In Fig. 8, we demonstrate that our framework is also capable of reconstructing 3D volumes. Here, we used the SSD distance together with curvature regularisation and the transport equation. We applied the 3D Radon transform to obtain ten measurements from angles within [0, 180]. The total computation time was roughly 800 seconds.
All in all, our results demonstrate that, given a suitable template image, very reasonable reconstructions can efficiently be obtained from only a few measurements, even in the presence of noise. Moreover, our examples show that the NCC-based distance adds robustness to the approach with regard to discrepancies in the image intensities.

Conclusions
Overall, our numerical examples show that our implementation yields good results, as long as the deformation between template and target is fairly regular. By using the NCCbased distance, robustness with respect to intensity differences between the template and the target image can be achieved. As already mentioned in the introduction, we do not follow the metamorphosis approach, since there is too much flexibility in the model and the source term is very likely to reproduce noise and artefacts if the data is too limited. It is left for further research to investigate possible adaptations of the model that allow for the appearance of new objects or structures in the reconstruction without reproducing noise or artefacts. Possibly, the results of our method can be used as better template for other algorithms that require template information. Finally, note that due to the great flexibility of the FAIR library, it is also possible to use a great variety of regularisation functionals for the velocities and other distances, see [40,Chaps. 7 and 8]. Additionally, our implementation is not necessarily restricted to the Radon transform and essentially every (continuous) operator can be used. The multilevel approach can be applied as long as a meaningful resampling procedure for the operator and the measured data can be provided.   [40] image using our approach, different regularisation functionals, and different PDE constraints. Here, ten measurements corresponding to ten angles equally distributed in the interval [0, 180] degrees were taken. The measured data was corrupted with five percent noise.   (c) Measured noisy Radon transform data. (d) Reconstruction using the SSD distance, curvature regularisation, and the transport equation, SSIM 0.887. Figure 8: Reconstructions of a 3D volume ('mice3D', see [40]) using our method. In Fig. 8a,  Fig. 8b, and Fig. 8d, slices (left to right, top to bottom) of each volume along the third coordinate direction are shown. In Fig. 8c, slices of the 3D Radon transform measurements are shown. Each slice corresponds to one measurements direction.
In total, only ten measurements were taken at angles equally distributed in [0, 180] degrees. As before, five percent noise was added.