Consistent Approximation of Interpolating Splines in Image Metamorphosis

This paper investigates a variational model for splines in the image metamorphosis model for the smooth interpolation of key frames in the space of images. The Riemannian manifold of images based on the metamorphosis model defines shortest geodesic paths interpolating two images as minimizers of the path energy which measures the viscous dissipation caused by the motion field and dissipation caused by the material derivative of the image intensity along motion paths. In this paper, we aim at smooth interpolation of multiple key frame images picking up the general observation of cubic splines in Euclidean space which minimize the squared acceleration along the interpolation path. To this end, we propose the spline functional which combines quadratic functionals of the Eulerian motion acceleration and of the second material derivative of the image intensity as the proper notion of image intensity acceleration. We propose a variational time discretization of this model and study the convergence to a suitably relaxed time continuous model via Γ\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\varGamma $$\end{document}-convergence methodology. As a byproduct, this also allows to establish the existence of metamorphosis splines for given key frame images as minimizers of the time continuous spline functional. The time discretization is complemented by effective spatial discretization based on finite differences and a stable B-spline interpolation of deformed quantities. A variety of numerical examples demonstrates the robustness and versatility of the proposed method in applications. For the minimization of the fully discrete energy functional, a variant of the iPALM algorithm is used.


Introduction
At first, we briefly review image metamorphosis as a flexible model for image morphing which generalizes the flow of diffeomorphism approach (cf. the textbook by Younes [You10]).It has been extensively investigated by Trouvé, Younes and coworkers [TY05b,TY05a].Image metamorphosis is a Riemannian manifold approach to the space of images and can be used to interpolate two images by a connecting geodesic path.Thereby, geodesic paths are minimizers of the associated path energy over all regular paths connecting a pair of input images.The underlying metric associates a cost both to the transport of image intensities via a viscous flow and to image intensity variations along motion paths.The path energy is given as time integral over the metric evaluated on a pair of transport velocities and material derivatives of image intensities.The metamorphosis model has been extended to discrete measures [RY13], to reproducing kernel Hilbert spaces [RY16], to functional shapes [CCT18], to images on Hadamard manifolds [ENR20] and to deep feature spaces [EKP + 21].Based on the time discrete metamorphosis model proposed in [BER15], Effland et al. introduced the Bézier curves interpolation in the space of images [ERS + 15], image geodesics for optical coherence tomography [BBE + 17] and discrete extrapolation [ERS18].The (time discrete) metamorphosis model was further used for regularization of several inverse problems in imaging: exemplar-based face colorization [PPS17], sparse and limited angle computerized tomography and super-resolution of images [NPS19], and joint (spatio-temporal) tomographic reconstruction and registration [GC Ö20].
Geodesic paths are the generalization of straight lines in Euclidean space to Riemannian manifolds.The need for an as smooth as possible interpolation of multiple data points is an extensively studied problem in mathematics.Cubic splines are a widespread and versatile solution for this problem.Due to the famous result by de Boor [dB63], in Euclidean space this approach amounts to finding a path t → u(t) which minimizes the integral of the squared acceleration 1 0 |ü(t)| 2 dt.In a Riemannian context, Noakes et al. [NHP89] introduced Riemannian cubic splines as stationary paths of the integrated squared covariant derivative of the velocity.In this paper, we aim at generalizing this ansatz to the Riemannian space of images equipped with the metamorphosis metric.Piecewise geodesic interpolation of images lacks smoothness at the interpolation data points leading to a jerking at the corresponding times when animating the resulting family of images.To overcome this shortcoming, we consider as the spline energy the sum of two integrals over squared acceleration quantities.The first quantity is the classical transport based Eulerian acceleration and the second quantity represents the second order material derivative of the image intensity.Minimizers of this spline energy are one way to naturally generalize Euclidean splines to the metamorphosis model on the space of images.In fact, our model separates in a physically intuitive way the Eulerian flow acceleration, and the second material derivative of the image intensity.This way differs from the spline energy introduced by Noakes with the spline energy being the integral over the squared covariant derivative of the path velocity (cf.also [NHP89,TV19,Via20] or [HRW18]).In this case the squared covariant derivative of the path velocity with respect to the Riemannian metric would lead to an interwoven model of the different types of acceleration.
One of the first applications of splines in computer vision was work by Mumford [Mum94] where splines appear as a maximum likelihood reconstruction of occluded edges due to the penalization of curvature of arc length parametrized curves.Today, there is a variety of spline approaches in nonlinear spaces and with applications to shape spaces.Trouvé and Vialard [TV12] studied a secondorder shape functional in landmark space based on an optimal control approach.Singh et al. [SVN15] introduced an optimal control method involving a functional which measures the motion acceleration in a flow of diffeomorphisms ansatz for image regression.Tahraoui and Vialard [TV19] consider a second-order variational model on the group of diffeomorphisms of the interval [0, 1].They proposed a relaxed model leading to a Fisher-Rao functional, as a convex functional on the space of measures.Vialard [Via20] showed the existence of a minimizer of the Riemannian acceleration energy on the group of diffeomorphisms endowed with a right-invariant Sobolev metric of high order.Benamou et al. [BGV19] and Chen et al. [CCG18] discussed spline interpolation in the space of probability measures endowed with the Wasserstein metric.Thereby, energy splines are defined as minimizers of the action functional on Wasserstein space which involves the acceleration of measure-valued paths, sharing similarities with the spline functional in the space of images introduced in this paper.The initial computational intractability of such approaches is tackled by a relaxation based on multimarginal optimal transport and entropic regularization.The transport problem this approach aims at solving might not have a Monge solution, which was remedied by a new method introduced by Chewi et al. [CCLG + 21] to construct measure-valued splines, dubbed transport splines.This method additionally enjoys substantial computational advantages.
In this paper, we will recall the image metamorphosis model and the corresponding path energy whose minimizers are geodesic paths.On that basis we then derive the novel spline energy functional.Given a set of key frames at disjoint times a spline curve is then given as a smooth in time minimizer of this spline energy respecting the key frame images as interpolation constraints.Furthermore, we will study a suitable time discrete variational model, which generalizes the time discrete metamorphosis model proposed in [BER15, EKP + 21] .The central contribution of this paper is the convergence of this time discrete model to the time continuous metamorphosis spline model in the sense of Mosco [Mos69].As a consequence, one obtains existence of metamorphosis spline paths.Finally, we discretize the model in space and we derive a numerical scheme to solve for fully discrete metamorphosis spline paths in the space of images.
Notation.Throughout this paper, we assume that the image domain Ω ⊂ R d for d ∈ {2, 3} is bounded and strongly Lipschitz.We use standard notation for Lebesgue and Sobolev spaces from the image domain Ω to a Banach space X, i.e.L p (Ω, X) and H m (Ω, X) and omit X if the space is clear from the context.The associated X norms are denoted by • L p (Ω) and • H m (Ω) , respectively, and the seminorm in for f ∈ H m (Ω).We use the notation C k,α (Ω, X) for Hölder spaces of order k ≥ 0 with Hölder regularity α ∈ (0, 1] for the k-th derivatives and the corresponding (semi)norm is The space AC p ([0, 1], X) is the space of absolutely continuous functions with the derivative in L p ((0, 1)).
The symmetric part of a matrix A ∈ R d,d is denoted by A sym , i.e.A sym = 1 2 (A+A ) and the symmetrized Jacobian of a differentiable function φ by ε[φ] = (Dφ) sym .We denote by 1 both the identity map and the identity matrix.Finally, f t refers to the evaluation of the function at time t, while ḟ refers to the temporal derivative of a differentiable function f .

