Image Extrapolation for the Time Discrete Metamorphosis Model - Existence and Applications

TL;DR

This paper develops a method for image extrapolation in the metamorphosis model using a discrete exponential map, enabling geodesic path computation and demonstrating its effectiveness through applications.

Contribution

It introduces a novel discrete exponential map for image extrapolation within the variational time discretization of the metamorphosis model, including existence and stability analysis.

Findings

Efficient and stable image extrapolation algorithm demonstrated.

Local existence and uniqueness proved for small time steps.

Practical algorithm implemented with Galerkin discretization.

Abstract

The space of images can be equipped with a Riemannian metric measuring both the cost of transport of image intensities and the variation of image intensities along motion lines. The resulting metamorphosis model was introduced and analyzed by Trouv\'e and Younes, and a variational time discretization for the geodesic interpolation was proposed by Berkels et al. In this paper, this time discrete model is expanded and an image extrapolation via a discrete exponential map is consistently derived for the variational time discretization. For a given weakly differentiable initial image and an initial image variation, the exponential map allows to compute a discrete geodesic extrapolation path in the space of images. It is shown that a time step of this shooting method can be formulated in the associated deformations only. For sufficiently small time steps local existence and uniqueness are proved using a suitable fixed point formulation and the implicit function theorem. A spatial Galerkin discretization with cubic splines on coarse meshes for the deformation and piecewise bilinear finite elements on fine meshes for the image intensities are used to derive a fully practical algorithm. Different applications underline the efficiency and stability of the proposed approach.

Figures

Equations

$$D^{m}f\cdot D^{m}g=\sum _{i_1,\dots ,i_m=1}^n\frac{{\partial }^{m}f}{{\partial }_{x_{i_1}}\cdots {\partial }_{x_{i_m}}}\cdot \frac{{\partial }^{m}g}{{\partial }_{x_{i_1}}\cdots {\partial }_{x_{i_m}}},\left|D^{m}f\right|={\left(D^{m}f\cdot D^{m}f\right)}^{\frac{1}{2}}.$$

$$\tilde{\mathcal{E}}[(\varphi (t){)}_{t\in [0,1]}]={\int }_0^1{\int }_{\Omega }L[v,v]\mathrm{d}x\mathrm{d}t.$$

$$\mathcal{E}[u]:={\int }_0^1\underset{(v,z)}{\inf }{\int }_{\Omega }L[v,v]+\frac{1}{\delta }{z}^2\mathrm{d}x\mathrm{d}t,$$

$$L[v,v]:=Dv:Dv+\gamma {\Delta }^{m}v\cdot {\Delta }^{m}v$$

$${\int }_0^1{\int }_{\Omega }\eta z\mathrm{d}x\mathrm{d}t=-{\int }_0^1{\int }_{\Omega }({\partial }_t\eta +\mathrm{div}(v\eta ))u\mathrm{d}x\mathrm{d}t$$

$$\mathcal{W}[u,\tilde{u}]:=\underset{\varphi \in \mathcal{A}}{\min }\left\{{\mathcal{W}}^D[u,\tilde{u},\varphi ]:={\int }_{\Omega }|D\varphi -\text{mathds}1{|}^2+\gamma |{\Delta }^m\varphi {|}^2+\frac{1}{\delta }(\tilde{u}\circ \varphi -u{)}^2\mathrm{d}x\right\},$$

$$\mathcal{A}:=\left\{\varphi \in H^{2m}(\Omega ,\Omega ):\varphi -\text{mathds}1\in H_0^{2m}(\Omega ,\Omega )\right\}.$$

$${\int }_{\Omega }|{\Delta }^m\psi {|}^2\mathrm{d}x={\int }_{\Omega }|D^{2m}\psi {|}^2\mathrm{d}x\text{for all }\psi \in H_0^{2m}(\Omega ,\Omega )$$

$${\int }_{\Omega }|\Delta \psi {|}^2\mathrm{d}x={\int }_{\Omega }\sum _{i,j=1}^n{\partial }_i^2\psi \cdot {\partial }_j^2\psi \mathrm{d}x={\int }_{\Omega }\sum _{i,j=1}^n{\partial }_i{\partial }_j\psi \cdot {\partial }_i{\partial }_j\psi \mathrm{d}x={\int }_{\Omega }|D^2\psi {|}^2\mathrm{d}x.$$

$$0\le \underset{\tilde{\varphi }\in \mathcal{A}}{\inf }{\mathcal{W}}^D[u,\tilde{u},\tilde{\varphi }]=\underset{j\to \infty }{\lim }{\mathcal{W}}^D[u,\tilde{u},{\varphi }^j]\le \overline{\mathbf{W}}:={\mathcal{W}}^D[u,\tilde{u},\text{mathds}1]=\frac{1}{\delta }\Vert \tilde{u}-u{\Vert }_{L^2(\Omega )}^2.$$

$${\int }_{\Omega }\gamma |D^{2m}({\varphi }^j-\text{mathds}1){|}^2\mathrm{d}x={\int }_{\Omega }\gamma |{\Delta }^m({\varphi }^j-\text{mathds}1){|}^2\mathrm{d}x={\int }_{\Omega }\gamma |{\Delta }^m{\varphi }^j{|}^2\mathrm{d}x\le \frac{1}{\delta }\Vert \tilde{u}-u{\Vert }_{L^2(\Omega )}^2\le \frac{C_{\mathcal{W}}^2}{\delta }.$$

$$\Vert {\varphi }^j-\text{mathds}1{\Vert }_{H^{2m}(\Omega )}\le C\Vert \tilde{u}-u{\Vert }_{L^2(\Omega )}\le C{C}_{\mathcal{W}}.$$

$$\Vert \det (D{\varphi }^j)-1{\Vert }_{L^{\infty }(\Omega )}\le C_d$$

$$\left|{\int }_{\Omega }|\tilde{u}\circ {\varphi }^j-u{|}^2-|\tilde{u}\circ \varphi -u{|}^2\mathrm{d}x\right|\le {\int }_{\Omega }(|\tilde{u}\circ {\varphi }^j-u|+|\tilde{u}\circ \varphi -u|)|\tilde{u}\circ {\varphi }^j-\tilde{u}\circ \varphi |\mathrm{d}x$$

$$\le \left(\Vert \tilde{u}\circ {\varphi }^j-u{\Vert }_{L^2(\Omega )}+\Vert \tilde{u}\circ \varphi -u{\Vert }_{L^2(\Omega )}\right)\Vert \tilde{u}\circ {\varphi }^j-\tilde{u}\circ \varphi {\Vert }_{L^2(\Omega )}\le 2\sqrt{\delta \overline{\mathbf{W}}}\Vert \tilde{u}\circ {\varphi }^j-\tilde{u}\circ \varphi {\Vert }_{L^2(\Omega )}.$$

$$\Vert \tilde{u}\circ {\varphi }^j-\tilde{u}\circ \varphi {\Vert }_{L^2(\Omega )}\le \Vert \tilde{u}\circ {\varphi }^j-{\tilde{u}}_i\circ {\varphi }^j{\Vert }_{L^2(\Omega )}+\Vert {\tilde{u}}_i\circ {\varphi }^j-{\tilde{u}}_i\circ \varphi {\Vert }_{L^2(\Omega )}+\Vert {\tilde{u}}_i\circ \varphi -\tilde{u}\circ \varphi {\Vert }_{L^2(\Omega )}.$$

$$\Vert \tilde{u}\circ {\varphi }^j-{\tilde{u}}_i\circ {\varphi }^j{\Vert }_{L^2(\Omega )}\le \Vert (\det (D{\varphi }^j)\circ ({\varphi }^j{)}^{-1}{)}^{-1}{\Vert }_{L^{\infty }(\Omega )}^{\frac{1}{2}}\Vert \tilde{u}-{\tilde{u}}_i{\Vert }_{L^2(\Omega )}\le \frac{1}{(1-C_d{)}^{\frac{1}{2}}}\Vert \tilde{u}-{\tilde{u}}_i{\Vert }_{L^2(\Omega )}.$$

$${\mathcal{W}}^D[u,\tilde{u},\varphi ]\le \underset{j\to \infty }{\mathrm{liminf}}{\mathcal{W}}^D[u,\tilde{u},{\varphi }^j],$$

$${\mathbf{E}}_K[u_0,\dots ,u_K]:=K\sum _{k=1}^K\mathcal{W}[u_{k-1},u_k].$$

$${\mathrm{EXP}}_{u_0}^k({\zeta }_1):=u_k$$

$${\mathrm{EXP}}_{u_0}^k({\zeta }_1)=u_k:={\mathrm{EXP}}_{u_{k-2}}^2({\zeta }_{k-1})$$

