# Differentiable Geodesic Distance for Intrinsic Minimization on Triangle
Meshes

Yue Li  ETH ZürichSwitzerland<EMAIL_ADDRESS>,  Logan Numerow  ETH
ZürichSwitzerland<EMAIL_ADDRESS>,  Bernhard Thomaszewski  ETH
ZürichSwitzerland<EMAIL_ADDRESS>and  Stelian Coros  ETH ZürichSwitzerland
<EMAIL_ADDRESS>
###### Abstract.

Computing intrinsic distances on discrete surfaces is at the heart of many
minimization problems in geometry processing and beyond. Solving these
problems is extremely challenging as it demands the computation of on-surface
distances along with their derivatives. We present a novel approach for
intrinsic minimization of distance-based objectives defined on triangle
meshes. Using a variational formulation of shortest-path geodesics, we compute
first and second-order distance derivatives based on the implicit function
theorem, thus opening the door to efficient Newton-type minimization solvers.
We demonstrate our differentiable geodesic distance framework on a wide range
of examples, including geodesic networks and membranes on surfaces of
arbitrary genus, two-way coupling between hosting surface and embedded system,
differentiable geodesic Voronoi diagrams, and efficient computation of Karcher
means on complex shapes. Our analysis shows that second-order descent methods
based on our differentiable geodesics outperform existing first-order and
quasi-Newton methods by large margins.

Embedded Elasticity, Intrinsic Minimization, Differentiable Simulation,
Geodesics, Differentiable Voronoi Diagram, Karcher Means

††submissionid: 133††journal: TOG††ccs: Computing methodologies Shape
modeling††ccs: Computing methodologies Physical simulation††copyright:
acmlicensed††journal: TOG††journalyear: 2024††journalvolume:
43††journalnumber: 4††article: 91††publicationmonth: 7††doi: 10.1145/3658122

##  1\. Introduction

Computing intrinsic distances on triangle meshes is a fundamental task in
geometry processing. From biological films on surfaces, to the fascia
enveloping our muscles, and to tight-fitting clothing—there are countless
examples of variational problems where the task is to minimize lengths on
discrete manifolds. Solving such problems on triangle meshes requires the
computation of geodesic distances on surfaces along with their derivatives.

Figure 1. Embedded two-way coupling. We simulate a geodesic elastic network
embedded in the surface of a deformable bunny, modeled with solid finite
elements. Our analytical geodesic distance derivatives allow us to use
Newton’s method for simulating the deformations induced by tightening the
network.

While gradients are readily computed, the convergence of first-order methods
is generally poor, in particularly for embedded elasticity problems with stiff
connections. Quasi-Newton methods such as L-BFGS may offer some acceleration,
but as we show in our analysis, they still converge slowly, especially for
stiff problems.

In this work, we present a novel approach for intrinsic minimization of
arbitrary objectives based on geodesic lengths on piece-wise linear surfaces.
Using a variational form of shortest-path geodesics, we obtain distance
derivatives in closed form using the implicit function theorem. We show that
first and second derivatives can be simplified to enable efficient
computation, allowing us to leverage powerful Newton-type methods for
simulation. Our approach is underpinned by the key observation that, although
geodesic paths on triangle meshes are generally not smooth functions of their
endpoints, geodesic distance is at least $C^{0}$-continuous everywhere—and
infinitely smooth if geodesics are unique. Our analysis reveals two types of
gradient discontinuities, one of which can be resolved through mollification
while the other one does not occur close to minima and can therefore be
ignored.

We demonstrate our differentiable geodesic distance framework on a large and
diverse set of minimization problems. Specifically, we construct embedded
elastic networks and membranes based on geodesic springs and triangle finite
elements, respectively. We present simulation examples on shapes with varying
topologies, including two-way coupling with the hosting surface. We
furthermore show how our formulation enables differentiable on-surface Voronoi
diagrams, whose site locations can be optimized for various objectives.
Finally, we demonstrate that our simulation framework can be used to compute
Karcher means on arbitrary surfaces in efficient and highly accurate ways.

##  2\. Related Work

#### Discrete Geodesic Distance

Geodesic distance has been extensively studied in the geometry processing
community (Crane et al., 2020; Peyré et al., 2010; Bose et al., 2011). Many
methods have been developed for calculating exact geodesic distances (Mitchell
et al., 1987; Chen and Han, 1990; Qin et al., 2016; Sharp and Crane, 2020b;
Liu et al., 2017a) as well as fast approximations (Crane et al., 2013; Belyaev
and Fayolle, 2015; Surazhsky et al., 2005; Edelstein et al., 2023; Belyaev and
Fayolle, 2020; Kimmel and Sethian, 1998; Ying et al., 2013; Trettner et al.,
2021; Pang et al., 2023). While approximations through the solution of smooth
energy functions offer advantages in terms of differentiability, accurately
tracing the geodesic path from the distance field is a challenging task in
itself (Crane et al., 2020). We use an exact geodesic distance computation
from which geodesic paths can be constructed. Although the numerical
properties of discrete geodesics have been studied intensively, we are not
aware of any second-order method for minimizing geodesics-based objectives. To
construct such a method, we propose a differentiable geodesic distance
formulation that is almost everywhere $C^{2}$-continuous and exhibits robust
convergence when integrated into Newton-type minimization algorithms.

#### Intrinsic Geometry Processing

The computer graphics community has made great strides in intrinsic geometry
processing (Liu et al., 2017b; Bobenko and Springborn, 2007; Xin et al., 2011;
Fisher et al., 2006; Sharp et al., 2019c; Sharp and Crane, 2020a; Gillespie et
al., 2021; Liu et al., 2023). Intrinsic triangulations facilitate tasks such
as Delaunay refinement (Sharp et al., 2019c), mesh simplification (Liu et al.,
2023), and construction of differential operators on non-manifold meshes
(Sharp and Crane, 2020a). Similar to intrinsic triangulations (Sharp et al.,
2019c), our method uses endpoints as only degrees of freedom—connecting
geodesic edges are defined implicitly and reconstructed on demand. Unlike
existing work, our approach enables the computation of first and second
derivatives of geodesic distances in closed form, thus opening the door to
efficient second-order minimization algorithms.

#### Structural Curve Networks

Structural curve networks have gained increasing attention in the fields of
computer graphics and robotics due to their aesthetic appeal and practical
utility. For instance, Schumacher  _et al_. (2018) characterize the mechanical
behavior of different families of tiling patterns, and Li  _et al_. (2022)
explore the direction-dependent stiffness of 3D-printed weave structures.
Another line of research focuses on the structural stability of curve networks
(Zehnder et al., 2016; Pérez et al., 2015; Miguel et al., 2016; Ren et al.,
2021; Neveu et al., 2022; Liu et al., 2021). Our work is similar in the sense
that we operate on curves embedded within surfaces. Instead of representing
these curves explicitly, however, they are represented implicitly as shortest
geodesic paths. Our differentiable geodesic distance formulation allows us to
solve intrinsic minimization problems defined on these implicit curve
networks.

#### Embedded Simulation

Many physical phenomena can be described through partial differential
equations on Riemannian manifolds. Examples include swirl dynamics on soap
bubbles (Yang et al., 2019; Huang et al., 2020), skin sliding (Li et al.,
2013), tight-fitting clothing (Montes et al., 2020), and elastic curve
networks embedded in curved surfaces (Zehnder et al., 2016). One approach to
this embedded elasticity problem is to use 3D spline curves whose control
vertices are constrained to lie on the surface (Lee and Lee, 2001; Hofer and
Pottmann, 2004; Wallner et al., 2005). Eulerian-on-Lagrangian methods are an
alternative representation for simulating constrained motion of deformable
bodies (Fan et al., 2013), rods (Sueda et al., 2011; Li et al., 2022), and
cloth (Cirio et al., 2014, 2015; Weidner et al., 2018). Within this context,
Li  _et al_. (2013) simulate skin sliding using texture-like material
coordinates as Eulerian degrees of freedom. In order for this approach to
work, however, a texture atlas of the underlying surface with quasi-isometric
charts is required. Montes  _et al_. (2020) propose a Lagrangian-on-Lagrangian
approach that combines subdivision surfaces with embedded triangle meshes to
allow for smooth sliding of skin-tight clothing across the underlying body.
However, they approximate on-surface lengths using Euclidean distances.
Consequently, surface and cloth discretizations must have adequate resolutions
to limit approximation error. Our simulation examples with embedded elastic
networks likewise use a Lagrangian-on-Lagrangian representation. However, we
avoid resolution dependence by computing exact geodesic distances on triangle
meshes.

##  3\. Background and Challenges

Our goal is to minimize distance-based functions on triangle meshes with
second-order optimizers. Naturally, this requires the distance metric to be
sufficiently smooth and differentiable. While Euclidean distance satisfies
these criteria, depending on the resolution of the embedded mesh, it can
deviate from the actual distance between two surface points, _i.e_., the
geodesic distance, to different degrees. This deviation can lead to undesired
local minima, hindering the optimization process (see Sec. 5.7). We therefore
opt to use geodesic distances for simulation. Before describing our method in
detail, we first briefly discuss the smoothness of geodesic distance.

###  3.1. Geodesic Distance

In the continuous setting, geodesics are locally shortest paths between two
points on a surface. These paths are also straightest, _i.e_., they exhibit
zero geodesic curvature. However, in the discrete setting with non-smooth
surfaces, these two definitions—straightest and shortest path—are not
equivalent and will generally lead to different geodesic paths. The shortest
path definition is a natural choice for boundary value problems with known
endpoints, while the straightest-path definition is best suited for initial
value problems with a given starting point and tangent direction. We refer to
Crane _et al_. (2020) for an in-depth discussion of this subject. In this
work, we use shortest-path geodesics, which can be computed by minimizing a
convex quadratic potential. To use this distance metric in a minimization
algorithm, we must first understand its smoothness properties.

###  3.2. Smoothness of Geodesic Distance

The geodesic distance between two points is a continuous function of the
points everywhere on the surface. The geodesic path, however, is not always
unique and can change abruptly for smooth motion of the endpoints. The
distance derivative is necessarily discontinuous in these singular
configurations. It is worth mentioning that this discontinuity exists even for
smooth surfaces and does not originate from discretization.

#### Continuous Setting.

To understand this discontinuity, consider a surface patch containing a bump
and a geodesic whose endpoints lie on opposite sides of the bump (see inset
figure). As we translate the upper endpoint to the right, the geodesic path
reaches a singular point where it suddenly flips to the right. At the
singularity, there exist two paths on different sides of the bump that have
the same geodesic distance. A local perturbation of the endpoint can easily
lead to flipping of the shortest geodesic. The collection of points

for which geodesics are not unique is referred to as the cut locus. The
shortest geodesic path from point $a$ to point $b$ is discontinuous as $b$
passes through the cut locus of $a$, and along the cut locus, the gradient of
the path is undefined. Although this discontinuity in the distance gradient
exists even in the continuous setting, it does not prevent gradient-based
minimization since points on the cut locus locally maximize the geodesic
distance from the source point.

#### Discrete Setting.

Geodesic distance remains continuous in the discrete setting, but geodesic
paths behave somewhat differently compared to the smooth setting. Special care
is required when geodesics pass close to mesh vertices, which we classify into
three categories depending on their discrete Gaussian curvature: spherical
vertices (positive curvature), hyperbolic vertices (negative curvature), and
planar vertices (zero curvature).

Figure 2. Behavior of geodesic paths passing over different types of vertices.
Fixing one endpoint of a geodesic, we translate the other such that the path
moves across the center vertex. Whereas the geodesic passes through the center
vertex in the planar and hyperbolic case, it jumps over the spherical vertex
to avoid the local distance maximum.

In the vicinity of planar vertices, the geodesic is a straight line and its
length reduces to the Euclidean distance. For spherical vertices, however, a
geodesic cannot pass through the vertex as there always exists a shorter path
going around it (see Fig. 2). Since geodesic paths vary discontinuously around
spherical vertices, geodesic distance is only $C^{0}$-continuous. Fortunately,
paths through spherical vertices are length-maximizing and are thus avoided
during distance minimization. Lastly, geodesic paths can continuously pass
through hyperbolic vertices, but the change in local tangent direction leads
to only $C^{1}$-continuity in these cases. In Sec. 4.1 we show that
$C^{2}$-continuity can be restored by adding a suitable mollifier.

###  3.3. Challenges

Our approach builds on the observation that, while geodesic paths are not
generally continuous functions of their endpoints, their lengths vary
continuously and can thus be used for gradient-based minimization. Computing
the required derivatives, however, is no trivial task. Geodesics on triangle
meshes generally span multiple faces and intersect with the edges of the
hosting surface. Explicitly tracking a varying number of intersection
variables is cumbersome. Using endpoints as degrees of freedom circumvents
this issue, but intersection points are then functions of the endpoints,
adding another layer of complexity to the computation of first and second
derivatives. In fact, we are not aware of any existing work that computes the
analytical derivatives of geodesic distance on triangle meshes. In the
following, we develop a differentiable geodesic distance formulation whose
first and second derivatives can be computed efficiently. Since our
formulation requires only intrinsic quantities, we refer to it as intrinsic
minimization.

##  4\. Intrinsic Minimization