Organization.
This paper is organized as follows.In Section 2 we review the most important properties of the flow of diffeomorphism model and the metamorphosis model as its extension.In Section 3 the time continuous spline energy is derived and the proper interplay between Lagrangian and Eulerian perspective is discussed.Then, in Section 4 a variational time discretization of the continuous spline energy is introduced and the existence of discrete splines is studied.Section 5 represents the time extension of the time discrete quantities, which allows study of convergence to the time continuous model, presented in Section 6. Section 7 explains the fully discrete scheme and Section 8 shows how to set up a suitable iPALM algorithm to numerically solve for a spline interpolation given a set of key frames.Finally, Section 9 experimentally demonstrates properties of the spline approach for image metamorphosis and shows applications of the proposed method.
In comparison with the conference proceeding [JRR21] we present a significantly more detailed study of both time discrete and time continuous metamorphosis splines.We prove existence of discrete metamorphosis splines, give a suitable continuous extension of discrete metamorphosis splines and show that the discrete spline functional converges in the sense of Mosco to the continuous spline functional.This in particular enables to establish the existence of continuous metamorphosis splines as minimizers of the continuous spline energy.Furthermore, we extended the applications section.

Review of Metamorphosis model
In this section, we briefly review the classical flow of diffeomorphism model and the metamorphosis model as its generalization.

Flow of Diffeomorphism
The flow of diffeomorphism model [DGM98, BMTY05, JM00, MTY02] is based on Arnold's paradigm [Arn66] to study of flow of ideal fluids by geodesics on the group of volume-preserving diffeomorphisms.To this end, given a domain Ω ⊂ R d , one considers a family of diffeomorphisms (ψ t ) t∈[0,1] : Ω → R d determined by its time-dependent Eulerian velocity v t = ψt • ψ −1 t .The Riemannian structure on this space is given by the metric and the path energy Here, L defines a quadratic form corresponding to a higher order elliptic operator.The specific choice, used throughout this paper is In this case the metric g ψt ( ψt , ψt ) describes the viscous dissipation in a multipolar fluid model as investigated by Nečas and Šilhavý [N Š91].The first term of the integrand represents the dissipation density in a simple Newtonian fluid and the second term can be regarded as a higher order measure of the fluid friction.Using that the metric g ψt is H m (Ω)-coercive [DGM98, Theorem 3.1] shows the existence of a flow of diffeomorphisms as a minimizer of the above energy.This minimizer represents a geodesic path connecting two fixed diffeomorphisms.
In the context of image morphing the flow of diffeomorphism is transporting image intensities along particle paths describing the temporal change of c-channel image intensities (u t ) t∈[0,1] : Ω → R c .This transport is given in terms of the equation u t (•) := u 0 • ψ −1 t (•) also known as brightness constancy assumption [HS81], which is equivalent to a vanishing material derivative D ∂t u := u + v • Du.Given two image intensity functions u A , u B , an associated geodesic path is a family of images subject to the constraint where the underlying family of diffeomorphisms (ψ t ) t∈[0,1] minimizes the path energy.

Metamorphosis Model
The metamorphosis approach originally proposed by Miller, Trouvé, Younes and coworkers in [MY01,TY05b,TY05a] generalizes the flow of diffeomorphism model by allowing for image intensity variations along motion paths and penalizing the squared material derivative in the metric.Under the assumption that the image path u is sufficiently smooth, the metric and the path energy read as Hence, the flow of diffeomorphism model is the limit case of the metamorphosis model for δ → 0. Since paths in the space of images do not exhibit any time nor space smoothness properties in general, the evaluation of the material derivative is not well-defined.As a solution to these problems, Trouvé and Younes [TY05a] proposed a nonlinear geometric structure in the space of images I := L 2 (Ω, R c ).In detail, for a given image path u ∈ L 2 ([0, 1], I) and an associated velocity field v ∈ L 2 ((0, 1), V), where V := H m (Ω, R d ) ∩ H 1 0 (Ω, R d ), we define the vector valued weak material derivative ẑ ∈ L 2 ((0, 1), L 2 (Ω, R c )) by With further consideration of equivalence classes of motion fields and material derivatives inducing the same temporal change of the image intensity we can consider (v, ẑ) as a tangent vector in the tangent space of I at the image u and write (v, ẑ) ∈ T u I.This gives rise to the notion H 1 ([0, 1], I) for regular paths in the space of images.The path energy in the metamorphosis model for a regular path u ∈ H 1 ([0, 1], I) is then defined as Image morphing of two input images u A , u B ∈ I then amounts to computing a shortest geodesic path u ∈ H 1 ([0, 1], I) in the metamorphosis model, which is defined as a minimizer of the path energy in the class of regular curves such that u(0) = u A and u(1) = u B .The infimum in (3) is attained, which is shown in [TY05a, Proposition 1 & Theorem 2], and the existence of a shortest geodesic is proven in [TY05a,Theorem 6].The Lagrangian formulation of variation of the image intensity along motion trajectories is This motivates a relaxed approach proposed in [ENR20, EKP + 21], where the material derivative quantity is retrieved from a variational inequality for a.e.x ∈ Ω and all 1 ≥ t > s ≥ 0. Formally, the scalar valued z = |ẑ| replaces the actually vector-valued material derivative.In fact, given a path u in L 2 ([0, 1], I) the inequality (5) defines a set C(u) ⊂ L 2 ((0, 1), V) × L 2 ((0, 1) × Ω) of admissible pairs (v, z) ∈ C(u) of the velocity and scalar weak material derivative fulfilling this inequality.Then the geodesic energy for a path u ∈ L 2 ([0, 1], I) is defined by The equivalence of the approaches (3) and (6) follows from [ENR20, Proposition 8], in the sense that for every z satisfying (5) there exists a ẑ satisfying (4) with |ẑ| ≤ z and for every ẑ satisfying (4) we have that z = |ẑ| satisfies (5).Furthermore, for every u ∈ L 2 ([0, 1], I) with non-empty C[u] one can show that u ∈ C 0 ([0, 1], I) (cf.[EKP + 21, Remark 1]), which allows point evaluation in time.In both [ENR20, EKP + 21] the existence of continuous time geodesic paths was shown, where the relaxed approach turned out to be very natural.