$$u_{k-1}=\underset{u\in H^1(\Omega )}{\mathrm{argmin}}\underset{{\varphi }_{k-1},{\varphi }_k\in \mathcal{A}}{\min }{\mathcal{W}}^D[u_{k-2},u,{\varphi }_{k-1}]+{\mathcal{W}}^D[u,u_k,{\varphi }_k].$$

$$\begin{array}{rl}{\partial }_{u_1}({\mathcal{W}}^D[u_0,u_1,{\varphi }_1]+{\mathcal{W}}^D[u_1,u_2,{\varphi }_2])(v)& =0,\\ {\partial }_{{\varphi }_1}{\mathcal{W}}^D[u_0,u_1,{\varphi }_1](\psi )& =0,\\ {\partial }_{{\varphi }_2}{\mathcal{W}}^D[u_1,u_2,{\varphi }_2](\psi )& =0,\end{array}$$

$${\int }_{\Omega }(u_1\circ {\varphi }_1-u_0)v\circ {\varphi }_1-(u_2\circ {\varphi }_2-u_1)v\mathrm{d}x$$

$${\int }_{\Omega }2D{\varphi }_1:D\psi +2\gamma {\Delta }^m{\varphi }_1\cdot {\Delta }^m\psi +\frac{2}{\delta }(u_1\circ {\varphi }_1-u_0)(\nabla u_1\circ {\varphi }_1)\cdot \psi \mathrm{d}x$$

$${\int }_{\Omega }2D{\varphi }_2:D\psi +2\gamma {\Delta }^m{\varphi }_2\cdot {\Delta }^m\psi +\frac{2}{\delta }(u_2\circ {\varphi }_2-u_1)(\nabla u_2\circ {\varphi }_2)\cdot \psi \mathrm{d}x$$

$$\Vert u_1-u_0{\Vert }_{L^2(\Omega )},\Vert u_2-u_1{\Vert }_{L^2(\Omega )}\le C_{\mathcal{W}},$$

$$\begin{array}{rl}{\int }_{\Omega }2\gamma {\Delta }^m{\varphi }_2\cdot {\Delta }^m\psi +2D{\varphi }_2:D\psi +\frac{2}{\delta }(u_1\circ {\varphi }_1-u_0)(\nabla u_1\cdot (D{\varphi }_2{)}^{-1}\psi )\circ {\varphi }_1& \\ +\frac{1}{\delta }\frac{(u_1\circ {\varphi }_1-u_0{)}^2}{\det D{\varphi }_1}\left((D{\varphi }_2{)}^{-T}:(D^2{\varphi }_2(D{\varphi }_2{)}^{-1}\psi )-(D{\varphi }_2{)}^{-T}:D\psi \right)\circ {\varphi }_1\mathrm{d}x& =0\end{array}$$

$$\begin{array}{rl}{\int }_{\Omega }& 2\gamma {\Delta }^m{\varphi }_2\cdot {\Delta }^m\psi +2D{\varphi }_2:D\psi \mathrm{d}x\\ ={\int }_{\Omega }& 2\gamma {\Delta }^m{\varphi }_1\cdot {\Delta }^m(((D{\varphi }_2{)}^{-1}\psi )\circ {\varphi }_1)+2D{\varphi }_1:D(((D{\varphi }_2{)}^{-1}\psi )\circ {\varphi }_1)\\ & -\frac{1}{\delta }\frac{(u_1\circ {\varphi }_1-u_0{)}^2}{\det D{\varphi }_1}\left((D{\varphi }_2{)}^{-T}:(D^2{\varphi }_2(D{\varphi }_2{)}^{-1}\psi )-(D{\varphi }_2{)}^{-T}:D\psi \right)\circ {\varphi }_1\mathrm{d}x.\end{array}$$

$${\mathcal{W}}^D[u_1,u_2,{\varphi }_2]={\int }_{\Omega }|D{\varphi }_2-\text{mathds}1{|}^2+\gamma |{\Delta }^m{\varphi }_2{|}^2+\frac{1}{\delta }\frac{(u_2-u_1\circ {\varphi }_2^{-1}{)}^2}{\det (D{\varphi }_2)\circ {\varphi }_2^{-1}}\mathrm{d}x,$$

$${\int }_{\Omega }2D{\varphi }_2:D\psi +2\gamma {\Delta }^m{\varphi }_2\cdot {\Delta }^m\psi +\frac{2}{\delta }(u_2-u_1\circ {\varphi }_2^{-1})\frac{(\nabla u_1\cdot (D{\varphi }_2{)}^{-1}\psi )\circ {\varphi }_2^{-1}}{\det (D{\varphi }_2)\circ {\varphi }_2^{-1}}$$

$$+\frac{1}{\delta }\frac{(u_2-u_1\circ {\varphi }_2^{-1}{)}^2}{(\det D{\varphi }_2{)}^2\circ {\varphi }_2^{-1}}\left(\mathrm{cof}D{\varphi }_2:(D^2{\varphi }_2(D{\varphi }_2{)}^{-1}\psi )-\mathrm{cof}D{\varphi }_2:D\psi \right)\circ {\varphi }_2^{-1}\mathrm{d}x$$

Full text

1 Introduction

Riemannian geometry has influenced imaging and computer vision tremendously in the past decades. In particular, many methods in image processing have benefited from concepts emerging from Riemannian geometry like geodesic curves, the logarithm, the exponential map, and parallel transport. For example, when considering the space of images as an infinite-dimensional Riemannian manifold, the exponential map of an input image w.r.t. an initial variation corresponds to an image extrapolation in the direction of this infinitesimal variation. In particular, the large deformation diffeomorphic metric mapping (LDDMM) framework proved to be a powerful tool underpinned with the rigorous mathematical theory of diffeomorphic flows. In fact, Dupuis et al. [6] showed that the resulting flow is actually a flow of diffeomorphism. Trouvé [22, 23] exploited Lie group methods to construct a distance in the space of deformations. In [14], Joshi and Miller applied this framework to (inexact and exact) landmark matching. Beg et al. [2] studied Euler–Lagrange equations for minimizing vector fields in the LDDMM framework and proposed an efficient algorithm incorporating a gradient descent scheme and a semi-Lagrangian method to integrate the velocity fields. Miller et al. [18] proved the conservation of the initial momentum in Lagrangian coordinates associated with a geodesic in the LDDMM framework, which allows for the stable computation of geodesic curves. Younes [29] used Jacobi fields in the flow of diffeomorphism approach to derive gradient descent methods for the path energy. In [11], Hart et al. exploited the optimal control perspective to the LDDMM model with the motion field as the underlying control. Vialard et al. [26, 27] studied methods from optimal control theory to accurately estimate this initial momentum and to relate it to the Hamiltonian formulation of the geodesic flow. Furthermore, they used the associated Karcher mean to compute intrinsic means of medical images. Vialard and Santambrogio investigated in [28] the flow of diffeomorphism approach for images in the space of functions of bounded variation. In particular, they were able to rigorously derive an Euler–Lagrange equation for the formulation with a matching energy. Lorenzi and Pennec [16] applied the LDDMM framework to compute geodesics and parallel transport using Lie group methods.

The metamorphosis model [19, 25] generalizes the flow of diffeomorphism approach allowing for intensity variations along transport paths and associated a corresponding cost functional with these variations. In [24], Trouvé and Younes rigorously analyzed the local geometry of the resulting Riemannian manifold and proved the existence of geodesic curves for square-integrable images and the (local) existence as well as the uniqueness of solutions of the initial value problem for the geodesic equation in the case of images with square-integrable weak derivatives. Holm et al. [12] studied a Lagrangian formulation for the metamorphosis model and proved existence for both the boundary value and the initial value problem in the case of measure-valued images. Hong et al. [13] proposed a metamorphic regression model and developed a shooting method to reliably recover initial momenta.

A comprehensive overview of most of the aforementioned topics is given in the book by Younes [30], for a historic account we additionally refer to [17].

In [4], a variational time discretization of the metamorphosis model based on a sequence of simple, elastic image matching problems was introduced and $\Gamma$-convergence to the time continuous metamorphosis model was proven. Furthermore, using a finite element discretization in space a robust algorithm was derived. Exploiting de Casteljau’s algorithm, this approach could also be used to compute discrete Riemannian Bézier curves in the space of images [8]. In [3], the geodesic interpolation proposed in [4] was employed to analyze the temporal evolution of a macular degeneration for medical images acquired with an optical coherence tomography device, where an efficient GPU implementation is used to speed up the registration subproblems.

