The Shapespace of Repulsive Shells

April 1, 2025

Figure 1: Computing the self-intersection free interpolation of two hand poses given the start and end pose

Abstract

This lab report discusses the computation of interpolations

between diﬀerent embeddings of a triangular mesh called

shapes. This is done by constructing a shapespace, which is

a Riemannian manifold, so that geodesics on it correspond to

"good" interpolations. This report introduces and discusses

the implementation of diﬀerent Riemannian metrics that

correspond to diﬀerent characteristics of the interpolation. In

particular, we discuss a simpliﬁed version of the framework

developed by Sassen et al. (2024), which leads to intersection-

free interpolations.

1 Introduction

In the ﬁeld of geometry processing, a shapespace is often

used to compute an interpolation between two isomorphic

shapes. With regard to this framework physically based

methods have been most extensively studied among the many

techniques developed. However, the methods extisting so far

have not been able to proactively prevent self-intersection.

Hence, Sassen et al. (2024) introduced a new method to

prevent this behaviour. In the lab report at hand we discuss

the theoretical foundations as well as our reformulation

and implementation of this new method. While we have

successfully reproduced the results for small examples, the

reduction of the required time per time step from

𝑂(𝑛

)

𝑂(𝑛 log 𝑛)

using an adaptive quadrature model is beyond

the scope of this report.

This lab report is structured as follows. Section 2, which

gives an overview of the literature, is followed by a brief

summary of the relevant topics of diﬀerential geometry,

namely geodesics on Riemannian manifolds. After that,

Section 4 discusses the construction of a shapespace, which

will be a special case of a Riemannian manifold. Moreover, it

shows how to interpolate two shapes by computing geodesics

on this shapespace. Section 5 gives an overview of our

implementation and Section 6 ﬁnally shows our results on a

variety of small examples.

2 Related Work

As the name suggests, shapespaces were originally used to

study the space of planar shapes. Beforehand, Robertson

(1977) utilised them for the classiﬁcation of triangles. Later,

Kendall (1984) deﬁned them more precisely. In the following

years they proved useful a variety of applications such as

the comparison (Klassen et al., 2004) and the matching

(Schmidt et al., 2007) of shapes.

These concepts were ﬁrst applied to triangle meshes

by Kilian et al. (2007), who showed how to use them to

interpolate between two shapes or extrapolate a series of

shapes. In the following, this framework was expanded by

Heeren (2017), who among other things showed that ﬁrstly

the discretization was well behaved and who secondly was

able to derive the shapespace from a physically sound model.

3 Fundamentals

This section is based on Lee (2018) and describes the

concepts underpinning the shapespace. We assume some

familarity with smooth manifolds (cf. Lee, 2012). More

speciﬁcally, we refer to chapter 3.1 and 3.2 of the previous

seminar report (Abt, 2024).

Throughout this section we make use of the Einstein

notation, which states that if an index appears exactly twice

in an expression - once as an upper index and once as a

lower index - we will assume a sum of all possible values of

that index. For example, given a basis

𝑒

, . . . 𝑒

𝑛

and a set of

coeﬃcients 𝑎

, . . . 𝑎

𝑛

, a vector 𝑣 can be expressed as

𝑣 = 𝑎

𝑖

𝑒

𝑖

𝑛



𝑖=0

𝑎

𝑖

𝑒

𝑖

. (1)

Throughout this section we assume that

𝑀

is an

𝑛−

dimensional smooth manifold and denote the space of

all real-valued functions from

𝑀

𝐶

∞

(𝑀)

. A coordinate

chart is a pair

(𝑈, 𝜑),

where

𝑈 ⊆ 𝑀

is an open set and

𝜑

𝑈 →

𝑈 ⊆ R

𝑛

is a homeomorphism. A tangent vector

𝑝

is a derivation, i.e. it is a linear map

𝑣

𝐶

∞

(𝑀) → R

satisfying

𝑣( 𝑓 · 𝑔) = 𝑓 (𝑝) · 𝑣(𝑔) + 𝑔(𝑝) · 𝑣 ( 𝑝).

The space

of all tangent vectors is an

𝑛

-dimensional vector space and

denoted by

𝑇

𝑝

𝑀

and disjoint union of all tangent spaces of

M is denoted by 𝑇 𝑀.

The coordinate map

𝜑

𝑀 → R

𝑛

can be decomposed into

the coordinate functions

(𝑥

. . . 𝑥

𝑛

)

, which together yield a

basis for the tangent space at each point 𝑝 ∈ 𝑀 by

𝜕

𝑖

𝑝

:= 𝑓 ↦→

𝜕

𝜕𝑥

𝑖



𝑝

( 𝑓 ) =

𝜕

𝜕𝑥

𝑖



𝜑 ( 𝑝)