Time Continuous Splines in Metamorphosis Model
In this section we study the spline interpolation for a set of given key frame images in the time continuous setting of metamorphosis model.In mathematical notation, given a set of J key frames u I j ∈ I we ask for spline interpolation (u t ) t∈[0,1] which satisfies the constraints To this end, we recall that cubic splines in Euclidean space minimize the integral over the squared motion acceleration subject to position constraints [dB63], whereas linear interpolation is associated with the minimization of the integral over the squared motion velocity.In our case, image morphing via minimization of the path energy (2) corresponds to this linear interpolation.Here, we introduce the acceleration quantities which will be penalized in the spline energy.The Eulerian flow acceleration is defined by and for image paths with enough smoothness the second order material derivative is given by As already mentioned in the introduction this splitting of acceleration term into flow acceleration and second order change of image intensity does not fully correspond to a Riemannian manifold approach since the stronger interference of the two entities is observed in the covariant derivative of the vector field (v, D ∂t u).The spline energy comprises of the integral over the two quadratic acceleration quantities: where for the simplicity we use the same elliptic operator (1).As in the case of geodesic path energy, we give rigorous formulations for general paths u ∈ L 2 ((0, 1), I).We observe the pairs (v, a) ∈ L 2 ((0, 1), V) × L 2 ((0, 1), V) which determine the system and the corresponding diffeomorphic flow ψ ∈ H 2 ((0, 1), H m (Ω, Ω)).Notice that from (8), one expects that the acceleration has one derivative less in comparison to the velocity.However, the approach we later take allows us to have the same number of derivatives.A similar gain of smoothness was noticed in [Via20].
Similar to (4) the Lagrangian formulation of the material acceleration ŵ ∈ L 2 ((0, 1), By further using (4) we have for all s, t ∈ (0, 1) and every τ , such that t + τ, s + τ ∈ [0, 1].Observe that for s = t − τ the right hand side of ( 10) is an integral version of the second order central difference.Because there is no differentiation involved in these definitions, they work for general image paths.Analogous to (5), we introduce the scalar quantity w ∈ L 2 ((0, 1) × Ω) as an relaxation of the weak second order material derivative for every s ≤ t ∈ (0, 1), τ ≥ 0, t + τ ≤ 1.This relaxed Lagrangian approach is substantially more handsome in comparison to the Eulerian approach (9) which will be exploited in the proof of consistency of continuous and time discrete approaches.We can show the equivalence of the two approaches corresponding to (10) and (11) (cf.[ENR20, Proposition 8]).
Proof.The first claim is obvious by the triangle inequality.To prove the converse, let z satisfy (5).We have from where we conclude that the function t )) which implies the a.e.differentiability and the existence of derivative z ∈ L 2 ((0, 1), for every s ≤ t ∈ [0, 1] and τ > 0 with t + τ ≤ 1, and the same holds for integration on interval [−τ, 0] for s − τ ≥ 0. Taking the limit as τ tends to zero and using Lebesgue's differentiation theorem [Fol99, Theorem 3.21] we conclude that for all s ≤ t ∈ [0, 1] and a.e.x ∈ Ω we have Thus, we can repeat the above procedure to conclude the existence of ŵ ∈ L 2 ((0, 1), from where the claim directly follows.
In [HRW18] a regularization of the spline path energy by addition of weighted geodesic path energy was necessary for existence and further analysis of the splines.We follow the analogous approach which seems to be unavoidable also in our model.
Definition 1 (Regularized spline energy).Let σ > 0, m > 1 + d 2 be an integer, and u ∈ L 2 ((0, 1), I) be an image curve.Then the regularized spline energy is defined by where Let us observe that similar to [EKP + 21, Remark 1] if C[u] is non-empty we have an additional regularity of the image curve with u ∈ C 1 ([0, 1], I).Motivated by [dB63], we now define the continuous time spline interpolation for given key frames.
If we do not impose any additional constraints, we say the continuous time spline interpolation has natural boundary constraints.Imposing periodic boundary condition is equivalent to defining the image curve u t on the sphere S 1 instead of the interval [0, 1].In the case of Hermite boundary conditions we additionally fix the values of velocity terms, for both the flow and the image intensity, at the initial and the final time point.Thus, we ask for v 0 = v A , v 1 = v B and, in light of Proposition 1 and a differentiation of the left hand side of (4), that ẑ0 = ẑA and ẑ1 = ẑB .Note that in the case of the Hermite boundary conditions (also called clamped boundary conditions), we implicitly require t 1 = 0 and t J = 1, so that u 0 and u 1 are also prescribed.

Variational Time Discretization
In this section we study the variational time discretization of the time continuous (regularized) spline energy.To this end, we pick up the approach of [BER15, EKP + 21] for the variational time discretization of geodesic energy.We consider a discrete image curve u = (u 0 , . . ., u K ) with u k ∈ I and define a set of admissible deformations for a fixed (small) > 0, which consists of Remark 1.The case = 0 is discussed in Remark 3.
Considering u ∈ I K+1 as time sampling at times k K , k = 0, . . ., K, of a smooth curve and ) and using forward finite difference approximations we obtain the discrete version of the Eulerian velocity v k := K(φ k − 1) and ẑk := K(u k • φ k − u k−1 ) for the discrete material derivative.Furthermore, by using central finite difference we define the discrete acceleration and as the discrete version of the second material derivative.Following [BER15, EKP + 21] we consider the discrete path energy where 2 is a simple elastic energy density.Then, the discrete counterpart of the spline energy is defined as with the energy density Finally, the regularized discrete spline energy is given by As in the continuous time model we have interpolation constraints.Let I K :=(i 1 , . . ., i J K ) be an index tuple with 2 ≤ J K ≤ K, i j ∈ {0, . . ., K} for j = 1, . . ., J K .We consider a J K -tuple I K f = (u I 1 , . . ., u I J K ) and define the set of admissible image vectors We are now in a position to define discrete splines.
Definition 3 (Discrete time regularized spline interpolation).Let u = (u 0 , . . ., u K ) ∈ I K adm .Then we set A discrete time regularized spline interpolation of {u I j } j=1,...,J K is a discrete (K + 1)-tuple that minimizes F σ,K over all discrete paths in I K adm .The presented discretization is valid in the case of natural boundary conditions, to which we restrict in further discussions.We remark that in the case of periodic boundary conditions we make an identification K ∧ = 0, K + 1 ∧ = 1 and the sum in (20) goes up to K.For the discrete version of Hermite boundary conditions we prescribe Next, we follow ideas from [EKP + 21] in order to prove the existence of discrete spline interpolations.The following lemma is the analogous result to [EKP + 21, Lemma 1] and we only state it for completeness.
Lemma 1.There exists a constant C which only depends on Ω, m, d, γ such that Proof.An application of the Gagliardo-Nirenberg inequality [Nir66] yields The last term in ( 23) is bounded by To estimate the lower order term on the right hand side we use Korn's and Poincare's inequality and write Thus, the lemma follows by combining ( 23), ( 24) and ( 25).
Remark 2. The analogous result holds for the boundedness of acceleration i.e.
Now, we show the well-posedness of (22).
Proposition 2. For every K ∈ N and every image vector u = (u 0 , . . ., u K ) ∈ I K adm there exists a deformation vector Φ = (φ 1 , . . ., φ K ) ∈ D K such that Proof.Let {Φ j } j∈N ∈ D K be a sequence for which it holds lim j→∞ F σ,K,D [u, By Lemma 1 we have Thus, {Φ j } j∈N is uniformly bounded in H m (Ω, Ω) K .By reflexivity of this space there exists a subsequence (with the same label) such that Φ j Φ in H m (Ω, Ω) K .By the compact Sobolev embedding, we have , which gives us that Φ ∈ D K .Analogously, we have boundedness of {a j } j∈N in H m (Ω, Ω) K−1 and thus a convergent subsequence satisfying a j a in H m (Ω, Ω) K−1 and a j → a in C 1,α (Ω, Ω) K−1 .Here, for every j ∈ N we have a j = (a j 1 , . . ., a j K−1 ) given by (18).From the strong convergence of deformations we have that for all k = 1, . . ., K − 1.Using weak lower semicontinuity of H m -seminorm and continuity of energy densities we have for all k, as j → ∞: To handle the rest of the terms we show that for all u ∈ I we have u To see this we approximate u by smooth functions {ũ i } i∈N ∈ C ∞ (Ω, Ω) with u − ũi I → 0.Then, using the transformation formula we obtain By first choosing i and then j we have that this expression converges to 0. Hence, for every k = 1, . . ., K we have . ., K − 1, which together with (26) finishes the proof.
In the next step, under suitable conditions, we prove that there exists a minimizing vector in I K adm (see ( 21)) for a fixed deformation vector Φ ∈ D K .Proposition 3. Let K ≥ 2, I K f and Φ ∈ D K be fixed.Assume that the deformations satisfy, for every x ∈ Ω, Then there exists a vector of images u ∈ I K adm such that Proof.Let {u j } j∈N ∈ I K adm be an approximation sequence such that Here, F σ,K,D := F σ,K,D [u, Φ] represents a finite upper bound for the energy with the vector of images where we used the transformation formula and (28).By further use of (28) we have from where we have, by induction, that {u j k } j∈N is uniformly bounded in I for every k = 0, . . ., K. By reflexivity, there exists a subsequence (labeled in the same way) such that u j k u k for some u ∈ I K adm .It remains to verify the weak lower semicontinuity of the matching functional, i.e. we have to show that for every k = 1, . . ., K and k = 1, . . ., K − 1, respectively.To this end, we first show which readily implies (30) and by applying the same technique in a nested fashion we get (31).Altogether, we have lim inf from where the optimality follows.
We are now in the position to show the existence of discrete time spline interpolations.
Theorem 1 (Existence of discrete time spline interpolations).Let K ≥ 2. Then for every Proof.Consider a sequence {u j } j∈N ∈ I K adm for which it holds lim j→∞ where the constant is independent of K. Furthermore, for every j let F σ,K [u j ] = F σ,K,D [u j , Φ j ], by Proposition 2. As in the proof of Proposition 2 we have weak convergence of {Φ j } j∈N in H m (Ω, Ω) K and strong in C 1,α (Ω, Ω) K to some Φ ∈ D K .Furthermore, we again have a j k a k in H m (Ω, Ω), and strongly in C 1,α (Ω, Ω), where a k = K 2 (φ k+1 • φ k − 2φ k + 1) and estimates from (26) are satisfied.By Proposition 3 we may replace u j by an energy optimal image vector.Keeping the same label and following the same arguments as in (29) we conclude that {u j } j∈N is uniformly bounded in I, weakly converging to u.We are left to verify the estimates for every k = 1, . . ., K and k = 1, . . ., K − 1, respectively.To that end, it is enough to show To see this, we first take into account the decomposition The second term is handled as in (27).It remains to consider the convergence properties of the first term.For a test function v ∈ I we obtain using the transformation rule The right hand side converges to 0 due to the convergence (det(Dφ This proves our claim and finally proves the theorem. Remark 3. The results of this section remain valid for any W D satisfying conditions (W 1) − (W 2) from [EKP + 21].Furthermore, let us observe the case when i j = K • t j and u ij = u I j for t j ∈ [0, 1] and u I j ∈ I for every j = 1, . . ., J (cf. (7) and (21)).Then for every large enough K (depending on u I j and t j ) one can show the existence of the discrete time spline interpolation even if = 0 in the definition of the admissible set (17).Namely, by (32) we have that F σ,K is a fixed finite upper bound for the discrete spline energy, independent of K.Then, using Lemma 1 and Sobolev embedding theorem we have max k=1,...,K By Lipschitz continuity of the determinant we have for large enough K, which proves min k=1,...,K det(Dφ k ) ≥ c det > 0 and we can proceed as before.

Temporal Extension Operators
In this section we define the suitable time extensions of the time discrete quantities from the previous section.This is necessary to compare discrete and continuous spline functional and to study convergence.
2 )τ , k = 0, 1, . . ., K. Furthermore, consider a vector of images u K = (u K 0 , . . ., u K K ) ∈ I K+1 and a vector of deformations Φ K = (φ K 1 , . . ., φ K K ) ∈ D K .We first define a discrete incremental transport path y K with y , 1], where This can be seen as the cubic Hermite interpolation on intervals (t , and an affine interpolation on [0, t K ] and (t with the corresponding slopes , respectively (cf. Figure 1 (left) for a sketch).We define the image extension operator and, for k = 1, . . ., K − 1 and t ∈ (t This can be seen as blending between the "half-way images" along the incremental transport path y K .This is depicted on Figure 1 right.The discrete velocity Left: Schematic drawing of the Hermite interpolation (blue) on the time interval ) t∈[0,1] from the left, plotted against time.Dots represent the values u K k , and crosses the "half-way" values 1 2 (u K k + u K k−1 ) along the discrete transport path.
field corresponding to the discrete incremental transport path y K t is given by ] and t ∈ (t K K− 1 2 , 1], respectively, and , 1], respectively, and for a small enough constant c > 0 (see Section 6, (39), (40)).We define the velocity and the acceleration along the incremental transport path by Now, the actual discrete flow given as the map (t, x) → ψ K t (x) is defined recursively by ], One can directly show that (37) is well-defined in the sense of equations ( 13) − (14), i.e.
Based on this, the first order scalar weak material derivative of u K can be defined as the absolute value of the material derivative along the paths t → ψ K t (x) with and For the second order scalar weak material derivative we have and w K t :=0 elsewhere, which is the absolute value of the second time derivative of U K [u K , Φ K ] along the path t → ψ K t (x).Indeed, one easily verifies (cf.[ENR20, Proposition 9]) that z K t and w K t are admissible in the sense of equations ( 15) and ( 16), i.e.
Remark 4. For periodic boundary conditions we make the cubic interpolation definition on t ∈ ].
Finally, we define the extension of the energy F σ,K to a functional F σ,K by if there exists such u K , Φ K and +∞ else.The existence of the infimum follows from the continuity of the constraint U K [u K , Φ K ] = u w.r.t. the weak convergence of u K and the strong convergence of Φ K (cf.[ENR20, Lemma 11]).

