# Discrete geodesic calculus in the space of viscous fluidic objects

Based on a local approximation of the Riemannian distance on a manifold by a computationally cheap dissimilarity measure, a time discrete geodesic calculus is developed, and applications to shape space are explored. The dissimilarity measure is derived from a deformation energy whose Hessian reproduces the underlying Riemannian metric, and it is used to define length and energy of discrete paths in shape space. The notion of discrete geodesics defined as energy minimizing paths gives rise to a discrete logarithmic map, a variational definition of a discrete exponential map, and a time discrete parallel transport. This new concept is applied to a shape space in which shapes are considered as boundary contours of physical objects consisting of viscous material. The flexibility and computational efficiency of the approach is demonstrated for topology preserving shape morphing, the representation of paths in shape space via local shape variations as path generators, shape extrapolation via discrete geodesic flow, and the transfer of geometric features.

## Authors

• 14 publications
• 10 publications
10/28/2019

### Image Morphing in Deep Feature Spaces: Theory and Applications

This paper combines image metamorphosis with deep features. To this end,...
12/16/2019

### Consistent Curvature Approximation on Riemannian Shape Spaces

We describe how to approximate the Riemann curvature tensor as well as s...
02/13/2016

### Manifolds of Projective Shapes

The projective shape of a configuration of k points or "landmarks" in RP...
08/21/2021

### ARAPReg: An As-Rigid-As Possible Regularization Loss for Learning Deformable Shape Generators

This paper introduces an unsupervised loss for training parametric defor...
02/16/2021

### ResNet-LDDMM: Advancing the LDDMM Framework Using Deep Residual Networks

In deformable registration, the geometric framework - large deformation ...
04/08/2016

### On the Hessian of Shape Matching Energy

In this technical report we derive the analytic form of the Hessian matr...
06/27/2019

### Geodesic analysis in Kendall's shape space with epidemiological applications

We analytically determine Jacobi fields and parallel transports and comp...
##### This week in AI

Get the week's most popular data science and artificial intelligence research sent straight to your inbox every Saturday.

## 1 Introduction

Geodesic paths in shape space allow to define smooth and in some sense geometrically or physically natural connecting paths , , between two given shapes , or they enable the extrapolation of a path from an initial shape and an initial shape variation

which encodes the path direction. Applications include shape modeling in computer vision

[17, 16], computational anatomy, where the morphing path establishes correspondences between a patient and a template [2, 26], shape clustering based on Riemannian distances [32], as well as shape statistics [9, 13], where geodesic paths in shape space transport information from the observed shapes into a common reference frame in which statistics can be performed.

As locally length minimizing paths, geodesic paths require to endow the space of shapes with a Riemannian metric which encodes the preferred shape variations. There is a rich diversity of Riemannian shape spaces in the literature. Kilian et al. compute isometry invariant geodesics between consistently triangulated surfaces [16], where the Riemannian metric measures stretching of triangle edges, while the metric by Liu et al. also takes into account directional changes of edges [22].

For planar curves, different Riemannian metrics have been devised, including the -metric on direction and curvature functions [17], the -metric on stretching and bending variations [31], as well as curvature-weighted - or Sobolev-type metrics [25, 34], some of which allow closed-form geodesics [37, 33]. A variational approach to the computation of geodesics in the space of planar Jordan curves has been proposed by Schmidt et al. in [30]. The extrapolation of geodesics in the space of curves incorporating translational, rotational, and scale invariance has been investigated by Mennucci et al. [24]. A Riemannian space of non-planar elastic curves has very recently been proposed by Srivastava et al. [32].

In the above approaches, the shape space is often identified with a so-called pre-shape space of curve parameterizations over a 1D domain (or special representations thereof) modulo the action of the reparameterization group. It is essential that the metric on the pre-shape space is invariant under reparameterization or equivalently that reparameterization represents an isometry in the pre-shape space so that the Riemannian metric can be inherited by the shape space. Such reparameterization-invariant metrics can also be defined on the space of parameterized 2D surfaces [1, 20]. For certain representations of the parameterization one is lead to a very simple form of the metric, e.g. an -type metric [19].

