Chong Chen
email: chench@lsec.cc.ac.cn

Barbara Gris

Ozan Öktem

Joint Image Reconstruction and Motion Estimation for Spatiotemporal Imaging *

Keywords: spatiotemporal imaging, image reconstruction, motion estimation, joint variational model, shape theory, large diffeomorphic deformations AMS subject classifications. 65F22, 65R32, 65R30, 65D18, 65J22, 65J20, 65L09, 68U10, 94A12, 94A08, 92C55, 54C56, 57N25, 47A52

We propose a variational model for joint image reconstruction and motion estimation applicable to spatiotemporal imaging. This model consists of two parts, one that conducts image reconstruction in a static setting and another that estimates the motion by solving a sequence of coupled indirect image registration problems, each formulated within the large deformation diffeomorphic metric mapping framework. The proposed model is compared against alternative approaches (optical flow based model and diffeomorphic motion models). Next, we derive efficient algorithms for a time-discretized setting and show that the optimal solution of the time-discretized formulation is consistent with that of the time-continuous one. The complexity of the algorithm is characterized and we conclude by giving some numerical examples in 2D space + time tomography with very sparse and/or highly noisy data.

two steps jointly [START_REF] Blume | Joint image and motion reconstruction for PET using a B-spline motion model[END_REF][START_REF] Gigengack | Motion Correction in Thoracic Position Emission Tomography[END_REF]. In the sequential approach, one starts with applying static image reconstruction on data from each of the gates, resulting in a series of (low-resolution) images.

Then, one selects a reconstructed image as target and registers the other reconstructed images against this target. The final reconstructed image is obtained as an average of all the registered images [START_REF] Bai | Regularized B-spline deformable registration for respiratory motion correction in PET images[END_REF][START_REF] Bai | Motion correction and attenuation correction for respiratory gated PET images[END_REF][START_REF] Dawood | Respiratory motion correction in 3-D PET data with advanced optical flow algorithms[END_REF][START_REF] Dawood | Lung motion correction on respiratory gated 3-d PET/CT images[END_REF][START_REF] Gigengack | Motion correction in dual gated cardiac PET using mass-preserving image registration[END_REF][START_REF] Gigengack | Motion Correction in Thoracic Position Emission Tomography[END_REF]. The other alternative that seeks to jointly perform image reconstruction and motion estimation is much more complex. Several approaches have been suggested for how to do this, such as [START_REF] Blume | Joint image and motion reconstruction for PET using a B-spline motion model[END_REF][START_REF] Schumacher | Combined reconstruction and motion correction in SPECT imaging[END_REF][START_REF] Taguchi | Motion Compensated Fan-Beam Reconstruction for Nonrigid Transformation[END_REF][START_REF] Gao | Robust principal component analysis-based four-dimensional computed tomography[END_REF][START_REF] Liu | 5D respiratory motion model based image reconstruction algorithm for 4D cone-beam computed tomography[END_REF][START_REF] Brehm | Self-adapting cyclic registration for motion-compensated cone-beam CT in image-guided radiation therapy[END_REF][START_REF] Brehm | Artifact-resistant motion estimation with a patient-specific artifact model for motion-compensated cone-beam CT[END_REF][START_REF] Ritschl | Iterative 4D cardiac micro-CT image reconstruction using an adaptive spatio-temporal sparsity prior[END_REF][START_REF] Jia | 4D computed tomography reconstruction from fewprojection data via temporal non-local regularization[END_REF][START_REF] Tang | A fully four-dimensional, iterative motion estimation and compensation method for cardiac CT[END_REF][START_REF] Bhagalia | Nonrigid registration-based coronary artery motion correction for cardiac computed tomography[END_REF][START_REF] Hinkle | 4D CT image reconstruction with diffeomorphic motion model[END_REF][START_REF] Hinkle | 4D CT image reconstruction with diffeomorphic motion model[END_REF][START_REF] Brehm | Cardiorespiratory motioncompensated micro-CT image reconstruction using an artifact model-based motion estimation[END_REF][START_REF] Burger | A Variational Model for Joint Motion Estimation and Image Reconstruction[END_REF].

The approach taken in this paper seeks to jointly perform image reconstruction and motion estimation. The motion is here modelled as deformations induced by diffeomorphisms that are provided by the large deformation diffeomorphic metric mapping (LDDMM) framework, which is a well-developed framework for diffeomorphic image registration [START_REF] Trouvé | Diffeomorphisms Groups and Pattern Matching in Imaging Analysis[END_REF][START_REF] Dupuis | Variational Problems on Flows of Diffeomorphisms for Imaging Matching[END_REF][START_REF] Miller | On the metrics and Euler-Lagrange equations of computational anatomy[END_REF][START_REF] Beg | Computing large deformation metric mappings via geodesic flow of diffeomorphisms[END_REF][START_REF] Grenander | Pattern Theory. From Representation to Inference[END_REF][START_REF] Younes | Shapes and Diffeomorphisms[END_REF][START_REF] Trouvé | Shape spaces[END_REF][START_REF] Bruveris | Geometry of image registration: The diffeomorphism group and momentum maps[END_REF][START_REF] Bruveris | On completeness of groups of diffeomorphisms[END_REF]. Joint image reconstruction and motion estimation can be based on the growth model in LDDMM [START_REF] Grenander | Pattern Theory. From Representation to Inference[END_REF] as shown in [START_REF] Hinkle | 4D CT image reconstruction with diffeomorphic motion model[END_REF] for 3D space + time computed tomography (CT).

It can also be used for indirect image registration as shown in [START_REF] Chen | Indirect Image Registration with Large Diffeomorphic Deformations[END_REF], which also proves the resulting model is stable. Specific contributions. The main contribution is a new variational model for joint image reconstruction and motion estimation in spatiotemporal imaging based on the LDDMM framework with deformable templates. As already mentioned, there are two components: one corresponding to modified static image reconstruction, and the other to estimating the motion. To handle the latter, we generalize the LDDMM framework to a setting that applies to a series of coupled sequential indirect image registration problems.

The mathematical properties of the proposed variational model is compared against the optical flow based model in [START_REF] Burger | A Variational Model for Joint Motion Estimation and Image Reconstruction[END_REF] and the diffeomorphic motion model in [START_REF] Hinkle | 4D CT image reconstruction with diffeomorphic motion model[END_REF]. The comparison shows that the proposed model has some desirable properties in terms of the optimal solution, e.g., guaranteeing elastically large diffeomorphic deformations, averagely distributed w.r.t. time t, and non-vanishing at both the initial and the end time points, etc. Moreover, a computationally efficient gradient-based iterative scheme is presented for a time-discretized formulation. More importantly, the optimal solution of the time-discretized problem is shown to be consistent with that of the time-continuous one. Most of the computationally demanding parts relate to computations involving linearized deformations [START_REF] Öktem | Shape based image reconstruction using linearised deformations[END_REF].

Outline. Section 2 presents a general variational model for joint image reconstruction and motion estimation. We also briefly review necessary parts of the LDDMM theory in subsection 3.1. Subsection 3.2 proposes the new variational model and subsection 3.3 makes the theoretical comparison between this ans other existing models. Section 4 gives the detailed numerical algorithms associated with the proposed model. Results of numerical experiments are presented in section 5 in the context of 2D space + time tomography. Finally, we conclude with section 6 that discusses future work.

2. A variational model for joint image reconstruction and motion estimation. Spatiotemporal imaging often leads to the task of recovering a spatially distributed quantity (image) that exhibits temporal variations given indirect time-dependent noisy observations (measured data). Hence, both the image and its motion are unknown.

General spatiotemporal inverse problem.

To formalise the notion of spatiotemporal imaging, let f : [t 0 , t 1 ] × Ω → R k denote the spatiotemporal image we need to recover. Here k This manuscript is for review purposes only.

is the number of channels/modalities (often k = 1) and Ω ⊂ R n is the spatial region containing the support of the image. Without loss of generality, the general (quasi-)time domain [t 0 , t 1 ] can be re-parameterized onto [0, 1].

The spatiotemporal inverse problem is now stated mathematically as the problem of reconstructing the spatiotemporal image f (t, • ) ∈ X given measured data g(t, • ) ∈ Y where (2.1) g(t, • ) = T t, f (t, • ) + g noise (t, • ) for t ∈ [0, 1].

Here, X (reconstruction space) is the vector space of all possible images on a fixed domain Ω, Y (data space) is the vector space of all possible data, and g noise (t, • ) ∈ Y is the observation noise in data. Furthermore, T (t, • ) : X → Y is a time-dependent forward operator, for short denoted by T t , that models how an image at time t gives rise to data in absence of noise or measurement errors (e.g., a stack of ray transforms with time varying data acquisition parameters in CT or the forward model in PET with time dependent attenuation, etc.) [START_REF] Natterer | The Mathematics of Computerized Tomography[END_REF].

To proceed, we further specify the form of the spatiotemporal image f (t, • ) by making use of deformable templates from shape theory.

2.2. Spatiotemporal inverse problem and shape theory. Shape theory seeks to develop quantitative tools to study shapes and their variability. The theory can be traced back to work by D'Arcy Thompson [START_REF] Thompson | On Growth and Form[END_REF]. The underlying idea is that shapes of objects are represented as a deformation that acts on a template [START_REF] Trouvé | Shape spaces[END_REF]. Hence, the template is the "shape invariant" part of the object whereas the set of deformations model the various shapes that can be generated from the template. Shape similarity between two objects can then be quantified as the "cost" of deforming one object into the other by means of a minimal deformation in the set of deformations. Further details are given in [START_REF] Younes | Shapes and Diffeomorphisms[END_REF][START_REF] Grenander | Pattern Theory. From Representation to Inference[END_REF][START_REF] Miller | Hamiltonian systems and optimal control in computational anatomy: 100 years since D'arcy Thompson[END_REF].