In this paper, we focus on the discrete exponential map associated with the time discrete metamorphosis model. The Euler–Lagrange equations of the discrete path energy proposed in [4] give rise to a set of equations, which characterize time steps of a discrete initial value problem for a given initial image and a given initial image variation. We study this time stepping problem both analytically and numerically. We will prove existence and uniqueness of solutions of the single time step problem, which is guaranteed to generate time discrete geodesics in the sense of the time discrete variational approach. A straightforward treatment, for instance via a Newton scheme, would lead to higher order derivatives of image functions concatenated with diffeomorphic deformations, which are both theoretically and numerically very difficult to treat. We show how to avoid these difficulties using a proper transformation of the defining Euler–Lagrange equations and reduce the number of unknowns. With respect to the existence, we apply a fixed point argument based on Banach’s fixed point theorem for images bounded in $H^{1}$ and an initial variation, which is supposed to be small in $L^{2}$. The uniqueness proof is based on an implicit function theorem argument for initial variations, which are small in $H^{1}$. Finally, the numerical algorithm picks up the fixed point approach for another variant of the Euler–Lagrange equations.

Compared to the proceedings paper [7], which introduced the discrete exponential map in the context of the time discrete metamorphosis model and the numerical optimization algorithm, we give in this paper the comprehensive derivation of the method, formulate and prove the existence of the discrete exponential map. Furthermore, two additional applications are presented.

The paper is organized as follows: In Section 2, we briefly recall the metamorphosis model in the time continuous and time discrete setting, respectively. Departing from the Euler–Lagrange equations of a time discrete geodesic, a single time step of the discrete exponential map is derived in Section 3. Then the discrete geodesic shooting relies on the iterative application of the one step extrapolation. In Section 4, local existence and local uniqueness of this discrete exponential map are proven based on a suitable combination of the implicit function theorem and Banach’s fixed point theorem. The fixed point formulation is also used in Section 5 to derive an efficient and stable algorithm. Finally, numerical results for different applications are presented in Section 6.

We use standard notation for Lebesgue and Sobolev spaces on the image domain $\Omega$, i.e. $L^{p}(\Omega)$ and $H^{m}(\Omega)=W^{m,2}(\Omega)$. The associated norms are denoted by $\|\cdot\|_{L^{p}(\Omega)}$ and $\|\cdot\|_{H^{m}(\Omega)}$, respectively, and the seminorm in $H^{m}(\Omega)$ is given by $|\cdot|_{H^{m}(\Omega)}$. Furthermore, $H^{m}_{0}(\Omega)$ is the closure of $C^{\infty}(\Omega)$ functions with compact support w.r.t. the norm $\|\cdot\|_{H^{m}(\Omega)}$ and its dual is denoted by $H^{-m}(\Omega)$. For any $f,g\in H^{m}(\Omega)$, $m\geq 1$, we set

[TABLE]

The polyharmonic operator is inductively defined by $\Delta^{m}f:=\Delta(\Delta^{m-1}f)$ for $f\in H^{2m}(\Omega)$ with $m\geq 2$. Depending on the context, ${\mathds{1}}$ denotes either the identity mapping or the identity matrix. For a matrix $A\in\mathbb{R}^{n,n}$, we refer to $A^{\mathrm{sym}}=\frac{1}{2}(A+A^{T})$ as the symmetric part of $A$. The symbol “:” indicates the sum over all pairwise products of two tensors. Finally, we denote the variational derivative of a functional $J$ at a point $A$ in a direction $B$ by \partial_{A}J[A](B)=\frac{\mathrm{d}}{\mathrm{d}\epsilon}J[A+\epsilon B]\big{|}_{\epsilon=0}.

2 Review of the metamorphosis model and its time discretization

In this section, we briefly recall in a non-rigorous fashion the Riemannian geometry of the space of images based on the flow of diffeomorphism and its extension, the metamorphosis model. For a detailed exposition of these models we refer to [6, 19, 25, 12, 24].

Throughout this paper, we suppose that the image domain $\Omega\subset\mathbb{R}^{n}$ for $n\in\{2,3\}$ has Lipschitz boundary. For a flow of diffeomorphisms $\phi(t):\overline{\Omega}\rightarrow\mathbb{R}^{n}$ for $t\in[0,1]$ driven by the Eulerian velocity $v(t)=\dot{\phi}(t)\circ\phi^{-1}(t)$ we take into account a quadratic form $L$ subjected to certain growth and consistency conditions, which can be considered as a Riemannian metric on the space of diffeomorphisms and thus on the space of diffeomorphic transformations $u(t)=u_{0}\circ\phi^{-1}(t)$ of a given reference image $u_{0}$. Based on these ingredients one can define the associated (continuous) path energy

[TABLE]

By construction, this model comes with the brightness constancy assumption in the sense that the material derivative $\frac{D}{\partial t}u=\dot{u}+v\cdot\nabla u$ vanishes along the motion paths. Contrary to this, the metamorphosis approach allows for image intensity variations along motion paths and penalizes the integral over the squared material derivative as an additional term in the metric. Hence, the path energy in the metamorphosis model for an image curve $u\in L^{2}((0,1),L^{2}(\Omega))$ and $\delta>0$ is defined as

[TABLE]

where the infimum is taken over all pairs $(v,z)$ which fulfill the transport equation $\frac{D}{\partial t}u=\dot{u}+v\cdot\nabla u=z$. Here, we consider

[TABLE]

with $\gamma>0$ and $2m>1+\frac{n}{2}$. To formulate this rigorously, one has to take into account the weak material derivative $z\in L^{2}((0,1),L^{2}(\Omega))$ defined via the equation

[TABLE]

for all $\eta\in C^{\infty}_{c}((0,1)\times\Omega)$. Geodesic curves are defined as minimizers of the path energy (2.1). Under suitable assumptions, one can prove the existence of a geodesic curve in the class of all regular curves with prescribed initial and end image. For the definition of regular curves and the existence proof we refer to [24].

In what follows, we consider the time discretization of the path energy (2.1) proposed in [4] adapted to the slightly simpler transport cost $L[\,\cdot\,,\,\cdot\,]$. To this end, we define for arbitrary images $u,\tilde{u}\in L^{2}(\Omega)$ the discrete matching energy

[TABLE]

which is composed of a rescaled thin plate regularization term (first two terms) and a quadratic $L^{2}(\Omega)$-mismatch measure (cf. [4, (6.2)]). The set of admissible deformations $\mathcal{A}$ is defined as

[TABLE]

Remark*.*

In [4], for every admissible deformation $\phi$ the weaker boundary condition $\phi={\mathds{1}}$ on $\partial\Omega$ instead of $\phi-{\mathds{1}}\in H_{0}^{2m}(\Omega,\Omega)$ was assumed. Here, this stronger condition is required for both a higher regularity result (cf. Proposition 3.2) and a higher order control of the deformations (cf. (4.1)). With these altered boundary conditions the equality

[TABLE]

holds true for all $m\geq 1$ (cf. [9, Section 2.2]). In fact, using integration by parts we exemplarily obtain for $m=1$

[TABLE]

In what follows, we need the existence of minimizers of this particular matching energy $\mathcal{W}$ for input images $u,\tilde{u}\in L^{2}(\Omega)$ with $\|\tilde{u}-u\|_{L^{2}(\Omega)}$ sufficiently small.

Proposition 2.1** (Existence of a minimizing deformation for $\mathcal{W}$).**

Let $u\in L^{2}(\Omega)$ and $2m-\frac{n}{2}>1$. Then there exists a constant $C_{\mathcal{W}}>0$ that solely depends on $\gamma$, $\delta$, $\Omega$ and $m$ such that for every $\tilde{u}\in L^{2}(\Omega)$ with $\|\tilde{u}-u\|_{L^{2}(\Omega)}\leq C_{\mathcal{W}}$ there is a minimizing deformation $\phi\in\mathcal{A}$ for $\mathcal{W}$, i.e. $\mathcal{W}[u,\tilde{u}]=\mathcal{W}^{D}[u,\tilde{u},\phi]$, and $\phi$ is a $C^{1}(\Omega)$-diffeomorphism.

Proof.

The proof is based on the direct method in the calculus of variations. Let $(\phi^{j})_{j\in\mathbb{N}}\subset\mathcal{A}$ be a minimizing sequence for $\mathcal{W}^{D}[u,\tilde{u},\phi^{j}]$ with monotonously decreasing energy and

[TABLE]

Since $\phi^{j}-{\mathds{1}}\in H^{2m}_{0}(\Omega)$, (2.3) and (2.4) imply

[TABLE]