The issue of reparameterization invariance does not occur when the mathematical description of the shape space is not based on parameterizations, which often simplifies the analysis (and is also the approach taken here). When warping objects in , a shape tube in is formed. Zolésio investigates geodesic in terms of shortest shape tubes [39]. The space of sufficiently smooth domains can be assigned a Riemannian metric by identifying the tangent space at with velocity fields and defining a metric on these. Dupuis et al. employ a metric

 G(v,v)=∫DLv⋅vdx

for a higher order elliptic operator on some computational domain [7], ensuring a diffeomorphism property of geodesic paths. A corresponding geodesic shooting method has been implemented in [3]. Fuchs et al. propose a viscous-fluid based Riemannian metric [12]. Fletcher and Whitaker employ a similar metric on pullbacks of velocity fields onto a reference shape [10]. Miller and Younes consider the space of registered images as the product space of the Lie group of diffeomorphisms and image maps. They define a Riemannian metric using sufficiently regular elliptic operators on the diffeomorphism-generating velocity fields, which may also depend on the current image [27]. A morphing approach based on the concept of optimal mass transport has been proposed by Haker et al. [14, 38]. An image or a shape is viewed as mass density, and for two such densities the Monge–Kantorovich functional

 ∫D|ψ(x)−x|2ρ0(x)dx

is minimized over all mass preserving mappings , i.e. mappings with . A morphing path then is given by for , . Like for our approach there is a continuum-mechanical interpretation of minimizing the action of an incompressible fluid flow [4], however, the flow typically does neither preserve local shape features or isometries nor the shape topology.

Very often, geodesics in shape space are approached via the underlying geodesic evolution equation, and geodesics between two shapes are computed by solving this ODE within a shooting method framework [17, 3, 1]. An alternative approach exploits the energy-minimizing property of geodesics: Schmidt et al. perform a Gauß-Seidel type fixed-point iteration which can be interpreted as a gradient descent on the path energy, and Srivastava et al. derive the equations of a gradient flow for the path energy which they then discretize [32]. In contrast, we employ an inherently variational formulation where geodesics are defined as minimizers of a time discrete path energy. Discrete geodesics are then defined consistently as minimizers of a corresponding discrete energy.

In this paper we start from this time discretization and consistently develop a time discrete geodesic calculus in shape space. The resulting variational discretization of the basic Riemannian calculus consists of an exponential map, a logarithmic map, parallel transport, and finally an underlying discrete connection. To this end, we replace the exact, computationally expensive Riemannian distance by a relatively cheap but consistent dissimilarity measure. Our choice of the dissimilarity measure not only ensures consistency for vanishing time step size but also a good representation of shape space geometry already for coarse time steps. For example, rigid body motion invariance is naturally incorporated in this approach. We illustrate this approach on a shape space consisting of homeomorphic viscous-fluid objects and a corresponding deformation-based dissimilarity measure.

Different from most approaches, which first discretize in space via the choice of a parameterization, a set of control points, or a mesh, and then solve the resulting transport equations by suitable solvers for ordinary differential equations (see the discussion above), our time discretization is defined on the usually infinite dimensional shape space. It results from a consistent transfer of time continuous to time discrete variational principles. Thereby, it leads to a collection of variational problems on the shape space, which in our concrete implementation of the proposed calculus consists of non-parameterized, volumetric objects.

Let us also already mention a further remarkable conceptual difference. The way the time discrete geodesic calculus is introduced differs substantially from the way the time continuous counterpart is usually developed. In classical Riemannian differential geometry one first defines a connection

for two vector fields

and on a manifold . With the connection at hand a tangent vector can be transported parallel along a path with motion field solving . Studying those paths where the motion field itself is transported parallel along the path (i.e. it solves the ODE ) one is led to geodesics. Next, the exponential map is introduced via the solution of the above ODE for varying initial velocity. Finally, the logarithm is obtained as the (local) inverse of the exponential map.

In the time discrete calculus we start with a time discrete formulation of path length and energy and then define discrete geodesics as minimizers of the discrete energy. Evaluating the initial step of a discrete geodesic path as the discrete counterpart of the initial velocity we are led to the discrete logarithm. Then, the discrete exponential map is defined as the inverse of the discrete logarithm. Next, discrete logarithm and discrete exponential allow to define a discrete parallel transport based on the construction of a sequence of approximate Riemannian parallelograms (commonly known as Schild’s ladder [8]). Finally, with the discrete parallel transport at hand, a time discrete connection can be defined.