We develop a formulation for differentiable geodesic distance on triangle
meshes that uses geodesic endpoints as the only explicit variables. The
intersection points on the geodesic path are implicitly defined through
equilibrium conditions of a shortest-path energy functional (Sec. 4.1). We
show that an intrinsic simulation framework driven by geodesic distance can be
formulated on this basis (Sec. 4.2) and the required derivatives can be
obtained in simple forms using sensitivity analysis and geometric insights. We
extend our distance-driven formulation from geodesic edges to triangles such
as to construct directionally-continuous deformation energies (Sec. 4.3). We
furthermore elaborate on two-way coupling effects between geodesic networks
and their embedding surfaces (Sec. 4.4). Finally, beyond simulating embedded
elasticity, we generalize the intrinsic minimization paradigm to
differentiable geodesic Voronoi diagrams (Sec. 4.5).

###  4.1. Differentiable Geodesic Distance

Our approach builds on the robust MMP algorithm (Mitchell et al., 1987) to
compute exact geodesic paths on triangle meshes. The endpoints of geodesic
paths are represented using barycentric coordinates $\mathbf{w}$ of the
hosting triangle mesh. These endpoints are the only simulation degrees of
freedom of a given geodesic path. The spatial coordinates of the endpoints
$\mathbf{c}$, and the intersections $\mathbf{x}$ between the path and edges of
the hosting mesh, are computed when either endpoint is moved (see inset). We
use a scalar parameter $t_{i}$ to denote the position of an intersection point
$\mathbf{x}_{i}$ along a corresponding mesh edge $e_{jk}$,

(1) $\mathbf{x}_{i}=\mathbf{v}_{j}+t_{i}(\mathbf{v}_{k}-\mathbf{v}_{j})\ ,$

where $\mathbf{v}$ are mesh vertices. We compute the geodesic distance $g$
between two points, $\mathbf{c}_{0}$ and $\mathbf{c}_{1}$, on the surface mesh
as the sum of the lengths of the line segments $l_{i}$ comprising the geodesic
path,

(2)
$g(\mathbf{c}_{0},\mathbf{c}_{1})=\sum_{i}l_{i}(\mathbf{c}_{0},\mathbf{c}_{1},\mathbf{x}(\mathbf{c}_{0},\mathbf{c}_{1}))\
.$

The relationship of the simulation variables $\mathbf{w}$ and geodesic
endpoints $\mathbf{c}$ are given explicitly by interpolating barycentric
coordinates. The relationship of the intersection variables $\mathbf{t}$ are
given implicitly by the algorithm for computing exact geodesics. All involved
quantities ultimately depend only on the simulation variables $\mathbf{w}$ and
we therefore write