Thus, the norm equivalence of $\|\cdot\|_{H^{2m}(\Omega)}$ and $|\cdot|_{H^{2m}(\Omega)}$ for the space $H^{2m}_{0}(\Omega)$, which follows by an iterative application of the Poincaré inequality (cf. [1, Corollary 6.31]), yields

[TABLE]

By taking into account the embedding $H^{2m}(\Omega)\hookrightarrow C^{1}(\overline{\Omega})$ and considering a smaller $C_{\mathcal{W}}$ if necessary we can assume

[TABLE]

for a constant $C_{d}\in(0,1)$, which implies that $\phi^{j}$ is $C^{1}(\Omega)$-diffeomorphism (see [5, Theorem 5.5-2]). Moreover, since $(\phi^{j})_{j\in\mathbb{N}}$ are uniformly bounded in $H^{2m}(\Omega)$ (cf. (2.5)), a subsequence (also denoted by $\phi^{j}$) converges weakly in $H^{2m}(\Omega)$ due to the reflexivity of this space and (strongly) in $C^{1,\alpha}(\overline{\Omega})$ for $\alpha\in(0,2m-1-\frac{n}{2})$ to a $C^{1}(\Omega)$-diffeomorphism $\phi\in\mathcal{A}$.

Next, we prove the convergence of the $L^{2}(\Omega)$-mismatch terms. To this end, we estimate

[TABLE]

Now, we approximate $\tilde{u}$ in $L^{2}(\Omega)$ by a sequence of smooth functions $(\tilde{u}_{i})_{i\in\mathbb{N}}$ such that $\|\tilde{u}-\tilde{u}_{i}\|_{L^{2}(\Omega)}\leq 2^{-i}$. Then,

[TABLE]

Next, applying the transformation formula yields

[TABLE]

Likewise, we can deduce $\|\tilde{u}\circ\phi-\tilde{u}_{i}\circ\phi\|_{L^{2}(\Omega)}\leq C\|\tilde{u}-\tilde{u}_{i}\|_{L^{2}(\Omega)}$. Furthermore, the middle term in (2.6) vanishes for fixed $i$ as $j\rightarrow\infty$. Finally, using the lower semicontinuity of the first two terms of the energy we get

[TABLE]

which proves this proposition. ∎

Following the general approach for the variational time discretization of geodesic calculus in [21] and the particular discretization of the metamorphosis model in [4], we define the discrete path energy $\mathbf{E}_{K}$ on a sequence of $K+1$ images $(u_{0},\ldots,u_{K})\in(L^{2}(\Omega))^{K+1}$ with $K\geq 2$ as the weighted sum of the discrete matching energy evaluated at consecutive images, i.e.

[TABLE]

A $(K+1)$-tuple $(u_{0},\ldots,u_{K})\in(L^{2}(\Omega))^{K+1}$ with given images $u_{0}$ and $u_{K}$ is defined to be a discrete geodesic curve connecting $u_{0}$ and $u_{K}$ if it minimizes $\mathbf{E}_{K}$ w.r.t. all other $(K+1)$-tuples with $u_{0}$ and $u_{K}$ fixed. For the proof of the existence of discrete geodesics we refer to [4]. It is also shown in [4] that a suitable extension of the discrete path energy $\mathbf{E}_{K}$ $\Gamma$-convergences to the continuous path energy $\mathcal{E}$. Let us finally mention that neither the matching deformation in (2.2) nor the discrete geodesic curve defined as the minimizer of (2.7) for given input images $u_{0}$ and $u_{K}$ are necessarily unique.

3 The time discrete exponential map

In this section, we define the discrete exponential map and derive optimality as well as regularity results, on which the study of existence and uniqueness in Section 4 and the algorithm introduced in Section 5 will be based.

Let us briefly recall the definition of the continuous exponential map on a Riemannian manifold. Let $y:[0,1]\rightarrow\mathcal{M}$ be the unique geodesic curve for a prescribed initial position $y(0)=y_{A}$ and an initial velocity $\dot{y}(0)=v$ on a Riemannian manifold $(\mathcal{M},g)$. The exponential map is then defined as $\exp_{y_{A}}(v)=y(1)$. Furthermore, one easily checks that $\exp_{y_{A}}(\frac{k}{K}v)=y(\frac{k}{K})$ for $0\leq k\leq K$. We refer to the textbook [15] for a detailed discussion of the (continuous) exponential map. Now, we ask for a time discrete counterpart of the exponential map in the metamorphosis model. To this end, we consider an image $u_{0}$ as the initial data and a second image $u_{1}$ such that $\zeta_{1}=u_{1}-u_{0}$ represents a small variation of the image $u_{0}$. This variation $\zeta_{1}$ is the discrete counterpart of the infinitesimal variation given by the velocity $v$ in the continuous case. For varying values of $K\geq 2$ we now ask for a discrete geodesic $(u_{0},u_{1},u_{2},\ldots,u_{K})$ described as the minimizer of the discrete path energy (2.7). Let us for the time being suppose that this geodesic is unique – a property to be verified later. Based on our above observation for the continuous exponential map we define $\mathrm{EXP}^{k}_{\ast}(\,\cdot\,)$ as the discrete counterpart of $\exp_{\ast}(\frac{k}{K}\,\cdot\,)$, i.e. we set

[TABLE]

for $k=1,\ldots,K$. The definition of the exponential map $\mathrm{EXP}^{k}_{u_{0}}(\zeta_{1})$ does not depend on the number of time steps $K$. Indeed, if $(u_{0},u_{1},u_{2},\ldots,u_{K})$ is a discrete geodesic, then $(u_{0},u_{1},u_{2},\ldots,u_{L})$ with $L\leq K$ is also a geodesic. Taking into account $k=2$ we immediately observe that the sequence of discrete exponential maps $(\mathrm{EXP}^{k}_{u_{0}}(\zeta_{1}))_{k=1,\ldots}$ can iteratively be defined as follows

[TABLE]

for $k\geq 2$, where $\zeta_{k-1}=u_{k-1}-u_{k-2}$, and for the sake of completeness we define $\mathrm{EXP}^{0}_{u_{0}}(\zeta_{1})=u_{0}$ and $\mathrm{EXP}^{1}_{u_{0}}(\zeta_{1})=u_{1}=u_{0}+\zeta_{1}$. Thus, it essentially remains to compute $\mathrm{EXP}^{2}$ for a given input image $u_{k-2}$ and an image variation $\zeta_{k-1}=u_{k-1}-u_{k-2}$ (see Figure 1). For a detailed discussion of the discrete exponential map in the simpler model of Hilbert manifolds we refer to [21]. The particular challenge here is that the matching energy $\mathcal{W}$ cannot be evaluated directly, but requires to solve the variational problem (2.2) for the matching deformation.

There are two major restrictions regarding the input images $u_{0}$ and $u_{1}$:

Firstly, the existence and uniqueness result for the discrete exponential map (cf. Section 4) will require weakly differentiable input images. These weak derivatives of images naturally arise in the Euler–Lagrange equations for $\mathrm{EXP}^{2}$ w.r.t. the deformations (see (3.5) and (3.6) below). Let us remark that the weak differentiability of the input data for the exponential map is also a crucial requirement in the initial value problem for the geodesic equation in [24]. Furthermore, the $H^{1}(\Omega)$-regularity property is inherited along discrete geodesics, i.e. $\mathrm{EXP}^{k}_{u_{0}}(u_{1}-u_{0})\in H^{1}(\Omega)$ for any $k\geq 1$ provided that $u_{0},u_{1}\in H^{1}(\Omega)$ (cf. [4, Remark 3.3 and Equation (3.2)]).

Secondly, the initial variation $\zeta_{1}=u_{1}-u_{0}$ is assumed to be sufficiently small in $L^{2}(\Omega)$ in order to ensure the existence of the initial deformation $\phi_{1}$ and guarantee the convergence of a suitable fixed point algorithm – a property which appears to be natural in light of the analogue assumption for the continuous exponential map [15]. We will also see that for fixed $K$ the variations $u_{k-1}-u_{k-2}$ for $k\leq K$ will remain small provided that $\zeta_{1}$ is small. Thus, for fixed $K$ the discrete exponential map $\mathrm{EXP}^{k}_{u_{0}}(\,\cdot\,)$ will be well-posed for a sufficiently small initial variation $\zeta_{1}$.

Hence, in what follows we consider images in $H^{1}(\Omega)$ and define $u_{k}:=\mathrm{EXP}^{2}_{u_{k-2}}(\zeta_{k-1})$ as the (unique) image in $H^{1}(\Omega)$ such that

[TABLE]