Convergence of Time Discrete Splines
In this section we study the convergence of extension of discrete regularized spline energy F σ,K to continuous counterpart F σ and the convergence of time discrete minimizers to a time continuous minimizer in the sense of Mosco.
In explicit (i) for every sequence {u K ∈ L 2 ((0, 1), I)} K∈N which converges weakly to u ∈ L 2 ((0, 1), (ii) for every u ∈ L 2 ((0, 1), I) there exists a sequence {u K } K∈N such that u K → u in L 2 ((0, 1), I) as Remark 5.The above result holds for any choice of W D satisfying assumption (W 1) − (W 3) from [EKP + 21], though we restrict ourselves to the special case Remark 6.In the case of periodic boundary conditions, the same applies in the topology L 2 (S 1 , I).
The proof requires minor alterations implied by already indicated changes in energy and interpolation.
Proof of the lim inf estimate.
Suppose we have a sequence {u K ∈ L 2 ((0, 1), I)} K∈N such that u K u in that space as K → ∞, and we suppose that F σ,K [u K ] < F < ∞.By definition for every K large enough an image vector u K ∈ I K+1 and a corresponding optimal vector of deformations Φ K ∈ D K exist, such that For the image vectors u K and the deformation vectors Φ K we define the discrete velocity and acceleration quantities as in Section 4 and their extensions as in Section 5. Using Lemma 1 we have max k=1,...,K which implies that y K k is invertible for any k = 0, . . ., K and any K large enough.Thus, we are able to define all the temporal extended quantities introduced in the previous section.Furthermore, for every t ∈ [0, 1] we have that y K k (t, •) is in C 1,α (Ω, Ω) and it is bounded uniformly in t, and by [BHS05, Theorem 2.1] the same holds for x K k (t, •) with an estimate In what follows, the quantities ẑK k and ŵK k are time discrete quantities from Section 4, while z K t and w K t are their time extensions from Section 5.The estimate (41) together with locally Lipschitz property of the determinant function gives Using the same ideas, together with )) on every subinterval we have lim This implies the uniform boundedness of {z K } K∈N and {w K } K∈N in L 2 ((0, 1) × Ω) and by reflexivity, the existence of weakly convergent subsequences (with the same labeling) to z and w, respectively.
Then by weak lower semicontinuity we have Using Korn's and Poincare's inequality we get The analogous estimates hold for v K and D m v K , where we additionally use 35)).Hence, we have that {v K } K∈N and {a K } K∈N are uniformly bounded in L 2 ((0, 1), V) and they have corresponding weak limits v and a in that space.
Furthermore, we compute the Taylor expansion of W A ((t For the remainder term we have r K a,k = O(K −6 |Da K k | 3 ) and by using Lemma 1 and (40) Then we use weak lower semicontinuity of the energy to write lim inf Analogous Taylor expansion arguments give lim inf Altogether, (42), ( 45) and (46) give lim inf We are left to show that the limit objects v, a, z, w are indeed corresponding quantities for the image curve u.First, let us observe that by using (41) (cf.also [Fio17, Proposition 1.2.4,1.2.7]) we have where ṽK t and ãK t are introduced in (36).This implies that we have the uniform boundedness of {ṽ K } K∈N in L 2 ((0, 1), V) and in L 2 ((0, 1), C 1,α (Ω, R d )) by ( 47).Then by [ENR20, Theorem 6] and Hölder's inequality we can infer that {ψ K } K∈N is uniformly bounded in C 0, 1 2 ([0, 1], C 1,α (Ω, Ω)).Furthermore, by compact embedding of Hölder spaces, we have for some To show that ψ is indeed the solution corresponding to v, we consider ψ v K , the solution corresponding to v K .By weak continuity of the solution operator mapping velocities to the flows [ENR20, Theorem 6] we have . Furthermore, by [ENR20, Remark 7] we have the Lipschitz continuity of the solution operator and using the spatial Lipschitz property of v K k together with (39) we have To show that the equation ψt = a t • ψ t is satisfied first observe that the equation ψK t = ãK t • ψ K t ensures the uniform boundedness of { ψK } K∈N in the space L 2 ((0, 1), C 1,α (Ω, Ω)).To this end we used the uniform boundedness of {ã K } K∈N in the space L 2 ((0, 1), C 1,α (Ω, Ω))(cf.(48)) and {ψ K } K∈N in C 0 ([0, 1], C 1,α (Ω, Ω)).This together with the previous paragraph implies that {ψ K } K∈N is uniformly bounded in H 2 ((0, 1), C 1,α (Ω, Ω)) and converges weakly to ψ in H 2 ((0, 1), C 1,β (Ω, Ω)) for 0 < β < 1 2 min{ 1 2 , α}.Altogether, we have that ψ ∈ H 2 ((0, 1), C 1,β (Ω, Ω)) and for all θ ∈ C ∞ c ((0, 1)×Ω).Hence, it sufficies to verify ãK •ψ K a•ψ in L 2 ((0, 1)×Ω).Since we already have a K • ψ K a • ψ in L 2 ((0, 1) × Ω), we conclude the proof by checking that ãK where we used the Lipschitz property of a K t , the transformation formula and finally (39).In order to show that z and w are indeed scalar weak material derivatives of u, first observe that from the weak convergence u K u and the strong convergence ψ K → ψ by similar approximation arguments as for (33) we obtain u where we used Hölder's inequality, the transformation formula and uniform boundedness of (ψ K ) −1 .Thus, for some finite L > 0.Then, by the weak closedness of A 1 2 ,L we obtain u • ψ ∈ A 1 2 ,L .Then for every Ω ⊂ Ω the functional u → Ω |u t (x) − u s (x)| dx is continuous because the point evaluation in time is continuous and since it is convex on A 1 2 ,L this implies weak lower semicontinuity and we obtain Since this holds for every Ω ⊂ Ω one obtains that z is the first order scalar material derivative for u.
We prove that w is the second scalar weak material derivative for u in an analogous way.This finally finishes the proof of the lim inf-inequality.
As a preparatory step for the proof of the lim sup-estimate we state a corollary of the preceding proof.
Proof.The functional F σ is coercive by Korn's inequality and Gagliardo-Nirenberg interpolation estimate and it is clearly weak lower semicontinuous.Since C[u] is a subset of a reflexive Banach space, then we just have to show weak closedness of the set.This is verified as above.
Proof of the lim sup estimate and the construction of a recovery sequence.
Consider an image curve u ∈ L 2 ([0, 1], I) with finite regularized spline energy.Then, the previous proposition guarantees the existence of an associated optimal velocity field, an acceleration field and the first and second order weak material derivatives, denoted by (v, a, z, w) ∈ C[u], respectively, i.e.
We define where ψ t is the flow associated with velocity v and ψ 0 = 1.We have max k=1,...,K by Lemma 1 and Cauchy's inequality.For the second inequality we used [BV17, Lemma 3.5] which states that Thereby, for K large enough we have Φ ∈ D K and we are in a position to define where the point time evaluation is possible since In what follows we present more detailed arguments for the second order terms, while the arguments for the first order terms follow as in [ENR20,Theorem 14].
First, we are able to relate the discrete and the continuous second order material derivative by Here, we first used (49) and the transformation formula for the first and the second equality, respectively, then the definition of the second order material derivative (11) for the first inequality, and the Cauchy-Schwarz inequality and (50) in the last two estimates.Summing the above expressions over k = 1, . . ., K − 1, we obtain Next, we express a K k in terms of a: where in the second equality we used (49), and in the last equality (14).Then, using the Cauchy-Schwarz inequality and (51) we obtain the following estimate The same Taylor expansion argument as in ( 43) and (44) now implies, together with (53) Applying Jensen's inequality twice on L, and taking into account the above relation between a K k and a gives We now estimate the summands of L individually.For the first term we use that | tr(AB)| ≤ | tr(A)|+ | trA(B − 1)|, (50) and transformation rule to get For the second term we use (51) and the fact for any 0 and iterating this argument and using the transformation formula we have Altogether with ( 53) and (54) we have which together with (52) gives This readily implies the lim sup-inequality for the pure spline part of the functional F σ,K .Following analogous steps as above, we can fully repeat the procedure from [ENR20, Theorem 14] to obtain the estimate for the lower order path energy which finally proves the lim sup-inequality.We are left to show that u K → u in L 2 ((0, 1), I) as K → ∞.To this end we introduce the piecewise constant interpolation Then, for t ∈ (t where we used the transformation formula, (50), ( 52) and an analogous estimate for the first order material derivative.For every t ∈ (0, 1) we can find a sequence {k(K)} K∈N such that for any ].We uniformly approximate the sequence {u K k(K) } K∈N by smooth functions (cf.( 27)) and use ( 50) and (53) to prove Plugging this back above we have that as we wanted to show.This concludes the proof of lim sup estimate and the construction of recovery and thus concludes the proof of Mosco convergence.
As a corollary of the previous theorem we are able to show the existence of the continuous time spline defined in Section 3. To this end, let J ≥ 2 and (t 1 , . . ., t J ) ⊂ [0, 1] ∩ Q be a sequence of fixed times.Then for infinitely many K ∈ N one can choose i K j := K • t j ∈ N for all j = 1, . . ., J. Let (u I j ) j=1,...,J ⊆ I be the set of constraint images at the corresponding constraint times.Theorem 3 (Convergence of discrete spline interpolations to continuous ones).For every K that satisfies the above condition, let u K ∈ L 2 ([0, 1], I) be a minimizer of F σ,K among the image curves satisfying u j for all j ∈ {1, . . ., J}.Then a subsequence of {u K } K∈N converges weakly in L 2 ([0, 1], I) to a minimizer of the continuous spline energy F σ as K → ∞.This minimizer satisfies u tj = u I j for all j ∈ {1, . . ., J} and the associated sequence of discrete energies converges to the minimal continuous spline energy.
Proof.For j = 1, . . ., J let η j : [0, 1] → R be smooth functions with η j (t i ) = δ ij .We define a smooth interpolating curve of the fixed images ũ(t) := J j=0 η j (t)u I j .Let us define the vector ũK := (ũ(t K 1 ), . . ., ũ(t K K )) and its time extension ũK := U K [ũ K , 1 K ].This image curve gives an admissible candidate for a minimizer of the functional F σ,K .Namely, where the upper bound is independent of K.As defined above, let {u K := U K [u K , Φ K ]} K∈N , where u K , Φ K are optimal pairs for the discrete spline (see Theorem 1).In particular, we have As in (29) we show uniform boundedness of u K k in L 2 (Ω).This further implies, together with boundedness of the deformations and the convergence of the discrete transport paths to the identity that u K t is uniformly bounded in L 2 (Ω).Therefore, {u K } K∈N is uniformly bounded in L ∞ ([0, 1], I) and a subsequence converges weakly to some u ∈ L 2 ([0, 1], I).Now, we follow the usual argument and assume that there exists an image path û ∈ L 2 ([0, 1], I) with finite energy such that F σ [û] < F σ [u].By the lim sup-part of Theorem 2 there exists a sequence {û K } K∈N ⊆ L 2 ((0, 1), I) such that lim sup K→∞ F σ,K [û K ] ≤ F σ [û].Now, we apply the lim inf-part of Theorem 2, thus obtaining which is a contradiction to the above assumption.Hence, u minimizes the continuous spline energy over all admissible image curves and discrete spline energies converge to the limiting spline energy along a subsequence, i.e. lim K→∞ F σ,K [u K ] = F σ [u], which follows from (56) by using û = u.Finally, we show that u tj = u I j .To this end recall that for the sequence of piecewise constant approximations {ū K } K∈N given by (55) we showed that u K t − ūK t → 0 in I, uniformly in t.Together with u K t u t in I for every t ∈ [0, 1] which is showed by a trace theorem type argument (see [BER15,Theorem 4.1, (iv)]) we have the needed result since ūK tj = u I j .
The analogous result for arbitrary (t 1 , . . ., t J ) ⊂ [0, 1] follows from density of Q in [0, 1].Let us remark that in light of Proposition 1 we have that for optimal scalar quantities z, w it holds z = |ẑ| and w = | ŵ|.

Fully Discrete Metamorphosis Splines
To numerically implement splines for image metamorphosis we have to further discretize the space.Here, we present a model for c image channels and a two dimensional image domain Ω := [0, 1] 2 and follow the fully discrete version of image metamorphosis introduced in [EKP + 21].Before presenting the details of the space discretization, to avoid double warping in the second material derivative term (19) and to further increase the robustness of the model we explicitly introduce a vector valued material derivative z ∈ L 2 ((0, 1), L 2 (Ω, R c )).This leads to a relaxation of (12): with a penalty on the misfit of the new variable z and the actual material derivative ẑ given by (4), while w is the material derivative of z.The time discrete counterpart F σ,K [u] is defined by where, for z = (z 1 , . . ., zK ) we write Here, wk = K(z k+1 • φ k − zk ) is the discrete material derivative of zk , while ẑk is the actual material derivative of u k , with the corresponding second order derivative ( ŵk ) K−1 k=1 given by (19).
For M, N ≥ 3 we define the computational domain (x,y)∈Ω M N u(x, y) p p .The discrete image space is I M N := {u : Ω M N → R c } and the set of admissible deformations is where the discrete Jacobian operator of φ at (x, y) ∈ Ω M N is defined as the forward finite difference operator with Neumann boundary conditions.Here and in the rest of the paper we used bold faced letters for fully discrete quantities.A spatial warping operator T that approximates the pullback of an image channel u j • φ at a point (x, y) ∈ Ω M N is defined by where s is the third order B-spline interpolation kernel.This form of warping is also used for composition of deformations, i.e we define the fully discrete acceleration as an approximation of (18) by ), j = 1, 2. In summary, the fully discrete spline energy in the metamorphosis model for a (K + 1)-tuple (u k ) K k=0 of discrete images, a K-tuple (z k ) K k=1 of discrete derivatives and a K-tuple (φ k ) K k=1 of discrete deformations reads as where .
While in the spatially continuous context the compactness induced by the H m -seminorm is indispensable, in this fully discrete model grid dependent regularity is ensured by the use of cubic B-splines.Thus, we dropped the higher order Sobolev norm terms in this fully discrete model.
To improve the robustness of the overall optimization, we take into account a multiresolution strategy.In detail, on the coarse computational domain of size M L × N L with M L = 2 −(L−1) M and N L = 2 −(L−1) N for a given L ≥ 1, a time discrete spline sequence (u k ) K k=0 is computed as minimizer of F σ,K M L N L subject to given fixed images u ij = u I j , j = 1, . . ., J.In subsequent prolongation steps, the width and the height of the computational domain are successively doubled and the initial deformations, images and derivatives are obtained via a bilinear interpolation of the preceding coarse scale solutions.

Numerical Optimization Using the iPALM Algorithm
In this section, we discuss the numerical solution of the above fully discrete variational problem based on the application of a variant of the inertial proximal alternating linearized minimization algorithm (iPALM, [PS16]).Using this algorithm effective optimization results were achieved for a wide range of non convex and non smooth problems.In particular, it was already used for numerical optimization in the context of the deep feature metamorphosis model [EKP + 21].Following [EKP + 21], to enhance the stability the warping operation is linearized with respect to the deformation at φ [β] ∈ D M N coming from the previous iteration which leads to the modified energies ] + ∇u j ).To further stabilize the computation, the Jacobian operator applied to the images is approximated using a Sobel filter.Here, •, • represent the pointwise product of the involved matrices.We use the proximal mapping of a functional f : D M N → (−∞, ∞] for τ > 0 given as Then, with the function values on ∂Ω M N remaining unchanged, the proximal operator we are interested in is given by prox k ).The first terms in both brackets are activated only for k < K.
The actual minimization F σ,K M N is performed by Algorithm 1.We used the following notation for the extrapolation with β > 0 of the k th path element in the i th iteration step while the acceleration a [β,i] is computed with correspondingly updated φ [β,i] values.Furthermore, I K is the set of fixed indices (cf.( 21)) and we denote by L[h] the Lipschitz constant of the gradient of the function h, which is determined by backtracking [BT09].The discrete deformations are initialized by the identity deformation, the discrete images by piecewise linear interpolation of the key frame images, while the discrete derivatives are initialized as differences of two consecutive images in the sequence.