Bearing in mind the above, we separate the spatial and temporal components of a spatiotemporal image as (2.2) f (t, • ) := W(φ t , I) for some φ t ∈ G and I ∈ X .

Here I : Ω → R (template) is the time-independent spatial component, W : G × X → X is a deformation operator the models how a deformation parameter φ t ∈ G deforms the an image, and φ t : Ω → Ω is the temporal evolution of the deformation parameter.

We will assume deformations are given by diffeomorphic group actions, so G is a suitable subgroup of the diffeomorphic group on Ω and W is given by the group action of G on X , i.e., W(φ t , I) := φ t . I. The spatiotemporal inverse problem can then be written as

(2.3) g(t, • ) = T t (φ t . I) + g noise (t, • ) for t ∈ [0, 1].
Notice that f (t, • ) = φ t . I is the spatiotemporal image at time t generated from the template I by the diffeomorphism φ t . Hence, the above inverse problem calls for jointly recovering the (time-independent) template I and the curve of diffeomorphisms t → φ t . (2.4)

φ t . I := I • φ -1 t .
This manuscript is for review purposes only.

An alternative group action corresponds to mass-preserving deformations [START_REF] Chen | Indirect Image Registration with Large Diffeomorphic Deformations[END_REF][START_REF] Younes | Shapes and Diffeomorphisms[END_REF]:

(2.5) φ t . I := D(φ -1 t ) I • φ -1 t .
Here, D(φ) denotes the determinant of the Jacobian of φ. This group action adjusts the intensity values but preserves the total mass. can be seen as a temporal version of the indirect image registration problem in [START_REF] Chen | Indirect Image Registration with Large Diffeomorphic Deformations[END_REF].

Inverse problem of the type (2.3) in imaging applications are often ill-posed. A flexible framework for regularizing many of these problems is through a variational formulation [START_REF] Scherzer | Variational Methods in Imaging[END_REF].

The idea is to add regularization functional that penalizes a maximum likelihood solution and thereby acts as a stabiliser by preventing overfitting. The variational formulation for (2.3) reads as

(2.6) min

I∈X φt∈G 1 0 D T t (φ t . I), g(t, • ) + µ 2 R 2 (φ t ) dt + µ 1 R 1 (I) for fixed µ 1 , µ 2 > 0.
In the above, µ 1 , µ 2 are the regularization parameters that must be chosen depending on the noise level in data. Next, D : Y × Y → R + above is the data discrepancy functional that quantifies the mismatch in data space Y . It is often taken as a suitable affine transform of the negative data likelihood for the data, so minimising it amounts to computing a maximim likelihood estimator. If one has data with additive Gaussian noise, then it is given by the squared L 2 -norm:

(2.7) D g 1 , g 2 := g 1 -g 2 2 2 for g 1 , g 2 ∈ Y .
Likewise, Poisson distributed data leads to the Kullback-Leibler (KL) divergence:

(2.8) D g 1 , g 2 := M g 1 (y) -g 2 (y) log g 1 (y) dy for g 1 , g 2 ∈ Y .
Moreover, the spatial regularization R 1 : X → R + introduces stability by encoding priori knowledge about the ground truth image. It is frequently based on a roughness prior given as an L p -norm of the gradient magnitude or L 1 -norm of a suitable sparse representation.

Typically, if X := H 1 (Ω) then one takes the squared L 2 -norm of the gradient magnitude:

(2.9)

R 1 (f ) := ∇f 2 2
This manuscript is for review purposes only.

This choice is known to produce smooth images whereas selecting the L 1 -norm of the gradient magnitude for X := BV (Ω) is better at preserving edges (total variation (TV) regularization)

[45]:

(2.10)

R 1 (f ) := ∇f 1
What remains is to describe how to generate the curve t → φ t and to select the shape regularization R 2 : G → R + . An overall difficulty is that diffeomorphisms G do not form a vector space. One option is to try re-parametrising the deformations using vector space elements. As we shall see next, this can be done within the LDDMM framework.

3. A variational model for joint image reconstruction and motion estimation. This section introduces a new variational model of the framework (2.6) based on LDDMM. We begin by recalling the basic principles of LDDMM.

3.1. The LDDMM framework. The LDDMM framework outlined here offers a powerful machinery for generating a flow of diffeomorphisms by means of a velocity field. One advantage is that the set of velocity fields form a vector space, so we have in this way re-parametrised deformations by vector space elements.

The idea in LDDMM is to consider a sequence of infinitesimally small vector field perturbations to the identity mapping. These vector fields can be seen as an instantaneous velocity field. Under certain regularity, the composition of such small deformations in the limit generates a flow of diffeomorphisms given as the solution to an ordinary differential equation (ODE) [START_REF] Younes | Shapes and Diffeomorphisms[END_REF]. More precisely, given a velocity field ν : [0, 1] × Ω → R n , a flow t → φ t is generated by the following ODE:

(3.1)    ∂ t φ t (x) = ν t, φ t (x) φ 0 (x) = x for x ∈ Ω and 0 ≤ t ≤ 1.
Note that φ 0 = Id, i.e., the flow starts at the identity mapping. If the velocity field ν is sufficiently regular (see Definition 3.1), then the above ODE has a well-defined solution that is a diffeomorphism at each time point, i.e., (3.1) defines a flow of diffeomorphisms.

Definition 3.1 (Admissible space [START_REF] Younes | Shapes and Diffeomorphisms[END_REF]). A Hilbert space

V ⊂ C 1 0 (Ω, R n ) is admissible if it is (canonically) embedded in C 1 0 (Ω, R n ) with the • 1,∞ norm, i.e.
, there exists a constant

C > 0 such that v 1,∞ ≤ C v V for all v ∈ V . In the above, v 1,∞ := v ∞ + Dv ∞ for v ∈ C 1 0 (Ω, R n ).
We now consider velocity fields that are L p in time and in an admissible space at every point in time, i.e.,

(3.2) L p ([0, 1], V ) := ν : ν(t, •) ∈ V and ν L p ([0,1],V ) < ∞ for 1 ≤ p ≤ ∞
with the associated norm

ν L p ([0,1],V ) := 1 0 ν(t, •) p V dt 1/p .
This manuscript is for review purposes only.

To simplify notation, let L p V (Ω) denote L p ([0, 1], V ) and note in particular that L 2 V (Ω) is a Hilbert space with inner product

ν, η L 2 V (Ω) = 1 0 ν(t, •), η(t, • ) V dt for ν, η ∈ L 2 V (Ω).
The following theorem states that admissible velocity field generates a flow of diffeomorphisms.

Theorem 3.2 ( [START_REF] Younes | Shapes and Diffeomorphisms[END_REF][START_REF] Bruveris | Geometry of image registration: The diffeomorphism group and momentum maps[END_REF]). Let V be an admissible Hilbert space and ν ∈ L 2 V (Ω) be a velocity field. Then the ODE in (3.1) admits a unique solution

φ ν ∈ C 1 0 ([0, 1] × Ω, Ω), such that for t ∈ [0, 1], the mapping φ ν t : Ω → Ω is a C 1 -diffeomorphism on Ω.
Hence, we define

(3.3) G V := φ : φ = φ ν 0,1 for some ν ∈ L 2 V (Ω) ,
which becomes a sub-group of the group of diffeomorphisms by Theorem 3.2 whenever V is admissible. Next, if φ ν t denotes the solution to the ODE in (3.1) with given ν ∈ L 2 V (Ω), then

(3.4) φ ν s,t := φ ν t • (φ ν s ) -1 for 0 ≤ t, s ≤ 1
Furthermore, φ ν 0 = Id and by (3.4) we also get

(3.5) φ ν t = φ ν 0,t , (φ ν t ) -1 = φ ν t,0 .
Next, several important properties about G V are summarised in the following theorem.

Theorem 3.3 ([60, 11]). Let V be an admissible Hilbert space, G V is defined as in (3.3),

and d G V : G V × G V → R + is defined as (3.6) d G V (φ, ψ) := inf ν∈L 2 V (Ω) ψ=φ•φ ν 0,1 ν L 2 V (Ω) for φ, ψ ∈ G V .
Then G V is a group for the composition of functions and G V is a complete metric space under the metric d G V . Furthermore, for each φ, ψ ∈ G V , there exists ν ∈ L 2 V (Ω) satisfying

ψ = φ • φ ν 0,1 , i.e., d G V (φ, ψ) = ν L 2 V (Ω) .
The metric d G V can now be used to define an energy term for regularizing the registration of a template image I 0 to a target image I 1 :

(3.7) inf φ∈G V φ . I 0 -I 1 2 L 2 (Ω) + µ d G V (Id, φ) 2
for fix regularization parameter µ > 0.

One can show that the infimum in (3.6) is reached, so it can be replaced by a minimum. Next one can also show [START_REF] Younes | Shapes and Diffeomorphisms[END_REF]Lemma 11.3] that (3.7) is equivalent to the following variational model that optimizes over velocity fields instead of diffeomorphisms:

(3.8) min

ν∈L 2 V (Ω) φ ν 0,1 .I 0 -I 1 2 L 2 (Ω) + µ 1 0 ν(t, •) 2 V dt with φ ν 0,1 ∈ G V as in (3.4).
This manuscript is for review purposes only.

In conclusion, the regularization term for image registration by LDDMM is