For the sake of simplicity, we restrict to the first step in the iterative computation of the discrete exponential map with $k=2$. Given $u_{0},u_{1}\in H^{1}(\Omega)$ the first order optimality conditions for (3.2) for $u_{2}\in H^{1}(\Omega)$ and $\phi_{1},\phi_{2}\in\mathcal{A}$ read as

[TABLE]

for all $v\in H^{1}(\Omega)$ and all $\psi\in H^{2m}_{0}(\Omega)$ for $2m-\frac{n}{2}>2$. The system (3.3) is equivalent to

[TABLE]

The subsequent lemma provides a reformulation of the above system of equations, in which the dependency of the unknown function $u_{2}$ in (3.6) is removed and in addition solely the function $u_{1}$ and no longer derivatives of $u_{1}$ appear.

Lemma 3.1** (Reformulation of the Euler–Lagrange equation for $\phi_{2}$).**

Let $u_{0},u_{1},u_{2}\in H^{1}(\Omega)$ such that (cf. Proposition 2.1)

[TABLE]

$2m-\frac{n}{2}>2$, and assume that (3.4) and (3.5) hold true. Let $\phi_{i}$ with $i=1,2$ be the minimizer of $\mathcal{W}^{D}[u_{i-1},u_{i},\,\cdot\,]$ on $\mathcal{A}$ according to Proposition 2.1.

1. (i)

Then (3.6) is equivalent to

[TABLE]

for all $\psi\in H^{2m}_{0}(\Omega)$. 2. (ii)

Under the additional assumptions that $\partial\Omega\in C^{4m}$ and $u_{0},u_{1},u_{2}\in L^{\infty}(\Omega)\cap H^{1}(\Omega)$ the equation (3.6) is equivalent to

[TABLE]

Here, the notation $(D^{2}\phi_{2}(D\phi_{2})^{-1}\psi)_{jk}=\sum_{i,l=1}^{n}\partial_{j}\partial_{k}\phi^{i}_{2}(D\phi_{2})^{-1}_{il}\psi_{l}$ is used.

Proof.

By using the transformation formula the energy $\mathcal{W}^{D}$ can be rewritten as follows

[TABLE]

since $\phi_{2}\in\mathcal{A}$ is a diffeomorphism (see Proposition 2.1). As a next step, we rewrite the Euler–Lagrange equation w.r.t. $\phi_{2}$ of $\mathcal{W}^{D}[u_{1},u_{2},\,\cdot\,]$. To this end, we use the identities $\partial_{\phi_{2}}\phi_{2}^{-1}(\psi)=-((D\phi_{2})^{-1}\psi)\circ\phi_{2}^{-1}$, which follows by differentiating $(\phi_{2}+\epsilon\psi)\circ(\phi_{2}+\epsilon\psi)^{-1}={\mathds{1}}$ w.r.t. $\epsilon$, and $\partial_{A}\det(A)(B)=\operatorname{cof}(A):B$ for $A\in GL(n)$ and $B\in\mathbb{R}^{n,n}$ with $\operatorname{cof}A=(\det A)A^{-T}$. Thus, we obtain

[TABLE]

A further application of the transformation formula w.r.t. $\phi_{2}$ yields

[TABLE]

To remove the dependency of the function $u_{2}$ above, we employ the pointwise condition

[TABLE]

for a.e. $x\in\Omega$, which follows directly from (3.4). Inserting this in (3.9) and using the integral transformation formula we achieve

[TABLE]

The identity $\operatorname{cof}(A)=\det(A)A^{-T}$ for $A\in GL(n)$ implies (i).

To show (ii), we take into account the test function $\zeta:=((D\phi_{2})^{-1}\psi)\circ\phi_{1}$ in (3.5). To justify this, we have to show that $\zeta\in H^{2m}_{0}(\Omega)$. To this end, we require $H^{2m+1}(\Omega)$-regularity of $\phi_{2}$, which will follow from Proposition 3.2, and classical differential calculus for Sobolev functions [1]. Inserting $\zeta$ into (3.5) we get

[TABLE]

By adding the above equation to (3.7) we have proven (ii). ∎

Proposition 3.2** (Maximal regularity of the deformations).**

Let $2m-\frac{n}{2}>2$ and $\partial\Omega\in C^{4m}$. Furthermore, let $u_{0},u_{1},u_{2}\in L^{\infty}(\Omega)\cap H^{1}(\Omega)$ and suppose that $\phi_{1},\phi_{2}\in\mathcal{A}$ are minimizers of $\mathcal{W}^{D}[u_{0},u_{1},\,\cdot\,]$ and $\mathcal{W}^{D}[u_{1},u_{2},\,\cdot\,]$, respectively. Then $\phi_{1},\phi_{2}\in\mathcal{A}\cap H^{4m}(\Omega)$.

Proof.

We only prove the result for $\phi_{2}$, for $\phi_{1}$ one proceeds analogously. Let $w$ be the displacement associated with $\phi_{2}$, i.e. $w=\phi_{2}-{\mathds{1}}\in H^{2m}_{0}(\Omega)$. Using integration by parts in (3.6) we obtain for a test function $\psi\in H^{2m}_{0}(\Omega)$

[TABLE]

with $f\in L^{2}(\Omega,\mathbb{R}^{n})$. Then, the assertion follows from the general $L^{2}$-regularity theory for polyharmonic equations as presented in [9, Section 2.5.2]. ∎

Remark*.*

Since $\phi_{2}$ is a diffeomorphism, (3.10) is equivalent to

[TABLE]

Here, the first summand reflects the intensity modulation along the geodesic, the second summand quantifies the contribution due to the transport.

We will use the first reformulation (3.7) (Lemma 3.1 (i)) of the Euler–Lagrange equation (3.6) with respect to $\phi_{2}$ to derive a fixed point iteration in the existence proof for the time discrete exponential map. The second reformulation (3.8) (Lemma 3.1 (ii)) will later be used in a modified and spatially discrete fixed point iteration in the numerical algorithm.

4 Local existence and uniqueness of the discrete exponential map

In this section, we prove local existence and local uniqueness for the discrete exponential map. At first, we make use of an argument based on Banach’s fixed point theorem applied to the reformulation of the Euler–Lagrange equation given in Lemma 3.1 (i) for a discrete geodesic $(u_{0},u_{1},u_{2})$ with deformations $\phi_{1}$ and $\phi_{2}$. For image pairs $(u_{0},u_{1})$ we establish the existence of a solution $(u_{2},\phi_{1},\phi_{2})$ to the system of equations (3.4), (3.5) and (3.6) provided that $u_{0}$ and $u_{1}$ are bounded in $H^{1}(\Omega)$ and close in $L^{2}(\Omega)$. This does not necessarily imply that for given $u_{0}$ and $u_{1}$ the resulting discrete path $(u_{0},u_{1},u_{2})$ is the unique discrete geodesic connecting $u_{0}$ and $u_{2}$. Thus, in a second step we will show that this indeed holds true if the images $u_{0}$ and $u_{1}$ are close in $H^{1}(\Omega)$. To this end, we apply an implicit function theorem argument (cf. the corresponding proof for the discrete exponential map on Hilbert manifolds given in [21]). Let us remark that this argument also allows to establish existence, but under the stronger assumption that the input images are close in $H^{1}(\Omega)$ compared to the requirement of closeness in $L^{2}(\Omega)$ and boundedness in $H^{1}(\Omega)$ for the existence proof via the fixed point theorem. Furthermore, the fixed point approach will be taken into account for the numerical approximation of the time discrete exponential map.

Theorem 4.1** (Existence of solutions of the Euler–Lagrange equations).**

Let $2m-\frac{n}{2}>2$, $u_{0}\in H^{1}(\Omega)$. Then there are constants $C_{u},c_{u}>0$ such that for every $u_{1}\in\left\{u\in H^{1}(\Omega):\ |u|_{H^{1}(\Omega)}\leq C_{u},\;\|u-u_{0}\|_{L^{2}(\Omega)}\leq c_{u}\right\}$ there exists a solution $(u_{2},\phi_{1},\phi_{2})\in H^{1}(\Omega)\times\mathcal{A}\times\mathcal{A}$ of (3.4), (3.5) and (3.6). In particular, the defining system of equations for $u_{2}=\mathrm{EXP}^{2}_{u_{0}}(u_{1}-u_{0})$ is solved.

Proof.

We begin with some preparatory considerations. Let $c_{u}\leq C_{\mathcal{W}}$ and $\phi_{1}\in\operatorname*{argmin}_{\phi\in\mathcal{A}}\mathcal{W}^{D}[u_{0},u_{1},\phi]$ be a minimizing deformation (cf. Proposition 2.1). Following the same line of arguments as for the estimate (2.5) in the proof of Proposition 2.1 we obtain

[TABLE]