Applications
In what follows, we investigate and discuss qualitative properties of the spline interpolation in the space of images, being aware that the superior temporal smoothness of this interpolation is difficult to show with series of still images.For all the examples we use L =5 levels in the multi-level approach and I = 250 iterations of iPALM algorithm on each level with the extrapolation parameter β = 1 √ 2 .For the first two examples M =N = 64, while for the others M =N = 128.Also, for the first two and the final example K = 8, while for the others K = 16.For plotting of images we crop the values to [0, 1], for the material derivative the values are scaled to the interval [0, 1] for plotting, while for the displacement and acceleration plots hue refers to the direction and the intensity is proportional to its norm as indicated by the color wheel.
Figure 2 shows a first test case.As key frames we consider three images showing two dimensional Gaussian distribution with small variance at different positions and of different mass.Spline interpolation is compared with piecewise geodesic interpolation.Furthermore, it is depicted that for the metamorphosis spline, the curve in (x, y, m)-space (position, mass) corresponds almost perfectly to the cubic spline interpolation of the parameters of the Gaussian distribution on the key frames.As a next step, we conceptually compare splines for metamorphosis and piecewise geodesic paths in Figure 3.For this specific example, we considered the image of a circle and two identical squares as key frames.The influence of the circle's curvature on the spline segment between the two identical squares is visible via concave 'edges', while for the piecewise geodesic interpolation any memory of the circle is naturally lost in the constant interpolation between the two squares.
Next, we consider spline interpolation between three human portraits on Figure 4.The plots of second material derivative and acceleration show a strong concentration around the key frames, where the spline is expected to be smooth and the piecewise geodesic path just Lipschitz.The analogous observations hold for Figure 5.
The impact of the coloring of the key frame images on the geometry of the spline interpolation and the interplay between the Eulerian flow acceleration and the second order material derivative along motion paths is depicted in Figure 6.Therein, we consider block-colored letters as key-frame images.In our first example, the coloring is consistent with the interpolating flow shown in the black and white spline interpolation in Figure 5. Hence, in the corresponding spline interpolation color is mainly passively transported along this flow rather than blending it (top row).In the second example, the red patch in the middle key-frame image is chosen to be located top right instead, further away from the red patches in the extremal key-frame images (bottom row).Hence, the coloring is no longer consistent to the flow in the black and white example.In fact, the transport in particular of the blue color is now strongly reconfigured as seen in the first interval between the 'P' and the 'A'.Furthermore, in the second interval between the 'A' and the 'Q' blending from blue to red and from red to blue occurs.
A comparison of splines and piecewise geodesic paths in shown in Figure 7.One particularly observes that for the faces the shown spline image is thicker and for the letters the spline image shows more round contours than the for the piecewise geodesic counterpart.This is again the non-local Next, we ask for a reconstruction of frames given certain frames at selected time stamps extracted from a video.Here, we compare the resulting spline interpolation and the piecewise geodesic interpolation directly with corresponding frames of the original video as a benchmark for both approaches.Indeed, Figure 8 shows this comparison of the original frames, spline interpolation and piecewise geodesic interpolation for the images extracted from a video made by David Rogers from Vanderbilt University in the 1950s2 , which shows the interaction between white blood cells and bacteria.The spline interpolation clearly shows less blending artifacts in comparison to the piecewise geodesic interpolation.