( 𝑓 ◦ 𝜑

−1

Expanding this tangent vector to the tangent bundle yields a

coordinate vector ﬁeld on U denoted by

𝜕

𝑖

. The set of all

coordinate vector ﬁelds form a basis of the space of all smooth

vector ﬁelds, which we denote by

𝔛( 𝑀)

. The diﬀerential

𝑑𝐹

𝐹

𝑀 → 𝑁

is a map

𝑑𝐹

𝑇

𝑝

𝑀 → 𝑇

𝐹 ( 𝑝)

𝑁

, so that it

satisﬁes

(𝑑𝐹

𝑝

(𝑣)) ( 𝑓 ) = 𝑣 ( 𝑓 ◦ 𝐹)

for every

𝑓 ∈ 𝐶

∞

(𝑀)

and

𝑣 ∈ 𝑇

𝑝

𝑀

. If

𝑀 = R

𝑛

and

𝑁 = R

𝑚

, it is the Jacobian matrix

𝑑𝐹 ∈ R

𝑚×𝑛

and generally it is represented in coordinate

bases by the Jacobian matrix of the coordinate representative

of 𝐹.

3.1 Riemannian Manifolds

In order to introduce geometry to a smooth manifold

𝑀

, we

make use of a Riemannian metric, which contrary to its

name it is not a metric, but instead an inner product deﬁned

on each tangent space, i.e. a map

𝑔

𝑝

𝑇

𝑝

𝑀 × 𝑇

𝑝

𝑀 ↦→ R

which is symmetric and positive deﬁnite. Thereby, it deﬁnes

a norm of any vector

𝑣 ∈ 𝑇

𝑝

𝑀

|𝑣|

𝑔



𝑔(𝑣, 𝑣)

at each

point. A Riemannian manifold is a pair

(𝑀, 𝑔)

consisting

of a smooth manifold

𝑀

and a Riemannian metric

𝑔

, which

is smooth over 𝑀.

The Riemannian metric can be written in a matrix form

where

𝑔

𝑖 𝑗

(𝑝) = 𝑔



𝜕

𝑖

𝑝

, 𝜕

𝑗



𝑝



This means that at a point

𝑝 ∈ 𝑀 with 𝑢, 𝑣 ∈ 𝑇

𝑝

𝑀,

𝑔

𝑝

(𝑢, 𝑣)⟨𝑢, 𝑣, ⟩

𝑔

= 𝑢

𝑇

𝐺𝑣,

where

𝐺 ∈ R

𝑛×𝑛

. Similarly to the Riemannian metric,

𝐺

symmetric and positive deﬁnite.

3.2 Geodesics

We want to expand the concepts of straight lines to Rieman-

nian manifolds, which are named geodesics in the general

case (cf. Figure 2. In the Euclidean case straight lines

are of minimal length and of zero acceleration. We will

generalize the last requirement, because although minimal

curves are geodesics, geodesics in general are not necessarily

of minimal length. Instead, they are locally minimizing,

which means that for every point on the curve there will be

a neighborhood so that the restriction of the curve to it is

minimizing.

Figure 2: Geodesics on a sphere are arcs of the great circles.

3.2.1 Vector Fields along Curves

A curve is a continuous map

𝛾

𝐼 → 𝑀,

where

𝐼 ⊆ R

Given a smooth curve

𝛾

, the velocity reﬂects the change of

curve. It is a vector ﬁeld 𝛾

′

: 𝐼 → 𝑇 𝑀, which is deﬁned as

𝛾

′

: 𝐼 → 𝑇 𝑀 := 𝑡 ↦→ (𝑑𝛾

𝑡

)



𝑇

𝑡

𝐼→𝑇

𝛾 (𝑡 ) 𝑀

( 𝑑/𝑑𝑡

𝑡

)

  

𝑇

𝑡

where

𝑑𝛾

𝑡

is the diﬀerential of

𝛾

and

( 𝑑/𝑑𝑡

𝑡

)

is the coor-

dinate basis of

𝑇

𝑡

, e.g. 1. The components of

𝛾

′

are the

ordinary derivatives of the component of 𝛾, i.e.

𝛾

′

(𝑡) = ¤𝛾

𝑖

(𝑡)𝜕

𝑖



𝛾 (𝑡 )

Notice that since

𝜕

𝑖

𝛾 (𝑡 )

∈ 𝑇

𝛾 (𝑡 )

𝑀

, the derivatives make

sense in the context of a Riemannian manifold and the

components of

𝛾

′

(𝑡)

are the ordinary

𝑡

-derivatives of the

components of

𝛾

with respect to the smooth local coordinates

𝑥

. . . 𝑥

𝑛

A curve is called a unit-speed curve if its velocity

|𝛾

′

(𝑡)|

𝑔

1 for all

𝑡 ∈ 𝐼 = [𝑎, 𝑏]

. Every smooth curve

has a unit-speed parameterization

˜𝛾

˜𝛾 = 𝛾 ◦ 𝑠

−1

where

𝑠(𝑡) =



𝑡

𝑎

𝛾

′

(𝑢)

𝑔

𝑑𝑢

. Thus, we lose no generality if we

require a curve to be unit-speed parameterized.

Next, the acceleration

𝛾

′′

reﬂects the change of the velocity

and thus is again a vector ﬁeld. In the Euclidean case the

acceleration ﬁeld is once more obtained by diﬀerentiating

the components of the acceleration by

𝛾

′′

(𝑡) = ¥𝛾

𝑖

(𝑡)𝜕

𝑖



𝛾 (𝑡 )

A unit-speed parameterized curve is a straight line if and

only if

𝛾

′′

≡

0. Diﬀerentiating the components is done

¥𝛾

𝑖

(𝑡) = lim

ℎ→0

¤𝛾

𝑖

(𝑡+ℎ) − ¤𝛾

𝑖

(𝑡 )

ℎ

which works out in the

Euclidean case because

¤𝛾(𝑡) ∈ 𝑇

𝛾 (𝑡 )

𝑛

) = 𝑇

𝛾 (𝑡+ℎ)

𝑛

)

However, on a general Riemannian manifold this equality

between diﬀerent tangent spaces is not given (see Abt (2024,

ch. 3.2.1) and Figure 3), which is why we need a few more

concepts to deﬁne geodesics.

Cartesian Coordinates Polar Coordinates

𝛾 (𝑡 )

𝛾 (𝑡 + ℎ)

𝛾

′

(𝑡 + ℎ)

𝛾

′

(𝑡)

𝛾

′′

(𝑡)

𝛾

′′

(𝑡 + ℎ)

𝛾 (𝑡 )

𝛾

′

(𝑡)

𝛾 (𝑡 + ℎ)

𝛾

′

(𝑡 + ℎ)

𝑥

𝑦

𝜙

𝑟

Figure 3:

𝛾

′

(𝑡)

and

𝛾

′

(𝑡 + ℎ)

lie in diﬀerent vector spaces

depending on the chart used to deﬁne the smooth manifold

3.2.2 Connection

A connection, named after its role to connect nearby tangent

spaces, generalizes the notion of directional derivatives of

vector ﬁelds to smooth manifolds. It is a map

∇ : 𝔛(𝑀)



direction

× 𝔛(𝑀)



function

→ 𝔛 (𝑀), (2)

written 𝑋, 𝑌 ↦→ ∇

𝑋

𝑌 satisfying the following properties.

1. ∇

𝑋

𝑌

is linear over

𝐶

∞

(𝑀)

𝑋

, because it is inde-

pendent of linear combinations of vector ﬁelds in the

ﬁrst parameter and their change: For

𝑓

, 𝑓

∈ 𝐶