Let us note that the approximation of parallel transport in shape space via Schild’s ladder has also been used in the context of the earlier mentioned flow of diffeomorphism approach [29, 23]. In our discrete framework, however, the notion of discrete parallel transport is directly derived from the parallelogram construction, consistently with the overall discrete approach to geodesics.

A related approach for time discrete geodesics has been presented in an earlier paper by Wirth et al. [36]. In contrast to [36], we here do not restrict ourselves to the computation of geodesic paths between two shapes but devise a full-fledged theory of discrete geodesic calculus (cf. Figure 1

). Furthermore, different from that approach we ensure topological consistency and describe shapes solely via deformations of reference objects instead of treating deformations and level set representations of shapes simultaneously as degrees of freedom, which in turn strongly simplifies the minimization procedure.

The paper is organized as follows. In Section 2 we introduce a special model for a shape space, the space of viscous fluidic objects, to which we restrict our exposition of the geodesic calculus. Here, in the light of the discrete shape calculus to be developed, we will review the notion of discrete path length and discrete path energy. After these preliminaries the actual time discrete calculus consisting of a discrete logarithm, a discrete exponential and a discrete parallel transport together with a discrete connection is introduced and discussed in Section 3. Then, Section 4

is devoted to the numerical discretization via characteristic functions and a parameterization via deformations over reference paths. Finally, we draw conclusions in Section

5.

## 2 A space of volumetric objects and an elastic dissimilarity measure

To keep the exposition focused we restrict ourselves to a specific shape model, where shapes are represented by volumetric objects which behave physically like viscous fluids. In fact, the scope of the variational discrete geodesic calculus extends beyond this concrete shape model. We refer to Section 5 for remarks on the application to more general shape spaces.

### 2.1 The space of viscous-fluid objects

Let us introduce the space of shapes as the set of all objects which are closed subsets of () and homeomorphic to a given regular reference object , i.e. for an orientation preserving homeomorphism . Furthermore, objects which coincide up to a rigid body motion are identified with each other. A smooth path in this shape space is associated with a smooth family of deformations. To measure the distance between two objects, a Riemannian metric is defined on variations of objects which reflects the internal fluid friction — called dissipation — that occurs during the shape variation. The local temporal rate of dissipation in a fluid depends on the symmetric part of the gradient of the fluid velocity (the antisymmetric remainder reflects infinitesimal rotations), and for an isotropic Newtonian fluid, we obtain the local rate of dissipation

 diss(∇v)=λ(trϵ[v])2+2μtr(ϵ[v]2), (1)

where are material-specific parameters. Given a family of deformations of the reference object , the change of shape along the path can be described by the (Lagrangian) temporal variation or the associated (Eulerian) velocity field

 v(t)=˙ϕ(t)∘ϕ−1(t)

on . Hence, the tangent space to at a shape can be identified with the space of initial velocities for deformation paths with . Here we identify those velocities which lead to the same effective shape variation, i.e. those with the same normal component on , where is the outer normal on . Now, integrating the local rate of dissipation for velocity fields on , we define the Riemannian metric on as the symmetric quadratic form with

 GO(v,v)=min{~v|~v⋅n=v⋅non∂O}∫Odiss(∇~v(x))dx. (2)

For the shape variation along a path described by the Eulerian motion field , path length and energy are defined as

 L[(O(t))t∈[0,1]] =∫10√GO(t)(v(t),v(t))dt, (3) E[(O(t))t∈[0,1]] =∫10GO(t)(v(t),v(t))dt. (4)

Paths which (locally) minimize the energy or equivalently the length are called geodesics (cf. Figure 2). A geodesic thus mimics the energetically optimal way to continuously deform a fluid volume.

### 2.2 Approximating the distance

The evaluation of the geodesic distance based on a direct space and time discretization of (2) and (3) turns out to be computationally very demanding (cf. for instance the approaches in [3, 7]). Hence, we use here an efficient and robust time discrete approximation based on an energy functional which locally behaves like the squared Riemannian distance (i.e. the squared length of a connecting geodesic):

Given two shapes and , we consider an approximation

 dist2(O,~O)≈WO[ψ], (5)