Conclusion
In this paper we introduced a novel model for the smooth interpolation of a set of given key frame images based on the image metamorphosis model.This is achieved by means of a minimization of the spline energy functional.Thereby, one penalizes acceleration reflected in the underlying transport and Eulerian acceleration defined in terms of the second order change of image intensity along the flow paths.Based on [BER15, EKP + 21] we proposed a variational time discretization by approximation of the aforementioned quantities by finite differences.We prove consistency of this time discrete model to the time continuous counterpart by the means of Mosco convergence.Space discretization is based on finite differences and approximation of warping operation by cubic B-splines, while for numerical optimization we use the iPALM algorithm.We present a range of numerical experiments where we compare the differences between the piecewise geodesic interpolation and our novel spline interpolation in the space image.The quantitative results show that the spline interpolation achieves significant improvement with respect to the smoothness in time, which is indicated in better control and regularity of acceleration quantities in the vicinity of the key frames.Furthermore, the influence of the preceding and/or succeeding key frame(s) is visible for the spline interpolation, which naturally relates to the larger support of the spline basis functions for Euclidean cubic splines.Still open is an investigation of the differences between our spline model based on the splitting of acceleration quantities and the model based on the Riemannian second order covariant derivative as relevant acceleration quantity.Furthermore, a challenge is the adaptation of our approach to textures and feature based image representations.M N , respectively.This is compared to the corresponding piecewise geodesic interpolation (green) (not visualized here, cf. Figure 7).