(3.9) R(φ) := d G V (Id, φ) 2 = 1 0 ν(t, •) 2 V dt, where ν above minimizes ν → d G V (Id, φ) = ν L 2 V (Ω)
where φ = Id •φ ν 0,1 (such a minimiser exists due to Theorem 3.3).

A final remark concerns the choice of V . If we choose it as a reproducing kernel Hilbert space (RKHS) with a symmetric and positive-definite reproducing kernel, then V becomes an admissible Hilbert space [START_REF] Bruveris | Geometry of image registration: The diffeomorphism group and momentum maps[END_REF]. Using such vector fields is also advantageous from a computational point of view as shown in [START_REF] Chen | Indirect Image Registration with Large Diffeomorphic Deformations[END_REF]. We will therefore assume that V is henceforth an admissible RKHS.

Spatiotemporal reconstruction with LDDMM.

The aim here is to formulate a special case of (2.6) for solving the spatiotemporal inverse problem in (2.3). Let φ t in (2.6) be generated by the flow equation (3.1) as in LDDMM, i.e., φ t = φ ν 0,t for some velocity field ν.

If ν ∈ L 2 V (Ω), then φ ν t is a diffeomorphism on Ω by Theorem 3.2. Consequently, combining
Theorem 3.3 with (3.9) implies that a shape regularizer R 2 for the temporal deformation given by φ ν 0,t in (2.6) can be defined as

(3.10) R 2 (φ ν 0,t ) := t 0 ν(τ, •) 2 V dτ.
The variational formulation in (2.6) for solving the spatiotemporal inverse problem in (2.3) now reads as min

I∈X ν∈L 2 V (Ω) 1 0 D T t (φ ν 0,t . I), g(t, • ) + µ 2 t 0 ν(τ, •) 2 V dτ dt + µ 1 R 1 (I)
s.t. φ ν 0,t solves (3.1).

(3.11)

Note that (3.11) is an ODE constrained optimization problem, which henceforth is referred to as time-continuous spatiotemporal reconstruction with LDDMM.

Remark 3.4. Changing the order of integration for the second term in (3.11) yields the following equivalent formulation:

min I∈X ν∈L 2 V (Ω) 1 0 D T t (φ ν 0,t . I), g(t, • ) + µ 2 (1 -t) ν(t, •) 2 V dt + µ 1 R 1 (I) s.t. φ ν 0,t solves (3.1).
(3.12)

3.2.1. Basic properties of optimal velocity fields. The aim here is to characterize temporal behaviour of velocity fields that solve (3.11) when V is a RKHS. For this purpose, we introduce some notation that will simplify the expressions:

(3.13) D gt f := D T t f , g(t, • ) for f ∈ X with given g(t, • ) ∈ Y .
This manuscript is for review purposes only.

By Theorem A.4 (in Appendix A), the optimal velocity field in (3.11) satisfies

(3.14) ν(t, • ) = 1 2µ 2 (1 -t) 1 t K ∇(φ ν 0,t . I) D(φ ν t,τ ) ∇D gτ (φ ν 0,τ . I)(φ ν t,τ ) dτ for 0 ≤ t < 1.
Here,

K(ϕ) = Ω K( • , y)ϕ(y) dy, with K : Ω × Ω → M n×n + denoting the reproducing kernel if V is a RKHS.
By the above, the optimal velocity field can be seen as a time-averaging. Next, ν(1, • ) is well-defined at t = 1 and reads as

(3.15) ν(1, •) = 1 2µ 2 K ∇(φ ν 0,1 . I) ∇D g 1 (φ ν 0,1 . I) .
Similarely, by (A.9) and (A.7) we get

(3.16) ν(0, •) = 1 2µ 2 K ∇I 1 0 D(φ ν 0,t ) ∇D gt (φ ν 0,t . I)(φ ν 0,t ) dt .
Hence, the optimal velocity field in (3.11) is averagely distributed w.r.t. time t, and nonvanishing at both the initial and the end time points.

Control theoretic formulation.

The aim here is to state an equivalent optimal control formulation of (3.11) in terms of a partial differential equation (PDE) constrained optimization problem. This makes it easier to compare our proposed approach in (3.11) against PDE based approaches for registration, as we do in subsection 3.3.

Theorem 3.5. Let X be a space of real valued functions that are sufficiently smooth, e.g., the space of functions with bounded variation given in subsection 3.3.1. Assume next that I ∈ X and define f : [0, 1] × Ω → R as (3.17) f (t, • ) := φ ν 0,t . I for 0 ≤ t ≤ 1 with φ ν 0,t given by (3.4).

Assume furthermore that f (t, • ) ∈ X for any 0 ≤ t ≤ 1. Then, (3.11) is under the group action in (2.4) (geometric deformation) equivalent to

min f (0, • )∈X ν∈L 2 V (Ω) 1 0 D T t f (t, • ) , g(t, • ) + µ 2 t 0 ν(τ, •) 2 V dτ dt + µ 1 R 1 f (0, • ) s.t. ∂ t f (t, • ) + ∇f (t, •), ν(t, •) R n = 0. (3.18)
Likewise, (3.11) is under the group action in (2.5) (mass-preserving deformation) equivalent to

min f (0, • )∈X ν∈L 2 V (Ω) 1 0 D T t f (t, • ) , g(t, • ) + µ 2 t 0 ν(τ, •) 2 V dτ dt + µ 1 R 1 f (0, • ) s.t. ∂ t f (t, • ) + ∇• f (t, • ) ν(t, • ) = 0. (3.19)
This manuscript is for review purposes only.

Proof. First we consider the geometric deformation in (2.4), so (3.17) reads as

(3.20) f (t, • ) = I • (φ ν 0,t ) -1 for 0 ≤ t ≤ 1. Obviously, f (0, • ) = I, f (1, • ) = I • (φ ν 0,1 ) -1
, and

(3.21) f t, φ ν 0,t = I for 0 ≤ t ≤ 1.
Differentiating (3.21) w.r.t. time t leads to

∂ t f t, φ ν 0,t + ∇f t, φ ν 0,t , ν t, φ ν 0,t R n = 0.
The above is the PDE constraint in (3.18), so a solution to (3.11) also solves (3.18).

We now consider the reverse implication, i.e., show that a solution to (3.18) also solves (3.11). Suppose that f and ν solve (3.18). Define the diffeomorphism ψ t that solves the ODE in (3.1) with the above given ν. Since f satisfies PDE constraint in (3.18), we get

d dt f (t, ψ t ) = ∂ t f (t, ψ t ) + ∇f (t, ψ t ), ν(t, ψ t ) R n = 0.
Hence, t → f t, ψ t is constant and in particular we have

f t, ψ t ≡ f 0, ψ 0 = f (0, • ). If f (0, • ) := I and ψ t := φ ν 0,t , then f t, • = I • (φ ν 0,t ) -1 .
Hence a solution to (3.18) also solves (3.11).

Using a mass-preserving deformation (2.5) as group action in (3.17) results in

(3.22) f (t, • ) = D (φ ν 0,t ) -1 I • (φ ν 0,t ) -1 for 0 ≤ t ≤ 1.
We then get that f (0,

• ) = I and f (1, • ) = D (φ ν 0,1 ) -1 I • (φ ν 0,1 ) -1 .
The symmetry of the mass-preserving property furthermore yields

(3.23) D(φ ν 0,t ) f (t, • ) • φ ν 0,t = I for 0 ≤ t ≤ 1.
Finally, differentiating (3.23) w.r.t. t leads to the constraint in (3.19). Hence, a minimizer of (3.11) with the group action given by (2.5) is also a minimizer of (3.19). Similar to the case of geometric deformation, it is not difficult to prove the reverse implication.

Comparison with existing approaches.

Here we compare (3.11) against several existing approaches (i.e., optical flow based model, diffeomorphic motion models). problem where the constraint is given by a brightness constancy equation. When applied to the spatiotemporal inverse problem in (2.3), it reads as min

f (t, • )∈X ν(t, • )∈BV (Ω) 1 0 D T t f (t, • ) , g(t, • ) + µ 1 R 1 f (t, • ) + µ 2 ν(t, • ) BV dt s.t. ∂ t f (t, • ) + ∇f (t, • ), ν(t, • ) R n = 0, (3.24)
This manuscript is for review purposes only.

where | • | BV is the semi-norm on the space of functions with bounded variation:

BV (Ω) := u ∈ L 1 (Ω) : |u| BV < ∞ with |u| BV := sup ϕ∈C 1 0 (Ω,R n ) ϕ ∞≤1 Ω u(x) div ϕ(x) dx.
In particular, ν(t, • ) BV denotes the sum of the semi-norm of all the elements in ν(t, • ) [START_REF] Aubert | Mathematical Problems in Image Processing: Partial Differential Equations and the Calculus of Variations[END_REF].

We will use the reformulation in (3.18) to compare the above optical flow based model from [START_REF] Burger | A Variational Model for Joint Motion Estimation and Image Reconstruction[END_REF] with our proposed model. The use of the brightness constancy equation points to using the geometric deformation in (2.4) as a group action in (3.17), i.e., we assume (3.20) holds. It is easy to see that the constraints in (3.18) and (3.24) are equivalent. Hence, the optical flow constraint given by the brightness constancy equation in [START_REF] Burger | A Variational Model for Joint Motion Estimation and Image Reconstruction[END_REF] is equivalent to using diffeomorphisms generated by LDDMM that act through geometric deformation.