where is the stored deformation energy of a deformation and is the minimizer of this energy over all such deformations with . Here, is a so-called hyperelastic energy density.

In correspondence to our assumption that objects are identical if they coincide up to a rigid body motion, we require to be rigid body motion invariant. Furthermore, we assume the objects to have no preferred material directions so that is in addition isotropic, which altogether leads to for all (cf. [5]). In the undeformed configuration for , energy and stresses (the first derivatives of ) are supposed to vanish so that we require , (where denotes the derivative with respect to the matrix argument). Furthermore, we need as to prohibit material self-penetration, which is linked to the preservation of topology. The approximation property (5) relies on a consistent choice of for the given metric which can be expressed by the relation

 12d2dt2WO[ψ(t)]∣∣∣t=0=12∫OD2W(1)(∇v,∇v)dx=∫Odiss(∇v)dx (6)

along any object path , , with and velocity field . Using the notion of the Hessian of a function on a manifold as the endomorphism representing its second variation in the metric, we can rephrase this approximation condition more geometrically as

 12HessMWO[id]=id

with the usual identification of objects and deformations . For the deformation energy density , this condition implies that its Hessian has to satisfy for all . A suitable example is

 W(A)=μ2tr(ATA)+λ4detA2−(μ+λ2)logdetA−dμ2−λ4.

Assume that the energy density satisfies the above-mentioned properties. We observe that the metric is the first non-vanishing term in the Taylor expansion of the squared length of a curve, i.e.

 (L[(O(t))t∈[0,T]])2=T2GO(0)(v,v)+O(T3)

with being the initial tangent vector along a smooth path . Thus, since the Hessian of the energy and the metric are related by (6), we obtain that the second order Taylor expansions of and in coincide and indeed

 dist2(O,~O)=min{ψ|ψ(O)=~O}WO[ψ]+O(dist3(O,~O)). (7)

Here, different from [36] we neither take into account mismatch penalties nor perimeter regularizing functionals for each object , .

### 2.3 Discrete length and discrete energy

Now, we are in a position to discretize length and energy of paths in shape space. To this end, we first sample the path at times for (), denote

, and obtain the estimates

 L[(O(t))t∈[0,1]] ≥ ∑Kk=1dist(Ok−1,Ok) E[(O(t))t∈[0,1]] ≥ 1τ∑Kk=1dist2(Ok−1,Ok)

for the length and the energy, where equality holds for geodesic paths. Indeed, the first estimate is straightforward, and the application of the Cauchy–Schwarz inequality leads to

 K∑k=1dist2(Ok−1,Ok) ≤ K∑k=1(∫kτ(k−1)τ√GO(t)(v(t),v(t))dt)2 ≤ K∑k=1τ∫kτ(k−1)τGO(t)(v(t),v(t))dt=τE[(O(t))t∈[0,1]]

which implies the second estimate.

Together with (7) this motivates the following definition of a discrete path energy and a discrete path length for a discrete path in shape space:

 L[(O0,…,OK)] =∑Kk=1√WOk−1[ψk], (8) E[(O0,…,OK)] =1τ∑Kk=1WOk−1[ψk], (9)

where (cf. also [36]). In fact, (8) and (9) can for general smooth paths even be proven to be first order consistent with the continuous length (3) and energy (4) as . For illustration, if is a two-dimensional manifold embedded in , we can interpret the terms as the stored elastic energies in springs which connect a sequence of points on the manifold through the ambient space. Then the discrete path energy is the total stored elastic energy in this chain of springs.

A discrete geodesic (of order ) is now defined as a minimizer of for fixed end points . The discrete geodesic is thus an energetically optimal sequence of deformations from into .

In the minimization algorithm to be discussed in Section 4.1 we do not explicitly minimize for the object contours as in [36] but instead for reference deformations defined on fixed reference objects. Figure 3 shows a discrete geodesic in the context of multicomponent objects, which is visually identical to that obtained by the more complex approach in [36]. Here, deformations are considered which map every component of a shape onto the corresponding component of the next shape in the discrete path as the obvious generalization of discrete geodesics between single component shapes.