∞

and

𝑋

, 𝑋

∈ 𝔛(𝑀),

∇

𝑓

𝑋

+ 𝑓

𝑋

𝑌 = 𝑓

∇

𝑋

𝑌 + 𝑓

∇

𝑋

𝑌 .

In contrast, the rate of change in the second parameter

inﬂuences the connection. Thus,

∇

𝑋

𝑌

is linear over

𝑌

(and not over

𝐶

∞

(𝑀)

): For

𝑎

, 𝑎

∈ R

and

𝑌

, 𝑌

∈ 𝔛,

∇

𝑋

(𝑎

𝑌

+ 𝑎

𝑌

) = 𝑎

∇

𝑋

𝑌

+ 𝑎

∇

𝑋

𝑌

Similarly to the Leibniz rule,

∇

satiﬁes the following

product rule: for 𝑓 ∈ 𝐶

∞

(𝑀),

∇

𝑋

( 𝑓 𝑌) = 𝑓 ∇

𝑋

𝑌 + (𝑋 𝑓 )𝑌.

Generalizing a connection to scalar ﬁelds, it turns out that

∇

is given by the ordinary diﬀerentiation of the function, i.e.

for 𝑓 ∈ 𝐶

∞

(𝑀),

∇

𝑋

𝑓 = 𝑋 𝑓 .

To give a connection that will become important later on,

consider the Euclidean connection for the special case if

𝑀 = R

𝑛

, which is deﬁned by

∇

𝑋

(𝑌) = 𝑋 (𝑌

)𝜕

+ . . . 𝑋 (𝑌

𝑛

)𝜕

𝑛

. (3)

Now, in the same way that we have generalized the Eu-

clidean Riemannian metric by allowing coeﬃcients with

non-matching indices to be multiplied together, we gener-

alize the Euclidean connection to be deﬁned by

𝑛

smooth

functions for each of the

𝑛

dimensions denoted by

𝑘

𝑖 𝑗

The vector ﬁeld

∇

𝜕

𝑖

𝜕

𝑗

in a smooth local coordinate frame

(𝜕

. . . , 𝜕

𝑛

) expands to

∇

𝜕

𝑖

𝜕

𝑗

= Γ

𝑘

𝑖 𝑗

𝜕

𝑘

(4)

and a connection of 𝑋 = 𝑋

𝑖

𝜕

𝑖

, 𝑌 = 𝑌

𝑗

𝜕

𝑗

expands to

∇

𝑋

(𝑌) = 𝑋 (𝑌

𝑘

)𝜕

𝑘

+ 𝑋

𝑖

𝑌

𝑗

𝑘

𝑖 𝑗

𝜕

𝑘

Thus, there is a correspondence between the choices of

𝑛

smooth functions and a connection, which we will limit

to deﬁne a unique connection reﬂecting the directional

derivative.

3.2.3 Riemannian connection

First, we restrict the possible connections with the help of the

Riemannian metric. Similarly to the Euclidean case, which

satiﬁes

∇

𝑋

⟨𝑌, 𝑍⟩ = ⟨

∇

𝑋

𝑌, 𝑍⟩ + ⟨𝑌 ,

∇

𝑋

𝑍⟩

, a connection is

compatible with g if

∇

𝑋

𝑔(𝑌 , 𝑍) = 𝑔(∇

𝑋

𝑌, 𝑍) + 𝑔(𝑌, ∇

𝑋

𝑍).

∇

𝑋

𝑔(𝑌 , 𝑍)

𝑔(∇

𝑋

𝑌, 𝑍)

𝑔(𝑌 , ∇

𝑋

𝑍)

Figure 4: For example, the Euclidean connection is a metric

connection, as visualized here for

𝑋 = (

)

(green),

𝑌 = (𝑥

, 𝑦)

(blue) and

𝑍 = (𝑥, 𝑦

)

(orange). In the right two

images, the entries of the function are visualized. In general,

the left size just means

𝑋 (𝑔(𝑌 , 𝑍))

viewing

𝑋

𝑝

𝐶

∞

→ R

for all

𝑝 ∈ 𝑀

. This ensures that if g is a geodesic,

|𝛾

′

(𝑡)|

𝑔

will be constant over all 𝑡.

However, this still allows for an inﬁnite amount of metrics.

Thus, we further restrict the connection to be symmetric,

i.e. to satisfy

𝑘

𝑖 𝑗

= Γ

𝑘

𝑗𝑖

which stems from the fact that

this equality holds for every connection of an embedded

manifold projected from R

𝑛

The fundamental theorem of Riemannian geometry states

that these two restrictions determine a unique connection

named the Riemannian connection. The connection co-

eﬃcients of Equation (4) in any smooth coordinate chart

𝜕

𝑖

. . . 𝜕

𝑛

are the so-called Christoﬀel symbols given by

𝑘

𝑖 𝑗

𝑔

𝑘𝑙