(3)
$g(\mathbf{c},\mathbf{x})=g(\mathbf{c}(\mathbf{w}),\mathbf{x}(\mathbf{t}(\mathbf{c}(\mathbf{w})))\
.$

To enable Newton-type solvers for intrinsic minimization, we must compute the
first and second derivatives of the above expression. The first derivative is
given by

(4) $\frac{\mathrm{d}g}{\mathrm{d}\mathbf{w}}=\frac{\partial
g}{\partial\mathbf{c}}^{\mathsf{T}}\frac{\partial\mathbf{c}}{\partial\mathbf{w}}+\frac{\partial
g}{\partial\mathbf{x}}^{\mathsf{T}}\frac{\partial\mathbf{x}}{\partial\mathbf{t}}\frac{\partial\mathbf{t}}{\partial\mathbf{c}}\frac{\partial\mathbf{c}}{\partial\mathbf{w}}\
,$

and the second derivative follows as

(5)
$\displaystyle\frac{\mathrm{d}^{2}g}{\mathrm{d}\mathbf{w}^{2}}=\frac{\partial\mathbf{c}}{\partial\mathbf{w}}^{\mathsf{T}}\frac{\partial^{2}g}{\partial\mathbf{c}^{2}}\frac{\partial\mathbf{c}}{\partial\mathbf{w}}+\frac{\partial\mathbf{c}}{\partial\mathbf{w}}^{\mathsf{T}}\frac{\partial^{2}g}{\partial\mathbf{c}\partial\mathbf{x}}\frac{\partial\mathbf{x}}{\partial\mathbf{t}}\frac{\partial\mathbf{t}}{\partial\mathbf{c}}\frac{\partial\mathbf{c}}{\partial\mathbf{w}}+\sum_{i}\frac{\partial
g}{\partial c_{i}}\frac{\partial^{2}c_{i}}{\partial\mathbf{w}^{2}}$
$\displaystyle+\left(\frac{\partial\mathbf{x}}{\partial\mathbf{t}}\frac{\partial\mathbf{t}}{\partial\mathbf{c}}\frac{\partial\mathbf{c}}{\partial\mathbf{w}}\right)^{\mathsf{T}}\left(\frac{\partial^{2}g}{\partial\mathbf{c}\partial\mathbf{x}}\frac{\partial\mathbf{c}}{\partial\mathbf{w}}+\frac{\partial^{2}g}{\partial\mathbf{x}^{2}}\frac{\partial\mathbf{x}}{\partial\mathbf{t}}\frac{\partial\mathbf{t}}{\partial\mathbf{c}}\frac{\partial\mathbf{c}}{\partial\mathbf{w}}\right)$
$\displaystyle+\left(\frac{\partial\mathbf{t}}{\partial\mathbf{c}}\frac{\partial\mathbf{c}}{\partial\mathbf{w}}\right)^{\mathsf{T}}\sum_{i}\frac{\partial
g}{\partial
x_{i}}\frac{\partial^{2}x_{i}}{\partial\mathbf{t}^{2}}\left(\frac{\partial\mathbf{t}}{\partial\mathbf{c}}\frac{\partial\mathbf{c}}{\partial\mathbf{w}}\right)$
$\displaystyle+\frac{\partial\mathbf{c}}{\partial\mathbf{w}}^{\mathsf{T}}\left(\sum_{i}\frac{\partial
g}{\partial\mathbf{x}}\frac{\partial\mathbf{x}}{\partial
t_{i}}\frac{\partial^{2}t_{i}}{\partial\mathbf{c}^{2}}\right)\frac{\partial\mathbf{c}}{\partial\mathbf{w}}+\sum_{i}\frac{\partial
g}{\partial\mathbf{x}}\frac{\partial\mathbf{x}}{\partial\mathbf{t}}\frac{\partial\mathbf{t}}{\partial
c_{i}}\frac{\partial^{2}c_{i}}{\partial\mathbf{w}^{2}}\ .$

While computing all of these terms is feasible, we show that both expressions
can be significantly simplified using geometric insight.

The shortest geodesic is a local minimizer of distance, giving rise to the
first-order optimality condition

(6) $\frac{\mathrm{d}g}{\mathrm{d}\mathbf{t}}=\frac{\partial
g}{\partial\mathbf{x}}^{\mathsf{T}}\frac{\partial\mathbf{x}}{\partial\mathbf{t}}=\mathbf{0}\
.$

If the interior vertices $\mathbf{x}$ satisfy this optimality condition,
perturbations $\delta\mathbf{t}$ to the intersection variables will not change
length to first order: moving along the path changes adjacent segment lengths,
but these changes sum to zero; moving orthogonal to the path simply rotates
the segments, which does not change their lengths to first order. We can
therefore simplify the total derivative to

(7) $\frac{\mathrm{d}g}{\mathrm{d}\mathbf{w}}=\frac{\partial
g}{\partial\mathbf{c}}^{\mathsf{T}}\frac{\partial\mathbf{c}}{\partial\mathbf{w}}\
.$

This expression is readily evaluated since it only requires the gradient of
length for the first and last segments and the (constant) Jacobian
$\frac{\partial\mathbf{c}}{\partial\mathbf{w}}$ of the barycentric
interpolation function.

Similarly, the Hessian can also be greatly simplified using implicit
differentiation. Note that all terms can be evaluated algebraically except for
$\frac{\partial\mathbf{t}}{\partial\mathbf{c}}$, _i.e_., the derivatives of
the edge intersection points with respect to the endpoints of the geodesic.
Computing this derivative explicitly requires differentiation through the
entire pipeline of the algorithm used for computing the exact geodesics. Apart
from the complexity of differentiating such a code, this strategy must fail
for cases in which geodesics are not unique and, hence, path derivatives do
not exist. Our key insight is that this problem can be entirely avoided by
implicit differentiation. We begin by differentiating both sides of (6) w.r.t.
$\mathbf{c}$,

(8)
$\frac{\partial\mathbf{x}}{\partial\mathbf{t}}^{\mathsf{T}}\left(\frac{\partial^{2}g}{\partial\mathbf{c}\partial\mathbf{x}}+\frac{\partial^{2}g}{\partial\mathbf{x}^{2}}\frac{\partial\mathbf{x}}{\partial\mathbf{t}}\frac{\partial\mathbf{t}}{\partial\mathbf{c}}\right)=\mathbf{0}\
.$

To obtain $\frac{\partial\mathbf{t}}{\partial\mathbf{c}}$, we rearrange and
solve the small linear system

(9)
$\left(\frac{\partial\mathbf{x}}{\partial\mathbf{t}}^{\mathsf{T}}\frac{\partial^{2}g}{\partial\mathbf{x}^{2}}\frac{\partial\mathbf{x}}{\partial\mathbf{t}}\right)\frac{\partial\mathbf{t}}{\partial\mathbf{c}}=-\frac{\partial\mathbf{x}}{\partial\mathbf{t}}^{\mathsf{T}}\frac{\partial^{2}g}{\partial\mathbf{c}\partial\mathbf{x}}\
.$

The effective size of this $m\times m$ system depends on the number $m$ of
intersections between a given geodesic and the hosting mesh. Multiplying (8)
by
$(\frac{\partial\mathbf{t}}{\partial\mathbf{c}}\frac{\partial\mathbf{c}}{\partial\mathbf{w}})^{\mathsf{T}}$
from the left and by $\frac{\partial\mathbf{c}}{\partial\mathbf{w}}$ from the
right gives

(10)
$\displaystyle\left(\frac{\partial\mathbf{x}}{\partial\mathbf{t}}\frac{\partial\mathbf{t}}{\partial\mathbf{c}}\frac{\partial\mathbf{c}}{\partial\mathbf{w}}\right)^{\mathsf{T}}\left(\frac{\partial^{2}g}{\partial\mathbf{c}\partial\mathbf{x}}\frac{\partial\mathbf{c}}{\partial\mathbf{w}}+\frac{\partial^{2}g}{\partial\mathbf{x}^{2}}\frac{\partial\mathbf{x}}{\partial\mathbf{t}}\frac{\partial\mathbf{t}}{\partial\mathbf{c}}\frac{\partial\mathbf{c}}{\partial\mathbf{w}}\right)$
$\displaystyle=\mathbf{0}\ .$

Substituting (6) and (10) into (5) and removing second derivatives of linear
terms, we obtain the substantially simplified expression

(11)
$\frac{\partial^{2}g}{\partial\mathbf{w}^{2}}=\frac{\partial\mathbf{c}}{\partial\mathbf{w}}^{\mathsf{T}}\frac{\partial^{2}g}{\partial\mathbf{c}^{2}}\frac{\partial\mathbf{c}}{\partial\mathbf{w}}+\frac{\partial\mathbf{c}}{\partial\mathbf{w}}^{\mathsf{T}}\frac{\partial^{2}g}{\partial\mathbf{c}\partial\mathbf{x}}\frac{\partial\mathbf{x}}{\partial\mathbf{t}}\frac{\partial\mathbf{t}}{\partial\mathbf{c}}\frac{\partial\mathbf{c}}{\partial\mathbf{w}}\
.$

#### Discontinuities & Mollification

With the expression for second derivatives in hand, we must still address
potential discontinuities when intersection points $\mathbf{x}(\mathbf{t})$
approach mesh vertices. Depending on to which side of a mesh vertex the
geodesic passes, it will intersect different mesh edges. The function
$\mathbf{x}(t)$ is therefore ill-defined when the geodesic path coincides with
a mesh vertex, _i.e_., when $t=0$ or $t=1$ for individual edges, and the
optimality condition (Eq. 6) does not hold. This discontinuity hinders
convergence, as shortest geodesics tend to pass through hyperbolic vertices
(Fig. 3a), and the total distance energy is not $C^{2}$-continuous. We resolve
this issue by adding a mollifier to $\mathbf{x}(t)$ such that
$\frac{d\mathbf{x}(t)}{dt}$ smoothly vanishes at vertex intersections.

A similar discontinuity appears when endpoints of geodesics coincide with
vertices of the hosting mesh. We therefore apply the same mollifier to the
barycentric coordinates of the endpoints.

Figure 3. Mollification for edge intersection points. We smoothly blend linear
edge intersection trajectories $\mathbf{x}(t)$ with cubic functions (b) to
achieve $C^{2}$-continuity when geodesics pass through mesh vertices (a).

#### Smooth Mollifier

We smoothly blend the linear intersection point parameterization (1) with two
cubic functions (see Fig. 3b) such that the derivatives vanish when the
geodesic curve passes through mesh vertices, _i.e_., when $t=0$ or $t=1$. The
mollifier function is defined as

(12)
$\hat{t}(t)=\begin{cases}-\frac{t^{3}}{\epsilon^{2}}+\frac{2t^{2}}{\epsilon}&0\leq
t<\epsilon,\\\ t&\epsilon\leq t\leq 1-\epsilon\\\
-\frac{(t-1)^{3}}{\epsilon^{2}}+\frac{2(t-1)^{2}}{\epsilon}+1&1-\epsilon<t\leq
1\ ,\end{cases}$

where $\epsilon=10^{-6}$ defines the smoothing length. This mollifier is
$C^{1}$-continuous along an edge, and as shown in Sec. 5.7, significantly
accelerates the convergence of our Newton solver.

###  4.2. Elastic Geodesic Networks

With our differentiable geodesic distance formulation established, we now turn
to energy minimization. As a first example, we consider elastic curve networks
comprised of geodesic springs. The total energy of this system is

(13)
$E_{\mathrm{network}}(\mathbf{w})=\sum_{i}(g_{i}(\mathbf{w})-\bar{g}_{i})^{2}\
,$

where $i$ runs over all geodesic springs, whose current and rest lengths we
denote as $g_{i}$ and $\bar{g}_{i}$, respectively. We minimize this energy
using Newton’s method with a standard backtracking line search. The search
direction per iteration is computed by solving the linear system

(14) $\mathbf{H}(\mathbf{w})\Delta\mathbf{w}=-\mathbf{y}(\mathbf{w})\ ,$

where $\mathbf{y}(\mathbf{w})$ and $\mathbf{H}(\mathbf{w})$ are the gradient
and Hessian of the objective function. We use adaptive diagonal regularization
(Nocedal and Wright, 1999) to ensure that the Hessian is positive definite.
Since all simulation variables are barycentric coordinates, the resulting
search direction yields endpoint updates expressed in the barycentric
coordinates of the corresponding mesh triangles. Inevitably, endpoints will
move across edges and vertices into neighboring triangles, which requires
updating the search direction to the new coordinate system. To this end, we
use the approach by Sharp  _et al_. (2019c) for tracing straightest geodesics
(Polthier and Schmies, 2006) and update local coordinates as

(15) $\displaystyle\mathbf{w}^{j+1}$
$\displaystyle=\mathbf{w}^{j}+\alpha\cdot\text{tracing}(\mathbf{w}^{j},\Delta\mathbf{w})\
,$

where $\alpha$ is a step size determined by line search to ensure monotonic
decrease in energy.

#### Karcher Means

Using a slightly different formulation, our method can be readily extended to
compute so-called Karcher means on triangle meshes. Consider a special case of
a geodesic spring network, where a point $\mathbf{p}$ is connected to several
anchored points $\mathbf{x}$ on a surface. The energy in this case amounts to
the variational formulation for the Karcher mean,

(16)
$E_{\mathrm{Karcher}}(\mathbf{w})=\frac{1}{2N}\sum_{i=1}^{N}g(\mathbf{p}(\mathbf{w}),\mathbf{x}_{i})^{2}\
.$

Sharp  _et al_. (2019d) showed that the gradient of this energy can be
approximated efficiently by solving a sparse linear system. While their
approach is robust and general, using only first-order derivatives precludes
second-order convergence. Our analytical second-order derivatives of geodesic
distance, in contrast, enable Newton-type minimization, which leads to
significantly improved convergence (see Sec. 5.2).

###  4.3. Elastic Geodesic Triangles

In the previous examples, we considered potentials defined on geodesic
networks. We now extend our formulation to model elastic membranes made of
geodesic triangles, _i.e_., finite element-like triangles with geodesic edges.

For linear triangle finite elements in Euclidean geometry—also known as
constant strain triangles—the deformation gradient maps edge vectors from the
undeformed configuration to corresponding deformed edges as

(17) $\mathbf{F}\bar{\mathbf{e}}_{ij}=\mathbf{e}_{ij}\ ,$

where $\bar{\mathbf{e}}_{ij}=(\bar{\mathbf{x}}_{j}-\bar{\mathbf{x}}_{i})$ and
${\mathbf{e}}_{ij}=({\mathbf{x}}_{j}-{\mathbf{x}}_{i})$ are undeformed and
deformed edges, respectively. Squaring both sides, we obtain

(18)
$\bar{\mathbf{e}}_{ij}^{\mathsf{T}}\mathbf{C}\bar{\mathbf{e}}_{ij}=|\mathbf{e}_{ij}|^{2}\
,\quad\mathbf{C}=\mathbf{F}^{\mathsf{T}}\mathbf{F}\ ,$

where $\mathbf{C}$ is the Cauchy-Green tensor. A direct translation of this
deformation measure to geodesic triangles would require geodesic edges,
expressed in a common coordinate system. This representation, however, is not
available since the vertices of a geodesic triangle generally fall into
different triangles of the hosting surface. Fortunately, we can construct the
Cauchy-Green tensor using only edge lengths (see inset figure). To this end,
we first note that

(19)
$\bar{\mathbf{e}}_{ij}^{\mathsf{T}}\mathbf{C}\bar{\mathbf{e}}_{ij}=g_{ij}^{2}\
,$

for all three edges of a geodesic triangle. Using the fact that $\mathbf{C}$
is a symmetric $2\times 2$ tensor, it is evident that its three unknown
entries can be computed by solving the three equations (19). With this
nonlinear deformation measurement at hand, we can then use a standard Neo-
Hookean constitutive model for membrane elasticity. The corresponding energy
density function is given as

(20)
$\Psi(\mathbf{C})=\frac{\mu}{2}(tr(\mathbf{C})-2-2\log(J))+\frac{\lambda}{2}\log(J)^{2}\
,$

where $J=\sqrt{\text{det}(\mathbf{C})}$. Since it only involves the lengths of
geodesic edges, this energy density is a continuous function of the simulation
variables $\mathbf{w}$.

#### Discussion

We assume that triangles are intrinsically flat when computing deformation
gradients from edge lengths. While this is an approximation, it allows us to
measure and penalize deformations in any direction, not just along edges. As
another point, although it is possible to compute the area enclosed by a
geodesic triangle, this

area can be discontinuous as the geodesic path flips in the vicinity of
spherical vertices. An example is shown in the inset figure, where a small
perturbation of one vertex changes the area of the two adjacent triangles
abruptly. Instead of directly computing the area of geodesic triangles, our
formulation captures changes in area through the determinant of the Cauchy-
Green tensor, $\text{det}(\mathbf{C})$. Since computing $\mathbf{C}$ requires
only the lengths of geodesic edges, not their paths, we avoid discontinuities
in area.

###  4.4. Two-way Coupling

The previous examples focused on geodesic elastic systems embedded in a rigid
hosting surface. However, our approach can also be extended to full two-way
coupling with deformable hosting objects. For these cases, the geodesic
distance function takes the form

(21)
$g(\mathbf{c}(\mathbf{w},\mathbf{v}),\mathbf{x}(\mathbf{t}(\mathbf{c}(\mathbf{w},\mathbf{v}),\mathbf{v})),\mathbf{v})\
,$

where $\mathbf{v}$ are the vertices of the hosting surface. The first and
second derivatives of this expression can again be obtained using sensitivity
analysis; see App. A. We describe the coupled system through its potential
energy,

(22)
$E_{\mathrm{coupling}}(\mathbf{q})=E_{\mathrm{network}}(\mathbf{q})+E_{\mathrm{host}}(\mathbf{q})\
,$

where $\mathbf{q}=(\mathbf{w},\mathbf{v})^{\mathsf{T}}$ concatenates the
barycentric coordinates of the embedded system and the vertices on the hosting
surface. The first term $E_{\mathrm{network}}(\mathbf{q})$ models the energy
of the embedded elastic system, whereas $E_{\mathrm{host}}(\mathbf{q})$ is the
elastic potential of the hosting object. We consider two types of models for
the hosting object: a discrete shell model (Grinspun et al., 2003) and a
standard volumetric finite element model (Kim and Eberle, 2020).

###  4.5. Differentiable Geodesic Voronoi Diagrams

Up to this point, we have focused on elastic geodesic systems embedded in
rigid and deformable hosting objects. We now extend our formulation to another
application of intrinsic distances, geodesic Voronoi diagrams.

#### Geodesic Voronoi Diagrams

For a given set of sample points (sites) on a surface, a surface Voronoi
diagram partitions the surface mesh $\mathcal{M}$ into cells. A cell $C_{i}$
is the set of all points $\mathbf{x}^{*}$ satisfying

(23) $C_{i}=\\{\mathbf{x}^{*}\in\mathcal{M}\ |\
g(\mathbf{x}^{*},\mathbf{s}_{i})<g(\mathbf{x}^{*},\mathbf{s}_{j}),\ \forall
j\neq i\\}\ ,$

where $\mathbf{s}_{i}$ and $\mathbf{s}_{j}$ are the barycentric coordinates of
two sites. The corresponding surface Voronoi diagram, using geodesic distance
as the metric, is referred to as a geodesic Voronoi diagram (GVD). To
represent the Voronoi diagram, we use only the site locations $\mathbf{s}$ as
explicit degrees of freedom. The cell boundary vertices $\tilde{\mathbf{x}}$,
which are either Voronoi vertices (with 3 or more adjacent cells) or
intersections between Voronoi edges and mesh edges (with 2 adjacent cells),
are given implicitly as the unique solution to the equidistance constraint

(24)
$g(\tilde{\mathbf{x}}_{i},\mathbf{s}_{0})-g(\tilde{\mathbf{x}}_{i},\mathbf{s}_{j})=0,j=1\dots
N_{i}-1$

where $N_{i}$ is the number of adjacent cells to cell boundary vertex
$\tilde{\mathbf{x}}_{i}$ and $\mathbf{s}_{0},\dots,\mathbf{s}_{N_{i}-1}$ are
the sites of adjacent cells. We compute the geodesic Voronoi diagram by
solving

(25)
$\tilde{\mathbf{x}}_{i}=\operatorname*{argmin}\limits_{\tilde{\mathbf{x}}^{*}}\>\sum_{j=1}^{N_{i}-1}(g(\tilde{\mathbf{x}}^{*},\mathbf{s}_{0})-g(\tilde{\mathbf{x}}^{*},\mathbf{s}_{j}))^{2}\
$

for each boundary vertex. Note that this requires prior knowledge of cell
connectivity, which we obtain using a fast GVD approximation (Xin et al.,
2022); see also Algorithm 1.

#### Bi-level Optimization

We now turn to optimizing geodesic Voronoi diagrams for higher-level objective
functions. Consider a constrained optimization problem of the form

(26)
$\min_{\tilde{\mathbf{x}},\mathbf{s}}O(\tilde{\mathbf{x}},\mathbf{s})\quad\text{s.t.}\quad\mathbf{f}(\tilde{\mathbf{x}},\
\mathbf{s})=\mathbf{0}\ ,$

where $O$ is an outer design objective function and $\mathbf{f}$ is the
vector-valued equidistance constraint function given by the gradient of the
inner objective (25). The states $\tilde{\mathbf{x}}$ are coupled with the
variables $\mathbf{s}$ through the constraint and we write
$\tilde{\mathbf{x}}=\tilde{\mathbf{x}}(\mathbf{s})$ to eliminate
$\tilde{\mathbf{x}}$ as explicit degrees of freedom. Optimizing this objective
with gradient-based methods requires the derivative
$\frac{\mathrm{d}\tilde{\mathbf{x}}}{\mathrm{d}\mathbf{s}}$, which we obtain
through sensitivity analysis. We refer to App. C for derivations.

The equidistance constraints derive from the definition of Voronoi diagrams,
_i.e_., a point on a Voronoi edge must have the same distance to all adjacent
sites. This includes Voronoi vertices as well as intersections between Voronoi
edges and mesh edges, for which $\mathbf{f}$ takes different forms. For
intersections between Voronoi edges and mesh edges, $\tilde{\mathbf{x}}_{i}$
is equidistant to two sites and $\mathbf{f}$ takes the form

(27)
$\mathbf{f}_{i}(\tilde{\mathbf{x}}_{i},\mathbf{s}_{0},\mathbf{s}_{1})=g(\tilde{\mathbf{x}}_{i}(\mathbf{s}_{0},\mathbf{s}_{1}),\mathbf{s}_{0})-g(\tilde{\mathbf{x}}_{i}(\mathbf{s}_{0},\mathbf{s}_{1}),\mathbf{s}_{1})\
.$

In this case, only one constraint equation is necessary to define a unique
point along the intersected mesh edge, parameterized by a scalar as in Eq.
(1). For Voronoi vertices that are equidistant to three or more sites (see
inset), we use two constraint equations

(28)
$\displaystyle\mathbf{f}_{i,0}(\tilde{\mathbf{x}}_{i},\mathbf{s}_{0},\mathbf{s}_{1})$
$\displaystyle=g(\tilde{\mathbf{x}}_{i}(\mathbf{s}),\mathbf{s}_{0})-g(\tilde{\mathbf{x}}_{i}(\mathbf{s}),\mathbf{s}_{1})\
,$
$\displaystyle\mathbf{f}_{i,1}(\tilde{\mathbf{x}}_{i},\mathbf{s}_{0},\mathbf{s}_{2})$
$\displaystyle=g(\tilde{\mathbf{x}}_{i}(\mathbf{s}),\mathbf{s}_{0})-g(\tilde{\mathbf{x}}_{i}(\mathbf{s}),\mathbf{s}_{2})\
.$

We minimize the inner objective (25) using Newton’s method and the outer
design objective (26) using quasi-Newton methods. A standard back-tracking
line search is applied to both optimization processes.

###  4.6. Implementation Details

Our code is implemented in C++ using Eigen (Guennebaud et al., 2010) for
primary data structure and Intel TBB for parallelization. Linear systems are
solved using CHOLMOD (Chen et al., 2008). We use the MMP algorithm (Mitchell
et al., 1987) implemented in the Geometry Central Library (Sharp et al.,
2019a) for computing exact geodesics. Finally, we credit Polyscope (Sharp et
al., 2019b) for producing figures. Our code is available through the
repository https://github.com/liyuesolo/DifferentiableGeodesics.

##  5\. Results

###  5.1. Energy-minimizing Geodesic Networks

In the first set of examples (Fig. 4), we simulate elastic curve networks
comprised of zero-length geodesic springs. We showcase our approach on two
challenging setups where the endpoints of the geodesic springs are positioned
close to mesh vertices or edges. Our approach converges robustly in both
scenarios. As can be seen in the close-up views, we arrive at energy minima
that differ significantly from the initial network.

Figure 4. Elastic geodesic spring networks. We initialize the nodes (shown in red) in close proximity to either the vertices of the hosting mesh or its edge midpoints. Our method converges robustly in both scenarios. Figure 5. Qualitative comparisons with the Vector Heat Method (Sharp et al., 2019d) for computing Karcher Means on a selection of meshes. The yellow curves show the optimization trajectories toward the Karcher mean of a given set of points (shown in red). The initial guesses (chosen randomly) are shown in white and the Karcher mean in blue. The color gradient indicates steps toward the solution. As can be seen from these examples, our second-order solver allows for larger steps toward the minima (col 1-3) and significantly fewer iterations to convergence (col 4-5). Table 1. Comparison to the Vector Heat Method (Sharp et al., 2019d) for computing Karcher Means. We report the runtimes and optimization statistics for the examples shown in Fig. 5. Heat* indicates using the convergence criterion in the official implementation (Sharp et al., 2019a), whereas Heat denotes using the same criterion as ours. While the Vector Heat Method converges to visually acceptable solutions in similar numbers of iterations, their gradient-based update is limited to linear convergence. Leveraging a Newton solver with analytical second derivatives, our approach converges quadratically. Examples |  Triangles |  Time[s] (Heat*/Heat/Ours) |  Iterations (Heat*/Heat/Ours) |  GradientNorm (Heat*/Heat/Ours) |  Objective (Heat*/Heat/Ours)
---|---|---|---|---|---
sphere |  1,280 |  0.0162 / 0.916 / 0.0284 |  2 / 200 / 2 |  4.35 $\times 10^{-4}$/ 3.04 $\times 10^{-5}$ / 2.31 $\times 10^{-11}$ |  0.0206 / 0.0204 / 0.00744
screwdriver |  6,786 |  0.119 / 17.5 / 0.291 |  3 / 200 / 3 |  3.04 $\times 10^{-4}$ / 2.67 $\times 10^{-4}$ / 1.57 $\times 10^{-11}$ |  0.00387 / 0.00383 / 0.00121
ear |  45,312 |  1.15 / 120 / 4.99 |  3 / 200 / 3 |  1.13 $\times 10^{-5}$ / 1.85 $\times 10^{-6}$ / 3.20 $\times 10^{-12}$ |  0.00974 / 0.00972 / 0.00225
protein |  60,820 |  17.0 / 120 / 12.8 |  25 / 200 / 8 |  1.81 $\times 10^{-1}$ / 9.64 $\times 10^{-2}$ / 3.05 $\times 10^{-7}$ |  0.970 / 0.991 / 0.208
spiral cup |  34,874 |  0.576 / 19.3 / 0.291 |  13 / 200 / 10 |  3.04 $\times 10^{-4}$ / 1.05 $\times 10^{-5}$ / 1.99 $\times 10^{-11}$ |  1.34 / 1.58 / 0.280

###  5.2. Karcher Means

Our approach enables second-order optimization for computing Karcher means on
triangle meshes. We compare our approach to the state-of-the-art Vector Heat
Method (Sharp et al., 2019d) on a set of examples with different geometries
and resolutions. While Sharp  _et al_. (2019d) propose an efficient method to
compute the energy gradient by solving sparse linear systems, their approach
is limited to linear convergence. As can be seen from the qualitative
comparisons in Fig. 5, our second-order solver enables larger steps and
converges to the solution in significantly fewer iterations. We use the open-
source reference implementation from the Geometry Central Library (Sharp et
al., 2019a), which adopts a relative convergence criterion based on the
current triangle area and step size. For our experiments, we impose an
absolute convergence criterion based on the gradient norm. Details on the
convergence criteria are given in App. B. We report statistics for both
convergence criteria in Table 1, showing that our approach converges to much
tighter tolerances with comparable performance. We set a maximum iteration of
200 for all approaches.

Another baseline for computing Karcher means is provided by Mancinelli and
Puppo (2023), which leverages approximate first and second derivatives of
geodesic distance. While these approximations enable fast computations, the
underlying convexity assumption for geodesic distance does not generally hold
for arbitrary triangle meshes. Unfortunately, the publicly available
algorithm111https://github.com/Claudiomancinelli90/RCM_on_meshes did not yield
successful outcomes for any of the problem instances shown in Fig. 5. We
therefore conduct a simpler test, computing the Karcher mean of five randomly
placed points on a spherical mesh (Fig. 5, column 1). We repeat this
experiment 100 times and report the average runtime, residual norm, and
success rate in Table 2. Despite being an order of magnitude slower, our
approach ensures robust convergence across all test cases.

Table 2. Quantitative comparison with Mancinelli and Puppo (2023) using 100 randomized problem instances on a spherical mesh. Averaged runtimes and gradient norms exclude failure cases. Method | Avg Time [s] | Avg —Grad— | Success Rate
---|---|---|---
MancinelliPuppo | $9.54\times 10^{-3}$ | $1.40\times 10^{-2}$ | 17%
Ours |  9.23 $\times 10^{-2}$ |  6.52 $\times 10^{-12}$ | 100%

###  5.3. Energy-minimizing Geodesic Triangles

We simulate membranes defined by elastic geodesic triangles. Inspired by
previous work from Li _et al_. (2013), we simulate an elastic membrane made of
geodesic triangles embedded in a rigid torus mesh (see Fig. 6). The torus mesh
is stretched along the horizontal axis of the image plane in a non-uniform
way. Simply keeping the barycentric coordinates of the embedded membrane
unchanged leads to large isolated distortions (Fig. 6, top row). Optimizing
these coordinates such as to minimize the membrane’s elastic energy yields
smoothly distributed deformations (Fig. 6, bottom row). Whereas simulating
this effect in a 2D parametric space requires UV unwrapping (Li et al., 2013),
our approach operates directly on the 3D surface mesh. Our Newton solver
converges in fewer than 5 iterations on average for this sequence.

Figure 6. Simulation of an elastic membrane using geodesic triangles. We
stretch the hosting torus mesh anisotropically to induce deformation of the
embedded triangles. Keeping the (barycentric) membrane coordinates unchanged
leads to large isolated distortions (top row). Optimizing coordinates such as
to minimize the membrane’s energy leads to smoothly distributed deformations
(bottom row).

###  5.4. Two-way Coupling

We now consider two-way coupling between a geodesic elastic network and its
hosting surface. In the first example, we simulate the tightening of an
elastic geodesic network embedded in an inflated spherical shell (see Fig. 7).
As the minimization proceeds, many regions of the shell bulge out in response
to the compression forces induce by the embedded network. The hosting surface
is simulated using a discrete shell model (Grinspun et al., 2003) augmented
with a volume preservation term. To avoid factorizing a dense Hessian matrix
resulting from the volume preservation term, we use the Sherman–Morrison
formula (1950) to compute the Newton step. In Fig. 1, we show that our
intrinsic minimization pipeline can be coupled with a volumetric mesh
discretized using tetrahedron finite elements. The hosting elastic object is
simulated using a Neo-Hookean constitutive model. The ears of the bunny bend
down and the back bulges in response to the tightening of the geodesic spring
network.

Figure 7. Two-way coupling. Simulation of a geodesic spring network on an
inflated spherical shell. Realistic bulging effects emerge upon tightening the
embedded curve network.

###  5.5. Optimization of Geodesic Voronoi Diagrams

We consider three examples with different design objectives defined on
geodesic Voronoi diagrams. Each example is formulated as a bi-level
optimization problem using sensitivity analysis as described in section 4.5.

#### Length similarity

In the first example, we optimize for uniform edge lengths such as to
facilitate potential manufacturing with equal-length rods. To this end, we
minimize the design objective

(29)
$\min_{\mathbf{s}}O(\tilde{\mathbf{x}}(\mathbf{s}),\mathbf{s})=\frac{1}{2}\sum_{(i,j)\in\mathcal{E}}(g(\tilde{\mathbf{x}}_{i},\tilde{\mathbf{x}}_{j})-L)^{2}\
,$

where $\mathcal{E}$ is the index set of all Voronoi edges and $L$ is the
target edge length, which we set to the average Voronoi edge length in the
initial configuration. As can be seen in Fig. 8, the optimized structures
exhibit much less variation in edge length compared to the initial
configurations.

Figure 8. Optimizing for uniform edge length. By minimizing the deviation of
all Voronoi edges from the target length, the optimized structure
significantly reduces the length variation compared to the initial
configuration.

#### Cell planarity

We consider the task of designing Voronoi diagrams with as-planar-as-possible
cells such that the 3D structure can be assembled from flat panels. We define
a corresponding objective as

(30)
$\min_{\mathbf{s}}O(\tilde{\mathbf{x}}(\mathbf{s}),\mathbf{s})=\frac{1}{2}\sum_{i=1}^{N}\sum_{j=1}^{N_{i}}\frac{1}{N_{i}}d(\tilde{\mathbf{x}}_{j},\mathcal{P}_{i}(\tilde{\mathbf{x}}))^{2}\
,$

where $N$ denotes the total number of Voronoi vertices and $N_{i}$ is the
number of vertices for Voronoi cell $i$. For each cell, we compute the least
squares-fitting plane $\mathcal{P}_{i}$ and penalize the corresponding point-
to-plane distance $d$ for each cell vertex. As can be seen from Fig. 9,
minimizing this planarity objective with our method leads to an order of
magnitude reduction in average and maximum vertex-to-plane distance.

Figure 9. Optimizing for planarity. We optimize site locations such that the
vertices of each Voronoi cell are as planar as possible. Our approach
successfully reduces the average and maximum vertex-to-plane distance by an
order of magnitude or more.

#### Cell regularity

Finally, we consider an objective to regularize cell shapes by minimizing the
squared distances between each site and its surrounding vertices,

(31)
$\min_{\mathbf{s}}O(\tilde{\mathbf{x}}(\mathbf{s}),\mathbf{s})=\frac{1}{2}\sum_{i=1}^{N}\sum_{j=1}^{N_{i}}g(\mathbf{s}_{i},\tilde{\mathbf{x}}_{j}(\mathbf{s}))^{2}\
.$

Fig. 10 illustrates the effect of this objective on a set of Voronoi diagrams
with irregular initial cell shapes. Despite the complexity of the hosting
surfaces, the resulting cell shape distributions show significant
improvements.

Figure 10. Optimizing cell regularity. Minimizing our regularity objective
leads to quasi-isotropic cells when starting from a poorly shaped initial
diagram.

###  5.6. Timings

We report detailed simulation statistics and averaged timings for all
examples. The timings listed in Tables 1—4 are measured on a workstation with
an AMD Ryzen Threadripper PRO 5995WX CPU.

Table 3. Simulation statistics and average timings. Example | Triangles | DoFs | Two-way | Exact Geodesics | Constructing Hessian | Solving Linear System | Newton Step
---|---|---|---|---|---|---|---
Fig. 1 | 9,396 | 20,838 | yes | 221.55 ms | 2.38 s | 690.8 ms | 3.64 s
Fig. 4 torus | 576 | 1,152 | no | 73.43 ms | 48.47 ms | 42.58 ms | 172.31 ms
Fig. 4 duck | 2,000 | 4,000 | no | 791.46 ms | 188.47 ms | 114.28 ms | 1.218 s
Fig. 5 sphere | 1,280 | 2 | no | 8.19 ms | 112.5 $\mu$s | 63.5 $\mu$s | 12.52 ms
Fig. 5 screwdriver | 6,786 | 2 | no | 63.20 ms | 12.24 ms | 43.67 $\mu$s | 94.47 ms
Fig. 5 ear | 45,312 | 2 | no | 1.18 s | 45.86 ms | 59.6 $\mu$s | 1.65 s
Fig. 5 protein | 60,820 | 2 | no | 710.72 ms | 948 ms | 40.86 $\mu$s | 1.668 s
Fig. 5 spiral cup | 34,874 | 2 | no | 214.67 ms | 285.3 ms | 40.4 $\mu$s | 555.91 ms
Fig. 6 | 2,304 | 282 | no | 89.36 ms | 5.91 ms | 596 $\mu$s | 96.75 ms
Fig. 7 | 1,280 | 30 | yes | 8.28 ms | 42.80 ms | 65.64 ms | 129.63 ms
Table 4. Statistics and average timing for energy-minimizing GVDs. Examples |  Construct GVD[ms] |  Time /Iter [s] |  #Tri |  #Sites
---|---|---|---|---
Fig. 8 sphere |  5.87 |  1.90 |  1,280 |  98
Fig. 8 duck |  9.97 |  4.39 |  2,000 |  125
Fig. 8 cactus |  9.71 |  4.03 |  2,000 |  107
Fig. 9 sphere |  14.37 |  0.474 |  1,280 |  396
Fig. 9 torus |  6.41 |  0.0733 |  576 |  155
Fig. 9 duck |  21.67 |  1.17 |  2,000 |  487
Fig. 10 fertility |  66.46 |  7.09 |  2,250 |  2,250
Fig. 10 cow |  93.12 |  39.3 |  5,856 |  2,928
Fig. 10 armadillo |  262.6 |  180 |  20,000 |  10,000

###  5.7. Ablation Study

#### Comparison with L-BFGS

We show that our analytical Hessian significantly improves convergence
compared to first-order or quasi-Newton methods without sacrificing
performance. We compare the convergence and performance of different methods
on two energy minimization problems involving geodesic spring networks. As can
be seen from Fig. 11 (left), while L-BFGS offers better convergence than
gradient descent (GD), Newton’s method outperforms both by a large margin. In
Fig. 11 (right), we show that our second-order approach also offers
significant performance gains. Thanks to the compact analytical expression of
the second derivative, computing the analytical Hessian of the geodesic
distance only takes a fraction of the overall computation time. The initial
and optimized configurations can be found in the inset figures. Dashed lines
in the figures indicate residual tolerances that lead to qualitatively on-par
results. We consider the results obtained from Newton and L-BFGS to be
visually indistinguishable when the maximum difference in nodal position is
smaller than $1\%$ of the scene bounding box diagonal.

Figure 11. Convergence comparison. We use our second-order approach, gradient
descent (GD), and L-BFGS to find equilibrium states of geodesic networks. Our
approach exhibits quadratic convergence and demonstrates significantly better
performance. Dashed lines indicate residual tolerances below which results
become visually indistinguishable. Inset figures show the system states
corresponding to the blue dots on the curves.

#### Euclidean Distance

Using Euclidean distance instead of geodesic distance greatly simplifies
computation but leads to additional and undesirable local minima in intrinsic
simulation. We illustrate this problem on an embedded curve consisting of
three vertices and two connecting segments modeled as zero-length springs
(shown in Fig. 12). The endpoints are anchored to the hosting surface while
the central vertex is free to move. We compare simulations using Euclidean and
geodesic distance metrics. In the left example, the hosting surface exhibits a
right angle, creating a local energy minimum for Euclidean-distance springs.
In the right example, multiple local minima exist for Euclidean springs and
the eventual equilibrium configuration depends on the initial position of the
central vertex. Using geodesic distance, however, the simulation converges to
the global optimum in both cases.

Figure 12. Spurious local minima when using Euclidean distance. In these two
examples, we consider two zero-length springs connected by a free vertex and
minimize the total elastic energy. Using Euclidean distance for the spring
energy leads to undesired local minima while using geodesic distance resolves
this problem. Geodesic paths are shown in all cases for visualization.

#### Mollification

To study the influence of our smooth mollifier, we set up a two-way coupling
example, where the geodesic path passes through a hyperbolic vertex in the
minimum-energy state (Fig. 13, right). In this example, a rectangle with an
inner crossbar is embedded in a deformable shell. The four corner vertices of
this rectangle are fixed to the surface and all edges are zero-length geodesic
springs. In this case, energy minimization results in the inner edge passing
through the hyperbolic vertex of the shell. As can be seen in the convergence
plot to the left, without the mollifier, Newton’s method struggles to converge
even after 200 iterations. When the energy is mollified into a
$C^{2}$-continuous function, Newton’s method exhibits quadratic convergence.

Figure 13. Mollification. In this two-way coupling example, the corners of the
orange rectangle are fixed while the two inner vertices are free to slide on
the surface. All embedded vertices are connected by geodesic springs with zero
rest length. To arrive at the energy minimum, the inner spring must pass
through the hyperbolic vertex of the hosting mesh. As can be seen from the
convergence plot, without the mollifier (orange), Newton’s method does not
achieve quadratic convergence due to the lack of $C^{2}$-continuity at the
minimum. With the mollifier, however, our approach converges quadratically
(blue).

##  6\. Conclusion & Future Work

We proposed a novel differentiable geodesic distance formulation for intrinsic
minimization on triangle meshes. By applying the implicit function theorem to
the variational formulation of shortest-path geodesics, we demonstrate that
closed-form first and second derivatives can be obtained and significantly
simplified, allowing for the use of Newton-type minimization solvers. We use
our method to minimize a diverse set of objectives, including two-way coupling
of embedded elastic structures with hosting geometries, Karcher means on
arbitrary triangle meshes, and differentiable geodesic Voronoi diagrams.

#### Limitations

While our intrinsic minimization paradigm enables robust second-order
optimization on many instances of embedded elasticity, in almost all cases,
the geodesics are under tension and the energy increases with geodesic length.
However, when a geodesic curve is under compression and increasing length is
energetically favorable, local distance maxima along the cut locus may lead to
convergence failure.

We use the open-source implementation of the MMP algorithm for computing
geodesic distance. While we perform this computation largely in parallel, it
remains the computational bottleneck of our method. While our triangle
elements based on geodesic edge lengths are promising for simulating embedded
elastic membranes, the approximation accuracy decreases as the curvature
within each geodesic triangle increases. Using denser geodesic triangulations
can alleviate this issue to some extent, but using internal quadrature points
and additional degrees of freedom could be a promising direction. Finally, we
only focused on triangle meshes in this work. Extending our intrinsic
minimization paradigm to polygonal meshes is another avenue for future
exploration.

###### Acknowledgements.

We thank Peiyuan Xie for his assistance with the visualization, and the
anonymous reviewers for their valuable feedback. This work was supported by
the European Research Council (ERC) under the European Union’s Horizon 2020
research and innovation program (grant agreement No. 866480), and the Swiss
National Science Foundation through SNF project grant 200021_200644.

## References

  * (1)
  * Belyaev and Fayolle (2020) Alexander Belyaev and Pierre-Alain Fayolle. 2020.  An ADMM-based scheme for distance function approximation.  Numerical Algorithms 84 (2020), 983–996.
  * Belyaev and Fayolle (2015) Alexander G Belyaev and Pierre-Alain Fayolle. 2015.  On variational and PDE-based distance function approximations. In Computer Graphics Forum, Vol. 34. Wiley Online Library, 104–118.
  * Bobenko and Springborn (2007) Alexander I Bobenko and Boris A Springborn. 2007.  A discrete Laplace–Beltrami operator for simplicial surfaces.  Discrete & Computational Geometry 38, 4 (2007), 740–756.
  * Bose et al. (2011) Prosenjit Bose, Anil Maheshwari, Chang Shu, and Stefanie Wuhrer. 2011.  A survey of geodesic paths on 3D surfaces.  Computational Geometry 44, 9 (2011), 486–498.
  * Chen and Han (1990) Jindong Chen and Yijie Han. 1990.  Shortest paths on a polyhedron. In Proceedings of the sixth annual symposium on Computational geometry. 360–369.
  * Chen et al. (2008) Yanqing Chen, Timothy A Davis, William W Hager, and Sivasankaran Rajamanickam. 2008\.  Algorithm 887: CHOLMOD, supernodal sparse Cholesky factorization and update/downdate.  ACM Transactions on Mathematical Software (TOMS) 35, 3 (2008), 1–14.
  * Cirio et al. (2014) Gabriel Cirio, Jorge Lopez-Moreno, David Miraut, and Miguel A Otaduy. 2014.  Yarn-level simulation of woven cloth.  ACM Transactions on Graphics (TOG) 33, 6 (2014), 1–11.
  * Cirio et al. (2015) Gabriel Cirio, Jorge Lopez-Moreno, and Miguel A Otaduy. 2015\.  Efficient simulation of knitted cloth using persistent contacts. In Proceedings of the 14th ACM SIGGRAPH/Eurographics Symposium on Computer Animation. 55–61.
  * Crane et al. (2020) Keenan Crane, Marco Livesu, Enrico Puppo, and Yipeng Qin. 2020.  A Survey of Algorithms for Geodesic Paths and Distances.  arXiv e-prints, Article arXiv:2007.10430 (July 2020), arXiv:2007.10430 pages.  arXiv:cs.GR/2007.10430
  * Crane et al. (2013) Keenan Crane, Clarisse Weischedel, and Max Wardetzky. 2013\.  Geodesics in heat: A new approach to computing distance based on heat flow.  ACM Transactions on Graphics (TOG) 32, 5 (2013), 1–11.
  * Edelstein et al. (2023) Michal Edelstein, Nestor Guillen, Justin Solomon, and Mirela Ben-Chen. 2023.  A Convex Optimization Framework for Regularized Geodesic Distances. In ACM SIGGRAPH 2023 Conference Proceedings (SIGGRAPH ’23). Association for Computing Machinery, New York, NY, USA, Article 2, 11 pages.  https://doi.org/10.1145/3588432.3591523
  * Fan et al. (2013) Ye Fan, Joshua Litven, David IW Levin, and Dinesh K Pai. 2013\.  Eulerian-on-lagrangian simulation.  ACM Transactions on Graphics (TOG) 32, 3 (2013), 1–9.
  * Fisher et al. (2006) Matthew Fisher, Boris Springborn, Alexander I Bobenko, and Peter Schroder. 2006.  An algorithm for the construction of intrinsic Delaunay triangulations with applications to digital geometry processing.  In ACM SIGGRAPH 2006 Courses. 69–74.
  * Gillespie et al. (2021) Mark Gillespie, Nicholas Sharp, and Keenan Crane. 2021\.  Integer coordinates for intrinsic geometry processing.  ACM Transactions on Graphics (TOG) 40, 6 (2021), 1–13.
  * Grinspun et al. (2003) Eitan Grinspun, Anil N Hirani, Mathieu Desbrun, and Peter Schröder. 2003.  Discrete shells. In Proceedings of the 2003 ACM SIGGRAPH/Eurographics symposium on Computer animation. Citeseer, 62–67.
  * Guennebaud et al. (2010) Gaël Guennebaud, Benoît Jacob, et al. 2010\.  Eigen v3.  http://eigen.tuxfamily.org. (2010).
  * Hofer and Pottmann (2004) Michael Hofer and Helmut Pottmann. 2004.  Energy-minimizing splines in manifolds.  In ACM SIGGRAPH 2004 Papers. 284–293.
  * Huang et al. (2020) Weizhen Huang, Julian Iseringhausen, Tom Kneiphof, Ziyin Qu, Chenfanfu Jiang, and Matthias B Hullin. 2020.  Chemomechanical simulation of soap film flow on spherical bubbles.  ACM Transactions on Graphics (TOG) 39, 4 (2020), 41–1.
  * Kim and Eberle (2020) Theodore Kim and David Eberle. 2020.  Dynamic deformables: implementation and production practicalities.  In ACM SIGGRAPH 2020 Courses. 1–182.
  * Kimmel and Sethian (1998) Ron Kimmel and James A Sethian. 1998.  Computing geodesic paths on manifolds.  Proceedings of the national academy of Sciences 95, 15 (1998), 8431–8435.
  * Lee and Lee (2001) Yun-Jin Lee and Seung-Yong Lee. 2001.  Geometric Snakes for Triangular Meshes.  Journal of the Korea Computer Graphics Society 7, 3 (2001), 9–18.
  * Li et al. (2013) Duo Li, Shinjiro Sueda, Debanga R Neog, and Dinesh K Pai. 2013\.  Thin skin elastodynamics.  ACM Transactions on Graphics (TOG) 32, 4 (2013), 1–10.
  * Li et al. (2022) Yue Li, Juan Montes, Bernhard Thomaszewski, and Stelian Coros. 2022.  Programmable digital weaves.  IEEE Robotics and Automation Letters 7, 2 (2022), 2891–2896.
  * Liu et al. (2017a) Bangquan Liu, Shuangmin Chen, Shi-Qing Xin, Ying He, Zhen Liu, and Jieyu Zhao. 2017a.  An optimization-driven approach for computing geodesic paths on triangle meshes.  Computer-Aided Design 90 (2017), 105–112.
  * Liu et al. (2023) Hsueh-Ti Derek Liu, Mark Gillespie, Benjamin Chislett, Nicholas Sharp, Alec Jacobson, and Keenan Crane. 2023.  Surface Simplification using Intrinsic Error Metrics.  ACM Transactions on Graphics (TOG) 42, 4 (2023), 1–17.
  * Liu et al. (2021) Jian Liu, Shiqing Xin, Xifeng Gao, Kaihang Gao, Kai Xu, Baoquan Chen, and Changhe Tu. 2021.  Computational Object-Wrapping Rope Nets.  ACM Transactions on Graphics (TOG) 41, 1 (2021), 1–16.
  * Liu et al. (2017b) Yong-Jin Liu, Dian Fan, Chun-Xu Xu, and Ying He. 2017b.  Constructing intrinsic Delaunay triangulations from the dual of geodesic Voronoi diagrams.  ACM Transactions on Graphics (TOG) 36, 2 (2017), 1–15.
  * Mancinelli and Puppo (2023) Claudio Mancinelli and Enrico Puppo. 2023.  Computing the Riemannian center of mass on meshes.  Computer Aided Geometric Design 103 (2023), 102203.
  * Miguel et al. (2016) Eder Miguel, Mathias Lepoutre, and Bernd Bickel. 2016\.  Computational design of stable planar-rod structures.  ACM Transactions on Graphics (TOG) 35, 4 (2016), 1–11.
  * Mitchell et al. (1987) Joseph SB Mitchell, David M Mount, and Christos H Papadimitriou. 1987.  The discrete geodesic problem.  SIAM J. Comput. 16, 4 (1987), 647–668.
  * Montes et al. (2020) Juan Montes, Bernhard Thomaszewski, Sudhir Mudur, and Tiberiu Popa. 2020.  Computational design of skintight clothing.  ACM Transactions on Graphics (TOG) 39, 4 (2020), 105–1.
  * Neveu et al. (2022) William Neveu, Ivan Puhachov, Bernhard Thomaszewski, and Mikhail Bessmeltsev. 2022.  Stability-aware simplification of curve networks. In ACM SIGGRAPH 2022 Conference Proceedings. 1–9.
  * Nocedal and Wright (1999) Jorge Nocedal and Stephen J Wright. 1999.  Numerical optimization.  Springer.
  * Pang et al. (2023) Bo Pang, Zhongtian Zheng, Guoping Wang, and Peng-Shuai Wang. 2023\.  Learning the geodesic embedding with graph neural networks.  ACM Transactions on Graphics (TOG) 42, 6 (2023), 1–12.
  * Pérez et al. (2015) Jesús Pérez, Bernhard Thomaszewski, Stelian Coros, Bernd Bickel, José A Canabal, Robert Sumner, and Miguel A Otaduy. 2015\.  Design and fabrication of flexible rod meshes.  ACM Transactions on Graphics (TOG) 34, 4 (2015), 1–12.
  * Peyré et al. (2010) Gabriel Peyré, Mickael Péchaud, Renaud Keriven, Laurent D Cohen, et al. 2010\.  Geodesic methods in computer vision and graphics.  Foundations and Trends® in Computer Graphics and Vision 5, 3–4 (2010), 197–397.
  * Polthier and Schmies (2006) Konrad Polthier and Markus Schmies. 2006.  Straightest geodesics on polyhedral surfaces.  In ACM SIGGRAPH 2006 Courses. 30–38.
  * Qin et al. (2016) Yipeng Qin, Xiaoguang Han, Hongchuan Yu, Yizhou Yu, and Jianjun Zhang. 2016.  Fast and exact discrete geodesic computation based on triangle-oriented wavefront propagation.  ACM Transactions on Graphics (TOG) 35, 4 (2016), 1–13.
  * Ren et al. (2021) Yingying Ren, Julian Panetta, Tian Chen, Florin Isvoranu, Samuel Poincloux, Christopher Brandt, Alison Martin, and Mark Pauly. 2021\.  3D weaving with curved ribbons.  ACM Transactions on Graphics 40, ARTICLE (2021), 127\.
  * Schumacher et al. (2018) Christian Schumacher, Steve Marschner, Markus Gross, and Bernhard Thomaszewski. 2018.  Mechanical characterization of structured sheet materials.  ACM Transactions on Graphics (TOG) 37, 4 (2018), 1–15.
  * Sharp et al. (2019b) Nicholas Sharp et al. 2019b.  Polyscope.  (2019).  www.polyscope.run.
  * Sharp and Crane (2020a) Nicholas Sharp and Keenan Crane. 2020a.  A laplacian for nonmanifold triangle meshes. In Computer Graphics Forum, Vol. 39. Wiley Online Library, 69–80.
  * Sharp and Crane (2020b) Nicholas Sharp and Keenan Crane. 2020b.  You can find geodesic paths in triangle meshes by just flipping edges.  ACM Transactions on Graphics (TOG) 39, 6 (2020), 1–15.
  * Sharp et al. (2019a) Nicholas Sharp, Keenan Crane, et al. 2019a.  GeometryCentral: A modern C++ library of data structures and algorithms for geometry processing.  https://geometry-central.net/.  (2019).
  * Sharp et al. (2019c) Nicholas Sharp, Yousuf Soliman, and Keenan Crane. 2019c.  Navigating Intrinsic Triangulations.  ACM Trans. Graph. 38, 4 (2019).
  * Sharp et al. (2019d) Nicholas Sharp, Yousuf Soliman, and Keenan Crane. 2019d.  The vector heat method.  ACM Transactions on Graphics (TOG) 38, 3 (2019), 1–19.
  * Sherman and Morrison (1950) Jack Sherman and Winifred J Morrison. 1950.  Adjustment of an inverse matrix corresponding to a change in one element of a given matrix.  The Annals of Mathematical Statistics 21, 1 (1950), 124–127.
  * Sueda et al. (2011) Shinjiro Sueda, Garrett L Jones, David IW Levin, and Dinesh K Pai. 2011.  Large-scale dynamic simulation of highly constrained strands.  In ACM SIGGRAPH 2011 papers. 1–10.
  * Surazhsky et al. (2005) Vitaly Surazhsky, Tatiana Surazhsky, Danil Kirsanov, Steven J Gortler, and Hugues Hoppe. 2005\.  Fast exact and approximate geodesics on meshes.  ACM transactions on graphics (TOG) 24, 3 (2005), 553–560.
  * Trettner et al. (2021) Philip Trettner, David Bommes, and Leif Kobbelt. 2021\.  Geodesic distance computation via virtual source propagation. In Computer graphics forum, Vol. 40. Wiley Online Library, 247–260.
  * Wallner et al. (2005) Johannes Wallner, Helmut Pottmann, and Michael Hofer. 2005\.  Fair curve networks in nonlinear geometries.  In ACM SIGGRAPH 2005 Sketches. 32–es.
  * Weidner et al. (2018) Nicholas J Weidner, Kyle Piddington, David IW Levin, and Shinjiro Sueda. 2018.  Eulerian-on-lagrangian cloth simulation.  ACM Transactions on Graphics (TOG) 37, 4 (2018), 1–11.
  * Xin et al. (2022) Shiqing Xin, Pengfei Wang, Rui Xu, Dongming Yan, Shuangmin Chen, Wenping Wang, Caiming Zhang, and Changhe Tu. 2022\.  SurfaceVoronoi: Efficiently Computing Voronoi Diagrams Over Mesh Surfaces with Arbitrary Distance Solvers.  ACM Transactions on Graphics (TOG) 41, 6 (2022), 1–12.
  * Xin et al. (2011) Shi-Qing Xin, Shuang-Min Chen, Ying He, Guo-Jin Wang, Xianfeng Gu, and Hong Qin. 2011\.  Isotropic mesh simplification by evolving the geodesic delaunay triangulation. In 2011 Eighth International Symposium on Voronoi Diagrams in Science and Engineering. IEEE, 39–47.
  * Yang et al. (2019) Bowen Yang, William Corse, Jiecong Lu, Joshuah Wolper, and Chen-Fanfu Jiang. 2019.  Real-time fluid simulation on the surface of a sphere.  Proceedings of the ACM on Computer Graphics and Interactive Techniques 2, 1 (2019), 1–17.
  * Ying et al. (2013) Xiang Ying, Xiaoning Wang, and Ying He. 2013.  Saddle vertex graph (SVG) a novel solution to the discrete geodesic problem.  ACM Transactions on Graphics (TOG) 32, 6 (2013), 1–12.
  * Zehnder et al. (2016) Jonas Zehnder, Stelian Coros, and Bernhard Thomaszewski. 2016\.  Designing structurally-sound ornamental curve networks.  ACM Transactions on Graphics (TOG) 35, 4 (2016), 1–10.

##  Appendix A First and Second Derivative of Geodesic Distance for Two-way
Coupling

The geodesic distance in this case is given by (21), which we also include
here for completeness

(32)
$g(\mathbf{c}(\mathbf{w},\mathbf{v}),\mathbf{x}(\mathbf{t}(\mathbf{c}(\mathbf{w},\mathbf{v}),\mathbf{v})),\mathbf{v}).$

The gradient of $g$ w.r.t. $\mathbf{w}$ is

(33) $\frac{\mathrm{d}g}{\mathrm{d}\mathbf{w}}=\frac{\partial
g}{\partial\mathbf{c}}^{\mathsf{T}}\frac{\partial\mathbf{c}}{\partial\mathbf{w}}+\frac{\partial
g}{\partial\mathbf{x}}^{\mathsf{T}}\frac{\partial\mathbf{x}}{\partial\mathbf{t}}\frac{\partial\mathbf{t}}{\partial\mathbf{c}}\frac{\partial\mathbf{c}}{\partial\mathbf{w}},$

and gradient of $g$ w.r.t. $\mathbf{v}$ is

(34) $\frac{\mathrm{d}g}{\mathrm{d}\mathbf{v}}=\frac{\partial
g}{\partial\mathbf{c}}^{\mathsf{T}}\frac{\partial\mathbf{c}}{\partial\mathbf{v}}+\frac{\partial
g}{\partial\mathbf{x}}^{\mathsf{T}}\frac{\partial\mathbf{x}}{\partial\mathbf{t}}(\frac{\partial\mathbf{t}}{\partial\mathbf{c}}\frac{\partial\mathbf{c}}{\partial\mathbf{v}}+\frac{\partial\mathbf{t}}{\partial\mathbf{v}})+\frac{\partial
g}{\partial\mathbf{x}}^{\mathsf{T}}\frac{\partial\mathbf{x}}{\partial\mathbf{v}}.$

Following the same procedure as detailed in Sec. 4.1, we again use the first-
order optimality condition,

(35) $\frac{\mathrm{d}g}{\mathrm{d}\mathbf{t}}:=\frac{\partial
g}{\partial\mathbf{x}}^{\mathsf{T}}\frac{\partial\mathbf{x}}{\partial\mathbf{t}}=\mathbf{0}$

to simplify the gradient expression to

(36) $\frac{\mathrm{d}g}{\mathrm{d}\mathbf{w}}=\frac{\partial
g}{\partial\mathbf{c}}^{\mathsf{T}}\frac{\partial\mathbf{c}}{\partial\mathbf{w}},\quad\frac{\mathrm{d}g}{\mathrm{d}\mathbf{v}}=\frac{\partial
g}{\partial\mathbf{c}}^{\mathsf{T}}\frac{\partial\mathbf{c}}{\partial\mathbf{v}}+\frac{\partial
g}{\partial\mathbf{x}}^{\mathsf{T}}\frac{\partial\mathbf{x}}{\partial\mathbf{v}}.$

The hessian is more complex due to the coupling. The complete expression
include these individual blocks,

(37)
$\displaystyle\frac{\partial^{2}g}{\partial\mathbf{w}^{2}}=\frac{\partial\mathbf{c}}{\partial\mathbf{w}}^{\mathsf{T}}\frac{\partial^{2}g}{\partial\mathbf{c}^{2}}\frac{\partial\mathbf{c}}{\partial\mathbf{w}}+\frac{\partial\mathbf{c}}{\partial\mathbf{w}}^{\mathsf{T}}\frac{\partial^{2}g}{\partial\mathbf{c}\partial\mathbf{x}}\frac{\partial\mathbf{x}}{\partial\mathbf{t}}\frac{\partial\mathbf{t}}{\partial\mathbf{c}}\frac{\partial\mathbf{c}}{\partial\mathbf{w}}+\sum_{i}\frac{\partial
g}{\partial c_{i}}\frac{\partial^{2}c_{i}}{\partial\mathbf{w}^{2}}$
$\displaystyle+\left(\frac{\partial\mathbf{x}}{\partial\mathbf{t}}\frac{\partial\mathbf{t}}{\partial\mathbf{c}}\frac{\partial\mathbf{c}}{\partial\mathbf{w}}\right)^{\mathsf{T}}\left(\frac{\partial^{2}g}{\partial\mathbf{c}\partial\mathbf{x}}\frac{\partial\mathbf{c}}{\partial\mathbf{w}}+\frac{\partial^{2}g}{\partial\mathbf{x}^{2}}\frac{\partial\mathbf{x}}{\partial\mathbf{t}}\frac{\partial\mathbf{t}}{\partial\mathbf{c}}\frac{\partial\mathbf{c}}{\partial\mathbf{w}}\right)$
$\displaystyle+\left(\frac{\partial\mathbf{t}}{\partial\mathbf{c}}\frac{\partial\mathbf{c}}{\partial\mathbf{w}}\right)^{\mathsf{T}}\sum_{i}\frac{\partial
g}{\partial
x_{i}}\frac{\partial^{2}x_{i}}{\partial\mathbf{t}^{2}}\left(\frac{\partial\mathbf{t}}{\partial\mathbf{c}}\frac{\partial\mathbf{c}}{\partial\mathbf{w}}\right)$
$\displaystyle+\frac{\partial\mathbf{c}}{\partial\mathbf{w}}^{\mathsf{T}}\left(\sum_{i}\frac{\partial
g}{\partial\mathbf{x}}\frac{\partial\mathbf{x}}{\partial
t_{i}}\frac{\partial^{2}t_{i}}{\partial\mathbf{c}^{2}}\right)\frac{\partial\mathbf{c}}{\partial\mathbf{w}}+\sum_{i}\frac{\partial
g}{\partial\mathbf{x}}\frac{\partial\mathbf{x}}{\partial\mathbf{t}}\frac{\partial\mathbf{t}}{\partial
c_{i}}\frac{\partial^{2}c_{i}}{\partial\mathbf{w}^{2}}.$ (38)
$\displaystyle\frac{\partial^{2}g}{\partial\mathbf{v}^{2}}=\frac{\partial\mathbf{c}}{\partial\mathbf{v}}^{\mathsf{T}}\frac{\partial^{2}g}{\partial\mathbf{c}^{2}}\frac{\partial\mathbf{c}}{\partial\mathbf{v}}+\frac{\partial\mathbf{c}}{\partial\mathbf{v}}^{\mathsf{T}}\frac{\partial^{2}g}{\partial\mathbf{c}\partial\mathbf{x}}\left(\frac{\partial\mathbf{x}}{\partial\mathbf{t}}\left(\frac{\partial\mathbf{t}}{\partial\mathbf{c}}\frac{\partial\mathbf{c}}{\partial\mathbf{v}}+\frac{\partial\mathbf{t}}{\partial\mathbf{v}}\right)+\frac{\partial\mathbf{x}}{\partial\mathbf{v}}\right)$
$\displaystyle+\sum_{i}\frac{\partial g}{\partial
c_{i}}\frac{\partial^{2}c_{i}}{\partial\mathbf{v}^{2}}+\left(\frac{\partial\mathbf{x}}{\partial\mathbf{t}}\left(\frac{\partial\mathbf{t}}{\partial\mathbf{c}}\frac{\partial\mathbf{c}}{\partial\mathbf{v}}+\frac{\partial\mathbf{t}}{\partial\mathbf{v}}\right)+\frac{\partial\mathbf{x}}{\partial\mathbf{v}}\right)^{\mathsf{T}}\frac{\partial^{2}g}{\partial\mathbf{c}\partial\mathbf{x}}^{\mathsf{T}}\frac{\partial\mathbf{c}}{\partial\mathbf{v}}$
$\displaystyle+\left(\frac{\partial\mathbf{x}}{\partial\mathbf{t}}\left(\frac{\partial\mathbf{t}}{\partial\mathbf{c}}\frac{\partial\mathbf{c}}{\partial\mathbf{v}}+\frac{\partial\mathbf{t}}{\partial\mathbf{v}}\right)+\frac{\partial\mathbf{x}}{\partial\mathbf{v}}\right)^{\mathsf{T}}\frac{\partial^{2}g}{\partial\mathbf{x}^{2}}\left(\frac{\partial\mathbf{x}}{\partial\mathbf{t}}\left(\frac{\partial\mathbf{t}}{\partial\mathbf{c}}\frac{\partial\mathbf{c}}{\partial\mathbf{v}}+\frac{\partial\mathbf{t}}{\partial\mathbf{v}}\right)+\frac{\partial\mathbf{x}}{\partial\mathbf{v}}\right)$
$\displaystyle+\sum_{i}(\frac{\partial\mathbf{t}}{\partial\mathbf{c}}\frac{\partial\mathbf{c}}{\partial\mathbf{v}}+\frac{\partial\mathbf{t}}{\partial\mathbf{v}})^{\mathsf{T}}\frac{\partial
g}{\partial
x_{i}}\left(\frac{\partial^{2}x_{i}}{\partial\mathbf{t}^{2}}\left(\frac{\partial\mathbf{t}}{\partial\mathbf{c}}\frac{\partial\mathbf{c}}{\partial\mathbf{v}}+\frac{\partial\mathbf{t}}{\partial\mathbf{v}}\right)+\frac{\partial^{2}x_{i}}{\partial\mathbf{t}\partial\mathbf{v}}\right)$
$\displaystyle+\frac{\partial
g}{\partial\mathbf{x}}^{\mathsf{T}}\sum_{i}\frac{\partial\mathbf{x}}{\partial
t_{i}}\left(\frac{\partial\mathbf{c}}{\partial\mathbf{v}}^{\mathsf{T}}\left(\frac{\partial^{2}t_{i}}{\partial\mathbf{c}^{2}}\frac{\partial\mathbf{c}}{\partial\mathbf{v}}+\frac{\partial^{2}t_{i}}{\partial\mathbf{c}\partial\mathbf{v}}\right)+\sum_{j}\frac{\partial
t_{i}}{\partial
c_{j}}\frac{\partial^{2}c_{j}}{\partial\mathbf{v}^{2}}+\frac{\partial^{2}t_{i}}{\partial\mathbf{c}\partial\mathbf{v}}\frac{\partial\mathbf{c}}{\partial\mathbf{v}}\right)$
$\displaystyle+\frac{\partial
g}{\partial\mathbf{x}}^{\mathsf{T}}\sum_{i}\frac{\partial\mathbf{x}}{\partial
t_{i}}\left(\frac{\partial^{2}t_{i}}{\partial\mathbf{v}^{2}}\right)+\sum_{i}\frac{\partial
g}{\partial
x_{i}}\left(\frac{\partial^{2}x_{i}}{\partial\mathbf{t}\partial\mathbf{v}}\left(\frac{\partial\mathbf{t}}{\partial\mathbf{c}}\frac{\partial\mathbf{c}}{\partial\mathbf{v}}+\frac{\partial\mathbf{t}}{\partial\mathbf{v}}\right)+\frac{\partial^{2}x_{i}}{\partial\mathbf{v}^{2}}\right).$
(39)
$\displaystyle\frac{\partial^{2}g}{\partial\mathbf{v}\partial\mathbf{w}}=\frac{\partial\mathbf{c}}{\partial\mathbf{w}}^{\mathsf{T}}\left(\frac{\partial^{2}g}{\partial\mathbf{c}^{2}}\frac{\partial\mathbf{c}}{\partial\mathbf{v}}+\frac{\partial^{2}g}{\partial\mathbf{c}\partial\mathbf{x}}\left(\frac{\partial\mathbf{x}}{\partial\mathbf{t}}\left(\frac{\partial\mathbf{t}}{\partial\mathbf{c}}\frac{\partial\mathbf{c}}{\partial\mathbf{v}}+\frac{\partial\mathbf{t}}{\partial\mathbf{v}}\right)+\frac{\partial\mathbf{x}}{\partial\mathbf{v}}\right)\right)$
$\displaystyle+\left(\frac{\partial\mathbf{x}}{\partial\mathbf{t}}\frac{\partial\mathbf{t}}{\partial\mathbf{c}}\frac{\partial\mathbf{c}}{\partial\mathbf{w}}\right)^{\mathsf{T}}\left(\frac{\partial^{2}g}{\partial\mathbf{c}\partial\mathbf{x}}\frac{\partial\mathbf{c}}{\partial\mathbf{v}}+\frac{\partial^{2}g}{\partial\mathbf{x}^{2}}\left(\frac{\partial\mathbf{x}}{\partial\mathbf{t}}\left(\frac{\partial\mathbf{t}}{\partial\mathbf{c}}\frac{\partial\mathbf{c}}{\partial\mathbf{v}}+\frac{\partial\mathbf{t}}{\partial\mathbf{v}}\right)+\frac{\partial\mathbf{x}}{\partial\mathbf{v}}\right)\right)$
$\displaystyle+\sum_{i}(\frac{\partial\mathbf{t}}{\partial\mathbf{c}}\frac{\partial\mathbf{c}}{\partial\mathbf{w}})^{\mathsf{T}}\frac{\partial
g}{\partial
x_{i}}\left(\frac{\partial^{2}x_{i}}{\partial\mathbf{t}^{2}}\left(\frac{\partial\mathbf{t}}{\partial\mathbf{c}}\frac{\partial\mathbf{c}}{\partial\mathbf{v}}+\frac{\partial\mathbf{t}}{\partial\mathbf{v}}\right)+\frac{\partial^{2}x_{i}}{\partial\mathbf{t}\partial\mathbf{v}}\right)$
$\displaystyle+\frac{\partial\mathbf{c}}{\partial\mathbf{w}}^{\mathsf{T}}\left(\frac{\partial
g}{\partial\mathbf{x}}\sum_{i}\frac{\partial\mathbf{x}}{\partial
t_{i}}\left(\frac{\partial^{2}t_{i}}{\partial\mathbf{c}^{2}}\frac{\partial\mathbf{c}}{\partial\mathbf{v}}+\frac{\partial^{2}t_{i}}{\partial\mathbf{c}\partial\mathbf{v}}\right)\right)+\sum_{i}\frac{\partial
g}{\partial\mathbf{x}}^{\mathsf{T}}\frac{\partial\mathbf{x}}{\partial\mathbf{t}}\frac{\partial\mathbf{t}}{\partial
c_{i}}\frac{\partial^{2}c_{i}}{\partial\mathbf{v}\partial\mathbf{w}}$
$\displaystyle+\sum_{i}\frac{\partial g}{\partial c_{i}}\frac{\partial
c_{i}^{2}}{\partial\mathbf{v}\partial\mathbf{w}}.$

We begin by simplifying the hessian expression using (35), and after removing
the second derivatives of linear functions we obtain:

(40) $\displaystyle\frac{\partial^{2}g}{\partial\mathbf{w}^{2}}$
$\displaystyle=\frac{\partial\mathbf{c}}{\partial\mathbf{w}}^{\mathsf{T}}\frac{\partial^{2}g}{\partial\mathbf{c}^{2}}\frac{\partial\mathbf{c}}{\partial\mathbf{w}}+\frac{\partial\mathbf{c}}{\partial\mathbf{w}}^{\mathsf{T}}\frac{\partial^{2}g}{\partial\mathbf{c}\partial\mathbf{x}}\frac{\partial\mathbf{x}}{\partial\mathbf{t}}\frac{\partial\mathbf{t}}{\partial\mathbf{c}}\frac{\partial\mathbf{c}}{\partial\mathbf{w}}$
$\displaystyle+\left(\frac{\partial\mathbf{x}}{\partial\mathbf{t}}\frac{\partial\mathbf{t}}{\partial\mathbf{c}}\frac{\partial\mathbf{c}}{\partial\mathbf{w}}\right)^{\mathsf{T}}\left(\frac{\partial^{2}g}{\partial\mathbf{c}\partial\mathbf{x}}\frac{\partial\mathbf{c}}{\partial\mathbf{w}}+\frac{\partial^{2}g}{\partial\mathbf{x}^{2}}\frac{\partial\mathbf{x}}{\partial\mathbf{t}}\frac{\partial\mathbf{t}}{\partial\mathbf{c}}\frac{\partial\mathbf{c}}{\partial\mathbf{w}}\right).$

(41)
$\displaystyle\frac{\partial^{2}g}{\partial\mathbf{v}^{2}}=\frac{\partial\mathbf{c}}{\partial\mathbf{v}}^{\mathsf{T}}\frac{\partial^{2}g}{\partial\mathbf{c}^{2}}\frac{\partial\mathbf{c}}{\partial\mathbf{v}}+\frac{\partial\mathbf{c}}{\partial\mathbf{v}}^{\mathsf{T}}\frac{\partial^{2}g}{\partial\mathbf{c}\partial\mathbf{x}}\left(\frac{\partial\mathbf{x}}{\partial\mathbf{t}}\left(\frac{\partial\mathbf{t}}{\partial\mathbf{c}}\frac{\partial\mathbf{c}}{\partial\mathbf{v}}+\frac{\partial\mathbf{t}}{\partial\mathbf{v}}\right)+\frac{\partial\mathbf{x}}{\partial\mathbf{v}}\right)$
$\displaystyle+\left(\frac{\partial\mathbf{x}}{\partial\mathbf{t}}\left(\frac{\partial\mathbf{t}}{\partial\mathbf{c}}\frac{\partial\mathbf{c}}{\partial\mathbf{v}}+\frac{\partial\mathbf{t}}{\partial\mathbf{v}}\right)+\frac{\partial\mathbf{x}}{\partial\mathbf{v}}\right)^{\mathsf{T}}\frac{\partial^{2}g}{\partial\mathbf{c}\partial\mathbf{x}}^{\mathsf{T}}\frac{\partial\mathbf{c}}{\partial\mathbf{v}}$
$\displaystyle+\left(\frac{\partial\mathbf{x}}{\partial\mathbf{t}}\left(\frac{\partial\mathbf{t}}{\partial\mathbf{c}}\frac{\partial\mathbf{c}}{\partial\mathbf{v}}+\frac{\partial\mathbf{t}}{\partial\mathbf{v}}\right)+\frac{\partial\mathbf{x}}{\partial\mathbf{v}}\right)^{\mathsf{T}}\frac{\partial^{2}g}{\partial\mathbf{x}^{2}}\left(\frac{\partial\mathbf{x}}{\partial\mathbf{t}}\left(\frac{\partial\mathbf{t}}{\partial\mathbf{c}}\frac{\partial\mathbf{c}}{\partial\mathbf{v}}+\frac{\partial\mathbf{t}}{\partial\mathbf{v}}\right)+\frac{\partial\mathbf{x}}{\partial\mathbf{v}}\right)$
$\displaystyle+\sum_{i}\frac{\partial g}{\partial
x_{i}}\left(\frac{\partial^{2}x_{i}}{\partial\mathbf{t}\partial\mathbf{v}}\left(\frac{\partial\mathbf{t}}{\partial\mathbf{c}}\frac{\partial\mathbf{c}}{\partial\mathbf{v}}+\frac{\partial\mathbf{t}}{\partial\mathbf{v}}\right)\right)+\sum_{i}\frac{\partial
g}{\partial
x_{i}}\left(\frac{\partial^{2}x_{i}}{\partial\mathbf{t}\partial\mathbf{v}}\left(\frac{\partial\mathbf{t}}{\partial\mathbf{c}}\frac{\partial\mathbf{c}}{\partial\mathbf{v}}+\frac{\partial\mathbf{t}}{\partial\mathbf{v}}\right)\right)^{\mathsf{T}}.$

The mixed second-order partial derivative of $g$ w.r.t. $\mathbf{w}$ and
$\mathbf{v}$ reads

(42)
$\displaystyle\frac{\partial^{2}g}{\partial\mathbf{v}\partial\mathbf{w}}=\frac{\partial\mathbf{c}}{\partial\mathbf{w}}^{\mathsf{T}}\left(\frac{\partial^{2}g}{\partial\mathbf{c}^{2}}\frac{\partial\mathbf{c}}{\partial\mathbf{v}}+\frac{\partial^{2}g}{\partial\mathbf{c}\partial\mathbf{x}}\left(\frac{\partial\mathbf{x}}{\partial\mathbf{t}}\left(\frac{\partial\mathbf{t}}{\partial\mathbf{c}}\frac{\partial\mathbf{c}}{\partial\mathbf{v}}+\frac{\partial\mathbf{t}}{\partial\mathbf{v}}\right)+\frac{\partial\mathbf{x}}{\partial\mathbf{v}}\right)\right)$
$\displaystyle+\left(\frac{\partial\mathbf{x}}{\partial\mathbf{t}}\frac{\partial\mathbf{t}}{\partial\mathbf{c}}\frac{\partial\mathbf{c}}{\partial\mathbf{w}}\right)^{\mathsf{T}}\left(\frac{\partial^{2}g}{\partial\mathbf{c}\partial\mathbf{x}}\frac{\partial\mathbf{c}}{\partial\mathbf{v}}+\frac{\partial^{2}g}{\partial\mathbf{x}^{2}}\left(\frac{\partial\mathbf{x}}{\partial\mathbf{t}}\left(\frac{\partial\mathbf{t}}{\partial\mathbf{c}}\frac{\partial\mathbf{c}}{\partial\mathbf{v}}+\frac{\partial\mathbf{t}}{\partial\mathbf{v}}\right)+\frac{\partial\mathbf{x}}{\partial\mathbf{v}}\right)\right)$
$\displaystyle+\sum_{i}(\frac{\partial\mathbf{t}}{\partial\mathbf{c}}\frac{\partial\mathbf{c}}{\partial\mathbf{w}})^{\mathsf{T}}\frac{\partial
g}{\partial
x_{i}}\frac{\partial^{2}x_{i}}{\partial\mathbf{t}\partial\mathbf{v}}+\sum_{i}\frac{\partial
g}{\partial c_{i}}\frac{\partial
c_{i}^{2}}{\partial\mathbf{v}\partial\mathbf{w}}.$

Now that all terms can be evaluated in close form except for
$\frac{\partial\mathbf{t}}{\partial\mathbf{c}}$ and
$\frac{\partial\mathbf{t}}{\partial\mathbf{v}}$ which can be computed using
sensitivity analysis. We begin by differentiate both sides of (6) w.r.t.
$\mathbf{c}$,

(43)
$\frac{\partial\mathbf{x}}{\partial\mathbf{t}}^{\mathsf{T}}\left(\frac{\partial^{2}g}{\partial\mathbf{c}\partial\mathbf{x}}+\frac{\partial^{2}g}{\partial\mathbf{x}^{2}}\frac{\partial\mathbf{x}}{\partial\mathbf{t}}\frac{\partial\mathbf{t}}{\partial\mathbf{c}}\right)=\mathbf{0}.$

To obtain $\frac{\partial\mathbf{t}}{\partial\mathbf{c}}$ we rearrange the
terms and solve for

(44) $\mathbf{A}\frac{\partial\mathbf{t}}{\partial\mathbf{c}}=\mathbf{B},$

where

(45)
$\displaystyle\mathbf{A}=\frac{\partial\mathbf{x}}{\partial\mathbf{t}}^{\mathsf{T}}\frac{\partial^{2}g}{\partial\mathbf{x}^{2}}\frac{\partial\mathbf{x}}{\partial\mathbf{t}},$
$\displaystyle\mathbf{B}=-\frac{\partial\mathbf{x}}{\partial\mathbf{t}}^{\mathsf{T}}\frac{\partial^{2}g}{\partial\mathbf{c}\partial\mathbf{x}}.$

Differentiating both sides of Eqn. 35 w.r.t. $\mathbf{v}$, we arrive at

(46)
$\displaystyle\frac{\partial\mathbf{x}}{\partial\mathbf{t}}^{\mathsf{T}}\left(\frac{\partial^{2}g}{\partial\mathbf{c}\partial\mathbf{x}}\frac{\partial\mathbf{c}}{\partial\mathbf{v}}+\frac{\partial^{2}g}{\partial\mathbf{x}^{2}}\left(\frac{\partial\mathbf{x}}{\partial\mathbf{t}}\left(\frac{\partial\mathbf{t}}{\partial\mathbf{c}}\frac{\partial\mathbf{c}}{\partial\mathbf{v}}+\frac{\partial\mathbf{t}}{\partial\mathbf{v}}\right)+\frac{\partial\mathbf{x}}{\partial\mathbf{v}}\right)\right)$
$\displaystyle+\sum_{i}\frac{\partial g}{\partial
x_{i}}\frac{\partial^{2}x_{i}}{\partial\mathbf{t}\partial\mathbf{v}}\left(\frac{\partial\mathbf{t}}{\partial\mathbf{c}}\frac{\partial\mathbf{c}}{\partial\mathbf{v}}+\frac{\partial\mathbf{t}}{\partial\mathbf{v}}\right)=\mathbf{0}.$

To obtain $\frac{\partial\mathbf{t}}{\partial\mathbf{v}}$ we again rearrange
the terms and solve for

(47) $\mathbf{C}\frac{\partial\mathbf{t}}{\partial\mathbf{v}}=\mathbf{D},$

where

(48)
$\displaystyle\mathbf{C}=\frac{\partial\mathbf{x}}{\partial\mathbf{t}}^{\mathsf{T}}\frac{\partial^{2}g}{\partial\mathbf{x}^{2}}\frac{\partial\mathbf{x}}{\partial\mathbf{t}},$
$\displaystyle\mathbf{D}=-\frac{\partial\mathbf{x}}{\partial\mathbf{t}}^{\mathsf{T}}\frac{\partial^{2}g}{\partial\mathbf{c}\partial\mathbf{x}}\frac{\partial\mathbf{c}}{\partial\mathbf{v}}-\frac{\partial\mathbf{x}}{\partial\mathbf{t}}^{\mathsf{T}}\frac{\partial^{2}g}{\partial\mathbf{x}^{2}}\left(\frac{\partial\mathbf{x}}{\partial\mathbf{t}}\frac{\partial\mathbf{t}}{\partial\mathbf{c}}\frac{\partial\mathbf{c}}{\partial\mathbf{v}}+\frac{\partial\mathbf{x}}{\partial\mathbf{v}}\right)-\sum_{i}\frac{\partial
g}{\partial
x_{i}}\frac{\partial^{2}x_{i}}{\partial\mathbf{t}\partial\mathbf{v}}.$

Multiply (43) by
$(\frac{\partial\mathbf{t}}{\partial\mathbf{c}}\frac{\partial\mathbf{c}}{\partial\mathbf{w}})^{\mathsf{T}}$
on the left and $\frac{\partial\mathbf{c}}{\partial\mathbf{w}}$ on the right
we have

(49)
$\displaystyle\left(\frac{\partial\mathbf{x}}{\partial\mathbf{t}}\frac{\partial\mathbf{t}}{\partial\mathbf{c}}\frac{\partial\mathbf{c}}{\partial\mathbf{w}}\right)^{\mathsf{T}}\left(\frac{\partial^{2}g}{\partial\mathbf{c}\partial\mathbf{x}}\frac{\partial\mathbf{c}}{\partial\mathbf{w}}+\frac{\partial^{2}g}{\partial\mathbf{x}^{2}}\frac{\partial\mathbf{x}}{\partial\mathbf{t}}\frac{\partial\mathbf{t}}{\partial\mathbf{c}}\frac{\partial\mathbf{c}}{\partial\mathbf{w}}\right)$
$\displaystyle=\mathbf{0}.$

Therefore, as in the one-way coupling case,
$\frac{\partial^{2}g}{\partial\mathbf{w}^{2}}$ can be further simplified to

(50) $\displaystyle\frac{\partial^{2}g}{\partial\mathbf{w}^{2}}$
$\displaystyle=\frac{\partial\mathbf{c}}{\partial\mathbf{w}}^{\mathsf{T}}\frac{\partial^{2}g}{\partial\mathbf{c}^{2}}\frac{\partial\mathbf{c}}{\partial\mathbf{w}}+\frac{\partial\mathbf{c}}{\partial\mathbf{w}}^{\mathsf{T}}\frac{\partial^{2}g}{\partial\mathbf{c}\partial\mathbf{x}}\frac{\partial\mathbf{x}}{\partial\mathbf{t}}\frac{\partial\mathbf{t}}{\partial\mathbf{c}}\frac{\partial\mathbf{c}}{\partial\mathbf{w}}.$

##  Appendix B Convergence Criteria

Here we provide the details for computing the convergence criteria shown in
Fig. 5. Let us use $\frac{\mathrm{d}O}{\mathrm{d}\mathbf{w}}^{i}$ to define
the gradient of the energy shown in (16) w.r.t. the barycentric coordinate
$\mathbf{w}$ of the Karcher mean point in iteration $i$, and the convergence
criterion defined in Geometry Central is given by

(51)
$\left|\alpha\cdot\frac{\mathrm{d}O}{\mathrm{d}\mathbf{w}}\right|<\frac{1}{3}\sqrt{\mathcal{A}^{i}},$

where $\mathcal{A}$ is the area of the triangle that the Karcher mean point
lies on in iteration $i$, and $\alpha$ is the line search scaling factor to
ensure energy decrease.

Our absolute convergence criterion is defined as

(52) $\left|\frac{\mathrm{d}O}{\mathrm{d}\mathbf{w}}\right|<10^{-6}.$

As a reference of scale, for a given input geometry, we perform a uniform
scaling such that the diagonal of its bounding box is unit length. While we
set the tolerance to be $10^{-6}$, our approach often arrives at values
orders-of-magnitude smaller thanks to the quadratic convergence.

##  Appendix C Geodesic Voronoi Diagram

In this section, we demonstrate how to leverage sensitivity analysis to
compute the gradient of the objective function defined on the GVDs, _i.e_.,
(26) in Sec. 4.5. The gradient of this objective is given by

(53) $\frac{\mathrm{d}O}{\mathrm{d}\mathbf{s}}=\frac{\partial
O}{\partial\mathbf{s}}+\frac{\partial
O}{\partial\tilde{\mathbf{x}}}^{\mathsf{T}}\frac{\mathrm{d}\tilde{\mathbf{x}}}{\mathrm{d}\mathbf{s}}.$

Whereas $\frac{\partial O}{\partial\mathbf{s}}$ and $\frac{\partial
O}{\partial\tilde{\mathbf{x}}}$ can be computed algebraically,
$\frac{\mathrm{d}\tilde{\mathbf{x}}}{\mathrm{d}\mathbf{s}}$ requires the
following sensitivity analysis

(54)
$\displaystyle\frac{\mathrm{d}\mathbf{f}}{\mathrm{d}\mathbf{s}}=\mathbf{0}$
$\displaystyle=\frac{\partial\mathbf{g}}{\partial\mathbf{s}}+\frac{\partial\mathbf{f}}{\partial\tilde{\mathbf{x}}}^{\mathsf{T}}\frac{\mathrm{d}\tilde{\mathbf{x}}}{\mathrm{d}\mathbf{s}},$
$\displaystyle\frac{\mathrm{d}\tilde{\mathbf{x}}}{\mathrm{d}\mathbf{s}}$
$\displaystyle=-\frac{\partial\mathbf{f}}{\partial\tilde{\mathbf{x}}}^{-1}\frac{\partial\mathbf{f}}{\partial\mathbf{s}}.$

#### Pseudocode

Alg. 1 provides pseudocode for computing the exact geodesic Voronoi diagram.
Alg. 2 gives the general recipe we follow for minimizing objectives defined on
GVDs.

Data: triangle mesh $\mathcal{M}$, number of Voronoi sites $n$, Voronoi site
location $\mathbf{s}$

Result: Voronoi boundary vertex position $\tilde{\mathbf{x}}$

$\tilde{\mathbf{x}}^{*}$ $\leftarrow$
ComputeGVDApproximation($\mathcal{M},n,\mathbf{s}$) // (Xin et al., 2022)

parallel for  _$i\leftarrow 1$ to len $(\tilde{\mathbf{x}}^{*})$_ do

$\tilde{\mathbf{x}}_{i}$ $\leftarrow$
ComputeExactGVD($\mathcal{M},n,\mathbf{s},\tilde{\mathbf{x}}^{*}$) // Eqn.(25)

end for

ALGORITHM 1 Compute Exact Geodesic Voronoi Diagram

Data: triangle mesh $\mathcal{M}$, number of Voronoi sites $n$, Voronoi site
location $\mathbf{s}$

Result: Voronoi site location $\mathbf{s}^{*}$

$\mathbf{s}^{*}\leftarrow\mathbf{s}$, $\epsilon\leftarrow 10^{-6}$ //
initialize

while  _True_ do

$\tilde{\mathbf{x}}$ $\leftarrow$
ComputeExactGVD($\mathcal{M},n,\mathbf{s}^{*}$) // Alg.(1)

$\nabla O\leftarrow$
computeGrad($\mathcal{M},n,\mathbf{s}^{*},\tilde{\mathbf{x}}$) // Eqn.(53)

if  _$|\nabla O| <\epsilon$_ then

return $s^{*}$;

end if

$O_{0}\leftarrow$
computeObj($\mathcal{M},n,\mathbf{s}^{*},\tilde{\mathbf{x}}$) // Eqn.(31)

$\nabla^{2}O\leftarrow$
computeHessApprox($\mathcal{M},n,\mathbf{s}^{*},\tilde{\mathbf{x}}$)

$\triangle\mathbf{s}\leftarrow-(\nabla^{2}O)^{-1}\nabla O$ // linear solve

$\alpha\leftarrow 1$ // line search step size

while  _True_ do

$\mathbf{s}^{*}\leftarrow\mathbf{s}^{*}+\alpha\cdot\triangle\mathbf{s}$ //
Eqn.(15)

$\tilde{\mathbf{x}}$ $\leftarrow$
ComputeExactGVD($\mathcal{M},n,\mathbf{s}^{*}$) // Alg.(1)

$O_{1}\leftarrow$
computeObj($\mathcal{M},n,\mathbf{s}^{*},\tilde{\mathbf{x}}$) // Eqn.(31)

if  _$O_{1} <O_{0}$_ then

break;

end if

$\alpha\leftarrow 0.5\cdot\alpha$

end while

end while

ALGORITHM 2 Energy-minimizing Geodesic Voronoi Diagram