While in the continuous case geodesic curves equally minimize length (3) and energy (4), minimizers of the discrete path length (9) are in general not related to discrete geodesics (and thus also not to continuous geodesics as ). Indeed, let us consider a two-dimensional manifold embedded in , paired with the deformation energy for a displacement vector in connecting points and on . Now take into account a continuous geodesic and a discrete path on where the end points are close to each other in the embedding space but far apart on the surface. Figure 4 depicts such a configuration with a discrete path which almost minimizes the discrete path length. A minimizer of the discrete path length will always jump through the protrusion and never approximate the continuous geodesic, whereas minimizers of the discrete path energy satisfy as and thus rule out such a short cut through the ambient space.

## 3 Time discrete geodesic calculus

With the notion of discrete geodesics at hand we will now derive a full-fledged discrete geodesic calculus based on a time discrete geometric logarithm and a time discrete exponential map, which then also give rise to a discrete parallel transport and a discrete Levi-Civita connection on shape space.

### 3.1 Discrete logarithm and shape variations

If is the unique geodesic on connecting and , the logarithm of with respect to is defined as the initial velocity of the geodesic path. In terms of Section 2.1 we have

 logO(~O)=v(0)

for , where defines the associated family of deformations. On a geodesically complete Riemannian manifold the logarithm exists as long as is sufficiently small. The associated logarithmic map represents (nonlinear) variations on the manifold as (linear) tangent vectors.

The initial velocity can be approximated by a difference quotient in time,

 v(0,x)=1τζ(x)+O(τ),

where denotes a displacement on the initial object . Thus, we obtain

 τlogO(~O)=ζ(x)+O(τ2).

This gives rise to a consistent definition of a time discrete logarithm. Let be a discrete geodesic between and with an associated sequence of matching deformations , then we consider for the displacement as an approximation of . Taking into account that we thus define the discrete -logarithm

 (1KLOG)O(~O):=ζ1. (10)

In the special case and a discrete geodesic we simply obtain

As in the continuous case the discrete logarithm can be considered as a representation of the nonlinear variation of in the (linear) tangent space of displacements on . On a sequence of successively refined discrete geodesics we expect

 (11)

for (cf. Figure 5 for an experimental validation of this convergence behaviour).

### 3.2 Discrete exponential and shape extrapolation

In the continuous setting, the exponential map maps tangent vectors onto the end point of a geodesic with and the given tangent vector at time . Using the notation from the previous section we have and, via a simple scaling argument, for . We now aim at defining a discrete power exponential map such that on a discrete geodesic of order with (the notation is motivated by the observation that on or more general matrix groups). Our definition will reflect the following recursive properties of the continuous exponential map,

 expO(1v)= (11logO)−1(v), expO(2v)= (12logO)−1(v), expO(kv)= expexpO((k−2)v)(2vk−1) for vk−1:=logexpO((k−2)v)expO((k−1)v).

Replacing by , by , and by we obtain the recursive definition

 EXP1O(ζ):= (11LOG)−1O(ζ), (12) EXP2O(ζ):= (12LOG)−1O(ζ), (13) EXPkO(ζ):= EXP2EXPk−2O(ζ)(ζk−1) (14)

It is straightforward to verify that as long as the discrete logarithm on the right is invertible. Equation (12) implies , and (13) in fact represents a variational constraint for a discrete geodesic flow of order  :

Given the object we consider discrete geodesic paths of order , where for any chosen the middle object is defined via minimization of (9) so that we may write . We now identify as the object for which , i.e. is the energetically optimal displacement from to and thus satisfies

 id+ζ=argmin{ψ1|ψ1(O)=~O1[~O2]}WO[ψ1] (15)

up to a rigid body motion.

Alternatively, the condition (15) can be phrased as

 id+ζ=argmin{ψ1}min{ψ2|(ψ2∘ψ1)(O)=~O2}(WO[ψ1]+Wψ1(O)[ψ2]). (16)

Figure 6 conceptually sketches the procedure to compute . For given initial object and initial displacement the discrete exponential is selected from a fan of discrete geodesics with varying as the terminal point of a discrete geodesic of order in such a way that (16) holds. To compute in the geodesic flow algorithm (13) and (14) we have to find the root of

 FO,ζ(~O2)=(12LOG)O(~O2)−ζ, (17)