By comparison, the primary distinction between (3.18) and (3.24) relates to the selection of the regularization term w.r.t. vector field ν(t, •). In the model (3.24), one uses the TV functional, so the space of vector fields is assumed to be in BV (Ω), which allows for a vector field that is a piecewise-constant vector-valued function distributed on Ω. In contrast, in the model (3.18), the space of vector fields reside in an admissible Hilbert space. Hence, the vector field is a sufficiently smooth vector-valued function distributed on Ω. This guarantees an elastic diffeomorphic deformation, which to some extent mimics the underlying physical mechanisms of organ motion [START_REF] Gigengack | Motion correction in dual gated cardiac PET using mass-preserving image registration[END_REF][START_REF] Burger | A Hyperelastic Regularization Energy for Image Registration[END_REF].

In addition to the above, both approaches also differ in the selection of regularization term R 1 . In (3.18) one only poses restriction on the initial image f (0, • ), whereas in (3.24) the whole time trajectory t → f (t, • ) is regularized. The treatment in (3.18) is reasonable to some extent, since regularity properties are assumed to preserved. More precisely, f (t, • ) is contained in the same space for all t. For example, if f (0, • ) ∈ BV (Ω), then there exists a unique weak solution f (t, Next, by Theorem A.4 it is easy to see that the L 2 V (Ω)-norm minimizer of (3.25) satisfies

• ) = f (0, • ) • (φ ν 0,t ) -1 to
(3.26) ν(t, •) = 1 2µ 2 1 t K ∇(φ ν 0,t . I) D(φ ν t,τ ) ∇D gτ (φ ν 0,τ . I)(φ ν t,τ ) dτ for 0 ≤ t ≤ 1.
It is not difficult to see from (3.26) that the optimal velocity field is not averagely distributed w.r.t. time t. In addition, the minimizer w.r.t. variations of the template satisfies (3.27)

1 0 D(φ ν 0,t ) ∇D gt (φ ν 0,t . I)(φ ν 0,t ) dt = 0.
Combining (3.26) and (3.27) yields ν(0, •) = ν(1, •) = 0, i.e., the optimal velocity field according to the model (3.25) vanishes at both the initial and end time points.

To summarise, both (3.11) and the approach taken in [START_REF] Hinkle | 4D CT image reconstruction with diffeomorphic motion model[END_REF] address joint reconstruction and motion estimation. Both approaches model the latter as diffeomorphic deformations generated by velocity fields within the LDDMM framework. A difference is that the regularization of the velocity field in [START_REF] Hinkle | 4D CT image reconstruction with diffeomorphic motion model[END_REF] is equally weighted over the entire time trajectory. Next, the optimal velocity field vanishes at both the initial and end time points. This is not the case with the model in (3.11).

A variant to (3.11) based on (2.6) would regularise the template I (image at t = 0) instead of its entire time evolution φ ν 0,t . I. This leads to min

I∈X ν∈L 2 V (Ω) 1 0 D T t (φ ν 0,t . I), g(t, • ) + µ 1 R 1 (φ ν 0,t . I) + µ 2 t 0 ν(τ, •) 2 V dτ dt s.t. φ ν 0,t solves (3.1).
(3.28)

By Theorem A.4, the optimal velocity field for 0 ≤ t < 1 satisfies

(3.29) ν(t, • ) = 1 2µ 2 (1 -t) 1 t K ∇(φ ν 0,t . I) D(φ ν t,τ ) ∇S gτ (φ ν 0,τ . I)(φ ν t,τ ) dτ
where

S gt f := D T t f , g(t, • ) + µ 1 R 1 f for fixed g(t, • ) ∈ Y and f ∈ X . The corre-
sponding optimal template satisfies (3.30)

1 0 D(φ ν 0,t ) ∇S gt φ ν 0,t . I φ ν 0,t dt = 0.
Evidently, the above optimal velocity field is also a time average of the integrand. Next, ν(1, • ) is well-defined at t = 1 and reads as

ν(1, •) = 1 2µ 2 K ∇(φ ν 0,1 . I) ∇S g 1 (φ ν 0,1 . I) .
This manuscript is for review purposes only.

However, by (3.29) and (3.30) at t = 0 we have ν(0, •) = 0. Hence, the optimal velocity field in (3.28) is averagely distributed w.r.t. time t, but vanishing at the initial time point. In addition, as stated in subsection 3.3.1, we only need to regularize I rather than φ ν 0,t . I to some extent, since I and its time evolution reside in the same space. This proposed model has a simpler structure to implement.

To summarise, the comparative analysis points to several advantages that come with using (3.11) over alternative approaches.

4. Numerical implementation. The first part of the numerical implementation is to derive a time-discretized formulation of (3.11).

4.1. Time-discretized formulation. Assume the time sampling of data is uniform, i.e., there is a uniform partition of [0, 1] and data is acquired at time points {t i } N i=0 with t i = i/N for 0 ≤ i ≤ 1. We refer to this as the gating grid and the time-discretized formulation of the general spatiotemporal inverse problem in (2.1) is the task to recover t → f (t, • ) from data

g(t i , • ) ∈ Y where (4.1) g(t i , • ) = T t i f (t i , • ) + g noise (t i , • ) for i = 1, . . . , N .
Then, a time-discretized version of (3.11) reads as min

I∈X ν∈L 2 V (Ω) 1 N N i=1 D T t i (φ ν 0,t i . I), g(t i , • ) + µ 2 t i 0 ν(τ, •) 2 V dτ + µ 1 R 1 (I) s.t. φ ν 0,t solves (3.1). (4.2)
Remark 4.1. The time-discretized version (4.2) can be also written such that the image in the first gate is the template:

min I∈X ν∈L 2 V (Ω) 1 N + 1 N i=0 D T t i (φ ν 0,t i . I), g(t i , • ) + µ 2 t i 0 ν(τ, •) 2 V dτ + µ 1 R 1 (I) s.t. φ ν 0,t solves (3.1).
Since (4.2) contains highly coupled arguments, it is difficult to jointly solve for the template I and the velocity field ν. An approach that is computationally more feasible is to compute I and ν through an intertwined iterative scheme. More precisely, if the velocity field ν is given then the spatiotemporal reconstruction problem (4.2) reduces to a static image reconstruction problem: Conversely, if the template I is fixed, then (4.2) reduces to a sequentially indirect image registration problem where we seek the velocity field ν from time-series data that are indirect observations of the target. Formalising the above, we try to solve (4.2) through the

following intertwined iterative scheme (I k , ν k ) ∈ X × L 2 V (Ω): (4.3)      I k ∈ arg min I∈X J ν k-1 (I) ν k ∈ arg min ν∈L 2 V (Ω) E I k (ν)
This manuscript is for review purposes only.

where the functionals J ν : X → R, for given ν ∈ L 2 V (Ω), and E I : L 2 V (Ω) → R, for given I ∈ X , are defined as

J ν (I) := 1 N N i=1 D T t i (φ ν k-1 0,t i . I), g(t i , • ) + µ 1 R 1 (I) (4.4) E I (ν) := 1 N N i=1 D T t i (φ ν 0,t i . I k ), g(t i , • ) + µ 2 t i 0 ν(τ, •) 2 V dτ (4.5) s.t. φ ν 0,t solves (3.1).
To proceed, we consider the special case when the data discrepancy term is the squared

L 2 -
I∈X ν∈L 2 V (Ω) 1 N N i=1 T t i φ ν 0,t i . I -g(t i , • ) 2 2 + µ 2 t i 0 ν(τ, •) 2 V dτ + µ 1 ∇I 1 s.t. φ ν 0,t solves (3.1). (4.6)
Correspondingly, (4.5) become

J ν (I) = 1 N N i=1 T t i φ ν 0,t i . I -g(t i , • ) 2 2 + µ 1 ∇I 1 (4.7) E I (ν) = 1 N N i=1 T t i φ ν 0,t i . I -g(t i , • ) 2 2 + µ 2 t i 0 ν(τ, •) 2 V dτ (4.8)
Inserting the above into (4.3) yields the following intertwined scheme for solving (4.6): (4.9)

     I k := arg min I∈X J ν k-1 (I) with J ν k-1
given by (4.7)

ν k := arg min ν∈L 2 V (Ω)
E I k (ν) with E I k given by (4.8).

Template reconstruction.

We henceforth consider geometric deformations where the deformation operator is given as φ 0,t . I := I • φ -1 0,t = I • φ t,0 . This is a common choice for image registration [START_REF] Blume | Joint image and motion reconstruction for PET using a B-spline motion model[END_REF][START_REF] Brehm | Self-adapting cyclic registration for motion-compensated cone-beam CT in image-guided radiation therapy[END_REF][START_REF] Liu | 5D respiratory motion model based image reconstruction algorithm for 4D cone-beam computed tomography[END_REF][START_REF] Burger | A Variational Model for Joint Motion Estimation and Image Reconstruction[END_REF].

Next, consider the problem of computing a minimiser to J ν k-1 (static reconstruction) in (4.7). This is a non-smooth TV-2 minimization problem that is convex when the forward operator T t i : X → Y is linear. Solving it by a first order method, like a gradient descent scheme as outlined in Algorithm 4.1, requires smoothing the TV-component, i.e., we seek to solve (4.10) min

I∈X 1 N N i=1 T t i I • φ ν t i ,0 -g(t i , • ) 2 2 + µ 1 Ω |∇I(x)| 2, dx ,
This manuscript is for review purposes only.

where |∇I(x)| 2, := i ∂ i I(x) 2 + with > 0 small, e.g., = 10 -12 . This is a frequently used modification of the TV functional in imaging applications [START_REF] Sidky | Accurate image reconstruction from few-views and limited-angle data in divergent-beam CT[END_REF][START_REF] Chen | Gradient-flow-based semi-implicit finite-element method and its convergence analysis for image reconstruction[END_REF][START_REF] Chen | Image reconstruction and scan configurations enabled by optimization-based algorithms in multispectral CT[END_REF]. A gradient descent scheme for (4.10) assuming a linear T t i reads as (4.11)