(𝜕

𝑖

𝑔

𝑗𝑙

+ 𝜕

𝑗

𝑔

𝑖 𝑗

+ 𝜕

𝑙

𝑔

𝑖 𝑗

3.2.4 Geodesic Equation

We are now able to deﬁne the Covariant Derivative Along

a Curve, which is the derivative of a vector ﬁeld along it.

For each curve, the Levi Civita connection determines a

unique operator

𝐷

𝑡

𝔛(𝛾) → 𝔛 (𝛾)

. In local coordinates

𝑉 ∈ 𝔛

is written as

𝑉 (𝑡) = 𝑉

𝑗

(𝑡) 𝜕

𝑗

𝛾 (𝑡 )

and

𝐷

𝑡

is deﬁned

by:

(𝐷

𝑡

𝑉)(𝑡) =



𝑉

𝑘

(𝑡) + ¤𝛾

𝑖

(𝑡)𝑉

𝑗

(𝑡)Γ

𝑘

𝑖 𝑗

(

𝛾(𝑡)

)



𝜕

𝑘

𝛾 (𝑡 )

. (5)

Finally Geodesics are now the curves with 0 acceleration,

i.e. 𝛾 is a geodesic if

𝐷

𝑡

𝛾

′

≡ 0. (6)

From Equation (5) and writing out the curve in smooth

coordinates coeﬃcients

𝛾(𝑡) = (𝑥

(𝑡), . . . , 𝑥

𝑛

(𝑡))

, we arrive

at the following geodesic equation:

¥𝑥

𝑘

(𝑡) + ¤𝑥

𝑖

(𝑡) ¤𝑥

𝑗

(𝑡)Γ

𝑘

𝑖 𝑗

(𝑥(𝑡)) = 0 (7)

3.3 Critical Points of 𝐿

𝑔

are Geodesics

Given a unit-speed curve

𝛾

𝐼 → 𝑀

, which maps from an

interval

𝐼 = [𝑎, 𝑏] ⊂ R

to a Riemannian manifold

𝑀

, let

𝐿

𝑔

(𝛾) =



𝑏

𝑎

|𝛾

′

(𝑡)|

𝑔

𝛾 (𝑡 )

𝑑𝑡

denote the length of this curve.

This section shows that the critical points of

𝐿

𝑔

which are

the curves where the derivative of

𝐿

𝑔

is zero are geodesics.

Since these are the curves which will be computed later on,

we go into more detail here and sketch the proof.

Families of Curves. To discuss this topic we must intro-

duce some terminology. Firstly, a family of curves is a

smooth map Γ : 𝐽 × 𝐼 → 𝑀. 2 Keeping the ﬁrst parameter

ﬁxed, the main curves are deﬁned by

𝑠

(𝑡) = Γ(𝑠, 𝑡)

for

𝑠 ∈ 𝐽

and keeping the seond parameter ﬁxed, the tranverse

curves are deﬁned by Γ

𝑡

(𝑠) = Γ(𝑠, 𝑡) for 𝑡 ∈ 𝐼.

Secondly, a vector ﬁeld along

is a continuous map

𝑉

𝐽 × 𝐼 → 𝑇 𝑀

so that

𝑉 (𝑠, 𝑡) ∈ 𝑇

Γ (𝑠,𝑡 )

𝑀

. Two

prominent examples are the velocity vectors of the main

Technically, throughout this section we do not require unit-speed curves,

but only piecewise continuous curves with

𝛾

′

(𝑡 ) ≠

0. We drop these terms

for curves and related concept by the reasons explained in Section 3.2.1.

A family of curves is actually allowed to be a continuous map, but this

makes the proof much harder, which is the reason why it is left out here.

and tranverse curves deﬁned by

𝜕

𝑡

Γ(𝑠, 𝑡) = (Γ

′

𝑠

)(𝑡)

and

𝜕

𝑠

Γ(𝑠, 𝑡) = (Γ

′

𝑡

)(𝑠) respectively.

Thirdly, given a unit-speed curve

𝛾

[𝑎, 𝑏] → 𝑀,

proper variation of

𝛾

is a family of curves

𝐽 × [𝑎, 𝑏] →

𝑀

so that

= 𝛾

as well as

𝑠

(𝑎) = 𝛾(𝑎)

and

𝑠

(𝑏) = 𝛾(𝑏)

for all

𝑠 ∈ 𝐽

. The special velocity ﬁeld of the transverse

curves at

𝛾

is called the variation ﬁeld of

deﬁned by

𝑉 (𝑡) = 𝜕

𝑠

Γ(

, 𝑡)

. It will be called proper if

𝑉 (𝑎) =

0 and

𝑉 (𝑏).

Fourthly, an important result is the fact that if

𝑉

is any

smooth vector ﬁeld along

𝛾

𝑉

will be the variation ﬁeld

of some variation of

𝛾

, which can be shown by smoothly

extending 𝑉 on the manifold.

Figure 5: every

proper vector ﬁeld

along

𝛾

is a proper

variation ﬁeld (Lee,

2018)

Figure 6: Because

𝑑

𝑑𝑠

𝐿

𝑔

for a

curve

𝛾

with

𝐷

𝑡

𝛾

′

≠

0 is neg-

ative, deforming

𝛾

towards its

acceleration decreases

𝐿

𝑔

(Lee,

2018).

Sketch of Proof. Finally, we want to show that the critical

points of

𝐿

𝑔

, i.e. the curves at which

𝑑

𝑑𝑠



𝑠=0

𝐿

𝑔

(Γ

𝑠

) =

are

geodesics for any proper variation Γ of 𝛾. Calculus tells us

that

𝑑

𝑑𝑠

𝐿

𝑔

(Γ

𝑠

) =



𝑏

𝑎

𝜕

𝜕𝑆



⟨𝜕

𝑡

Γ, 𝜕

𝑡

Γ⟩

𝑔

𝑑𝑡



𝑏

𝑎

⟨𝜕

𝑡

Γ, 𝜕

𝑡

Γ⟩

−1/2

𝑔

· 2⟨𝐷

𝑠

(𝜕

𝑡

Γ), 𝜕

𝑡

Γ⟩

𝑔

𝑑𝑡.

Because of the symmetric property of the Riemannian con-

nection described in Section 3.2.3, it is possible to show

that

𝐷

𝑠

(𝜕

𝑡

Γ) = 𝐷

𝑡

(𝜕

𝑠

Γ),

which allows us to rewrite the

integrand to

𝑑

𝑑𝑠

𝐿

𝑔

(Γ

𝑠

) =



𝑏

𝑎

|𝜕

𝑡

Γ|

⟨𝐷

𝑡

(𝜕

𝑠

Γ), 𝜕

𝑡

Γ⟩

𝑔

𝑑𝑡.

Finally, setting

𝑠

to zero and noting that

𝛾

is unit-speed

parameterized, the last equation simpliﬁes to

𝑑

𝑑𝑠



𝑠=0

𝐿

𝑔

(Γ

𝑠

) = −



𝑏

𝑎

⟨𝜕

𝑠

(𝑡), 𝐷

𝑡

𝛾

′

⟩

𝑔

𝑑𝑡.

Now consider

𝜕

𝑠

(𝑡)

as a smooth vector ﬁeld

𝑉

along

𝛾

Since the last equation must hold for all proper variations,

𝑉

can be set to an arbitrary smooth proper vector ﬁeld. Thus,

we can set 𝑉 = 𝐷

𝑡

𝛾

′

to obtain

0 = −



𝑏

𝑎

𝐷

𝑡

𝛾

′

𝑔

𝑑𝑡. (8)

Since the integrand is non-negative,

𝐷

𝑡

𝛾

′

0 must hold for

all 𝑡 ∈ [𝑎, 𝑏].

4 Methodology

Next, we discuss the construction of a shapespace which will

allow us to compute "good" interpolations between shapes.

4.1 Shape Space

In the continuous case a shapespace is the space of all

immersions of a smooth manifold

𝑁

into

𝑑

of a given

smooth manifold, i.e.

𝑀

𝐶

∞

(𝑁, 𝑂 ↩→ R

𝑑

)

. A classical

result from global analysis is that this space is itself a smooth

manifold. 3

Any open curve between two embeddings is a smooth

interpolation between them. Next, a Riemannian metric is

constructed in a way that the geodesics are "good" curves, i.e.

that the corresponding interpolations are visually pleasing,

physically correct or without any self-intersections.

As to the discrete case, we approximate the smooth man-

ifold by a simplicial complex and an embedding of it by a

mapping of all vertices into

𝑛

, i.e. the shapespace of a

simplicial complex is

𝑀

= R

𝑛

where

𝑛

is the multiplication

of the number of vertices by three.

In order to construct the Riemannian metric at any point

𝑝 ∈ 𝑀

, we divide it into a sum of

𝑚

Riemannian metrics,

which each encode a property which we wish to see fullﬁlled

by the geodesics. The Riemannian metric

ˆ𝑔

is thereby

deﬁned by

ˆ𝑔(𝑢, 𝑣) = 𝑢

𝑇

𝐺𝑣 = 𝑢

𝑇



𝑚



𝑖=0

𝐺

𝑖



𝑣, (9)

where

𝑢, 𝑣 ∈ 𝑇

𝑝

𝑀 = 𝑇

𝑝

𝑛

= R

𝑛

and

𝐺

𝑖

∈ R

𝑛×𝑛

for all

metrics.

We postpone the description of the implemented metrics

until Section 4.4 and ﬁrst describe how they are used to

compute geodesics.

4.2 Path energy

Recall that in Section 3.3 critical points of

𝐿

𝑔

have proved

to be geodesics. Thus, we can ﬁnd a geodesic by ﬁnding the

critical point the path energy deﬁned by

E(𝛾) :=



𝑏

𝑎

∥𝛾

′

(𝑡)∥

𝑔

𝛾 (𝑡 )

𝑑𝑡, (10)

where we can integrate over the square norm for eﬃciency,

because quadratic optimization yields the same minima.

As to the discrete case we approximate the curve

ˆ𝛾

a piecewise linear curve of

𝑘

segments, in which each

Because

𝑀

is of inﬁnite dimensions its construction is Frechet manifold,

i.e. it is modeled on the Frechet space instead on

𝑛

. For a proof that the

M is indeed smooth see Kriegl and Michor (1997, sec. 42.4)

vertex reﬂects a shape

𝑠

𝑖

∈ R

𝑛

. The tangent vector at each

vertex

𝑖

towards the next vertex can thus be computed by

𝑡

𝑖

𝑛

→ R

𝑛

= 𝑠

𝑖

↦→ 𝑠

𝑖+1

− 𝑠

𝑖

and the Riemannian metric

can be computed by

𝑔

𝑖

𝑛

→ R

= 𝑣 ↦→ 𝑔

𝑠

𝑖

(𝑣, 𝑣),

where

𝑔

𝑠

𝑖

is given by Equation (9). Thus, the discrete path energy

E =

𝑘



𝑖=1

(𝑔

𝑖

◦ 𝑡

𝑖

)(𝑆

𝑖

). (11)

Equation (11) is optimized using a quasi-Newton method,

which requires computing the derivative

∇

which is given

∇

E( ˆ𝛾) =

𝑘



𝑖=1

∇

𝑡

𝑖

(𝑠

𝑖

)

(𝑔

𝑖

) ◦ ∇

𝑠

𝑖

(𝑡

𝑖

) (12)

𝑘



𝑖=1

−2𝐺

𝑖

𝑠

𝑖

+ 2𝐺

𝑖

𝑠

𝑖+1

. (13)

However, because the start and end shape are not allowed to

change, we ﬁx the gradient for 𝑠

and 𝑠

𝑘

at 0.

4.3 Reﬁnement

In order to improve eﬃciency we ﬁrst start with a coarse

discretization of the geodesic and ﬁnd the local minima of

Next, we allow for reﬁnement of the curve by inserting new

vertices of the curves, i.e. new shapes, which can then again

be optimized to minimize

. This reﬁnement is known as

time reﬁnement and can be repeated as often as required.

Kilian et al. (2007) reﬁne the geodesic with respect to time.

They insert a new vertex by lineary blending the existing

vertices of the curve. Instead, we construct a cubic spline

for each vertex of the mesh and resample this spline. In our

tests, the resulting curves proved to be closer to the ﬁnal

geodesic.

Similarly, it is possible to ﬁrst simplify the mesh, optimize

this coarse representation and gradually reﬁne it back until

to its original shape. However, this is out of scope for this

lab.

Even though both of these approaches introduce additional

steps, the overall approach proved to be much faster, because

the optimization on coarse mesh is much more performant

while subsequent optimization only needs to ﬁnetune the

geodesic.

4.4 Metrics

This section gives an overview of the diﬀerent metrics while

their implementations can be found in Section 5. Let

𝑀

a (discrete) shapespace. Throughout this section we will

denote a point in the shapespace

𝑀

𝑠 ∈ 𝑀

, its set of

triangles of a shape by

𝑇

𝑠

, its set of edges by

𝐸

𝑠

and its set

of vertices by

𝑉

𝑠

. The position of a vertex is denoted by

𝑠

𝑣

∈ R

for

𝑣 ∈ 𝑉

𝑠

, the tangent vector by

𝑥, 𝑦 ∈ 𝑇

𝑠

𝑀

and

the change of the position of a vertex by

𝑥

𝑣

∈ R

for

𝑣 ∈ 𝑉

𝑠

4.4.1 𝐿

Metric

Most other metrics will turn out to be only semi-Riemannian

metrics, i.e. they will be allowed to be zero for non-zero

tangent vectors. Thus, a regularization term is added, which

is given by

𝑔

𝐿

(𝑥, 𝑦) =



𝑣∈𝑉

𝑠

⟨𝑥

𝑣

, 𝑦

𝑣

⟩ · 𝐴

𝑣

where

𝐴

𝑣

is one-third of the area of the surrounding triangles

of 𝑣 (Kilian et al., 2007).

4.4.2 Isometric Metric

A deformation of a mesh is called isometric if distances on

the mesh are preserved. In order to interpolate shapes as

isometrically as possible we employ the isometric metric,

which increases the geodesic length depending on the defor-

mations of the the edges (Kilian et al., 2007). At each shape

𝑠 ∈ 𝑀, the isometric metric is deﬁned as

𝑔

𝐼

(𝑥, 𝑦) =



(𝑢,𝑣) ∈𝐸

𝑠

(⟨𝑥

𝑢

− 𝑥

𝑣

, 𝑠

𝑢

− 𝑠

𝑣

⟩ · ⟨𝑦

𝑢

− 𝑦

𝑣

, 𝑠

𝑢

− 𝑠

𝑣

⟩).

(14)

This expression can be reformed in terms of each component

𝑢, 𝑣

(cf. Section 5), which can then be used to construct

4.4.3 Membrane Metric

The membrane metric

𝑔

𝑚

originates from the membrane

energy used in continuum mechanics. Together with the

bending energy (cf. Section 4.4.4), it models the energy

required to deform a thin deformable material called a shell.

We denote the original shell by

𝑠

and the deformed shell and

derived concepts by

˜𝑠

. The membrane energy accounts for

the tangential stretching of the surface and, in the continuous

case, is the smooth function

𝒲

𝑚

: 𝑀 × 𝑀 → R := 𝑠, ˜𝑠 ↦→



𝑠

𝑊 (I

−1

𝑠



𝑠

) 𝑑𝐴, (15)

where

denotes the ﬁrst fundamental form in its matrix form

and W describes the response of the material by

𝑊 :R

2×2

→ R := 𝐴 ↦→

𝜇

tr 𝐴 +

𝜆

tr 𝐴 − (

𝜇

𝜆

) log(det 𝐴) − 𝜇 −

𝜆

As to the discrete case we approximate the integral by a sum

over all triangles to get



𝒲

𝑚

:= 𝑠, ˜𝑠 ↦→



𝑡 ∈𝑇

𝑠

𝑎

𝑡

𝑊 (I

−1

𝑡



𝑡

where I

𝑡

= 𝐽

𝑇

𝑡

𝐽

𝑡

can be computed using

𝐽

𝑡



𝑥

𝑗

− 𝑥

𝑖

𝑥

𝑘

− 𝑥

𝑖



∈ R

3×2

Next, Heeren (2017, p. 75) use the Hessian as the Rieman-

nian metric, that is they use the function

∇

𝑋

∇

𝑌

𝑢

. However,

this requires itself a connection and thus if we use the Rie-

mannian connection a Riemannian metric, but luckily it

turns out that at a critical point of function the Hessian is

independent of the choice of metric. The properties of the

resulting Riemannian metric can be easily veriﬁed to deﬁne

a riemmanian metric by

𝑔

𝑚

= 𝑑

˜𝑠

𝒲

𝑚

(𝑠, ˜𝑠)



𝑠= ˜𝑠

4.4.4 Bending Metric

Similarly, the bending metric is deﬁned by

𝑔

𝑏

= 𝑑

˜𝑠

𝒲

𝑏

(𝑠, ˜𝑠)



𝑠= ˜𝑠

where

𝒲

𝑏

denotes the bending energy. The latter energy

models the energy required to bend a material in its normal

direction and in the continuous case is deﬁned by

𝒲

𝑏

: 𝑀 × 𝑀 → R := 𝑠, ˜𝑠 ↦→ 𝛿



𝑠

∥



II − II∥

𝐹

𝑑𝐴,

where II denotes the second fundamental form.

As to the discrete case we approximate the integral by the

discrete bending area, which is deﬁned by



𝒲

𝑏

:= 𝑠, ˜𝑠 ↦→ 𝛿



𝑒∈𝐸

𝑠

𝑒 is interior



2tan





𝜃

𝑒



− 2tan

(

𝜃

𝑒

)



𝑙

𝑒

/𝑎

𝑒

where

𝜃

𝑒

denotes the angle between the triangle normals

and

𝑎

𝑒

= (𝑎

𝑡

+ 𝑎

𝑡

3 denotes the edge area depending on

the two neighbouring triangles.

4.4.5 Repulsive Metric

In order to disallow self-intersections, we make use of a

repulsive potential which is based on the tangent point

energy. In the continuous case the latter one is deﬁned

𝜙(Σ) =



𝐾 (𝑥, 𝑦) 𝑑𝑥 𝑑𝑦,

where

𝐾 (𝑥, 𝑦)

reﬂects the

inverse of the radius of the smallest sphere through

𝑥

and

𝑦

tangent to Σ at 𝑥. It is deﬁned as

𝐾 (𝑥, 𝑦) :=

⟨𝑛

𝑥

, 𝑥 − 𝑦⟩

𝛼

𝑥 − 𝑦

, (16)

where

𝑛

𝑥

denotes the unit normal at

𝑥

. In the discrete case,

we compute 𝐾 for each pair of triangles, i.e.



𝜙(x) :=



𝑡

∈𝑇



𝑡

∈𝑇

𝑎

𝑡

𝑎

𝑡

𝐾 (𝑐

𝑡

, 𝑐

𝑡

, 𝑛

𝑡

where 𝑎

, 𝑎

denote the areas of the triangles.

The repulsive metric

𝑔

reﬂects the change of repulsive

potential

𝑔

(𝑢, 𝑣) = ∇

𝑢

𝜙 · ∇

𝑣

𝜙. (17)

In case that the shapespace is

𝑛

, the connection

∇

is the

directional derivative and can be expressed as

∇

𝑢

= ⟨𝑑𝜙, 𝑢⟩,

where 𝑑𝜙 is the gradient. Thus, Equation (17) becomes

𝑔

(𝑢, 𝑣) = ⟨𝑑𝜙, 𝑢⟩ · ⟨𝑑𝜙, 𝑢⟩ = ⟨𝑢, 𝑑𝜙⟩ · ⟨𝑑𝜙, 𝑢⟩.

As to the discrete case,

ˆ𝑔

is due to the symmetry of the inner

product

ˆ𝑔

(𝑢, 𝑣) = (𝑢

𝑇

𝑑𝜙)(𝑑𝜙

𝑇

𝑣) = 𝑢

𝑇

𝐺

𝑣, (18)

where the associated matrix 𝐺

= ((𝑑𝜙)(𝑑𝜙)

𝑇

The repulsive metric is special in the sense that it is

everywhere non-zero except the the diagonal. Thus, unlike

the previous metrics it is not eﬃcienctly computed using

a sparse matrix. In addition, its construction is in

𝑂(𝑛

)

That is why Sassen et al. (2024) employ a hierarchical fast

multipole method, which only computes the repulsive metric

and its gradient where needed. In addition, their method

subdivides close triangles to avoid any irregularities by the

discretization stemming from the fact the discretization only

takes the centers of the triangles into account. As this is

incompatible with our framework, it is out of scope of this

lab.

5 Implementation

The implementation named Khaki is written in c++ and

divided into 3 major parts:

The application Khaki-app implements a GUI, which

allows conﬁguring the shapespace and visualizing the inter-

polations.

The library Khaki-lib mainly implements the shapespace

as well as all metrics discussed in Section 4.4. In addition,

it implements the framework of the exercises in the class

that this lab is based on. This framework provides the mesh

class based on OpenMesh (Bischoﬀ et al., 2002) and the

visualization based on the Hazel game engine (Chernikov,

2018). Moreover, Khaki-lib is being tested against unit tests

implemented in Khaki-test.

The optimization is done using Qiu (2023) with Equa-

tion (11) and Equation (13) being computed in parallel using

a threadpool. Every Riemannian metric described in Sec-

tion 4.4 must return the matrix form of

𝑔

provided a shape

𝑠.

5.1 Computation of Riemannian Metrics

Next we describe how we have derived the methods to

compute the matrix forms of the supported metrics. The

implementation of 𝑔

𝐿

is straight forward (Algorithm 1).

Algorithm 1 Computation of Matrix of 𝑔

𝐿

𝐺 ← 0

𝑛×𝑛

for 𝑣 ∈ 𝑉

𝑠

𝑎 ← 0

for 𝑓 in neighbouring faces of 𝑣 do

𝑎 ← 𝑎 + area( 𝑓 )/3

end for

for 𝑑 = 0 to 3 do

𝐺 [idx(𝑣, 𝑑), idx(𝑣, 𝑑)] ← 𝑎

end for

return 𝐺

Similarly, we compute the matrix form of the isometric

metric

𝑔

𝐼

by adding all summands. However, the transfor-

mation of Equation (14) into the form

𝑥

𝑇

𝐺𝑦

for

𝑥, 𝑦 ∈ 𝑇

𝑠

𝑀

is quite convoluted. Therefore, we have used the computer

algebra program Sympy (Meurer et al., 2017) to rewrite the

summands in terms of one sample edge, which are then used

to set the entries of the matrix.

We have tried to use the same approach to compute the

entries of the repulsive metric

𝑔

𝑟

. However, its expression

expanded so largely that taking the derivative was infeasible.

Thus, we switched to GiNaC (Bauer et al., 2002), which

proved to be more performant and ergonomic, because it

supports operator overloading and seemingly integrated into

our codebase (e.g. listing 1). In order to compute the

entries of

𝑔

𝑟

we make use of the functionality of GiNaC to

compute the derivatives with respect to the vertex positions.

Afterwards, these expressions are converted to c functions,

which are linked to our program.

Listing 1: constructing expression of Equation (16)

GiNaC : : ex K(

c o n s t GiNaC : : vec3& x ,

c o n s t GiNaC : : vec3& y ,

c o n s t GiNaC : : vec3& n ,

c o n s t GiNaC : : r e a l s y m b o l a l p h a

){

r e t u r n GiNaC : : pow (

GiNaC : : ab s ( n . d o t ( x−y ) ) ,

a l p ha

) / GiNaC : : pow (

( x−y ) . l2norm ( ) ,

2∗ a lpha

) ;

}

In order to further increase the performance of the deriva-

tives it will be more eﬃcient to optimize the expression

before and after computing the derivatives. This is achieved

by Enzyme (Moses and Churavy, 2020), which provides a

plugin for Clang. We saw a speedup of about 10x compared

with the derivatives generated by GiNaC. However, Enzyme

requires relatively simple functions and cannot handle func-

tions involving vectors. Thus, we use a combined approach

by ﬁrst collecting the expression using GiNaC and writing it

to a c function before using Enzyme to generate the gradient

and Hessian.

6 Evaluation

We evaluated our framework on a variety of small examples.

All tests were performed on a workstation with an AMD

Ryzen 7 2700X eight-core processor and 32 GB of RAM.

For all test cases, the start and end positions were given and

ﬁxed. However, it should be noted that the parameterization

of the diﬀerent metrics was very important as otherwise

certain metrics would dominate all the other ones.

The runtime depends on the metrics used as well as on the

model. Except for the repulsive metric, our implementation

computed geodesics of small examples almost instanta-

neously and of large examples in a few minutes. Since we

did not include the hierarchical fast multipole method, one

cannot expect good performance for large meshes. The

largest mesh was the hand in Figure 12 consisting of 600

triangles, which was computed in about 20 minutes.

As described by Sassen et al. (2024), we often need

a center surface so that the initial guess is without any

self-intersections. Otherwise, the time reﬁnement leads to

self-intersections of the models.

Figure 7: Geodesic computed only using the

𝐿

-metric

results in to linear interpolations of the start and endshape

Figure 8: Geodesic computed only using 0

001

· 𝑔

𝐿

, 1

· 𝑔

𝐼

stretches the geodesic horizontally

Figure 9: Geodesic ﬂipping the position of two quads,

computed using 0.001𝑔

𝐿

, 1𝑔

𝐼

and 1𝑔

𝑟

with 𝛼 = 2

Figure 10: Geodesic ﬂipping the position of two grids,

computed using 0.001𝑔

𝐿

, 1𝑔

𝐼

and 1𝑔

𝑟

with 𝛼 = 2

Figure 11: Geodesic computed using 0

001

𝑔

𝐿

, 1

𝑔

𝐼

𝑔

𝑚

and 0.1𝑔

𝑟

, with 𝛼 = 2

Figure 12: Geodesic computed using 0

001

𝑔

𝐿

, 1

𝑔

𝐼

and

1.0𝑔

𝑟

, with 𝛼 = 2

Figure 13: Geodesic computed using 0

001

𝑔

𝐿

, 1

𝑔

𝐼

, with

𝛼 = 2

7 Future Work

So far, all discretizations of a shapespace have only deﬁned a

shape through the embeddings of all vertices in

. However,

we theoretically only require smooth maps between all

embeddings along a curve. This could prove useful by

allowing the shape to be only dynamically reﬁned when

needed, similarly to Section 4.4.5. However, this requires a

diﬀerent discretization of the shapespace.

Instead of solving the integrals numerically, we propose

solving them by using Monte Carlo Sampling. For example,

in order to compute the repulsive metric with respect to one

point

𝑥 ∈ R

, one would sample

𝑛

random other points

𝑦

𝑖

which are sampled according to a probability distribution

function 𝑝, and estimate the value by

𝑛



𝑛

𝑖=1

𝐾 ( 𝑥, 𝑦,𝑛

𝑥

)

𝑝 ( 𝑦

𝑖

)

The hierarchical methods in Section 4.4.5 would need to

be represented by

𝑝

to require as few samples as possible.

However, we believe that this approach would allow for

further optimization as all techniques developed for Monte

Carlo sampling could be applied.

8 Conclusion

We haved discussed one possible implementation of shapes-

paces and the repulsive metric introduced by Sassen et al.

(2024) along with the necessary foundations of Riemannian

manifolds.

While we were not able to reproduce the performance

of their implementation, our framework has proved to be

very general and allows computing geodesics for any given

metrics. This was done by a rigorous approach to the

shapespace from a mathematical background to derive the

equations for the derivative of the discrete path energy.

In addition, we have computed the Riemannian metrics

by employing computer algebra systems and automatic

diﬀerentiation, which is an approach we have not seen so

far.

References

Abt, B. (2024). “Eﬃcient Computation of Minimal Surfaces

Using Exterior Calculus, Hodge Theory and Geometric

Measure Theory”. In: url:

https://spooky.moe/d00.

html.

Bauer, C. W., A. Frink, and R. B. Kreckel (2002). “Introduc-

tion to the GiNaC Framework for Symbolic Computation

within the C++ Programming Language”. In: J. Symb.

Comput. 33.1, pp. 1–12. arXiv: cs/0004015 [cs.SC].

Bischoﬀ, B. S., M. Botsch, S. Steinberg, S. Bischoﬀ, L.

Kobbelt, and R. Aachen (2002). “OpenMesh–a generic

and eﬃcient polygon mesh data structure”. In: In openSG

symposium. Vol. 18.

Chernikov, Y. (2018). Hazel.

https:/ /github.com /

TheCherno/Hazel.

Heeren, B. (2017). “Numerical methods in shape spaces and

optimal branching patterns”. PhD thesis. Universitäts-und

Landesbibliothek Bonn.

Kendall, D. G. (1984). “Shape manifolds, procrustean

metrics, and complex projective spaces”. In: Bulletin of

the London mathematical society 16.2, pp. 81–121.

Kilian, M., N. J. Mitra, and H. Pottmann (2007). “Geometric

modeling in shape space”. In: ACM SIGGRAPH 2007

papers, 64–es.

Klassen, E., A. Srivastava, M. Mio, and S. H. Joshi (2004).

“Analysis of planar shapes using geodesic paths on shape

spaces”. In: IEEE transactions on pattern analysis and

machine intelligence 26.3, pp. 372–383.

Kriegl, A. and P. W. Michor (1997). The convenient setting

of global analysis. Vol. 53. American Mathematical Soc.

Lee, J. M. (2012). Smooth manifolds. Springer.

–

(2018). Introduction to Riemannian manifolds. Vol. 2.

Springer.

Meurer, A., C. P. Smith, M. Paprocki, O. Čertík, S. B.

Kirpichev, M. Rocklin, A. Kumar, S. Ivanov, J. K. Moore,

S. Singh, T. Rathnayake, S. Vig, B. E. Granger, R. P.

Muller, F. Bonazzi, H. Gupta, S. Vats, F. Johansson, F.

Pedregosa, M. J. Curry, A. R. Terrel, Š. Roučka, A. Saboo,

I. Fernando, S. Kulal, R. Cimrman, and A. Scopatz (Jan.

2017). “SymPy: symbolic computing in Python”. In:

PeerJ Computer Science 3, e103. url:

https://doi.

org/10.7717/peerj-cs.103.

Moses, W. and V. Churavy (2020). “Instead of Rewrit-

ing Foreign Code for Machine Learning, Automatically

Synthesize Fast Gradients”. In: Advances in Neural In-

formation Processing Systems. Ed. by H. Larochelle, M.

Ranzato, R. Hadsell, M. F. Balcan, and H. Lin. Vol. 33.

Curran Associates, Inc., pp. 12472–12485. url:

https:

//proceedings.neurips .cc/ paper/2020/file/

9332c513ef44b682e9347822c2e457ac-Paper.pdf.

Qiu, Y. (2023). LBFGS++.

https : / / github . com /

yixuan/LBFGSpp/.

Robertson, S. (1977). “Classifying triangles and quadri-

laterals”. In: The Mathematical Gazette 61.415, pp. 38–

49.

Sassen, J., H. Schumacher, M. Rumpf, and K. Crane (2024).

“Repulsive shells”. In: ACM Transactions on Graphics

43.4, pp. 1–22.

Schmidt, F. R., D. Farin, and D. Cremers (2007). “Fast

matching of planar shapes in sub-cubic runtime”. In:

2007 IEEE 11th International Conference on Computer

Vision. IEEE, pp. 1–6.