Furthermore, taking into account $\mathcal{W}^{D}[u_{0},u_{1},\phi_{1}]\leq\mathcal{W}^{D}[u_{0},u_{1},{\mathds{1}}]$ we infer

[TABLE]

Now, we define the fixed point iteration and prove the contraction property in several steps:

(i) Defining the fixed point mapping $\mathcal{F}$. Using Lemma 3.1 (i) we define for a fixed deformation $\phi_{1}$ the operators $\mathcal{T},\mathcal{R}:\mathcal{A}\to H^{-2m}(\Omega)$ as

[TABLE]

for a diffeomorphism $\phi\in\mathcal{A}$ and all $\psi\in H^{2m}_{0}(\Omega)$, which allows us to reformulate the Euler–Lagrange equation w.r.t. the deformation $\phi_{2}$ in (3.7) as

[TABLE]

Next, we will study the invertibility of the linear operator $\mathcal{R}$ and the Lipschitz continuity of $\mathcal{T}$ with a Lipschitz constant which depends monotonically on $c_{u}$ and vanishes for $c_{u}\searrow 0$. This will imply that $\mathcal{F}:=\mathcal{R}^{-1}\circ\mathcal{T}$ is a contraction for sufficiently small $c_{u}$. The fixed point iteration to compute the unknown deformation $\phi_{2}$ reads as $\phi^{j+1}=\mathcal{F}[\phi^{j}]$ for $j\in\mathbb{N}$ and $\phi^{0}={\mathds{1}}$.

(ii) Lipschitz continuity of $\mathcal{T}$. In what follows, we assume that

[TABLE]

for a sufficiently small $\epsilon>0$, the dependency of $\epsilon$ on $C_{u}$ and $c_{u}$ is discussed below. By the embedding $H^{2m}(\Omega)\hookrightarrow C^{2}(\overline{\Omega})$ and for $\epsilon$ sufficiently small we may assume that

[TABLE]

for all deformations $\phi$ considered. Since $\det(D\phi(x))\geq\frac{1}{2}$ for all $x\in\Omega$, this ensures that such deformations are in $\mathcal{A}$ and $C^{1}(\Omega)$-diffeomorphisms (see [5, Theorem 5.5-2]). Furthermore, we obtain

[TABLE]

and deduce from $(D\phi)^{-1}=(\det(D\phi))^{-1}\operatorname{cof}(D\phi)^{T}$

[TABLE]

for deformations $\phi,\tilde{\phi}\in B_{\epsilon}({\mathds{1}})$. Thus, using the Cauchy-Schwarz inequality, the transformation formula, (4.2), (4.3) and (4.4) we achieve the following estimate corresponding to the first term of $\mathcal{T}$:

[TABLE]

Likewise, for the second term of $\mathcal{T}$ we obtain by the transformation formula and by the embedding $H^{2m}(\Omega)\hookrightarrow C^{2}(\overline{\Omega})$

[TABLE]

To conclude, for $\|u_{1}-u_{0}\|_{L^{2}(\Omega)}\leq c_{u}$ and $|u_{1}|_{H^{1}(\Omega)}\leq C_{u}$ the mapping $\mathcal{T}$ is indeed Lipschitz continuous on $B_{\epsilon}({\mathds{1}})\subset\mathcal{A}$ and the Lipschitz constant is bounded by $C(C_{u}c_{u}+c_{u}^{2})$.

(iii) Invertibility of $\mathcal{R}$. The bilinear form

[TABLE]

is bounded in $H^{2m}_{0}(\Omega)$. Furthermore, $\widetilde{\mathcal{R}}$ is coercive since for any $\psi\in H^{2m}_{0}(\Omega)$ we obtain

[TABLE]

due to (2.3) and the iterative application of the Poincaré inequality (cf. [1, Corollary 6.31]). Hence, by the Lax-Milgram Theorem (cf. [10]) there exists for each $z\in H^{-2m}(\Omega)$ a unique $\zeta\in H^{2m}_{0}(\Omega)$ such that $\widetilde{\mathcal{R}}[\zeta](\psi)=z(\psi)$ and $\widetilde{\mathcal{R}}^{-1}:H^{-2m}(\Omega)\rightarrow H^{2m}_{0}(\Omega)$ is a bounded operator. Finally, since $\mathcal{R}[\phi]=\widetilde{\mathcal{R}}[\phi-{\mathds{1}}]$ we can infer that $\mathcal{R}$ is a bounded and invertible operator with inverse $\mathcal{R}^{-1}[z]={\mathds{1}}+\widetilde{\mathcal{R}}^{-1}[z]$.

(iv) Contraction property of $\mathcal{F}$. Using the boundedness of $\mathcal{R}^{-1}$ and the Lipschitz-continuity of $\mathcal{T}$ we obtain for $\mathcal{F}:=\mathcal{R}^{-1}\circ\mathcal{T}$

[TABLE]

for $\phi,\tilde{\phi}\in B_{\epsilon}({\mathds{1}})$, which proves that $\mathcal{F}$ is contractive for sufficiently small $C_{u}$, $c_{u}$ and $\epsilon$.

Next, we prove $\mathcal{F}:B_{\epsilon}({\mathds{1}})\rightarrow B_{\epsilon}({\mathds{1}})$ for a proper choice of $C_{u}$, $c_{u}$ and $\epsilon$. By using the boundedness of $\mathcal{R}^{-1}$ and $\mathcal{R}[{\mathds{1}}]=0$ one can infer

[TABLE]

Thus, for any $\phi\in B_{\epsilon}({\mathds{1}})$ one gets

[TABLE]

Now, choosing $C_{u}$, $c_{u}$ small enough and $\epsilon$ such that the conditions in (4.3) are satisfied for any $\phi\in B_{\epsilon}({\mathds{1}})$ and for $\phi_{1}$, $\mathcal{F}$ maps $B_{\epsilon}({\mathds{1}})$ onto $B_{\epsilon}({\mathds{1}})$.

Hence, the application of Banach’s fixed point theorem proves the existence of a unique deformation $\phi_{2}$ in $B_{\epsilon}({\mathds{1}})\subset\mathcal{A}$ solving (3.7). Then, the unique image $u_{2}$ associated with $(\phi_{1},\phi_{2})$ can be computed using the formula (3.11). Thus, there exists a solution $(u_{2},\phi_{1},\phi_{2})\in H^{1}(\Omega)\times\mathcal{A}\times\mathcal{A}$ of (3.4), (3.5) and (3.6), and this solution is unique in a small neighborhood around $(u_{0},{\mathds{1}},{\mathds{1}})$. ∎

Theorem 4.2** (Local uniqueness and well-posedness of the discrete exponential map).**

Let $2m-\frac{n}{2}>2$ and $u_{0}\in H^{1}(\Omega)$. Then there exist neighborhoods ${\bf{U}}\subset H^{1}(\Omega)$ of $u_{0}$ and ${\bf{D}}\subset\mathcal{A}$ of ${\mathds{1}}$ such that for every $u_{2}\in{\bf{U}}$ there exists at most one solution $(u_{1},\phi_{1},\phi_{2})\in{\bf{U}}\times{\bf{D}}\times{\bf{D}}$ of the equations (3.4)-(3.6). In particular, the discrete exponential map is locally well-posed and

[TABLE]

Proof.

At first, we get rid of the unknown image $u_{1}$. To this end, we observe that the sum of the two matching terms in $\mathcal{W}^{D}[u_{0},u_{1},\phi_{1}]+\mathcal{W}^{D}[u_{1},u_{2},\phi_{2}]$ can be rewritten as follows:

[TABLE]

Therefore, the image $u_{1}$ minimizing the above integral is characterized pointwise a.e. on $\Omega$ by

[TABLE]

and can thus be written as a function $u_{1}(u_{2},\phi_{1},\phi_{2})$ of the image $u_{2}$ and the deformations $\phi_{1}$ and $\phi_{2}$ (we omit the dependence on the image $u_{0}$). Hence, the Euler–Lagrange equations (3.4)-(3.6) can be reformulated as

[TABLE]

where $\mathcal{J}$ is a functional on $H^{1}(\Omega)\times\mathcal{A}\times\mathcal{A}$ with

[TABLE]