I k+1 = I k -α k 2 N N i=1 D φ ν 0,t i T * t i T t i (I k • φ ν t i ,0 ) -g(t i , • ) (φ ν 0,t i ) + µ 1 ∇ * ∇I k |∇I k | 2, .
Here, α k is the stepsize for the k-th iteration and The time interval [0, 1] is subdivided uniformly into M N parts thereby forming a discretized time grid that is given as τ j = j/(M N ) for j = 0, 1, . . . , M N . Evidently, τ iM = t i for i = 0, 1, . . . , N , so each subinterval [t i , t i+1 ] is segmented into M even parts. The M is named as a factor of discretized time degree. If M = 1, then τ i = t i , implying that the discretized time grid is consistent with the gating grid. Note also that the different subintervals of the gating grid can be discretized adaptively according to the degree of motions.

Within a short-time interval one can approximate the diffeomorphic deformation with linearized deformations [START_REF] Öktem | Shape based image reconstruction using linearised deformations[END_REF]. More precisely, let s = τ j , t = τ j-1 and τ j+1 in (4.13). Then, the small deformations φ ν τ i ,τ i-1 and φ ν τ i ,τ i+1 can be approximated by

(4.14) φ ν τ j ,τ j-1 ≈ Id - 1 M N ν(τ j , • ), and φ ν τ j ,τ j+1 ≈ Id + 1 M N ν(τ j , • ).
This manuscript is for review purposes only.

Moreover, by (3.4) we get φ ν τ j ,0 = φ ν τ j-1 ,0 • φ ν τ j ,τ j-1 , which combined with (4.14) yields (4.15)

φ ν τ j ,0 ≈ φ ν τ j-1 ,0 • Id - 1 M N ν(τ j , • ) for j = 1, 2, . . . , M N .
This yields the following estimate for I • φ ν τ j ,0 : (4.16)

I • φ ν τ j ,0 ≈ I • φ ν τ j-1 ,0 • Id - 1 M N ν(τ j , • ) for j = 1, 2, . . . , M N
and I • φ ν τ 0 ,0 = I. Similarly, (3.4) also implies φ ν τ j ,t i = φ ν τ j+1 ,t i • φ ν τ j ,τ j+1 for i ≥ 1, which combined with (4.14) gives the following approximation:

(4.17)

φ ν τ j ,t i ≈ φ ν τ j+1 ,t i • Id + 1 M N ν(τ j , • ) for j = iM -1, iM -2, . . . , 0
and φ ν t i ,t i = Id. To summarize, the deformation between two images at adjacent points in the discretized time grid is approximately represented by a linearized deformation.

Computing mass-preserving deformations.

Expression for the gradient of the data discrepancy term for the mass-preserving group action in (2.5) involves terms of the type

D φ ν 0,t i T * t i T t i I • φ ν t i ,0 -g(t i , • ) • φ ν 0,t i for i ≥ 1.
Starting with the Jacobian determinant, by (4.17) we get

(4.18) D(φ ν τ j ,t i ) ≈ 1 + 1 M N div ν(τ j , • ) D(φ ν τ j+1 ,t i ) • Id + 1 M N ν(τ j , • )
for j = iM -1, iM -2, . . . , 0 and where D(φ ν t i ,t i ) = 1. Next, (4.17) also yields the following approximation:

(4.19) T * t i T t i I • φ ν t i ,0 -g(t i , • ) • φ ν τ j ,t i ≈ T * t i T t i I • φ ν t i ,0 -g(t i , • ) • φ ν τ j+1 ,t i • Id + 1 M N ν(τ j , • ) for j = iM -1, iM -2, . . . , 0. For simplicity, let (4.20) η I,ν τ,t = D φ ν τ,t T * t T t I • φ ν t,0 -g(t, • ) • φ ν τ,t .
Then multiplying (4.18) by (4.19), and using (4.20), η I,ν 0,t i for i ≥ 1 is computed by

(4.21) η I,ν τ j ,t i ≈ 1 + 1 M N div ν(τ j , • ) η I,ν τ j+1 ,t i • Id + 1 M N ν(τ j , • ) for j = iM -1, iM -2, . . . , 0 with η I,ν t i ,t i = T * t i T t i I • φ ν t i ,0 -g(t i , • ) .
Based on the above derivations, the concrete implementation is given as the gradient descent scheme in Algorithm 4.1.

This manuscript is for review purposes only.

Algorithm 4.1 Gradient descent scheme for minimizing I → J ν (I) in (4.7)

1: Initialize: 2: k ← 0. 3: t i ← i N for i = 0, 1, . . . , N . 4: τ j ← j
M N for j = 0, 1, . . . , M N . 5: Given ν. 6: I k ← I 0 . Here I 0 is a given initial template. 7: Spatial regularization parameter µ 1 > 0. 8: Error tolerance I > 0, stepsize α k = α > 0, and iteration number K I > 0. 9: Loop:

10: Compute I k • φ ν τ j ,0 for 1 ≤ j ≤ M N by I k • φ ν τ j ,0 ← I k • φ ν τ j-1 ,0 • Id - 1 M N ν(τ j , • ) with I k • φ ν 0,0 = I k . 11: Update η I k ,ν t i ,t i for 1 ≤ i ≤ N by η I k ,ν t i ,t i ← T * t i T t i (I k • φ ν t i ,0 ) -g(t i , • )
.

12:

Compute η I k ,ν 0,t i for 1 ≤ i ≤ N by η I k ,ν τ j ,t i ← 1 + 1 M N div ν(τ j , • ) η I k ,ν τ j+1 ,t i • Id + 1 M N ν(τ j , • )
for j = iM -1, iM -2, . . . , 0.

13:

Evaluate I k+1 by

I k+1 ← I k -α 2 N N i=1 η I k ,ν 0,t i + µ 1 ∇ * ∇I k |∇I k | 2,
.

14:

If I k+1 -I k > I and k < K I , then k ← k + 1, goto Loop. 15: Output I k+1 .

Velocity field estimation.

The aim here is to provide an algorithm for minimising E I in (4.8), which amounts to sequential indirect image registration. We will for this purpose use a gradient descent scheme of the form

(4.22) ν k+1 = ν k -β k ∇E I (ν k ).
Here, β k is the step-size in the k-th iteration, and

∇E I (ν k ) ∈ L 2 V (Ω) is calculated by (B.1).
The challenge here lies in the computation of this gradient (subsection 4.3.1) and the final algorithm for the gradient descent scheme in (4.22) is given in Algorithm 4.2.

This manuscript is for review purposes only.

4.3.1.

Computing the gradient ∇E I . We begin by introduce notations

(4.23) h I,ν τ,t :=    η I,ν τ,t 0 ≤ τ ≤ t ≤ 1 0 t < τ,
and ν τ,t :=

   ν(τ, • ) 0 ≤ τ ≤ t ≤ 1, 0 t < τ.
Theorem B.1 gives an expression for ∇E I where the kernel function K : Ω × Ω → M n×n + is evaluated on points that do not move as iteration proceeds. By choosing a translation invariant kernel and points on a regular grid in Ω, we can use FFT-based convolution scheme to efficiently evaluate the velocity field at each iteration. This is computationally more feasible than letting the kernel depend on points that move in time as in the shooting method [START_REF] Miller | Geodesic shooting for computational anatomy[END_REF][START_REF] Vialard | Diffeomorphic 3D image registration via geodesic shooting using an efficient adjoint calculation[END_REF].

In what follows, we write out the explicit derivations for computing ∇E I . As derived in subsection 4.2.1, I • φ ν τ j ,0 can be approximated by (4.16). The key step is now to update h I,ν τ j ,t i for {i : t i ≥ τ j } in (B.1). We know, by (4.23), (4.24) h I,ν τ j ,t i = η I,ν τ j ,t i for t i ≥ τ j .

By using (4.21) for 1 ≤ i ≤ N , we can compute h I,ν τ j ,t i as

(4.25) h I,ν τ j ,t i ≈ 1 + 1 M N div ν(τ j , • ) h I,ν τ j+1 ,t i • Id + 1 M N ν(τ j , • ) for j = iM -1, iM -2, . . . , 0 and with h I,ν t i ,t i = T * t i T t i (I • φ ν t i ,0 ) -g(t i , • ) . Hence, at t = τ j , we get from (B.1) that (4.26) ∇E I (ν)(τ j , x) = - 2 N {i≥1:t i ≥τ j } Ω K(x, y)∇ I • φ ν τ j ,0 (y)h I,ν τ j ,t i (y) dy -µ 2 ν(τ j , x)
for 0 ≤ j ≤ M N and x ∈ Ω. In particular, for j = M N (i.e., τ j = 1) we have

∇E I (ν)(1, x) = - 2 N Ω K(x, y)∇ I • φ ν 1,0 (y)h I,ν 1,1 (y) dy -µ 2 ν(1, x) .
Remark 4.3. It is easy to verify that the optimal solution of the time-discretized version of the proposed model is consistent with that of the time-continuous one. Here the concept of consistency is that the time-discretized solution converges to the time-continuous solution for increasingly fine time discretization. This is however not the case for the diffeomorphic motion model (3.25) in [START_REF] Hinkle | 4D CT image reconstruction with diffeomorphic motion model[END_REF]. As an example, at τ j = 1, the optimal velocity field of the time-discretized problem in [START_REF] Hinkle | 4D CT image reconstruction with diffeomorphic motion model[END_REF] satisfies