implicitly assuming that is small enough so that discrete geodesics are unique (cf. Section 4.2 for the algorithmic realization based on a representation of the unknown domain via a deformation).

Equation (14) describes the recursion to compute based on the above single step scheme: For given and one first retrieves from the previous step, where

 ψk−1=argmin{ψ|ψ(Ok−2)=Ok−1}WOk−2[ψ].

Then (13) is applied to compute from and as the root of .

For sufficiently small we expect to be well-defined. Since by definition, every triplet of the sequence is a geodesic of order 2 and minimizes , the resulting family indeed is a discrete geodesic of order .

In fact, discrete geodesics that are variationally described as discrete energy minimizing paths between two given objects can be reproduced via the discrete geodesic flow associated with the discrete exponential map (cf. Figure 7).

As for the discrete logarithm we experimentally observe convergence of the discrete exponential map in the sense

 EXPkO(1kζ)→expO(ζ)for k→∞ (18)

as shown in Figure 5. An example of geodesic shape extrapolation for multicomponent objects is depicted in Figure 8.

### 3.3 Discrete parallel transport and detail transfer

Parallel transport allows to translate a vector (which is considered as the variation of an object ) along a curve in shape space. The resulting changes as little as possible while keeping the angle between and the path velocity fixed. Using the Levi-Civita connection this can be phrased as . There is a well-known first-order approximation of parallel transport called Schild’s ladder [8, 15], which is based on the construction of a sequence of geodesic parallelograms, sketched in Figure 9, where the two diagonal geodesics always meet at their midpoints. Given a curve and a tangent vector , the approximation of the parallel transported vector via a geodesic parallelogram can be expressed as

 Opk−1 =expO((k−1)τ)ζk−1, O×k =expOpk−112logOpk−1O(kτ), Opk =expO((k−1)τ)2logO((k−1)τ)O×k, ζk =logO(kτ)Opk.

Here, is the midpoint of the two diagonals of the geodesic parallogramm with vertices , , , and . This scheme can be easily transferred to discrete curves in shape space based on the discrete logarithm and the discrete exponential introduced above. In the th step of the discrete transport we start with a displacement on and compute

 Opk−1 =EXP1Ok−1ζk−1, O×k =EXP1Opk−1((12LOG)Opk−1(Ok)), Opk =EXP2Ok−1((11LOG)Ok−1(O×k)), ζk =(11LOG)Ok(Opk),

where is the transported displacement on . Here, is the midpoint of the two discrete geodesics of order with end points , and , , respectively. Since the last of the above steps is the inverse of the first step in the subsequent iteration, these steps need to be performed only for . We will denote the resulting transport operator by . Figure 10 shows examples of discrete parallel transport for feature transfer along curves in shape space.

Remark: As in the continuous case, the discrete parallel transport can be used to define a discrete Levi-Civita connection. For and for a vector field in the tangent bundle one computes and then defines

 ∇τξη:=1τ(PO,Oτη(Oτ)−η(O))

as the time discrete connection with time step size .

## 4 Numerical discretization

The proposed discrete geodesic calculus requires an effective and efficient spatial discretization of

• volumetric objects in the underlying shape space,

• of nonlinear deformations to encode matching correspondences,

• and of linear displacements as approximate tangent vectors.

We restrict ourselves here to the case of objects . To this end we consider the space of piecewise affine finite element functions on a regular simplicial mesh over a rectangular computational domain . Here indicates the grid size, where in our applications ranges from a coarse grid size to a fine grid size . Then, deformations and displacements are considered as functions in . Objects , the original degrees of freedom in our geometric calculus, will be represented via deformations over reference objects (e.g. ), i.e. . These reference objects are encoded by approximate characteristic functions and the deformations are considered as injective deformations and discretized as elements in .

### 4.1 Parameterization of discrete geodesics

To compute a discrete geodesic — different from [36] — we now replace the objects as arguments of the energy (9) by associated deformations over a set of reference domains as described above. By this technique, instead of deformations and domain descriptions (e.g. via level sets) as in [36] we will be able to consider solely parameterizing deformations, which turns out to be a significant computational advantage. Next, we assume that reference matching deformations are given with (cf. Figure 11).

Now, we express the matching deformations over which we minimize in (9) in terms of the parameterizing deformations and the reference matching deformations and set

 ψk