for ${h}(A)=\frac{\det A}{1+\det A}$ and $\mathcal{K}:H^{1}(\Omega)\times\mathcal{A}\times\mathcal{A}\to(H^{2m}_{0}(\Omega)\times H^{2m}_{0}(\Omega))^{\prime}$. Now, we will show that in a neighborhood of $(u_{0},{\mathds{1}},{\mathds{1}})$ one obtains an explicit representation $(\phi_{1},\phi_{2})[u_{2}]$ for the implicit equation $0=\mathcal{K}[u_{2},\phi_{1},\phi_{2}]$ via the implicit function theorem. Hence, for every $u_{2}$, which is close to $u_{0}$ in $H^{1}(\Omega)$, there exists in a small neighborhood of $\mathcal{A}\times\mathcal{A}$ a unique tuple $(\phi_{1},\phi_{2})$, which solves the above implicit equation. This indeed proves the claim. To apply the implicit function theorem, we have to show that $\partial_{(\phi_{1},\phi_{2})}\mathcal{K}[u_{0},{\mathds{1}},{\mathds{1}}]=\partial^{2}_{(\phi_{1},\phi_{2})}\mathcal{J}[u_{0},{\mathds{1}},{\mathds{1}}]$ is invertible with bounded inverse. At first, we compute the different components of $\partial^{2}_{(\phi_{1},\phi_{2})}\mathcal{J}[u_{2},\phi_{1},\phi_{2}]$. For this reason, we focus here on the variation of

[TABLE]

the derivatives of the other components of $\mathcal{J}$ are straightforward. We use integration by parts to avoid derivatives of the involved image intensities. For the first variation with respect to $\phi_{1}$ and $\phi_{2}$ we obtain

[TABLE]

using the following different versions of the chain rule:

[TABLE]

Then, for the second order variations one gets

[TABLE]

Evaluating the second order variational derivatives at the point $(u_{0},{\mathds{1}},{\mathds{1}})$ yields

[TABLE]

Here, we have used the following identities, which rely on integration by parts,

[TABLE]

for any vector field $v\in H^{1}_{0}(\Omega,\mathbb{R}^{n})$. Altogether, taking also into account the second order variation of the remaining terms of $\mathcal{J}$ we obtain

[TABLE]

It is straightforward to verify that $\partial^{2}_{(\phi_{1},\phi_{2})}\mathcal{J}[u_{0},{\mathds{1}},{\mathds{1}}]$ is a continuous bilinear form on $H^{2m}_{0}(\Omega)\times H^{2m}_{0}(\Omega)$ by taking into account the estimate

[TABLE]

the coercivity follows by analogous arguments as in the proof of the coercivity of $\mathcal{R}$ (cf. Theorem 4.1). Thus, the Lax-Milgram Theorem ensures the required invertibility. ∎

5 Spatial discretization and fixed point algorithm

In what follows, we introduce a spatial discretization scheme as well as an algorithm to compute the discrete exponential map based on the time discrete operator $\mathrm{EXP}^{2}_{u_{0}}(u_{1}-u_{0})$ for given images $u_{0}$ and $u_{1}$. Let us recall that the computation of $\mathrm{EXP}^{k}$ for $k>2$ requires the iterative application of $\mathrm{EXP}^{2}$ as defined in (3.1). In explicit, we ask for a numerical approximation of the matching deformations $\phi_{1}$, $\phi_{2}$ and the actual succeeding image $u_{2}=\mathrm{EXP}^{2}_{u_{0}}(u_{1}-u_{0})$ along the shot discrete path. Here, we restrict to two dimensional images and for the sake of simplicity we assume that the image domain is the unit square, i.e. $\Omega=(0,1)^{2}$. Conceptually, the generalization to three dimensions is straightforward. As a simplification for the numerical implementation, we restrict to the case $m=1$ despite the theoretical requirement that $m>1+\frac{n}{4}=\frac{3}{2}$. Below, we will introduce the space of tensor product cubic splines for the discretization of deformations. For such discrete deformations the reformulation in Lemma 3.1 (ii) holds true (the regularity result in Proposition 3.2 is only required for the reformulation in the spatially continuous case). We experimentally observed that the spatially discretized model ensures sufficient regularity of the deformations to reliably solve the Euler–Lagrange equations numerically. To sum up, the discrete energy density that we will employ in all numerical computations is given by

[TABLE]

for $u,\tilde{u}\in H^{1}(\Omega)$ and $\phi\in\mathcal{A}$.

The algorithm to compute $\phi_{2}$ is based on a spatially discrete fixed point iteration similar to the one used in proof of Theorem 4.1. In explicit, we follow the derivation of the fixed point mapping in this proof using now the reformulation (ii) instead of (i) in Lemma 3.1 as a starting point and define

[TABLE]

for all $\psi\in H^{2m}_{0}(\Omega)$. Here, we used integration by parts to get the second equality. This ansatz is numerically beneficial because it avoids the evaluation of gradients of image intensities. In fact, we experimentally observed that the evaluation of the expression $\int_{\Omega}(u_{1}\circ\phi_{1}-u_{0})(\nabla u_{1}\cdot(D\phi)^{-1}\psi)\circ\phi_{1}{\,\mathrm{d}}x$ appearing in the definition of $\mathcal{T}$ in proof of Theorem 4.1 suffers from accuracy problems in the proximity of interfaces of $u_{1}$ due to the approximate numerical quadrature. To further improve the stability of the numerical algorithm with respect to the evaluation of the first integrand, we additionally rewrite this expression by making use of $A:B=\mathrm{tr}(A^{T}B)$ as follows

[TABLE]

The second operator $\mathcal{R}$ is chosen identically to the one in the proof of Theorem 4.1. Then, taking into account the identity $\widetilde{\mathcal{T}}[\phi_{2}](\psi)=\mathcal{R}[\phi_{2}](\psi)$ for all test functions $\psi\in H^{2m}_{0}(\Omega)$ the modified fixed point equation based on (3.8) reads as

[TABLE]

for $j\in\mathbb{N}$.

We use different discrete ansatz spaces for the deformations and the images. As the discrete ansatz space for deformations we choose the conforming space of cubic splines ${\mathcal{S}_{H}}\subset C^{2}(\Omega)$. Here, $H=2^{-N}$ with $N\in\mathbb{N}$ denotes the grid size of the underlying uniform and rectangular mesh, and the basis functions are vector-valued B-splines. Moreover, we only impose the Dirichlet boundary condition $\mathbf{\Phi}={\mathds{1}}$ on $\partial\Omega$ instead of the stronger boundary conditions $\mathbf{\Phi}-{\mathds{1}}\in H^{2}_{0}(\Omega)$ for the discrete deformations $\mathbf{\Phi}\in{\mathcal{S}_{H}}$. Indeed, we experimentally observed that these Dirichlet boundary conditions allow to reliably compute proper deformations. The gray value images are approximated with finite element functions in the space ${\mathcal{V}_{h}}$ of piecewise bilinear and globally continuous functions on $\Omega$ with input intensities in the range $[0,1]$. The underlying grid consists of uniform and quadratic cells with mesh size $h=2^{-M}$ with $M>N$, the index set of all grid nodes is denoted by ${I_{\mathcal{V}_{h}}^{N}}$. We take into account the usual Lagrange basis functions $\{\Theta^{i}\}_{i\in{I_{\mathcal{V}_{h}}^{N}}}$ to represent image intensities $\mathbf{U}\in{\mathcal{V}_{h}}$. In our numerical experiments we set $M=N+1$.

Now, we are in the position to define spatially discrete counterparts of the energy and the operators involved in the fixed point iteration. We apply a Gaussian quadrature of order $5$ on both meshes. The discrete energy for $\mathbf{U},\tilde{\mathbf{U}}\in{\mathcal{V}_{h}}$ and $\mathbf{\Phi}\in{\mathcal{S}_{H}}$ is defined as (cf. (2.2))

[TABLE]

where we sum over all grid cells $c_{H}$ of the spline mesh and all local quadrature points within these cells indexed by $q_{H}$ with respect to the deformation energy and over all grid cells $c_{h}$ of the finer finite element mesh and all local quadrature points within these cells indexed by $q_{h}$. Here, $(\omega^{c_{H}}_{q_{H}},\mathbf{x}^{c_{H}}_{q_{H}})$ and $(\omega^{c_{h}}_{q_{h}},\mathbf{x}^{c_{h}}_{q_{h}})$ are the pairs of quadrature weights and points on the spline mesh and the finite element mesh, respectively. For the fully discrete counterparts of the operators $\widetilde{\mathcal{T}}$ and $\mathcal{R}$ one gets

[TABLE]

for $\mathbf{\Psi}\in{\mathcal{S}_{H}}$ with $\mathbf{\Psi}=0$ on $\partial\Omega$. Finally, one obtains the following fixed point iteration to compute the spatially discrete $\mathbf{\Phi}_{2}$

[TABLE]

for all $j\geq 0$ and initial data $\mathbf{\Phi}^{0}={\mathds{1}}$. The application of $\mathbf{R}^{-1}$ requires the solution of the associated linear system of equations.