ν(1, x) = 1 µ 2 Ω K(x, y)∇ I • φ ν 1,0 (y)h I,ν 1,1 (y) dy.
On the other hand, as derived in subsection 3.3.2, the optimal velocity field at t = 1 of its time-continuous problem satisfies ν(1, x) = 0. This obviously causes inconsistencies and our consistent approach is an advantage compared to the approach in [START_REF] Hinkle | 4D CT image reconstruction with diffeomorphic motion model[END_REF].

Finally, Algorithm 4.2 outlines the procedure for computing the gradient descent scheme (4.22) that makes use of the above derivations.

This manuscript is for review purposes only.

Algorithm 4.2 Gradient descent scheme for minimizing ν → E I (ν) in (4.8)

1: Initialize: 2: k ← 0. 3: t i ← i N for i = 0, 1, . . . , N . 4: τ j ← j
M N for j = 0, 1, . . . , M N . 5: Fixed I. 6: ν k (τ i ) ← ν 0 (τ i ), where ν 0 is a given initial velocity field. 7: Fixed kernel function K( • , • ). 8: Shape regularization parameter µ 2 > 0. 9: Error tolerance ν > 0, stepsize β k = β > 0, and maximum iterations K ν > 0. 10: Loop:

11: Update I • φ ν k τ j ,0 for 1 ≤ j ≤ M N by I • φ ν k τ j ,0 ← I • φ ν k τ j-1 ,0 • Id - 1 N ν k (τ j , • ) with I • φ ν k 0,0 = I.
12:

Update h I,ν k t i ,t i for 1 ≤ i ≤ N by h I,ν k t i ,t i ← T * t i T t i (I • φ ν k t i ,0 ) -g(t i , • ) . 13: Compute h I,ν k τ j ,t i for 1 ≤ i ≤ N by h I,ν k τ j ,t i ← 1 + 1 M N div ν(τ j , • ) h I,ν k τ j+1 ,t i • Id + 1 M N ν k (τ j , • )
for j = iM -1, iM -2, . . . , 0.

14:

Compute ∇ ν E I (ν k )(τ j , • ) (using FFT for convolutions) by

∇E I (ν k )(τ j , x) ← - 2 N {i≥1:t i ≥τ j } Ω K(x, y) ∇ I • φ ν k τ j ,0 (y)h I,ν k τ j ,t i (y) dy -µ 2 ν k (τ j , x)
for 0 ≤ j ≤ M N .

15:

Update ν k (τ j , • ) for 0 ≤ j ≤ M N by:

ν k+1 (τ j , • ) ← ν k (τ j , • ) -β ∇E I (ν k )(τ j , • ).
16:

If ν k+1 -ν k > ν and k < K ν , then k ← k + 1, goto Loop. 17: Output ν k+1 .
This manuscript is for review purposes only.

4.4. Alternating template reconstruction and velocity field estimation. As described in the beginning of section 4, we aim to solve (4.6) by an iterative scheme where iterates for template image I and velocity field ν are updated in an alternating manner as in (4.9). Hence, at each iterative step requires solving two sub-problems, one for updating I given ν (Algorithm 4.1) and the other for updating ν given I (Algorithm 4.2).

The algorithms for solving the two sub-problems are iterative, so there are inner iterations for each outer iterative step that update the template and velocity field. In order to have a computationally feasible approach, we limit the number of inner iterations to one. The final algorithm for recovering the template and velocity field is presented in Algorithm 4.3. The proposed model is non-convex so a convergence analysis of the algorithm towards a global optima is currently too difficult to obtain.

1: Initialize: 2: Given M, N . 3: k ← 0. 4: t i ← i
N for i = 0, 1, . . . , N . This subdivides the time [0, 1] uniformly into N parts. 5: τ j ← j M N for j = 0, 1, . . . , M N . This subdivides the time interval [0, 1] uniformly into M N parts. 6: Fixed kernel function K( • , • ). 7: Given regularization parameters µ 1 , µ 2 > 0. 8: I k ← I 0 , where the template is initialized. 9: ν k (τ i ) ← 0, where the velocity field is initialized to a zero velocity field. 10: Error tolerances I , ν > 0, stepsizes α k = α > 0, β k = β > 0, and maximum iteration number K > 0. 11: Loop:

If ν k+1 -ν k > ν or I k+1 -I k > I , and k < K, then k ← k + 1, goto Loop. 15: Output I k+1 , ν k+1 .
Complexity analysis. Here we consider the complexity, including computational cost and memory footprint, of Algorithm 4.3. The main computationally demanding steps at each iteration is located on lines 12-13 (actually lines 10-13 in Algorithm 4.1 and lines 11-15 in Algorithm 4.2), so we restrict our complexity analysis to these parts. For ease of description, we assume Ω ⊂ R 2 and the size of the image to be reconstructed is n × n pixels. On line 10 of Algorithm 4.1 and line 11 of Algorithm 4.2, we need to update I • φ ν τ j ,0 for j = 1, . . . , M N . Moreover, each of them should be used to compute the gradient of the objective functional on line 14 of Algorithm 4.2, so they need to be stored at hand. Hence, in these two steps, the computational cost is O(n 2 M N ) and the memory footprint is O(n 2 M N ).

For line 11 of Algorithm 4.1 and line 12 of Algorithm 4.2, the η I,ν t i ,t i (i.e., h I,ν t i ,t i ) need to be updated and then stored for i = 1, . . . , N . The computational cost is

O(n 2 N d N ), where N d
This manuscript is for review purposes only.

is the number of data points. Actually, the N d is at least proportional to the size of n, which is often √ 2nN v with N v denoting the number of views. Hence, the computational cost scales as O(n 3 N N v ). Since the calculation for the forward and backward projections is on the fly, the required space is not too demanding. Furthermore, on line 12 of Algorithm 4.1 and line 13 of Algorithm 4.2, for i = 1, . . . , N , the η I,ν τ j ,t i (i.e., h I,ν τ j ,t i ) need to be updated and stored for j from iM -1 to 0, then are used to compute the gradient of the objective functional for each time point on line 14 of Algorithm 4.2.

Therefore, the computational cost is O(n 2 M N 2 ). For lines 11-12 of Algorithm 4.1 and lines 12-13 of Algorithm 4.2, the memory footprint is O(n 2 M N 2 ).

For line 13 of Algorithm 4.1, we need to update I once, the computational cost is O(n 2 N ) and the memory footprint is O(n 2 N ). At each time point, the FFT is used to compute the gradient of the objective functional on line 14 of Algorithm 4.2. Hence the computational cost for this line is O(M N 2 n 2 log n). For line 15 of Algorithm 4.2, we need to update a vector field at each time point. Since a vector field would take twice more memory than a scalar field on 2D domain, we spend twice more computational cost to update that. Even so, the computational cost is O(n 2 M N ) and the memory footprint is O(n 2 M N ).

In summary, for Algorithm 4.3, the computational cost scales as O(n 3 N ) and the memory footprint scales as O(n 2 M N 2 ).

Numerical experiments.