Figure 2 :
Figure 2: Left: Time discrete spline with framed key frame images (first row), color-coded displacement field (second row), discrete second order material derivative (third row) and color-coded discrete acceleration field (fourth row) for the Gaussians example and values of the parameters δ = 5 • 10 −3 , σ = 1, θ = 5 • 10 −5 .The colors and their intensities indicate the direction and the intensity of the field, as indicated by the color wheel on the left.Right: Euclidean splines in (x, y, m) coordinates for the input parameters versus splines for metamorphosis in (x, y, m) extracted from the numerical results in post-processing, with (x, y) denoting the center of mass and m the mass of the distribution.

Figure 3 :
Figure 3: Left: Time discrete spline (top row) and piecewise geodesic (middle row) interpolation with framed key frames.The bottom row shows the difference in intensity between the different interpolations, using the color map −0.35 0.35.Right: Width of the interpolated shape measured at the horizontal axis of symmetry (in number of pixels) for a spline interpolation (orange) and piecewise geodesic interpolation (green) showing the concavities (δ = 5 • 10 −3 , σ = 1, θ = 5 • 10 −4 ).

Figure 4 :
Figure 4: Time discrete spline with framed fixed images (first and second row), first order material derivative slack variable z (third and fourth row), second order material derivative with energies comparison (fifth and sixth row) and color-coded acceleration field with energies comparison (seventh and eight row) for values of the parameters δ = 2 • 10 −2 , σ = 2, θ = 8 • 10 −4 .The graphics on the right in row four and six show for the spline (orange) time plots of the L 2 -norm of the actual second order material derivative ŵ and the dissipation energy density reflecting motion acceleration W A (∇ M N a k ) L 1M N , respectively.This is compared to the corresponding piecewise geodesic interpolation (green) (not visualized here, cf.Figure7).

Figure 7 :
Figure 7: Top row: Image u 11 for the human face example, second order material derivative ŵk and acceleration field a k for k = 8, for the time discrete piecewise geodesic (left image of each panel pair) and spline (right image of each panel).The pairs of material derivatives and the acceleration fields are jointly scaled to reflect the differences in intensities.Bottom row: Same visualization of image u 4 from the letter example.