In a preparatory step, the deformation $\mathbf{\Phi}_{1}\in\operatorname*{argmin}_{\mathbf{\Phi}\in{\mathcal{S}_{H}}}\mathbf{W}^{D}[\mathbf{U}_{0},\mathbf{U}_{1},\mathbf{\Phi}]$, which is used in the first step of a time discrete geodesic shooting, is calculated using a Fletcher-Reeves nonlinear conjugate gradient descent multilevel scheme with an Armijo step size control.

Then, the deformation in the current step is computed using the fixed point iteration (5.1), which is stopped if the $L^{\infty}$-difference of the deformations in two consecutive iterations is below the threshold value $\mathrm{THRESHOLD}=10^{-12}$. To compute $\mathbf{U}_{2}$, we employ the spatially discrete analog of the update formula (3.11)

[TABLE]

Here, we evaluate (5.2) at all grid nodes of the finite element grid. To compute approximate inverse deformations $\mathbf{\Phi}_{i}^{-1}\in{\mathcal{S}_{H}}$, $i\in\{1,2\}$, all cells of the grid associated with ${\mathcal{S}_{H}}$ are traversed and the deformed positions $\mathbf{\Phi}_{i}(\mathbf{x}_{j})$ for all vertices $\mathbf{x}_{j}$, $j\in\{1,\ldots,4\}$, of the current element are computed. Then, we use a bilinear interpolation of these deformed positions to define an approximation of $\mathbf{\Phi}_{i}^{-1}(x)$ for $x\in\Omega$. Furthermore, we explicitly ensure the boundary condition $\mathbf{\Phi}_{i}^{-1}(x)=x$ for $x\in\partial\Omega$.

In our numerical experiments on real image data, we observed slight local oscillations emerging from the inexact evaluation of the expression $\mathbf{J}_{k}=\mathbf{U}_{k}(\mathbf{x})-\mathbf{U}_{k-1}\circ\mathbf{\Phi}_{k}^{-1}(\mathbf{x})$ in the quadrature of the intensity modulation. Since the calculation of $\mathrm{EXP}^{k}$ requires a recursive application of $\mathrm{EXP}^{2}$, these oscillations turn out to be sensitive to error propagation, and it is advantageous to apply in a post-processing step one iteration of the anisotropic diffusion filter $(\mathbf{M}_{h}+\tau\mathbf{S}_{h}[\mathbf{J}_{k},\lambda])^{-1}\mathbf{M}_{h}$ to $\mathbf{J}_{k}$ for weight parameters $\tau,\lambda>0$ (see [20]). Here, $\mathbf{M}_{h}$ is the usual mass matrix and $\mathbf{S}_{h}[\mathbf{J}_{k},\lambda]$ the anisotropic stiffness matrix associated with ${\mathcal{V}_{h}}$, i.e. $(\mathbf{S}_{h}[\mathbf{J}_{k},\lambda])_{i,j}=\sum_{c_{h},q_{h}}\omega^{c_{h}}_{q_{h}}(1+\lambda^{-2}\|\nabla\mathbf{J}_{k}(\mathbf{x}^{c_{h}}_{q_{h}})\|^{2})^{-1}\nabla\Theta^{i}(\mathbf{x}^{c_{h}}_{q_{h}})\cdot\nabla\Theta^{j}(\mathbf{x}^{c_{h}}_{q_{h}})$ for $i,j\in{I_{\mathcal{V}_{h}}^{N}}$. Furthermore, in all following applications except the first test case (Figure 2) we choose $\tau=\beta^{k-2}\cdot 10^{-3}$ as the exponentially decaying time step size ($k$ denoting the index of the image in the sequence and $\beta=0.8$) and $\lambda=0.5$ as the smoothing parameter along the discrete geodesic. The impact of this filtering can be seen for instance in Figure 3.

6 Numerical results

In this section, we present applications of the fully discrete exponential map. In all computations, we use the parameters $\gamma=10^{-4}$ and $\delta=10^{-2}$.

As a first example, we investigate an artificial test case consisting of an input image $\mathbf{U}_{0}$ with three ellipses of different intensities and an associated variation $\mathbf{U}_{1}-\mathbf{U}_{0}$. The first row in Figure 2 depicts distinct images of the image sequence $\mathrm{EXP}^{k}_{\mathbf{U}_{0}}(\mathbf{U}_{1}-\mathbf{U}_{0})$ for time steps $k=0,1,2,3,6,9$, the input images $\mathbf{U}_{0}$ and $\mathbf{U}_{1}$ with resolution $257\times 257$ are framed in red. In the initial variation $\mathbf{U}_{1}-\mathbf{U}_{0}$ underlying the exponential shooting, the upper left ellipse is slightly translated to the bottom and simultaneously expanded. The upper right ellipse is undergoing a small rotation and the third one is also slightly translated with some modulation of the shading. The initial variation encoded in the image pair $(\mathbf{U}_{0},\mathbf{U}_{1})$ is prolongated along the sequence generated by an iterative application of the discrete exponential map $\mathrm{EXP}^{2}$. In the second row, for each $k>0$ the discrete intensity modulations $\mathbf{I}_{k}=\mathbf{U}_{k}\circ\mathbf{\Phi}_{k}-\mathbf{U}_{k-1}$ are visualized. Here, on the left the color bar with bounds coinciding with the extremal values for this image sequence is displayed. The third row depicts the discrete velocity fields $\mathbf{v}_{k}=\frac{1}{\tau}(\mathbf{\Phi}_{k}-{\mathds{1}})$ for each $k$, where $\tau=\frac{1}{K}$ is the associated time step size. Here, the hue refers to the direction and the color intensity is proportional to the local norm of $\mathbf{v}_{k}$ as represented by the leftmost color wheel. In particular, one observes that the resulting underlying velocity field $\mathbf{v}_{k}$ is not constant in time.

In [4, Figure 6.2], a geodesic sequence between two female portrait paintings222first painting by A. Kauffmann (public domain, see http://commons.wikimedia.org/wiki/File:Angelika_Kauffmann_-_Self_Portrait_-_1784.jpg), second painting by R. Peale (GFDL, see http://en.wikipedia.org/wiki/File:Mary_Denison.jpg) was computed using the finite element discretization for both the images and the deformations on the same grid. The image resolution is $257\times 257$ ($M=8$). We recomputed this geodesic sequence $(\tilde{\mathbf{U}}_{0},\tilde{\mathbf{U}}_{1},\ldots,\tilde{\mathbf{U}}_{16})$ with $K=16$ for the discrete function spaces ${\mathcal{V}_{h}}$ and ${\mathcal{S}_{H}}$ with $N=7$, the resulting sequence is shown in the first row of Figure 3 with framed input images $\tilde{\mathbf{U}}_{0}$ and $\tilde{\mathbf{U}}_{16}$. This is compared with the discrete exponential shooting for the initial image pair $\mathbf{U}_{0},\,\mathbf{U}_{1}$ taken from this geodesic sequence. One observes that the discrete exponential map is capable to recover the original geodesic sequence for small time steps $k$. Only in late stages visible differences become apparent. Again, we highlight that the discrete motion fields significantly alter in time.

Figure 4 depicts a picture details for the time steps $k=0,1,2,4,8,16$ of the discrete exponential map applied to two different pairs of photos. These photos show human faces and small variations of them and the resolution of the underlying full images is $1025\times 1025$. The red boxes indicate these input images (first and third row), which are consecutive photos of a series at $5$ and $7$ fps, respectively, taken with a digital camera. We observe that small initial variations result in a nonlinear deformation of the lips (first row) and of the lips, the cheeks and the eyes (third row), respectively. Furthermore, the textures are transported along the sequence. The second and fourth row depict the color coded time varying velocity fields.

Figure 5 shows the discrete exponential map for $k=0,\ldots,7$ applied to a pair of images of a dog for a resolution of $1025\times 1025$. Again, the input pictures are consecutive photos of a series with $7$ fps taken with a digital camera. The initial image pair shows a slight rotation of the dog’s head and a small opening of its eyes. The proposed algorithm generates a extrapolation of this movement. As a consequence in particular of the rotation, the method fills in reasonable image features below the mouth and right to the ear which correspond to hidden object regions in $\mathbf{U}_{0}$ and $\mathbf{U}_{1}$.

Supplementary material.

The supplementary material of this publication includes video sequences with animations of the discrete exponential map shown in Figure 3 (together with the discrete geodesic interpolation), Figure 4 and Figure 5.

Acknowledgements.

A. Effland and M. Rumpf acknowledge support of the Hausdorff Center for Mathematics and the Collaborative Research Center 1060 funded by the German Research Foundation.