In this section, the proposed method for joint image reconstruction and motion estimation is applied to parallel beam tomography with very sparse or highly noisy data in spatiotemporal (e.g., 2D + time) imaging. We use the intensity-preserving group action to consider the involved deformations. The algorithms were implemented in Python using the Operator Discretization Library (http://github.com/odlgroup/odl), which offers GPU parallellized routines for evaluating the ray transform and it adjoint. The source code is available from https://github.com/chongchenmath/SpatiotemporalImaging. The numerical experiments ran on ThinkStation with Xeon E5-2620 v4 2.10 GHz CPU, 64Gb RAM, TITAN Xp GPU, and Ubuntu 16.04 OS. Although not a complete evaluation, the experiments illustrate the performance of the proposed method.

5.1. Spatiotemporal 2D CT. Image reconstruction in classical static 2D CT amounts to inverting the 2D ray transform. In our spatiotemporal setting, we assume that the temporal variation comes from the object itself that undergoes motion. Hence, the forward operator T t : X → Y , which is the 2D ray transform, is not dependent on time t:

T (f )(ω, x) := R f (x + sω) ds for ω ∈ S 1 and x ∈ ω ⊥ .
Here, S 1 is the unit circle and (ω, x) encodes a line on R 2 with direction ω through x.

Moreover, consider V as the space of vector fields that is a RKHS with a reproducing The images of all gates are supported on Ω. At each gate, the noise-free data per view is

This manuscript is for review purposes only.

measured by evaluating the 2D parallel beam scanning geometry. Then the additive Gaussian white noise at varying levels is added onto the noise-free data, which leads to the noise data.

As in [START_REF] Chen | Indirect Image Registration with Large Diffeomorphic Deformations[END_REF], the noise level in data is quantified in terms of signal-to-noise ratio (SNR) defined in logarithmic decibel (dB).

5.2.

Test suites and results. The test suites seek to assess the performance against different noise levels, and the sensitivity against various selections of regularization parameters µ 1 , µ 2 , and kernel width σ. We also compare the proposed method to TV-based static reconstruction method.

5.2.1. Test suite 1: Overview performance. Here we consider a test for evaluating the overview performance. This test uses a multi-object phantom with five gates (i.e., N = 5).

The used phantom is shown in the last row of Figure 2, which is taken from [START_REF] Chen | Indirect Image Registration with Large Diffeomorphic Deformations[END_REF].

The image at each gate is consisting of six separately star-like objects with grey-values over [0, 1], which is digitized using 438 × 438 pixels. The images of all gates are supported on The regularization parameters (µ 1 , µ 2 ) are selected as (0.05, 10 -7 ) for data noise level

This manuscript is for review purposes only. Apart from the visual perception, the reconstruction is quantitatively compared using structural similarity (SSIM) and peak signal-to-noise ratio (PSNR), which is frequently used to evaluate image quality [START_REF] Wang | Image quality assessment: From error visibility to structural similarity[END_REF]. The SSIM and PSNR values are tabulated in 2.

Comparison against static TV-regularized reconstruction. It is well-known that tomographic reconstruction by TV-regularization outperforms other methods, such as filtered back projection (FBP), when the gradient of the image is sparse. This is furthermore especially notable when data is under-sampled. In our tests we use a phantom (ground truth image) that has sparse gradient, so comparing against static TV-regularized reconstruction pitches our approach against one of the best static reconstruction methods.

For static TV-regularized reconstruction we disregard any temporal evolution, which is equivalent to simplify the spatiotemporal problem into a static problem. The whole tomographic data set will then have 60 projection views. The regularization parameter for static TV-regularized reconstruction is selected depending on the SNR of data in the same way as for spatiotemporal reconstruction. In addition, the corresponding SSIM and PSNR values are listed in Table 2 and compared to Table 1, the values of SSIM and PSNR for static TV-regularization are lower than those with the proposed method.

This manuscript is for review purposes only. 3 for the images. Each entry has two values, where the upper is the value of SSIM and the bottom is the value of PSNR.

Test suite 2: Sensitivity against selections of regularization parameters.

To solve the proposed model, three regularization parameters µ 1 , µ 2 and σ need to be selected.

Hence the sensitivity test should be concerned against the selections of these parameters.

As shown in the last row of Figure 5, a heart phantom with four gates (i.e., N = 4) is used in this test, which is originated from [START_REF] Grenander | Pattern Theory. From Representation to Inference[END_REF]. The image at each gate is consisting of a heart-like object with grey-values in [0, 1], which is digitized using 120 × 120 pixels For comparison, we also present the result for static TV-regularized reconstruction in

This manuscript is for review purposes only. The last row shows the ground truth for each gate. The regions of interest in the small boxes (column 3, rows 6-7) are enlarged in Figure 6. Furthermore, SSIM and PSNR values are reported in Table 3. As can be seen, the SSIM and PSNR values for the proposed method are better than those obtained by static TVregularized reconstruction. Furthermore, they are quite similar for different parameter choices.

As shown in Table 3, these values are a little bit decreased when the value of kernel parameter σ is changed from 1.0 to 0.5 with fixed µ 1 and µ 2 , as compared the values between row 1 and row 3, also row 2 and row 4, for instance. Therefore, this test demonstrates that to some extent the proposed method is not sensitive to the precise selection of the regularization parameters under the visual perception and the quantitative comparison (SSIM and PSNR).

However, those values are selected too big or too small, which would causes over-or underregularized results.

As indicated in subsection 3.3, the optimal velocity field of the proposed model is nonvanishing at both the initial and the end time points. To verify this standpoint in numerical way, we display the computed optimal velocity field at times t = 0 and t = 1 in Figure 7 for the example with parameter pair (0.005, 10 -7 , 0.5) on Row 4 of Figure 5. problem. Based on the equivalency, we made a mathematical comparison against the joint TV-TV optical flow based model [START_REF] Burger | A Variational Model for Joint Motion Estimation and Image Reconstruction[END_REF], which showed that our method can guarantee elastically diffeomorphic deformations, and is of benefit to the practical computation additionally.

This manuscript is for review purposes only. Figure 7: The computed optimal velocity fields at times t = 0 (left) and t = 1 (right) for the example with parameter pair (0.005, 10 -7 , 0.5) on row 4 in Figure 5.

Furthermore, the theoretical comparison was also performed between the proposed model and other diffeomorphic motion models, which demonstrated that the optimal velocity field of our model is distributed w.r.t. time t averagely, and non-vanishing at both the initial and the end time points. We also presented an efficient computational method for the time-discretized formulation and showed that its optimal solution is consistent with that of the time-continuous one. This is not the case for the diffeomorphic motion model in [START_REF] Hinkle | 4D CT image reconstruction with diffeomorphic motion model[END_REF].

This manuscript is for review purposes only.

3. 3 . 1 .

 31 Optical flow based models. Recently, [13] proposes an optical flow based variational model (joint TV-TV optical flow) for joint motion estimation and image reconstruction in spatiotemporal imaging. The model is formulated as a PDE-constrained optimal control

Algorithm 4 . 3

 43 Alternately minimizing model(4.8)

12 : 13 :

 1213 Let ν = ν k . Perform lines 10-13 in Algorithm 4.1. Output I k+1 . Let I = I k+1 . Perform lines 11-15 in Algorithm 4.2. Output ν k+1 . 14:

  kernel represented by symmetric and positive definite Gaussian function K : Ω × Ω → M y ∈ R 2 and σ > 0 (kernel width).

  Figure 1.

Figure 1 :

 1 Figure 1: Data at the first view for Gate 1. The left, middle, and right figures show data at the first view for different noise levels 4.71dB, 7.7dB, and 14.67dB, respectively. The blue smooth curve is noise-free data, and the red jagged curve is noisy data.

4 .Figure 2 :

 42 Figure 2: Multi-object phantom. Columns represent the gates and the first three rows are reconstructed spatiotemporal images for the data with noise levels 4.71dB, 7.7dB, and 14.67dB, respectively. The last row shows the ground truth for each gate. The regions of interest (the small boxes in images at the 3rd column in rows 3 and 4) is enlarged in Figure 4.

Figure 3 :Figure 3

 33 Figure 3: Static TV-regularized reconstructions for the measured data with different noise levels 4.71dB (left), 7.7dB (middle), and 14.67dB (right), respectively. The region of interest in the small box (right) is enlarged in Figure 4.

Figure 4 :

 4 Figure 4: Regions of interest for the images in column 3 in Figure 2. Left images shows the region of interest in row 4 (ground truth), middle is for row 3 (the proposed method)) and right is from Figure 3 (static TV-regularization).

  . The images of all gates are supported on a fixed rectangular domain Ω = [-4.5, 4.5] × [-4.5, 4.5]. At each gate, the noise-free data per view is measured by evaluating the 2D parallel beam scanning geometry with uniform 170 bins, which is supported on the range of [-6.4, 6.4]. Then the additive Gaussian white noise is added onto the noise-free data. The resulting SNR is about 14.9dB. For gate i (1 ≤ i ≤ N ), the scanning views are distributed on [(i-1)π/5, π+(i-1)π/5] evenly, which totally has 20 views. The factor of discretized time degree is M = 2. The gradient stepsizes are set as α = 0.01 and β = 0.05, respectively. We first employ Algorithm 4.1 to gain an initial template after 50 iterations, then use Algorithm 4.3 to solve the proposed model. With selecting different values for regularization parameters, after 200 iterations, the reconstructed results are obtained, as shown in Figure 5. The runtime for each example is about 3 minutes. The detailed selections of varying parameters can be referred to the caption of Figure 5.

Figure 5 :

 5 Figure 5: Heart phantom. The columns are the 4 gates and the first 6 rows are reconstructed spatiotemporal images with parameter pairs (µ 1,µ 2 , σ) chosen as (0.01, 10 -7 , 1.0), (0.01, 10 -6 , 1.0), (0.01, 10 -7 , 0.5), (0.005, 10 -7 , 0.5), (0.01, 10 -6 , 0.5), and (0.005, 10 -6 , 0.5). The last row shows the ground truth for each gate. The regions of interest in the small boxes (column 3, rows 6-7) are enlarged in Figure6.

Figure 6

 6 Figure6as we did in the first test. As shown in Figure5, the related reconstructed results are almost the same and close to the counterpart ground truth. However, the reconstructed result by static TV-regularization in Figure6is severely degraded. To illustrate this, we enlarge the regions of interest in the small boxes of Figure5(column 3, rows 6-7) and Figure6for comparison, which are shown in Figure6.

Figure 6 :

 6 Figure6: Leftmost image shows reconstruction using static TV-regularization from the measured data with comparable noise level with the region of interest shown separately in the image ot the right. The next image to the right is the same region taken from column 3, row 6 in Figure5(the proposed method). Finally, the rightmost image is from column 3, row 7 in Figure5(ground truth).

6 .

 6 Conclusions and the future work. A general framework of variational model has been investigated for joint image reconstruction and motion estimation in spatiotemporal imaging, which is based on the deformable templates from shape theory. Along this framework, we proposed a new variational model for solving the above joint problem using the principle of LDDMM. The proposed model is equivalent to a PDE-constrained optimal control

  , and the kernel width σ. Each table entry has two values, the top being the SSIM and the bottom being the PSNR. The image is shown in the corresponding row/gate in Figure5. The bottom row reports SSIM and PSNR values of TV-regjularized reconstructions compared to each ground truth from gates 1-4 by the measured data with comparable noise level.

  2.3. Joint image reconstruction and motion estimation. It is not difficult to observe that the inverse problem in (2.3) decomposes into two sub-problems. Given a diffeomorphism φ t , the original problem reduces to a static image reconstruction problem where one seeks to recover the template "I" from noisy measured data. The data represents indirect noisy measurements of the template that has undergone a known diffeomorphic deformation. Conversely, given the template I, the original problem reduces to an indirect registration problem at each point in time t. More precisely, the task is to recover a curve of diffeomorphisms "t → φ t " that registers the template at time t against a target that is known indirectly

through noisy measured data g(t, • ). Such a "sequential indirect image registration problem"

  Remark 3.6. Note that in(3.11), we regularize the velocity field more at the beginning, and which is relevant because the template is selected at the initial time t = 0, then the This manuscript is for review purposes only. velocity field at time t i has influence on the geodesic trajectory at t ≥ t i . Equivalently, the geodesic trajectory at time t i depends on the velocity field at t ≤ t i . The earlier on we are at the velocity field, the more influence one has on the entire geodesic trajectory. Since the optimal velocity field (3.14) is averagely distributed w.r.t. time t (subsection 3.2.1), the motion of the object is close to a uniform speed. This can also result in a velocity field that is non-vanishing at the end time points (see(3.15) and(3.16)).

	(3.25)	I∈X ν∈L 2 V (Ω)	0	1	D T 2 V dt
		s.t. φ ν 0,t solves (3.1).
					2
	V term. In contrast, (3.11) uses a non-uniformly weighting (see (3.12) for more clarity) that
	puts more weight on the previous time instance.

the PDE constraint of

(3.18)

is L ∞ in time and in BV (Ω) at any t [21,

Theorem 4]

. Hence, (3.18) has a simpler structure which is also beneficial in implementation.

3.3.2. Diffeomorphic motion models.

A diffeomorphic motion model for 4D CT image reconstruction was proposed in

[START_REF] Hinkle | 4D CT image reconstruction with diffeomorphic motion model[END_REF]

. It is based on the LDDMM growth model

[START_REF] Grenander | Pattern Theory. From Representation to Inference[END_REF]

with a time-continuous model that reads as min t (φ ν 0,t . I), g(t,

• ) + µ 2 ν(t, •)

Compared to

(3.11)

, the above has no regularization term R 1 . Another difference relates to the selection on the shape regularization R 2 . In (3.25) one uniformly weights the ν(t, •)

  norm as in (2.7) and the spatial regularization is the TV functional as in (2.10). Then

	(4.2) reads as
	min

  T * t i is the adjoint of T t i . The iterative scheme (4.11) is used for updating the template in the intertwined scheme (4.3) for solving (4.6).

	Remark 4.2. There are several optimization techniques for solving convex non-smooth
	problems, like minimising I → J (4.13) φ ν s,t = Id +

ν (I) in

(4.7)

. Proximal gradient methods

[START_REF] Sidky | Convex optimization problem prototyping for image reconstruction in computed tomography with the Chambolle-Pock algorithm[END_REF][START_REF] Chen | The Linearized Split Bregman Iterative Algorithm and Its Convergence Analysis for Robust Tomographic Image Reconstruction[END_REF][START_REF] Nien | Fast X-ray CT image reconstruction using a linearized augmented Lagrangian method with ordered subsets[END_REF][START_REF] Chen | A new linearized split Bregman iterative algorithm for image reconstruction in sparse-view X-ray computed tomography[END_REF][START_REF] Bauschke | Convex Analysis and Monotone Operator Theory in Hilbert Spaces[END_REF]

are an important class of methods that are suitable for solving non-differentiable convex optimization problems, so they can be applied for solving (4.7). However, evaluating the proximal operator of a function itself involves solving a small convex optimization problem. For this reason, these more advanced methods introduce further auxiliary variables/parameters and they come with slower convergence rates. 4.2.1. Computing diffeomorphic deformations. Updating the template requires access to diffeomorphic deformations φ ν t i ,0 and φ ν 0,t i for 1 ≤ i ≤ N . By definition, φ ν s,t solves the flow equation (4.12)    ∂ t ϕ(t, x) = ν t, ϕ(t, x) ϕ(s, x) = x for x ∈ Ω and 0 ≤ s, t ≤ 1, where s is a fixed time point. Integrating w.r.t. time t in (4.12) yields t s ν τ, φ ν s,τ dτ for 0 ≤ t ≤ 1.

Table 1 .

 1 As listed in the above table, the corresponding SSIM and PSNR values are depended on SNR of the data. The higher SNR, the larger values of SSIM and PSNR.

	Figure 2 Gate 1 Gate 2 Gate 3 Gate 4 Gate 5
	Row 1	0.4069 0.4208 0.4273 0.4305 0.4337 22.10 23.02 23.27 23.40 23.64
	Row 2	0.5934 0.6086 0.6131 0.6149 0.6156 25.36 27.22 27.37 27.66 27.86
	Row 3	0.8411 0.8523 0.8564 0.8576 0.8587 28.30 31.49 32.48 32.65 32.76

This manuscript is for review purposes only.

Table 1 :

 1 SSIM and PSNR values of reconstructed spatiotemporal images for data with varying noise levels, see Figure2for the images. Comparisons are made against ground truth. Each table entry has two values, the top being the SSIM and the bottom being the PSNR. The image is shown in the corresponding row/gate in Figure

Table 2 :

 2 SSIM and PSNR values of TV-regularized reconstructions compared to each ground truth. Data is from Gates 1-5 with varying noise levels, see Figure

	Figure 3 Gate 1 Gate 2 Gate 3 Gate 4 Gate 5
	Left	0.3012 0.3163 0.3202 0.3146 0.3030 18.57 19.94 20.42 19.98 18.80
	Middle	0.4673 0.4867 0.4910 0.4840 0.4694 20.44 22.82 23.76 22.90 20.79
	Right	0.6004 0.6239 0.6291 0.6212 0.6029 21.42 24.71 26.40 25.00 21.95

Table 3 :

 3 Gate 1 Gate 2 Gate 3 Gate 4 Rows 1-6: SSIM and PSNR values of reconstructed spatiotemporal images for varying values of the regularization parameters µ 1 , µ 2

	Row 1	0.8928 0.9382 0.9340 0.9235 24.25 28.44 27.64 26.28
	Row 2	0.8960 0.9415 0.9346 0.9234 24.30 28.47 27.67 26.37
	Row 3	0.9103 0.9497 0.9459 0.9343 25.33 29.41 28.97 27.78
	Row 4	0.8940 0.9368 0.9361 0.9291 25.13 29.06 28.60 27.65
	Row 5	0.9087 0.9472 0.9462 0.9336 25.14 29.30 28.83 27.68
	Row 6	0.8884 0.9339 0.9358 0.9295 25.23 29.06 28.65 27.74
	Static TV	0.5641 0.7310 0.7458 0.5969 14.09 19.09 18.96 14.01

This manuscript is for review purposes only.

Acknowledgments. The authors would like to thank Alain Trouvé for his helpful discussions. The authors also would like to express their gratitude to the anonymous reviewers for their comments, which helped to improve the presentation of this paper.

* Submitted to the editors 2018-12-18. The work of all authors was supported by the Swedish Foundation for Strategic Research grant No. AM13-0049. Chen was in addition supported by the Beijing Natural Science Foundation grant No. Z180002.

An alternating gradient descent algorithm was designed to solve the time-discretized proposed model, where the main calculations were only based on the easily-implemented linearized deformations. For spatiotemporal (2D space + time) parallel beam tomographic imaging, the computational cost of the algorithm is then O(n 3 N N v ) and its memory requirement scales as O(n 2 M N 2 ). With Algorithm 4.3, we have evaluated the performance of the proposed model in dealing with the 2D space + time tomography in the case of very sparse and/or highly noisy data. As shown in these visual and quantitative results, the new method can yield reconstructed spatiotemporal images of high quality for the above difficult problems.

The future work will focus on the theoretical analysis of the proposed model, such as the existence and uniqueness of the solution, the convergence analysis of the proposed algorithm, and its extensions and applications to more complicated modalities in spatiotemporal imaging.

Appendix A. Optimality conditions. To derive optimality conditions for (3.11), we begin with the following result, which is proved in [START_REF] Younes | Shapes and Diffeomorphisms[END_REF].

Lemma A.1 can be used to prove the following statement:

Lemma A.2. Let the assumptions in Lemma A.1 hold and assume

Then, under the group action in (2.4) we have

Proof. By the chain rule we get

Inserting (A.4) into (A.3), we immediately prove (A.2).

The following result is a direct consequence of the above definition and Lemma A.2.

Lemma A.3. Let the assumptions in Lemma A.2 hold and define the data discrepancy functional D gt : X → R as in (3.13). If D gt is differentiable, then

This manuscript is for review purposes only.

We are now ready to characterize optimality conditions for (3.11).

Assuming that the regularization term R 1 is differentiable, and V is a RKHS with a reproducing kernel K : Ω × Ω → M n×n + . Then, the optimality conditions for (3.11) read as

Here, the L 2 V (Ω)-gradient w.r.t. the velocity field ν of L(I, • ) :

for 0 ≤ t ≤ 1 with K(ϕ) := Ω K( • , y)ϕ(y) dy. Furthermore, the gradient w.r.t. the template

where ∇R 1 denotes the gradient of R 1 : X → R.

Proof. Applying the result in Lemma A.3, we immediately have

Changing the order of integration in the above equation gives

As V is a RKHS with a reproducing kernel represented by

Combining (A.10) with (A.11) proves (A.8). Finally, the results in (A.9) and (A.7) are rather straightforward to obtain, so we omit their proofs.

This manuscript is for review purposes only.

Appendix B. First-order variation of E I .

Theorem B.1. Let the assumptions in Lemma A.2 hold and E I : L 2 V (Ω) → R is given as in (4.8). Assume furthermore that V is a RKHS with a reproducing kernel K : Ω × Ω → M n×n + .

Then, the

for 0 ≤ t ≤ 1 and x ∈ Ω.

Proof. From Lemma A.2 it is not difficult to derive

h I,ν τ,t i ∇(I • φ ν τ,0 ), η(τ, • )