# Biharmonic Coordinates and their Derivatives for Triangular 3D Cages

Jean-Marc Thiery, Élie Michel, Jiong Chen

► To cite this version:

Jean-Marc Thiery, Élie Michel, Jiong Chen. Biharmonic Coordinates and their Derivatives for Triangular 3D Cages. ACM Transactions on Graphics, 2024, 43 (138), pp.17. 10.1145/3658208. hal-04812859

HAL Id: hal-04812859

https://hal.science/hal-04812859v1

Submitted on 1 Dec 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Biharmonic Coordinates and their Derivatives for Triangular 3D Cages

JEAN-MARC THIERY, Adobe Research, France

ÉLIE MICHEL, Adobe Research, France

JIONG CHEN, Inria, France

<figure><img src="image_2.png" /><figcaption>Fig. 1. While Mean-value coordinates interpolate the Dirichlet boundary condition exactly (i.e., the deformation is interpolatory on the cage surface) and allow artists to control the shape precisely near the cage, they lead to poor deformations, as the cage deformation is by nature continuous only across its edges. Being based on harmonic functions, that by definition can fit only one type of boundary condition at a time, Green coordinates follow only loosely the cage deformation, even on rather simple cage edits. We introduce biharmonic coordinates for triangular cages in 3D, that allow obtaining biharmonic 3D deformations that conform better to both Dirichlet (position) and Neumann (normal derivative) boundary conditions. By constraining the normal derivative stretch in this example, the deformed shape remains in the center of the deformed cage.</figcaption></figure>

As a natural extension to the harmonic coordinates, the biharmonic coordinates have been found superior for planar shape and image manipulation with an enriched deformation space. However, the 3D biharmonic coordinates and their derivatives have remained unexplored. In this work, we derive closed-form expressions for biharmonic coordinates and their derivatives for 3D triangular cages. The core of our derivation lies in computing the closed-form expressions for the integral of the Euclidean distance over a triangle *and* its derivatives. The derived 3D biharmonic coordinates not only fill a missing component in methods of generalized barycentric coordinates but also pave the way for various interesting applications in practice, including producing a family of biharmonic deformations, solving variational shape deformations, and even unlocking the closed-form expressions for recently-introduced Somigliana coordinates for both fast and accurate evaluations.

**CCS Concepts: • Computing methodologies → Shape modeling; Parametric curve and surface models; Volumetric models.**

**Additional Key Words and Phrases:** cage-based modeling, 3D shape deformation, biharmonic functions, biharmonic coordinates

**ACM Reference Format:**

Jean-Marc Thiery, Élie Michel, and Jiong Chen. 2024. Biharmonic Coordinates and their Derivatives for Triangular 3D Cages. *ACM Trans. Graph.* 43, 4, Article 138 (July 2024), 17 pages. https://doi.org/10.1145/3658208

## 1 INTRODUCTION

Cage deformers are a powerful tool for deforming complex shapes in computer graphics. They allow artists to manipulate a high-resolution mesh through a simpler enclosing cage, providing a simple and efficient approach to generating smooth, intuitive volumetric deformations. They are highly complementary of point-based (brushes) and skeleton-based deformers both in nature (points vs. regions vs. line segments vs. surfaces) and in deformation capabilities, and are now ubiquitous in modeling applications as well as in the animation industry.

Early works designed cage deformers based on generalized barycentric coordinates, such as the Mean Value Coordinates (MVC) [Floater 2003; Ju et al. 2005; Lipman et al. 2007] or Harmonic Coordinates (HC) [Joshi et al. 2007], which ensure exact interpolation of the deformation on the cage surface. One clear benefit of these approaches is that they provide precise and intuitive shape control to artists, as the deformed shape follows the cage tightly near its boundary.

However, the main issue observed with MVC is their non-positivity for non-convex cages: when two different parts of the cage are close in the rest pose (e.g., limbs of a character), and spread apart in the deformation, some cage coordinates become negative, which results in unnatural deformation behaviors (e.g., translating parts of the cage vertices in a given direction results in parts of the shape being moved in the opposite direction). Joshi et al. [2007] partly addressed this issue by computing non-negative barycentric coordinates that are a solution to the harmonic equation. Unfortunately, their coordinates have no closed-form expressions, requiring a global optimization inside the cage per coordinate, and their evaluation is as a result computationally expensive and error-prone. Lipman et al. [2007] used hardware accelerated rasterization to estimate Positive MVC,

Authors' Contact Information: Jean-Marc Thiery, Adobe Research, France, jthiery@adobe.com; Élie Michel, Adobe Research, France, emichel@adobe.com; Jiong Chen, Inria, France, jiong.chen@inria.fr.

© 2024 Copyright held by the owner/author(s). Publication rights licensed to ACM. This is the author's version of the work. It is posted here for your personal use. Not for redistribution. The definitive Version of Record was published in *ACM Transactions on Graphics*, https://doi.org/10.1145/3658208.
by simply computing the MVC w.r.t. the portion of the cage that is
visible from the point of view of the encoded mesh vertex (this point
being encoded w.r.t. to back-facing triangles only, all resulting coor-
dinates are positive). Their solution appears simpler to implement
and more robust than HC, but their coordinates do not come with a
closed-form expression, and result in continuous-only deformations
due to the use of the (discontinuous) visibility term.

While all those coordinates offer precise control near the cage,
since they ensure interpolation on its surface, they unfortunately
lead to visible artifacts as a direct result. This issue does not come
from the way those methods spread the deformation from the cage
to its interior, it rather comes from the nature of the cage defor-
mation itself: the cage surface being a polygonal mesh (and not a
smooth surface), the resulting deformations are continuous only
at the cage edges, and the cage's polygonal discretization leads
to visible deformation artifacts near the cage surface (see Fig. 1,
Mean-value coordinates).

The Green coordinates (GC) introduced by [Lipman et al. 2008]
instead strike for a better balance between the demands of shape-
preserving deformation and low-cost coordinate evaluation. Unlike
generalized barycentric coordinates like MVC, GC leverage both
cage vertices and cage face normals to control the interior's defor-
mation, and can thus infer local shape rotations from cage vertex
translations. Being based on Green's third identity to diffuse both the
Dirichlet (position) and Neumann (normal derivative) boundary con-
ditions, they yield conformal mappings in 2D and quasi-conformal
mappings in 3D. However, as they are harmonic functions – which
are uniquely defined by one boundary condition (most commonly,
Dirichlet), they can only loosely approximate the two combined
boundary conditions, and the deformed shapes may largely deviate
from the deformed cage, even for simple deformations (see Fig. 1).

The inability of harmonic functions to obey many boundary con-
ditions at once has motivated the use of higher-order basis func-
tions in many works on variational surface modeling and optimiza-
tion [Botsch and Kobbelt 2004; Jacobson et al. 2010; Stein et al.
2019], or skeleton-based deformations [Jacobson et al. 2011, 2012b].
Building on those observations, Weber et al. [2012] derived the
BiHarmonic Coordinates (BiHC) for 2D planar deformations as a
natural extension to GC. BiHC allow controlling not only the bound-
ary shapes but also boundary derivatives, enriching the deformation
space and enabling better alignment between the boundary cage and
its interior geometry. While the resulting coordinates have shown
clear advantages over MVC and GC in 2D, their 3D counterparts, to
the best of our knowledge, have remained unexplored so far.

We introduce in this paper this missing component in the family
of 3D cage deformers, and derive the closed-form expressions for
the biharmonic coordinates and their derivatives for triangular 3D
cages. We demonstrate that the foundation of our derivation lies in
computing the integral of the Euclidean distance function over a
3D triangle and its derivatives. Upon analytically evaluating these
values, we show that the 3D biharmonic coordinates, for both cage
faces and vertices, can be constructed by composing this distance
integral, its derivatives, and Green coordinates. This opens up av-
enues for applications such as direct cage-based modeling (where
an artist interacts with the cage vertices directly) as well as 3D

shape variational deformations (where the cage allows defining a
well-behaved subspace for non-linear shape optimization).

We make the following technical contributions:

* We derive closed-form expressions for the integrals of the biharmonic kernels on 3D triangles, from which we obtain a Boundary Element Basis of 3D biharmonic functions. Their gradient and Hessian are also derived;

* Based on the above, we devise 3D biharmonic coordinates for cage-based deformations. Inspired by Weber et al. [2012], we introduce a family of such biharmonic coordinates, which enables powerful boundary control;

* We investigate variational biharmonic deformations;

* Finally, we discover the relationship between the our coordi-
nates and the recently-introduced elasticity-based Somigliana
coordinates (SC) [Chen et al. 2023]. As a result, we obtain
closed-form expressions for SC, which significantly outper-
form numerical integration in both efficiency and accuracy.

1.1 Related work

*Generalized barycentric coordinates.* Generalized barycentric co-
ordinates were first introduced by Möbius [1827]. Since then, many
variants have been proposed and widely applied in graphics appli-
cations for real-time shape editing. The popular Mean-Value Co-
ordinates (MVC) were first proposed by Floater [2003] in 2D, and
later successively extended to 3D triangular cages [Ju et al. 2005],
3D cages made of planar n-gons [Langer et al. 2006], and triquad
cages in 3D [Thiery et al. 2018]. The MVC are derived from the
mean-value property of harmonic functions which can be computed
using closed-form expressions, and they are positive for convex
cages. However, MVC can have negative values for non-convex
cages, which may result in unintuitive deformations and undesir-
able artifacts. Both Positive MVC [Lipman et al. 2007] and Harmonic
Coordinates [Joshi et al. 2007] address the issue of non-positivity, but
these coordinates cannot be evaluated with closed-form expressions,
requiring expensive and approximate computations.

Maximum entropy coordinates [Hormann and Sukumar 2008]
(MEC) differ from most coordinates, as they use input prior mass
functions {$m_i$} guiding their definition, which is based on Infor-
mation Theory. MEC are the valid barycentric coordinates {$\lambda_i(\eta)$}
(i.e., for {$v_i$} the cage vertices and $\eta \in \mathbb{R}^3$ the evaluation point:
$\sum_i \lambda_i(\eta)(v_i; 1) = (\eta; 1)$) best-fitting {$m_i$} in that they maximize the
corresponding cross-entropy: $H(\{\lambda_i\}, \{m_i\}) := - \sum_i \lambda_i \log(\lambda_i/m_i)$.
Setting positive {$m_i$} ensures positive {$\lambda_i$} ($\lambda_i$ has the same sign
as $m_i$ by design, as $\lambda_i \log(\lambda_i/m_i)$ is otherwise ill-defined). Their
evaluation requires however a non-linear optimization scheme and
no solution may exist if {$m_i$} are defined carelessly (e.g., providing
positive {$m_i$} for points $\eta$ lying outside the convex hull of {$v_i$} will
fail). The geometric properties of MEC heavily depend on {$m_i$}, and
the regularizing property of MEC has made them a useful tool in
several works [Corda et al. 2020; Thiery et al. 2018]. Linked to Infor-
mation Theory as well, Maximum Likelihood Coordinates [Chang
et al. 2023] define coordinates as the ones maximizing the likelihood
$\mathcal{L}(\{\lambda_i\}) := \sum_i \log(\lambda_i)$, eliminating the priors from their definition.

Though difficult to define rigorously, *locality* is a property that
may be important in specific scenarios. Inspired by earlier works
in the Signal Processing community, Local Barycentric Coordinates [Zhang et al. 2014] (LBC) were defined as the minimizer of the total variation. Variational Barycentric Coordinates [Dodik et al. 2023] (VBC) later extended LBC, as they allow for minimizing various objective energies (which include the coordinates' total variation, but also their Dirichlet energy – minimized by harmonic functions [2007]) using a novel dedicated neural network architecture. Like MEC and MLC, LBC and VBC require an optimization process, and do not admit closed-form expressions.

**Cage coordinates with boundary derivative control.** Cage coordinates offering control over not only boundary values but also boundary derivatives (i.e., normals) have been introduced more recently. Based on Green's third identity for harmonic functions, Green Coordinates were introduced by Lipman et al. [2008], providing smooth and shape preserving deformations that achieve conformal mapping in 2D and quasi-conformal mappings in 3D. While those coordinates are not guaranteed to be positive, their use of the cage face normals make the deformations cage-aware in practice, and GC provide more intuitive deformations than MVC for cages featuring close limbs meant to be separated. Since then, GC have been extended several times. Weber et al. [2009] have shown that 2D GC are a special variant of their 2D complex-valued barycentric coordinates, which are derived from the Cauchy integral formula for holomorphic functions. GC for triquad cages [Thiery and Boubekeur 2022] in 3D and its polynomial extension in 2D [Michel and Thiery 2023] have also been investigated, further broadening its potential applications. Weber et al. [2012] demonstrated that, by using biharmonic deformations that expose an explicit control over boundary derivatives, the interior deformation can better respect the boundary cage. This key property enables producing intuitive smooth deformations, but unfortunately authors provided BiHC only for the 2D case, and 3D biharmonic cage-based deformations have remained unexplored.

Recently, Somigliana coordinates [Chen et al. 2023] were derived from the boundary integral formulation of linear elasticity, similarly to how GC had been derived from Green's third identity for harmonic functions. SC enable volume control of cage deformation, addressing potential undesirable volume scaling produced by GC. However, the authors did not derive closed-form expressions at the time; instead, they rely on quadrature rules to compute their coordinates, which can be inefficient and inaccurate in practice.

**Variational surface deformation.** Using cage-based deformations as a compact smooth deformation subspace enables the direct As-Rigid-As-Possible (ARAP) manipulation of a shape without requiring an explicit discretization of the domain [Ben-Chen et al. 2009] or deformation transfer [Chen et al. 2010], as the cost for minimization is much lower on the coarse cage than directly on the high-resolution deformed geometry. Using cage coordinates in this context amounts to using a well-behaved subspace for optimization. Designing subspaces that allow for rich deformations while minimizing deformation artifacts becomes critical, which has been the topic of interest of many research works [Brandt et al. 2018; Jacobson et al. 2012a; von Radziewsky et al. 2016; Von-Tycowicz et al. 2015; Wang et al. 2015]. For 2D cage-based deformations, [Weber et al. 2012] have demonstrated that biharmonic coordinates enrich harmonic functions in a way that they offer additional flexibility

and deformation capabilities for variational deformations, while remaining smooth and well-behaved. We march in their steps, and investigate the use of our 3D biharmonic coordinates in this context as well in Section 4.2.

## 2 BIHARMONIC COORDINATES FOR TRIANGULAR CAGES IN 3D

We note $G_k(\xi, \eta)$ the fundamental solution to the $k$-harmonic equation in $\mathbb{R}^3$, i.e., $\Delta_{\xi}^k G_k(\xi, \eta) = \delta_0(\|\xi - \eta\|)$, $G_1(\xi, \eta) := G(\xi, \eta)$.

We can express a biharmonic function $f$ in a bounded 3D domain $\Omega$ from its boundary conditions as

$$f(\eta) = \int_{\xi \in \partial\Omega} f(\xi) \frac{\partial_1 G_1}{\partial n}(\xi, \eta) d\xi - \int_{\xi \in \partial\Omega} G_1(\xi, \eta) \frac{\partial f}{\partial n}(\xi) d\xi \quad (1) \\ + \int_{\xi \in \partial\Omega} \Delta f(\xi) \frac{\partial_1 G_2}{\partial n}(\xi, \eta) d\xi - \int_{\xi \in \partial\Omega} G_2(\xi, \eta) \frac{\partial \Delta f}{\partial n}(\xi) d\xi,$$

where $n$ is the normal of the cage at $\xi$. More generally, any $k$-harmonic function can be expressed as

$$f(\eta) = \sum_{j=1}^{k} \int_{\xi \in \partial\Omega} \Delta^{j-1} f(\xi) \frac{\partial_1 G_j}{\partial n}(\xi, \eta) d\xi - \int_{\xi \in \partial\Omega} G_j(\xi, \eta) \frac{\partial \Delta^{j-1} f}{\partial n}(\xi) d\xi,$$

where $\Delta^0 f$ is simply $f$. It is easy to see that $G_k$ takes the following form: $G_k(\xi, \eta) = \lambda_k \|\xi - \eta\|^{2k-3}$, and by recurrence, we obtain

$$\begin{aligned} \Delta(G_{k+1}(\xi, \eta)) &= \lambda_{k+1} \Delta (\|\xi - \eta\|^{2k-1}) \\ &= \lambda_{k+1} (2k)(2k-1) \|\xi - \eta\|^{2k-3} \\ &= G_k(\xi, \eta) = \lambda_k \|\xi - \eta\|^{2k-3}, \end{aligned}$$

leading to $\lambda_{k+1} = \lambda_k / ((2k)(2k-1))$ and $\lambda_{k+1} = \frac{\lambda_1}{(2k)!}$. For example, $\lambda_1 = -1/(4\pi)$, $\lambda_2 = -1/(8\pi)$, and

$$\begin{aligned} G_1(\xi, \eta) &= \frac{-1}{4\pi \|\xi - \eta\|}, \\ G_2(\xi, \eta) &= \frac{-\|\xi - \eta\|}{8\pi}. \end{aligned} \quad (2)$$

### 2.1 Boundary-driven biharmonic functions

A biharmonic function ($\Delta^2 = 0$) inside a bounded domain $\Omega \subset \mathbb{R}^3$ is uniquely defined by the combination of the Dirichlet and Neumann boundary conditions. Deriving biharmonic functions using

$$f(\eta) = \int_{\xi \in \partial\Omega} f^a(\xi) \frac{\partial_1 G_1}{\partial n}(\xi, \eta) d\xi - \int_{\xi \in \partial\Omega} G_1(\xi, \eta) f^b(\xi) d\xi \quad (3) \\ + \int_{\xi \in \partial\Omega} f^c(\xi) \frac{\partial_1 G_2}{\partial n}(\xi, \eta) d\xi - \int_{\xi \in \partial\Omega} G_2(\xi, \eta) f^d(\xi) d\xi$$

results in a function $f$ that effectively fits the four (redundant) boundary conditions ($f = f^a, \partial f/\partial n = f^b, \Delta f = f^c, \partial \Delta f/\partial n = f^d$ on the boundary $\partial\Omega$ if and only if they are compatible with each other (i.e., there effectively exists a biharmonic function $f$ that fits the four conditions). In general, setting four arbitrary decorrelated conditions will result in a function $f$ that fits none exactly, but that follows them loosely only. However, any function $f$ computed from Eq. (3) is
biharmonic, as $f$ is expressed as a convolution with the biharmonic
kernels (one can differentiate Eq. (3) w.r.t. $\eta$ beneath the integral sign,
following Lebesgue's theorem). In the following, we will use Eq. (3)
to derive biharmonic coordinates for spatial deformation. The re-
sulting deformations $f$ will fit none of the four redundant boundary
conditions we will set up. Fitting those is however not our objec-
tive: we focus on providing efficient and high-quality biharmonic
deformations. Note that a similar observation was already made
in [Thiery and Boubekeur 2022] for GC, and that the same exact
paradigm was used as well for 2D biharmonic deformations [Weber
et al. 2012] or more recently Somigliana coordinates [Chen et al.
2023] and extended 2D GC [Michel and Thiery 2023].

## 2.2 Boundary finite element (BFE) basis functions

We note $F_1(i)$ the set of faces adjacent to vertex $i$ and $\Gamma_i$ the piece-
wise linear "hat function" that takes value 1 on vertex $i$ and 0
on the others. By setting the boundary conditions (see Eq. (3))
$f^a \approx f, f^c \approx \Delta f$ to be piecewise affine (affine on each trian-
gle), and $\{f^b \approx \frac{\partial f}{\partial n}, f^d \approx \frac{\partial \Delta f}{\partial n}\}$ to be piecewise constant (constant
within each triangle), we obtain a boundary finite elements basis to
represent biharmonic functions $f: \Omega \subset \mathbb{R}^3 \mapsto \mathbb{R}^n$ as

$$
f(\eta) = \sum_{i \in \mathcal{V}} \phi_i(\eta) a_i + \sum_{j \in \mathcal{T}} \psi_j(\eta) b_j + \sum_{i \in \mathcal{V}} \bar{\phi}_i(\eta) c_i + \sum_{j \in \mathcal{T}} \bar{\psi}_j(\eta) d_j, \quad (4)
$$

with $(a_i, b_j, c_i, d_j) \in \mathbb{R}^n$ and

$$
\phi_i(\eta) := \int_{\xi \in F_1(i)} \Gamma_i(\xi) \frac{\partial_1 G_1}{\partial n_{\xi}}(\xi, \eta) d\xi, \quad (5)
$$

$$
\psi_j(\eta) := \int_{\xi \in t_j} -G_1(\xi, \eta)d\xi,
$$

$$
\bar{\phi}_i(\eta) := \int_{\xi \in F_1(i)} \Gamma_i(\xi) \frac{\partial_1 G_2}{\partial n_{\xi}}(\xi, \eta) d\xi,
$$

$$
\bar{\psi}_j(\eta) := \int_{\xi \in t_j} -G_2(\xi, \eta)d\xi.
$$

We note $\bar{\phi}_i^t(\eta)$ (resp. $\phi_i^t(\eta)$) the contribution of triangle $t$ to $\bar{\phi}_i(\eta)$ (resp. $\phi_i(\eta)$), for $t \in F_1(i)$. We obtain trivially that

$$
\phi_i(\eta) = \sum_{t \in F_1(i)} \phi_i^t(\eta), \quad (6)
$$

$$
\bar{\phi}_i(\eta) = \sum_{t \in F_1(i)} \bar{\phi}_i^t(\eta).
$$

Biharmonic coordinates. The functions {$\phi_i(\eta), \psi_j(\eta), \bar{\phi}_i(\eta), \bar{\psi}_j(\eta)$} are generalized barycentric coordinates of $\eta$, as Eq. (4) allows reproducing the identity mapping, with $a_i$ set to the position $v_i$ of the vertex $i$, $b_j$ set to the normal $n_j$ of triangle $j$, and $c_i = d_j = 0$:

$$
\eta = \sum_{i \in \mathcal{V}} \phi_i(\eta) v_i + \sum_{j \in \mathcal{T}} \psi_j(\eta) n_j, \quad (7)
$$

leading to the known Green coordinates [Lipman et al. 2008].

# 3 COMPUTATIONS OF BFE BASIS FUNCTIONS

As a disclaimer, this section presents mathematical expressions that are relatively involved and hard to digest, but not needed for the reader mostly interested in reviewing applications of biharmonic functions, who may go to Section 4 directly. For readability, we only present the expressions here, and defer the mathematical proofs to the appendix.

Notations for affine basis functions in-
side triangles. We note in the following
$t = (t^0, t^1, t^2)$ a triangle which we con-
sider oriented counterclock-wise, with
normal $n_t := N_t / \|N_t\|$, $N_t := (t^1 - t^0) \times (t^2 - t^0)$ (see inset).

<figure><img src="image_3.png" /></figure>

For a point $\xi \in t$, its barycentric coordinate $\Gamma_{ti}^t(\xi)$ with respect to vertex $t^i$ (restricted to $t$) is: $\Gamma_{ti}^t(\xi) := \frac{(\xi - t^{i+1}) \cdot r^i}{h_i}$, all indices being taken modulo 3, with $r^i := n_t \times (t^{i+2} - t^{i+1}) / \|t^{i+2} - t^{i+1}\|$ the unit vector orthogonal to the edge opposite to $v_i$ inside $t$, and $h_i = (t^i - t^{i+1}) \cdot r^i$ being the distance from $t^i$ to its opposite edge. For simplicity, we note $g_i^t$ the (constant) gradient of the barycentric coordinate function $\Gamma_{ti}^t(\xi)$, i.e., $g_i^t := \nabla_{\xi} \Gamma_{ti}^t(\xi) = \frac{r^i}{h_i}$. Using this notation, one can write $\Gamma_{ti}^t(\xi) = (\xi - t^{i+1}) \cdot g_i^t$.

We also note $d_t := (\eta - t^0) \cdot n_t$ the signed distance of $\eta$ to the tangent plane of the triangle $t$.

Green coordinates – Urago’s formulation. For compactness, we refer to the annex of [Ben-Chen et al. 2009] for the formulas of the {$\phi, \psi$} harmonic coordinates and their gradients and Hessians.

Relating $\bar{\phi}_i^t(\eta)$ to $\bar{\psi}_t(\eta)$. Looking at $\bar{\phi}_{tk}^t(\eta)$ (the contribution of triangle $t$ in $\bar{\phi}_{tk}(\eta)$), we observe that:

$$
\begin{align*}
\bar{\phi}_{tk}^t(\eta) &:= \int_{\xi \in t} \Gamma_{tk}^t(\xi) \frac{\partial G_2}{\partial n}(\xi, \eta) d\xi = \lambda_2 \int_{\xi \in t} \Gamma_{tk}^t(\xi) \frac{(\xi - \eta) \cdot n_t}{\|\xi - \eta\|} d\xi \\
&= -\lambda_2 d_t \int_{\xi \in t} \frac{\Gamma_{tk}^t(\xi)}{\|\xi - \eta\|} d\xi = -\lambda_2 d_t \int_{\xi \in t} \frac{(\xi - t^{k+1}) \cdot g_k^t}{\|\xi - \eta\|} d\xi \\
&= -\lambda_2 d_t \int_{\xi \in t} \frac{(\xi - \eta) \cdot g_k^t}{\|\xi - \eta\|} d\xi - \lambda_2 d_t \int_{\xi \in t} \frac{(\eta - t^{k+1}) \cdot g_k^t}{\|\xi - \eta\|} d\xi \\
&= \lambda_2 d_t \nabla_{\eta} \left( \int_{\xi \in t} \|\xi - \eta\| d\xi \right) \cdot g_k^t - \lambda_2 d_t \int_{\xi \in t} \frac{(\eta - t^{k+1}) \cdot g_k^t}{\|\xi - \eta\|} d\xi \\
&= -d_t(\eta) \nabla_{\eta}(\bar{\psi}_t) \cdot g_k^t + \frac{\lambda_2}{\lambda_1} d_t \Gamma_{tk}^t(\eta) \int_{\xi \in t} \frac{-\lambda_1}{\|\xi - \eta\|} d\xi \\
&= -d_t \nabla_{\eta}(\bar{\psi}_t) \cdot g_k^t + \frac{1}{2} d_t \Gamma_{tk}^t(\eta) \psi_t(\eta).
\end{align*}
$$

We can therefore relate the function $\bar{\phi}_{tk}^t(\eta)$ to the gradient of $\bar{\psi}_t(\eta)$ as well as the harmonic basis function $\bar{\psi}_t(\eta)$. We already have closed-form expressions for the latter, as well as for its derivatives, and we now propose a closed-form expression for $\bar{\psi}_t(\eta)$ as it is thus key to deriving all other coordinates and their derivatives.
## 3.1 Closed-form expression of $\tilde{\psi}_t(\eta)$ and its derivatives

We provide our derivation in Appendix A. We denote $\omega_t(\eta)$ the signed solid angle of the triangle $t$ seen from $\eta$, and we consider the three consistently oriented edges $e$ of $t$, noting for each edge $e_0$ (resp. $e_1$) its origin (resp. destination). In the remainder of this section, we will note many expressions (Eqs. (10),(12), (14)) of the form

$$ \int_{\xi \in t} fd\xi = A(t) + \sum_{e \in t} B(e), \qquad (8) $$

$B(e)$ being a per-edge $e$ quantity. Some of the terms in $B(e)$ are ill-defined for $\eta$ perfectly aligned with the edge $e$ ($\eta \in (e_0e_1)$), but all those converge to 0 in this case (e.g., $\lim_{x \to 0^+} x \log(x) = 0$). We color in blue the terms that remain in $B(e)$ for $\eta \in (e_0e_1)$; all other terms of $B(e)$ are safely replaced by 0 in this case.

Using these notations, $\tilde{\psi}_t(\eta)$ is given by:

$$ \tilde{\psi}_t(\eta) := -\lambda_2 \int_{\xi \in t} \|\xi - \eta\| d\xi, \qquad (9) $$

with

$$ \int_{\xi \in t} \|\xi - \eta\| d\xi = \frac{d_t^3 \omega_t(\eta)}{3} \qquad (10) \\ + \sum_{e \in t} \frac{a_e^t}{6} \left( (2d_t^2 + D_e^2) \log\left(\frac{l_{e1} - \zeta_{e1}}{l_{e0} - \zeta_{e0}}\right) - l_{e1}\zeta_{e1} + l_{e0}\zeta_{e0} \right), $$

and

$$ \begin{aligned} d_t &= (\eta - e_0) \cdot n_t \\ u_e &:= \vec{e} / \|\vec{e}\|; \quad r_e^t := n_t \times u_e \\ \vec{l}_{e_i} &:= \eta - e_i; \quad l_{e_i} = \|\vec{l}_{e_i}\|; \quad \zeta_{e_i} = \vec{l}_{e_i} \cdot u_e \\ a_e^t &= r_e^t \cdot \vec{l}_{e_k}, \quad k = 0 \text{ or } 1 \\ D_e^2 &= \|(I_3 - u_e \cdot u_e^T) \cdot \vec{l}_{e_k}\|^2 =: \|P_e \cdot \vec{l}_{e_k}\|^2 =: \|\vec{D}_e\|^2, \quad k = 0 \text{ or } 1 \end{aligned} $$

where $I_3$ is the identity matrix of $3 \times 3$.

One can already note that, in case $\eta \in (e_0e_1)$, $a_e^t = D_e = d_t = 0$, which we use to simplify our expressions in this case, as we extend by continuation $x \log(x)$ to 0 at $x = 0$.

By direct differentiation of Eq. (10), we obtain its gradient as:

$$ \nabla_{\eta}(\tilde{\psi}_{t}) := -\lambda_{2}\nabla_{\eta}\left(\int_{\xi \in t} \|\xi - \eta\| d\xi\right), \qquad (11) $$

with

$$ \begin{aligned} \nabla_{\eta} \left( \int_{\xi \in t} \|\xi - \eta\| d\xi \right) = & d_t^2 n_t \omega_t + \frac{d_t^3}{3} \nabla_{\eta}(\omega_t) \qquad (12) \\ & + \sum_{e \in t} \frac{r_e^t}{6} \left( (2d_t^2 + D_e^2) \log\left(\frac{l_{e1} - \zeta_{e1}}{l_{e0} - \zeta_{e0}}\right) + l_{e0}\zeta_{e0} - l_{e1}\zeta_{e1} \right) \\ & + \sum_{e \in t} \frac{a_e^t}{6} \left( (4d_t n_t + 2\vec{D}_e) \log\left(\frac{l_{e1} - \zeta_{e1}}{l_{e0} - \zeta_{e0}}\right) \right) \\ & + \sum_{e \in t} \frac{a_e^t}{6} \left( (2d_t^2 + D_e^2) [m_{e1} - m_{e0}] \right) \\ & + \sum_{e \in t} \frac{a_e^t}{6} \left( -u_{e1}\zeta_{e1} + u_{e0}\zeta_{e0} + (l_{e0} - l_{e1})u_e \right), \end{aligned} $$

and

$$ \left\{ \begin{aligned} u_{e_i} &:= \nabla_{\eta}(l_{e_i}) = \frac{\eta - e_i}{\|\eta - e_i\|} = \frac{\vec{l}_{e_i}}{l_{e_i}} \\ m_{e_i} &:= \frac{u_{e_i} - u_e}{l_{e_i} - \zeta_{e_i}} \\ \nabla_{\eta}(\omega_t) &:= \sum_{e \in t} \frac{2(l_{e0} + l_{e1})}{(l_{e0} + l_{e1})^2 - \|e\|^2} u_{e0} \times u_{e1} \end{aligned} \right. . $$

By differentiating Eq. (12) again, and denoting $\mathcal{T}(u; v) := u \cdot \vec{v}^T + v \cdot u^T$ the (symmetric) tensor matrix obtained from the 3D vectors $u$ and $v$, we obtain its Hessian matrix as

$$ H_{\eta}(\tilde{\psi}_{t})(\eta) := -\lambda_{2}H_{\eta}\left(\int_{\xi \in t} \|\xi - \eta\| d\xi\right), \qquad (13) $$

with

$$ \begin{aligned} H_{\eta}\left(\int_{\xi \in t} \|\xi - \eta\| d\xi\right) = & d_t \omega_t \mathcal{T}(n_t; n_t) + d_t^2 \mathcal{T}(n_t; \nabla_{\eta}(\omega_t)) + \frac{d_t^3}{3} H_{\eta}(\omega_t) \\ & + \sum_{e \in t} \frac{1}{6} \log \left( \frac{l_{e1} - \zeta_{e1}}{l_{e0} - \zeta_{e0}} \right) \mathcal{T}(r_e^t; 4d_t n_t + 2\vec{D}_e) \\ & + \sum_{e \in t} \frac{1}{6} \mathcal{T}(r_e^t; u_{e0}\zeta_{e0} - u_{e1}\zeta_{e1} + (l_{e0} - l_{e1})u_e) \\ & + \sum_{e \in t} \frac{a_e^t}{6} \log \left( \frac{l_{e1} - \zeta_{e1}}{l_{e0} - \zeta_{e0}} \right) (4n_t n_t^T + 2P_e) \\ & + \sum_{e \in t} \frac{1}{6} \mathcal{T}(6a_e^t d_t n_t + 3D_e^2 r_e^t; m_{e1} - m_{e0}) \\ & + \sum_{e \in t} \frac{a_e^t}{6} (\zeta_{e0} U_{e0} - \zeta_{e1} U_{e1} + \mathcal{T}(u_{e0} - u_{e1}; u_e)) \\ & + \sum_{e \in t} \frac{a_e^t}{6} (2d_t^2 + D_e^2) (M_{e1} - M_{e0}), \end{aligned} \qquad (14) $$
and

$$
\begin{align*}
U_{e_i} &:= \nabla_{\eta}(u_{e_i}) = \frac{l_{e_i}^2 I_3 - \vec{l}_{e_i} \cdot \vec{l}_{e_i}^T}{l_{e_i}^3} \\
M_{e_i} &:= \nabla_{\eta}(m_{e_i}) = \frac{U_{e_i}}{l_{e_i} - \zeta_{e_i}} - m_{e_i} \cdot m_{e_i}^T \\
H_{\eta}(\omega_t) &= \sum_{e \in t} \frac{(l_{e_0} + l_{e_1}) \mathcal{T}(u_{e_1} \times u_{e_0}; l_{e_0} u_{e_1} + l_{e_1} u_{e_0})}{[(l_{e_0} + l_{e_1})^2 - \|e\|^2] l_{e_0} l_{e_1}} \\
&\quad + \sum_{e \in t} \left( \frac{1/2}{(l_{e_0} + l_{e_1} - \|e\|)^2} + \frac{1/2}{(l_{e_0} + l_{e_1} + \|e\|)^2} \right) \mathcal{T}(u_{e_1} \times u_{e_0}; u_{e_1} + u_{e_0}).
\end{align*}
$$

3.2 Closed-form expression of $\bar{\phi}_{t_i}^t(\eta)$ and its derivatives

Following our earlier remark linking both types of coordinates,
$\bar{\phi}_{i}^{t}(\eta)$ can then be obtained as

$$
\bar{\phi}_{t_i}^t(\eta) = -dt \nabla_{\eta}(\bar{\psi}_t) \cdot \mathbf{g}_i^t + \frac{1}{2} dt \Gamma_{t_i}^t(\eta) \psi_t(\eta). \quad (15)
$$

Its gradient can thus be computed as

$$
\nabla_{\eta} \bar{\phi}_{t_i}^t (\eta) = - (\nabla_{\eta}(\bar{\psi}_t) \cdot \mathbf{g}_i^t) n_t - dt H_{\eta}(\bar{\psi}_t) \mathbf{g}_i^t \\
+ \frac{1}{2} \Gamma_{t_i}^t (\eta) \psi_t (\eta) n_t + \frac{1}{2} dt \psi_t (\eta) \mathbf{g}_i^t + \frac{1}{2} dt \Gamma_{t_i}^t (\eta) \nabla_{\eta} \psi_t (\eta). \quad (16)
$$

Finally, its Hessian matrix can be computed as

$$
\begin{align}
H_{\eta} \bar{\phi}_{t_i}^t (\eta) = & - (n_t \mathbf{g}_i^T \cdot H_{\eta}(\bar{\psi}_t) + H_{\eta}(\bar{\psi}_t) \mathbf{g}_i^T n_t^T) \tag{17} \\
& - dt \nabla_{\eta}^T (H_{\eta}(\bar{\psi}_t) \mathbf{g}_i^T) \nonumber \\
& + \frac{1}{2} \psi_t(\eta) \mathcal{T} (n_t; \mathbf{g}_i^T) + \frac{1}{2} dt \mathcal{T} (\nabla_{\eta}(\psi_t); \mathbf{g}_i^T) \nonumber \\
& + \frac{1}{2} \Gamma_{t_i}^T(\eta) \mathcal{T} (\nabla_{\eta}(\psi_t); n_t) + \frac{1}{2} dt \Gamma_{t_i}^T(\eta) H_{\eta}(\psi_t). \nonumber
\end{align}
$$

All the necessary terms but $\nabla_{\eta}^T (H_{\eta}(\bar{\psi}_t) \mathbf{g}_i^T)$ can be computed using the expressions provided in Section 3.1.

We sketch in Appendix B the derivation of $\nabla_{\eta}^{T}(H_{\eta}(\bar{\psi}_{t}) \mathbf{g}_{i}^{T})$, and give here its formula, concluding the exposition of our closed-form expressions.

We note that for 3D vectors $p, q, r$ ($r$ not depending on $\eta$):

$$
\nabla_{\eta}^{T} (\mathcal{T}(p; q) \cdot r) = (q^{T} \cdot r I_{3} + q \cdot r^{T}) J p + (p^{T} \cdot r I_{3} + p \cdot r^{T}) J q.
$$

Using $\mathcal{T}(M; L) := M \cdot L^T + L \cdot M^T$ for matrices as well, we obtain:

$$
\nabla_{\eta}^{T} (H_{\eta}(\bar{\psi}_{t}) \mathbf{g}_{i}^{t}) = -\lambda_{2} \nabla_{\eta}^{T} \left( H_{\eta} \left( \int_{\xi \in t} \|\xi - \eta\| d\xi \right) \mathbf{g}_{i}^{t} \right), \quad (18)
$$

with

$$
\begin{align*}
\nabla_{\eta}^{T} \left( H_{\eta} \left( \int_{\xi \in t} \|\xi - \eta\| d\xi \right) r \right) &= 4d_{t}^{2}An_{t}n_{t}^{T} + \frac{4d_{t}^{3}}{3} \mathcal{T}(\nabla_{\eta} A; n_{t}) + \frac{d_{t}^{4}}{3} H_{\eta}(A) \\
&\quad + \sum_{e \in t} \frac{1}{6} \mathcal{T}(r_{e}^{T} r^{T}; \mathcal{T}(u_{e0} - u_{e1}; u_{e}) + \zeta_{e0} U_{e0} - \zeta_{e1} U_{e1}) \\
&\quad + \sum_{e \in t} \frac{D_{e}^{2}}{2} \mathcal{T}(r_{e}^{T} r^{T}; M_{e1} - M_{e0}) + \sum_{e \in t} a_{e}^{T} dt \mathcal{T}(n_{t} r^{T}; M_{e1} - M_{e0}) \\
&\quad + \sum_{e \in t} \mathcal{T}(\vec{D}_{e} r^{T}; (m_{e1} - m_{e0}) r_{e}^{T}^{T}) + \sum_{e \in t} dt ((m_{e1} - m_{e0}) \cdot r) \mathcal{T}(n_{t}, r_{e}^{T}) \\
&\quad + \sum_{e \in t} \frac{a_{e}^{T}}{6} \mathcal{T}(u_{e} r^{T}; U_{e0} - U_{e1}) + \sum_{e \in t} a_{e}^{T} ((m_{e1} - m_{e0}) \cdot r) (r_{e}^{T} r_{e}^{T} + n_{t} n_{t}^{T}) \\
&\quad + \sum_{e \in t} \log \left( \frac{l_{e1} - \zeta_{e1}}{l_{e0} - \zeta_{e0}} \right) (r_{e}^{T} \cdot r) P_{e} + \sum_{e \in t} \frac{D_{e}^{2}}{2} (r_{e}^{T} \cdot r) (M_{e1} - M_{e0}) \\
&\quad + \sum_{e \in t} \frac{r_{e}^{T} \cdot r}{6} (\mathcal{T}(u_{e0} - u_{e1}; u_{e}) + \zeta_{e0} U_{e0} - \zeta_{e1} U_{e1}) + \sum_{e \in t} \frac{a_{e}^{T}(u_{e} \cdot r)}{6} (U_{e0} - U_{e1}) \\
&\quad + \sum_{e \in t} \frac{a_{e}^{T}(2d_{t}^{2} + D_{e}^{2})}{6} (M_{e1}^{r} - M_{e0}^{r}) + \sum_{e \in t} \frac{a_{e}^{T}}{6} (\zeta_{e0} u_{e0}^{r} - \zeta_{e1} u_{e1}^{r})
\end{align*}
$$

for any $r \in \mathbb{R}^3$ tangent to $t$ (i.e., such that $r \cdot n_t = 0$), with

$$
\mathcal{U}_{e_i}^r := \nabla_{\eta}^T (U_{e_i} r) = \frac{-\mathcal{T}(r; u_{e_i})}{l_{e_i}^2} - \frac{(u_{e_i}^T \cdot r) U_{e_i}}{l_{e_i}} + 2 \frac{(u_{e_i}^T \cdot r) u_{e_i} u_{e_i}^T}{l_{e_i}^2},
$$

$$
\mathcal{M}_{e_i}^r := \nabla_{\eta}^T (M_{e_i} r) = \frac{\mathcal{U}_{e_i}(r)}{l_{e_i} - \zeta_{e_i}} - \frac{\mathcal{T}((u_{e_i} - u_e)r^T; U_{e_i})}{(l_{e_i} - \zeta_{e_i})^2} \\
- \frac{(u_{e_i} - u_e)^T \cdot r U_{e_i}}{(l_{e_i} - \zeta_{e_i})^2} + \frac{(u_{e_i} - u_e)^T \cdot r \mathcal{T}(u_{e_i} - u_e; u_{e_i} - u_e)}{(l_{e_i} - \zeta_{e_i})^3},
$$

and

$$
\begin{align*}
A &:= \sum_{e \in t} \frac{-2 \|e\| (l_{e0} + l_{e1}) r_e^{T} r}{\lambda_e}, \\
\lambda_e &:= ((l_{e0} + l_{e1})^2 - \|e\|^2) l_{e0} l_{e1}, \\
\nabla_{\eta} A &:= \sum_{e \in t} \frac{-2 \|e\| r_e^{T} \cdot r (u_{e0} + u_{e1})}{\lambda_e} + \frac{2 \|e\| (l_{e0} + l_{e1}) r_e^{T} \cdot r}{\lambda_e^2} \Lambda_e, \\
\Lambda_e &:= \nabla_{\eta} \lambda_e = 2(l_{e0} + l_{e1}) l_{e0} l_{e1} (u_{e0} + u_{e1}) \\
&\quad + ((l_{e0} + l_{e1})^2 - \|e\|^2) (l_{e1} u_{e0} + l_{e0} u_{e1}),
\end{align*}
$$

$$
\begin{align*}
H_{\eta} \lambda_e ={}& 2 l_{e0} l_{e1} (u_{e0} + u_{e1}) (u_{e0} + u_{e1})^T \\
& + 2(l_{e0} + l_{e1}) \mathcal{T}(u_{e0} + u_{e1}; l_{e1} u_{e0} + l_{e0} u_{e1}) \\
& + ((l_{e0} + l_{e1})^2 - \|e\|^2)(l_{e1} U_{e0} + l_{e0} U_{e1} + \mathcal{T}(u_{e0}; u_{e1})) \\
& + 2(l_{e0} + l_{e1}) l_{e0} l_{e1} (U_{e0} + U_{e1}),
\end{align*}
$$

$$
H_{\eta} A := \sum_{e \in t} \frac{-2 \|e\| r_e^T r (U_{e0} + U_{e1})}{\lambda_e} + \frac{2 \|e\| r_e^T r}{\lambda_e^2} \mathcal{T}(u_{e0} + u_{e1}; \Lambda_e) \\
+ \sum_{e \in t} \frac{-4 \|e\| (l_{e0} + l_{e1}) r_e^T r}{\lambda_e^3} \Lambda_e \Lambda_e^T + \frac{2 \|e\| (l_{e0} + l_{e1}) r_e^T r}{\lambda_e^2} H_{\eta} \lambda_e.
$$
# 4 APPLICATIONS AND RESULTS

## 4.1 Biharmonic deformations

In this section, we show what new deformation controls are offered by the setting of higher-order boundary constraints. In this setup, a point $η$ is transformed using the following expression:

$$f(η) = \sum_{i \in V} \phi_i(η)v'_i + \sum_{j \in T} \psi_j(η)n'_j + \sum_{i \in V} \bar{\phi}_i(η)L'_i + \sum_{j \in T} \bar{\psi}_j(η)D'_j, \quad (19)$$

$(L'_i, D'_j)$ denoting the Laplacian of the deformation function at vertex $i$ and its normal derivative on triangle $j$ respectively. Note again that the rest pose shape is obtained by setting the following boundary conditions: $v'_i = v_i, n'_j = n_j; L'_i = D'_j = 0$.

While manually playing with $(L', D')$ allows deforming shapes in novel manners, the resulting deformations might exhibit an unpredictable behavior. We present here a strategy inspired by the work of Weber et al. [2012], allowing expressing $f$ as

$$f(η) = \sum_{i \in V} \alpha_i(η)v'_i + \sum_{j \in T} \beta_j(η)n'_j, \quad (20)$$

$(\alpha_i, \beta_j)$ being biharmonic coordinates of $η$, with respect to cage vertices and triangle normals only. The derivation relies on the evaluation of $f$ and its derivatives on locations sampled on the cage, and setting $(L'_i, D'_j)$ for those quantities to be enforced at sampled locations. This allows expressing $(L'_i, D'_j)$ in terms of the main deformation variables (new vertex locations $V'$, new triangle normal constraints $N'$). Note that, while the derivation is inspired by Weber et al. [2012], their strategy cannot be adapted from the 2D case to 3D directly, as it relies on the fact that 2D cages have the same number of vertices and faces (edges). This is obviously not true at all for triangle cages in 3D – thus the differences that both strategies exhibit.

As a reminder, while our deformation function $f$ is derived from affine $1^{st}/3^{rd}$ and per-triangle-constant $2^{nd}/4^{th}$ order boundary conditions using Eq. (3), $f$ fits none of our four redundant conditions exactly (see Section 2.1). In the following, we call $(k, m)$-regularized biharmonic coordinates those obtained by enforcing the fitting of the $k^{th}$ and $m^{th}$ boundary conditions at point locations on the cage.

### 4.1.1 (1, 3)-regularized biharmonic deformations. Since

$$\Delta f(η) = \sum_{i \in V} \phi_i(η)L'_i + \sum_{j \in T} \psi_j(η)D'_j, \quad (21)$$

one can evaluate this expression at points $p_k \in \partial\Omega$, and derive a first constraint by enforcing linear interpolation of $L'$ at $p_k$:

$$\sum_{i} \Gamma_{i}(p_{k}) L'_{i} = \Delta f(p_{k}) = \sum_{i \in V} \phi_{i}(p_{k}) L'_{i} + \sum_{j \in T} \psi_{j}(p_{k}) D'_{j}.$$

Considering this constraint for a collection of samples $\{p_k\}$ leads to the following matrix constraint:

$$(\Phi_L - M_L)L' + \Psi_L D' = 0, \quad (22)$$

$M_L$ denoting a mass matrix that averages, for each line $k$ corresponding to the $k^{th}$ sample, the quantities defined on the vertices (i.e., each row of $M_L$ has three $\Gamma_i(p_k)$ entries, that sum to 1, at the columns of the vertices indexed by the triangle on which $p_k$ lies).

Considering similarly the interpolation of the Dirichlet constraints at sparse locations $q_k \in \partial\Omega$ leads similarly to

$$\sum_{i} \Gamma_{i}(q_{k}) v'_{i} = \sum_{i \in V} \phi_{i}(q_{k}) v'_{i} + \sum_{j \in T} \psi_{j}(q_{k}) n'_{j} + \sum_{i \in V} \bar{\phi}_{i}(q_{k}) L'_{i} + \sum_{j \in T} \bar{\psi}_{j}(q_{k}) D'_{j}.$$

Stacking all constraints for a collection of sampled constraints $\{q_k\}$ leads to the following matrix constraint:

$$(\Phi_V - M_V)V' + \Psi_V N' + \bar{\Phi}_V L' + \bar{\Psi}_V D' = 0. \quad (23)$$

Solving Eqs. (22) and (23) in the least-squares sense leads to the following minimization problem:

$$\begin{pmatrix} L' \\ D' \end{pmatrix} = \underset{X}{\operatorname{argmin}} \|AX - b\|^2 + \|BX\|^2, \quad \text{with} \quad (24)$$

$$\begin{cases} A := (\bar{\Phi}_V; \bar{\Psi}_V) \\ B := (\Phi_L - M_L; \Psi_L) \\ b := (M_V - \Phi_V)V' - \Psi_V N' \end{cases}$$

leading to, with $\begin{pmatrix} C_L \\ C_D \end{pmatrix} := (A^T A + B^T B)^{-1} A^T:$

$$\begin{pmatrix} L' \\ D' \end{pmatrix} = \begin{pmatrix} C_L \\ C_D \end{pmatrix} \left( (M_V - \Phi_V)V' - \Psi_V N' \right) \quad (25)$$

Considering an input mesh $\mathcal{M}$ encoded using our unconstrained biharmonic coordinates (see Eq. (19)) as (in matrix form):

$$\mathcal{M}' = \Phi_M V' + \Psi_M N' + \bar{\Phi}_M L' + \bar{\Psi}_M D', \quad (26)$$

plugging in our matrix constraints (Eq. (25)) leads to the following expression for our regularized biharmonic coordinates (see Eq. (20)):

$$\begin{align} \mathcal{M}' &= [\Phi_M + (\bar{\Phi}_M C_L + \bar{\Psi}_M C_D)(M_V - \Phi_V)] V' \nonumber \\ &\quad + [\Psi_M - (\bar{\Phi}_M C_L + \bar{\Psi}_M C_D)\Psi_V] N' \nonumber \\ &=: [\bar{\Phi}_M + \bar{\Phi}_M C_{11} + \bar{\Psi}_M C_{12}] V' + [\Psi_M + \bar{\Phi}_M C_{21} + \bar{\Psi}_M C_{22}] N' \nonumber \\ &=: \alpha_M^{(1,3)} V' + \beta_M^{(1,3)} N'. \tag{27} \end{align}$$

Evaluating those expressions requires the computation of the biharmonic coordinates on the cage itself, which can be done by taking limits of the expressions. We give those in Appendix C. Fig. 2 illustrates the behavior of our coordinates, and compares them with GC and MVC. In this figure, we set the Neumann condition to be the deformed triangle normal (i.e., $\frac{\partial f}{\partial n} = n'_i$ on the triangle $t$).

**Overconstrained regularization.** Note that we described our construction by considering an arbitrarily large number of Laplacian and Dirichlet interpolation constraints, leading to an overconstrained system minimized in Eq. (24). We found that sampling more constraints than needed improved the quality of our deformations. Note that harmonic coordinates (surprisingly, not the biharmonic parts) are discontinuous at the cage vertices. This is because $\phi$ integrates the normal derivative of the Green kernel, which is singular on $\partial\Omega$.

This leads to a discontinuity as the cage normal field is discontinuous at the cage vertices and edges. We thus avoid sampling cage constraints on the cage vertices and edges, and sample constraints strictly inside triangles only (we use a simple regular sampling in barycentric space, corresponding to a recursive 1-4 triangle split – see inset).

<figure><img src="image_4.png" /></figure>
<figure><img src="image_5.png" /><figcaption>Fig. 2. GC-based deformations are quasi-conformal, but largely ignore the prescribed boundary conditions, and do not follow the cage precisely (Spiky-Box, FireHydrant, Snake examples). MVC are (Dirichlet-) *interpolant* and follow the cage exactly, but may exhibit severe deformation artifacts (Spiky-Box, Cactus, Snake examples). Our (1, 3)–regularized biharmonic coordinates are *pseudo-interpolant* (i.e., interpolation is loosely enforced at sparse locations only, instead of everywhere on $\partial\Omega$). In these examples, we used a unit normal stretch (i.e., $\frac{\partial f}{\partial n} = n'_t$ on $t$). Like GC, our deformations are cage-aware, and do not result in the visible artifacts of MVC on the Cactus/Snake examples (due to some MVC being negative here).</figcaption></figure>

*Comparison with alternate constructions.* We evaluate in Fig. 3 two alternative versions of regularized biharmonic coordinates on simple examples, and give in Appendix D their derivation for completeness.

<figure><img src="image_6.png" /><figcaption>Fig. 3. Our (1, 3)–regularized biharmonic deformations vs. alternatives.</figcaption></figure>

One remarkable point is that, while biharmonic deformations are generally used to enforce both Dirichlet and Neumann boundary conditions (e.g., in Finite Elements Methods), using (1, 2)–regularized biharmonic deformations leads to undesirable effects in practice. The reason for this is simple: while we can indeed enforce the Dirichlet and Neumann boundary conditions on the cage (loosely if the number of constraints is larger than the number of DoFs $\{L'_i, D'_j\}$), those joint conditions are too agressive to be fit with biharmonic functions defined by such a small subspace. To fit those constraints, high-amplitude values of $(L', D')$ are set as $3^{rd}/4^{th}$ order boundary conditions, leading to large oscillations inside $\Omega$.

In light of this, the choice of our (1, 3)–regularized biharmonic coordinates is motivated by the following remarks:

* Fitting both Dirichlet and Neumann constraints requires strong oscillations/variations of the deformation near the cage, as its deformation is $C^0$ only at the cage edges and vertices. In this context, as artists deform cages by moving their vertices, it appears more natural to regularize the Dirichlet condition rather than the Neumann one;

* Enforcing loose-only Neumann boundary condition allows deriving smooth deformations in this context;

* Ensuring interpolation of $\Delta f$ along the cage faces enforces deformations with minimal oscillations in practice (as $\Delta f$ across the triangles is then heavily constrained). Since $\Delta f$ is harmonic (as $f$ is biharmonic), it verifies the maximum property, and controlling $\Delta f$ on $\partial\Omega$ ensures a controlled behavior of $\Delta f$ inside $\Omega$, thus ensuring minimal oscillations inside $\Omega$.

Note that there are minor differences in spirit between our (1, 3)–coordinates and the ones presented by Weber et al. [2012] in 2D. In their optimization, the Laplacian regularization constraints are first solved by optimizing for $\{D'_j\}$ (see Eq. (63)); only then $\{L'_i\}$ are optimized to fit the Dirichlet regularization constraints (see Eq. (64)). In our setup, we jointly optimize both sets of variables to best fit both constraints. In a sense, we provide more degrees of freedom ($L'$ and $D'$, not just $D'$) to regularize $\Delta f$ on the cage. We observe that, as a result, our coordinates are consistently better-behaved and do not exhibit similar oscillations (see Fig. 3). We stress that, while we named the construction presented in Appendix D.2 “Weber’s formulation” to acknowledge the source of inspiration, their exact setup cannot be reproduced in 3D.
### 4.1.2 Balancing regularization energies.
We can balance both constraint types (e.g., to obtain "better-behaved" deformations at the cost of following the Dirichlet constraints more loosely) by minimizing

$$ \begin{pmatrix} L' \\ D' \end{pmatrix} = \underset{X}{\operatorname{argmin}} \gamma_D \|AX - b\|^2 + \|BX\|^2 \quad (28) $$

in place of Eq. (24), leading to

$$ \begin{pmatrix} C_L \\ C_D \end{pmatrix} := \underbrace{\gamma_D (\gamma_D A^T A + B^T B)^{-1} A^T}_{\text{used for } \gamma_D < 1} = \underbrace{(A^T A + \gamma_D^{-1} B^T B)^{-1} A^T}_{\text{used for } \gamma_D \ge 1} \quad (29) $$

$$ \begin{pmatrix} L' \\ D' \end{pmatrix} = \begin{pmatrix} C_L \\ C_D \end{pmatrix} \left( (M_V - \Phi_V)V' - \Psi_V N' \right). \quad (30) $$

Interestingly, setting $\gamma_D = 0$ results in the Green coordinates, as only purely harmonic functions (i.e., $L'_i = D'_j = 0$) can exactly lead to a strict interpolation of their Laplacian across the cage faces with the subspace we use. If, on the contrary, one sets large values for $\gamma_D$, the position boundary constraint will be better fit. We show in Fig. 4 the influence of $\gamma_D$ on challenging deformations. Apart from this figure, all other results are given for $\gamma_D = 1$.

### 4.1.3 Variational Neumann conditions.
While enforcing a straightforward unit Neumann condition can lead to satisfactory results on many examples, it can also produce undesirable deformation behavior: the deformations may contain foldovers or high oscillations. The main reason is that the Neumann condition is a local quantity,

<figure><img src="image_7.png" /><figcaption>Fig. 4. Influence of $\gamma_D$ on our (1, 3)-regularized biharmonic coordinates. We highlight unwanted wavy artifacts due to poor Laplacian regularization. A larger $\gamma_D$ will result in deformations that will follow the cage more closely, at the cost of eventually featuring such deformation artifacts. For reference, we show the construction adapted from the work of Weber and colleagues [Weber et al. 2012] on 2D biharmonic deformations.</figcaption></figure>

and avoiding the mentioned artifacts require global knowledge of the deformation. For example, enforcing a unit normal derivative may result in foldovers when the shape's thickness is drastically diminished (e.g., when opposite structures are brought closer to each other – see the 2^nd^ example of Fig. 5).

Such difficulties were already observed by Weber et al. [2012] for 2D deformations, and motivated more complex Neumann condition optimization strategies. Marching in their steps, we show how to derive the Neumann condition automatically from the Dirichlet condition, in order to satisfy quality-preserving criteria.

Next, we introduce a set of deformation energies formulated as

$$ \mathcal{E}_{\text{REG}} := \|\Lambda_V V' + \Lambda_N N'\|^2, \quad (31) $$

$\Lambda_V$ and $\Lambda_N$ being matrices, which we detail for the various energies we introduce next (see 4.1.3: Deformation energies).

**Normal alignment constraints.** We adapt the normal alignment strategy of Weber and colleagues, and enforce the normal derivative to be aligned with the deformed triangle normal, while adjusting its scale to minimize $\mathcal{E}_{\text{REG}}$ (Eq. (31)).

Considering $\mathcal{E}_{\text{REG}}$ in serialized form (i.e., stacking all variables: $V' \in \mathbb{R}^{3V}$, $\notin \mathbb{R}^{V \times 3}$), we express our constrained $N' \in \mathbb{R}^{3T}$ as

$$ N' = \begin{pmatrix} \sigma_0 n'_0 \\ \sigma_1 n'_1 \\ \vdots \end{pmatrix} = \underbrace{\begin{pmatrix} n'_0 & 0 & 0 & \cdots \\ 0 & n'_1 & 0 & \cdots \\ \vdots & \vdots & \vdots & \vdots \\ \vdots & \vdots & \vdots & \vdots \end{pmatrix}}_{=:M_{N'} \in \mathbb{R}^{3T \times T}} \underbrace{\begin{pmatrix} \sigma_0 \\ \sigma_1 \\ \vdots \end{pmatrix}}_{=:\sigma \in \mathbb{R}^T} \quad (32) $$

Minimizing $\mathcal{E}_{\text{REG}}$ in this setup amounts to minimizing

$$ \mathcal{E}_{\text{REG}}^{\text{align}}(N'; V') := \|\Lambda_N M_{N'} \sigma + \Lambda_V V'\|^2. \quad (33) $$

**Soft alignment constraints.** Rather than ensuring explicit alignment between $N'$ and the normals of the deformed cage (given vertex positions $V'$), we can enforce it by adding a soft penalty term:

$$ \mathcal{E}_{\text{REG}}^{\text{soft align}}(N'; V'; \gamma_a) := \|\Lambda_N N' + \Lambda_V V'\|^2 + \gamma_a \|\text{Al}(V')N'\|^2, \quad (34) $$

each $k^{th}$ row of $\text{Al}(V')$ being a $3 \times 3$ block row where $(n'_k)_{[\times]}$ appears at column $k$ $((a)_{[\times]}$ denoting the $3 \times 3$ matrix such that $(a)_{[\times]}u = a \times u, \forall u \in \mathbb{R}^3)$. Note that, by setting $\gamma_a \mapsto \infty$, we recover the hard constraints case described earlier. On the contrary, setting $\gamma_a = 0$ results in dropping the normal alignment constraint entirely.

**Implementation notes.** The left-hand sides of Eqs. (33) and (34) depend on the cage vertex positions. Those cannot be factorized once and for all, and a complete minimization is thus required each time a new deformation is performed by the artist.

As advocated in [Weber et al. 2012], in order to encourage non-overlapping local deformations, one may constrain $\sigma$ to be strictly positive. While it may help locally, it does not prevent global intersections nor non-overlapping deformations far from the cage. In our implementation, we ignore this constraint, and perform a simple unconstrained least-squares minimization of $\mathcal{E}_{\text{REG}}$ (Eqs. (33) and (34)).
<figure><img src="image_8.png" /><figcaption>Fig. 5. Comparison of the Neumann unit stretch condition with optimized Neumann conditions (i.e., minimizing $\mathcal{E}_{\text{REG}}$). $\gamma_a = +\infty$ (resp. = 0) implies hard (resp. no) alignment constraint.</figcaption></figure>

*Deformation energies.* We sample point constraints $C = \{q_k\}$ on $\partial\Omega$, to discretize our three different deformation energies:

* As-Harmonic-As-Possible (AHAP), minimizing $\sum_k ||\Delta f(q_k)||^2$;

* As-Affine-As-Possible (AAAP), minimizing $\sum_k ||Hf(q_k)||^2$;

* Thickness-preserving (TP), minimizing $\sum_k \left\|\frac{\partial^2 f}{\partial n_k^2}(q_k)\right\|^2$.

To compute those energies, we compute the biharmonic coordinates and derivatives at sample points $C$ and stack those to build the following matrices:

* $\Phi_C := \Delta\bar{\Phi}_C \in \mathbb{R}^{C \times V}$ (resp. $\Psi_C \in \mathbb{R}^{C \times T}$): row $k$ being the $\phi$ (resp. $\psi$) coordinates of sample $q_k$;

* $\partial_{n^2}\Phi_C \in \mathbb{R}^{C \times V}$: row $k$ being the second normal derivative of $\phi$ at $q_k$. Noting $n_k$ the cage's unit normal at $q_k$, we have $\partial_{n^2}\Phi_C(k, v) = n_k^T H \phi_v(q_k) n_k$;

* $\partial_{n^2}\Psi_C, \partial_{n^2}\bar{\Phi}_C, \partial_{n^2}\bar{\Psi}_C$ built similarly;

* $h\Phi_C \in \mathbb{R}^{9 \times V}$: rows 9k to 9k + 8 being the serialized 9D vectors of $H\phi$ at $q_k$;

* $h\Psi_C, h\bar{\Phi}_C, h\bar{\Psi}_C$ built similarly.

Equipped with those matrices, and using the notations of Eq. (27), our deformation energies read:

$$\mathcal{E}_{\text{AHAP}} := \| \Delta \alpha_C^{(1,3)} V' + \Delta \beta_C^{(1,3)} N' \|^2,$$

$$\mathcal{E}_{\text{AAAP}} := \| h \alpha_C^{(1,3)} V' + h \beta_C^{(1,3)} N' \|^2,$$

$$\mathcal{E}_{\text{TP}} := \| \partial_{n^2} \alpha_C^{(1,3)} V' + \partial_{n^2} \beta_C^{(1,3)} N' \|^2,$$

with, for each linear operator $O \in \{\Delta; h; \partial_{n^2}\}$:

$$O \alpha_C^{(1,3)} = O \Phi_C + O \bar{\Phi}_C C_{11} + O \bar{\Psi}_C C_{12}, \quad (35)$$

$$O \beta_C^{(1,3)} = O \Psi_C + O \bar{\Phi}_C C_{21} + O \bar{\Psi}_C C_{22}. \quad (36)$$

We compare in Fig. 5 our (1, 3)-constrained biharmonic deformations with $\sigma = 1$ (unit normal derivative constraint) and $\sigma = \arg\min \mathcal{E}_{\text{REG}}$ for our three energies.

While Weber and colleagues advocate the use of AHAP deformations in 2D, we observe that they do not always produce satisfactory deformations in 3D (especially when enforcing the alignment onto the deformed normals – see $\gamma_a = +\infty$ on the last example of Fig. 5). Our understanding is that harmonic deformations are much more permissive in 3D than in 2D, while boundary conditions are more intricate. The AAAP and TP energies (the thickness-preserving energy being introduced in [Weber et al. 2012]) always produce viable deformations, and are sometimes preferable to unit-stretch biharmonic deformations – depending on the application scenario.

Focusing on the Gear example in Fig. 5, we observe that interestingly, our AAAP optimizer recovers indeed *exactly affine deformations* when the cage is deformed accordingly (here, the cage has been stretched uniformly along the z-axis), when no alignment constraints are set ($\gamma_a = 0$). If the Neumann condition is constrained to be aligned onto the deformed normals ($\gamma_a = +\infty$), an unavoidable local rotation of the shape is observed, as stretching the cage vertex positions uniformly results in rotation of the normals.

**Published additional material.** We provide as additional material:

* a video comparing our (1, 3)-regularized coordinates with MVC, GC, and the two alternatives to our regularized biharmonic coordinates, as well as showcasing the impact of $\gamma_D$ on several examples;

* C++ code for the computation of our coordinates (Eq. (19)), their gradients and Hessians, their special case computations, and the matrices needed for computation of our regularized coordinates (Eq. (27)). The code is accessible at https://gitlab.com/jthiery/biharmonic3dcordinates_siggraph2024.

## 4.2 Variational shape deformation

As demonstrated by Ben-Chen et al [2009], Green coordinates are an efficient subspace basis for variational deformations, as $\{\phi_i, \psi_j\}$ are harmonic (smooth) and they allow deforming 3D shapes in a *cage-aware, quasi-conformal* manner. Ben-Chen et al. used those as the backbone of their variational as-rigid-as-possible (ARAP) optimization method. Using as input a triangular manifold mesh acting as the cage, they define a harmonic function everywhere inside its volume as

$$f(\eta) = \sum_{i \in V} \phi_i(\eta) a_i + \sum_{j \in T} \psi_j(\eta) b_j, \quad (37)$$
$\{a_i, b_j\}$ being arbitrary 3D parameters associated with cage vertices and triangles (note that $b_j$ is not a 3D normal constrained to be orthogonal to the triangle $j$ defined by its $\{a_i\}$ corner values).

The gradient and Hessian of $f$ are obtained at any point using the gradients and Hessians of the coordinates, as

$$Jf(\eta) = \sum_{i \in V} a_i \cdot \nabla \phi_i(\eta)^T + \sum_{j \in T} b_j \cdot \nabla \psi_j(\eta)^T, \quad (38)$$

$$Hf(\eta) = \sum_{i \in V} a_i \cdot H\phi_i(\eta) + \sum_{j \in T} b_j \cdot H\psi_j(\eta). \quad (39)$$

To achieve real-time deformations, the authors proposed to couple the complexity of the deformed mesh from the complexity of the expected deformation, by setting appropriate deformation constraints on the ambient space rather than on the mesh directly. In particular, as volumetric ARAP deformations reproduce empirically local rotations near the medial axis of the deformed shapes, they propose to sample sparse rigidity constraints on points $\mathcal{R}_S$ located on the medial axis. To enforce smoothness and minimize deformation undulations, they additionally sample points $\mathcal{H}_S$ near the surface of the cage, where the Hessian of the deformation is minimized.

Considering additional user-specified positional $\mathcal{P}_C$ and orientation constraints $\mathcal{J}_C$, the cage geometric parameters $\{a_i, b_j\}$ are optimized to minimize the final deformation energy:

$$\mathcal{E} = \sum_{p \in \mathcal{P}_C} \|f(p) - \bar{p}\|^2 + \alpha \sum_{p \in \mathcal{J}_C} \|Jf(p) - \bar{J}_p\|^2 + \quad (40)$$

$$\beta \sum_{p_j \in \mathcal{R}_S} \|Jf(p_j) - R_j\|^2 + \gamma \sum_{p_h \in \mathcal{H}_S} \|Hf(p_h)\|^2.$$

This energy is minimized using a global/local approach, by alternating optimization of auxiliary unknown rotations $\{R_j\}$ (local solve) and optimization of the main cage parameters $\{a_i, b_j\}$ (global solve), which is done simply by inverting a linear system in the least-squares sense. The result of the optimization is the set of parameters $\{a_i, b_j\}$ that best fits the deformation constraints, and that can be used to extend the deformation to any 3D shape embedded in the cage at binding time. Note that optimizing for as-similar-as-possible deformations can be done by constraining each Jacobian at $\mathcal{R}_S$ to be a similarity matrix instead.

Since the energy of Eq. (40) depends only on the few sampled constraints as well as the cage complexity, it can be performed in real time, regardless of the size of the input shape to deform. Later transferring the deformation from the cage to the shape involves matrix-vector multiplications only, and can be performed efficiently on the GPU as well, as long as the model fits the video card memory.

We illustrate in Fig. 6 the use of our biharmonic coordinates as an alternate subspace for variational deformations on a toy example (top). As can be seen, even using a simple box cage split in two (left column, V=12), our biharmonic subspace is rich enough to allow for smooth rotations of the shape. In the middle, while increasing the number of degrees of freedom allows for richer deformations using a harmonic subspace, the cage discretization is still revealed in the shape deformation. Further increasing the number of DoFs (right column), this effect is reduced but still present. Without claiming

<figure><img src="image_9.png" /><figcaption>Fig. 6. **Top:** Comparing harmonic and biharmonic subspaces for variational deformation on a toy example. **Bottom:** Rich smooth biharmonic variational deformation using a simple cage. We sampled Hessian constraints near the cage surface, as well as rotation constraints on the shape's 1D medial axis (i.e., mimicking the setup presented in [Ben-Chen et al. 2009]). We use our unconstrained biharmonic subspace (Eq. (19)) for those examples.</figcaption></figure>

any related contribution, we believe that this effect is in fact inherent to harmonic deformations, and that it remains regardless of the cage local density. In spirit, ARAP deformations require indeed fitting both positional constraints (akin to Dirichlet constraints) and gradient constraints (akin to Neumann constraints). Fitting both constraint types while allowing for smooth deformations require using biharmonic (or higher-order) functions, as harmonic functions are already uniquely defined from either a Dirichlet or a Neumann boundary condition. The bottom row shows a smooth ARAP biharmonic deformation obtained with a simple cage structure.

## 4.3 Affine Neumann conditions and their link with elasticity-based Somigliana coordinates

We finally present in this section formulas for substituting an affine Neumann condition to a constant one per cage triangle $t = (t^0, t^1, t^2)$:

$$\frac{\partial f}{\partial n}(\xi) = \sum_{k=0}^{2} \Gamma_{tk}^{t}(\xi) n_{t}^{k} \quad \forall \xi \in t, \quad (41)$$

$n_t^k$ being the 3D normal on triangle $t$ at its $k^{th}$ corner. This results in the following Neumann contribution:

$$f_N(\eta) = \sum_{j \in T} \sum_{k=0}^{2} \psi_t^k(\eta) n_t^k, \quad (42)$$

$$\psi_t^k(\eta) := \int_{\xi \in t} -\Gamma_{tk}^t(\xi)G(\xi, \eta)d\xi = -\lambda_1 \int_{\xi \in t} \frac{\Gamma_{tk}^t(\xi)}{\|\xi - \eta\|} d\xi. \quad (43)$$
We observe that

$$ 
\begin{aligned}
\psi_t^k(\eta) &= \frac{\lambda_1}{\lambda_2 d_t} \bar{\phi}_{t^k}^t(\eta) \\
&= -2 \nabla_{\eta} (\bar{\psi}_t)^T \cdot g_{t^k}^t + \Gamma_{t^k}^t(\eta) \psi_t(\eta).
\end{aligned} 
\quad (44) $$

It follows that

$$ \nabla_{\eta} \psi_t^k(\eta) = -2H_{\eta}(\bar{\psi}_t) \cdot g_{t^k}^t + \Gamma_{t^k}^t(\eta) \nabla_{\eta} \psi_t(\eta) + \psi_t(\eta) g_{t^k}^t, \quad (45) $$

$$ H_{\eta}(\psi_t^k)(\eta) = -2\nabla_{\eta}^T (H_{\eta}(\bar{\psi}_t) \cdot g_{t^k}^t) + \mathcal{T}(\nabla_{\eta}\psi_t(\eta); g_{t^k}^t) + \Gamma_{t^k}^t(\eta)H_{\eta}(\psi_t)(\eta). \quad (46) $$

While playing with this novel boundary condition offers little new expressive power to artists (in our early experiments), it is mathematically required to derive closed-form expressions for Somigliana coordinates, which we demonstrate next.

**Closed-form expressions for Somigliana coordinates.** Somigliana coordinates in 3D, as introduced in [Chen et al. 2023], are defined as

$$ 
\begin{aligned}
K_t(\eta) &:= \int_{\xi \in t} \mathcal{K}(\eta, \xi) d\xi, \\
T_{t^i}^t(\eta) &:= \int_{\xi \in t} \Gamma_{t^i}^t(\xi) \mathcal{T}(\eta, \xi) d\xi.
\end{aligned} 
\quad (47) $$

where $\mathcal{K}$ and $\mathcal{T}$ are the fundamental solutions of displacement and traction respectively, written as

$$ 
\begin{aligned}
\mathcal{K}(\xi, \eta) &= \frac{(a-b)}{\|r\|} I_3 + \frac{b}{\|r\|^3} rr^T, \\
\mathcal{T}(\xi, \eta) &= \frac{\mu(a-2b)}{\|r\|^3} [(n_t^T r)I_3 + n_t r^T - rn_t^T] + \frac{6\mu b}{\|r\|^5} (n_t^T r)rr^T,
\end{aligned} 
\quad (48) $$

where $r = \xi - \eta$, $a = 1/4\pi\mu$, $b = a/4(1 - \nu)$ with $\mu$ being the Lamé constant (shear modulus) and $\nu$ the Poisson ratio, describing the stiffness and incompressibility of the elastic material, respectively.

We note that $n_t^T r = n_t^T (\xi - \eta) = -d_t \forall \xi \in t$, and differentiate twice Eqs. (9) and (43) directly w.r.t. $\eta$ to obtain:

$$ 
\begin{aligned}
H_{\eta}(\bar{\psi}_t)(\eta) &= -\lambda_2 \int_{\xi \in t} \frac{I_3}{\|\xi - \eta\|} d\xi + \lambda_2 \int_{\xi \in t} \frac{(\xi - \eta)(\xi - \eta)^T}{\|\xi - \eta\|^3} d\xi, \\
\nabla_{\eta}(\psi_t^i)(\eta) &= \lambda_1 \int_{\xi \in t} \frac{\Gamma_{t^i}^t(\xi)(\eta - \xi)}{\|\xi - \eta\|^3} d\xi, \\
H_{\eta}(\psi_t^i)(\eta) &= \lambda_1 \int_{\xi \in t} \frac{\Gamma_{t^i}^T(\xi)I_3}{\|\xi - \eta\|^3} d\xi - 3\lambda_1 \int_{\xi \in t} \frac{\Gamma_{t^i}^t(\xi)(\eta - \xi)(\eta - \xi)^T}{\|\xi - \eta\|^5} d\xi.
\end{aligned} 
$$

By identification, we relate Somigliana coordinates to ours as:

$$ K_t(\eta) = \underbrace{\frac{b}{\lambda_2} H_{\eta}(\bar{\psi}_t)(\eta) - \frac{a}{\lambda_1} \psi_t(\eta) I_3}_{\text{Eq. (13)}} \quad (49) $$

$$ T_{t^i}^t(\eta) = \underbrace{\frac{2bd_t}{\lambda_1} H_{\eta}(\psi_t^i)(\eta) - \frac{a}{\lambda_1} \phi_{t^i}^t(\eta) I_3 - \frac{a-2b}{\lambda_1} \left[ n_t \nabla_{\eta}^T(\psi_t^i) - \nabla_{\eta}(\psi_t^i) n_t^T \right]}_{\text{Eq. (46)}} \quad (49) $$

The numerical integration adopted in [Chen et al. 2023] can be inaccurate when a query point $\eta$ is too close to the cage surface,

<figure><img src="image_10.png" /><figcaption>Fig. 7. We compared the averaged time cost of integration over a single cage facet for one point using the closed-form expressions and quadrature rule (top). Insufficient quadratures lead to error in satisfying the partition of unity, which may cause notable artifacts (bottom). Using more quadratures alleviates the issue, but the integration is slower by orders of magnitude than evaluating the closed-form expressions.</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>Fig. 8. The quadrature rule for computing Somigliana coordinates can introduce noticeable artifacts in deformation due to its limited accuracy. In contrast, our closed-form expressions are both accurate and efficient.</figcaption></figure>

because of the near singular integration caused by the Green's function. This inaccuracy can result in visible artifacts in deformation, such as artificial spikes on the deformed model in Fig. 8. Therefore, obtaining accurate results may require a large number of quadrature points, which makes the coordinate computation much less efficient, especially with limited multithreading. For instance, the time cost of computing Somigliana coordinates for one query point per cage face increases linearly with the number of quadratures. Fig. 7 shows that using a few quadratures does not provide enough accuracy for even linear reproduction. When the error of partition of unity is reduced to a level where no significant artifact is present (e.g.,
<figure><img src="image_12.png" /><figcaption>Fig. 9. Like all coordinates designed for triangle cages, such as MVC and GC, ours suffer from asymmetric artifacts, when cage quads are arbitrarily cut into two triangles.</figcaption></figure>

around $10^{-3}$), the numerical integration becomes much slower than our analytical integration by about two orders of magnitude.

## 5 DISCUSSION

Our technique suffers from a number of limitations. The main one (in our opinion) is that our biharmonic coordinates are compatible with *triangle* cages only (i.e., cages made of triangle facets). As a result, all our results suffer from asymmetric artifacts when cage quads are arbitrarily cut into two triangles (see Fig. 9). Extending our work to handle triquad cages, in the same way QMVC [Thiery et al. 2018] (resp. QGC [Thiery and Boubekeur 2022]) have extended MVC (resp. GC), is the future work we consider most worthy of pursuing.

As we have demonstrated that our coordinates and Somigliana coordinates are related, we envision that it may be feasible to obtain closed-form expressions for the derivatives of those Somigliana coordinates using the formulas we derived. As our coordinates are well-behaved and their singular cases are well identified (restricted to points lying on the lines supporting the cage edges) and easy to handle, it may be sufficient to obtain those derivatives using automatic differentiation tools; this has to be checked however. This may open interesting avenues for variational elastic deformations based on Somigliana coordinates.

While we have experimented with variational biharmonic deformations, we believe that we have only scratched the surface, as we have considered those for well-known deformation energies only. Designing richer 3D deformation energies may offer additional flexibility and expressiveness to 3D modelers.

Finally, even if we have obtained closed-form expressions, those are rather involved and computationally expensive. Simplifying those may be possible, while deriving those from physically-inspired principles rather than mathematical ones may shed some light on possible usage opportunities and links to other works than Somigliana-based deformations.

## ACKNOWLEDGMENTS

We would like to thank reviewers for their valuable comments on improving this work. J. Chen was partially funded by a Choose France Inria chair (Desbrun/Geomerix).

## REFERENCES

* Mirela Ben-Chen, Ofir Weber, and Craig Gotsman. 2009. Variational harmonic maps for space deformation. *ACM ToG* 28, 3 (2009), 1–11.

* Mario Botsch and Leif Kobbelt. 2004. An intuitive framework for real-time freeform modeling. *ACM Transactions on Graphics* (TOG) 23, 3 (2004), 630–634.

* Christopher Brandt, Elmar Eisemann, and Klaus Hildebrandt. 2018. Hyper-reduced projective dynamics. *ACM Transactions on Graphics* (TOG) 37, 4 (2018), 1–13.

* Qingjun Chang, Chongyang Deng, and Kai Hormann. 2023. Maximum Likelihood Coordinates. In *Computer Graphics Forum*, Vol. 42. Wiley Online Library, e14908.

* Jiong Chen, Fernando de Goes, and Mathieu Desbrun. 2023. Somigliana Coordinates: an elasticity-derived approach for cage deformation. In *Special Interest Group on Computer Graphics and Interactive Techniques Conference Conference Proceedings (SIGGRAPH '23 Conference Proceedings)*. 8 pages. https://doi.org/10.1145/3588432.3591519

* Lu Chen, Jin Huang, Hanqiu Sun, and Hujun Bao. 2010. Cage-based deformation transfer. *Computers & Graphics* 34, 2 (2010), 107–118.

* Fabrizio Corda, Jean Marc Thiery, Marco Livesu, Enrico Puppo, Tamy Boubekeur, and Riccardo Scateni. 2020. Real-Time Deformation with Coupled Cages and Skeletons. *Computer Graphics Forum* (to appear) (2020). https://doi.org/10.1111/cgf.13900

* Ana Dodik, Oded Stein, Vincent Sitzmann, and Justin Solomon. 2023. Variational Barycentric Coordinates. *ACM Transactions on Graphics* (TOG) 42, 6 (2023), 1–16.

* Michael S Floater. 2003. Mean value coordinates. *Computer aided geometric design* 20, 1 (2003), 19–27.

* Kai Hormann and Natarajan Sukumar. 2008. Maximum entropy coordinates for arbitrary polytopes. In *Computer Graphics Forum*, Vol. 27. Wiley Online Library, 1513–1520.

* Alec Jacobson, Ilya Baran, Ladislav Kavan, Jovan Popović, and Olga Sorkine. 2012a. Fast automatic skinning transformations. *ACM Transactions on Graphics* (TOG) 31, 4 (2012), 1–10.

* Alec Jacobson, Ilya Baran, Jovan Popović, and Olga Sorkine. 2011. Bounded Biharmonic Weights for Real-Time Deformation. *ACM Transactions on Graphics* (proceedings of *ACM SIGGRAPH*) 30, 4 (2011), 78:1–78:8.

* Alec Jacobson, Elif Tosun, Olga Sorkine, and Denis Zorin. 2010. Mixed Finite Elements for Variational Surface Modeling. *Computer Graphics Forum* (proceedings of *EUROGRAPHICS/ACM SIGGRAPH Symposium on Geometry Processing*) 29, 5 (2010), 1565–1574.

* Alec Jacobson, Tino Weinkauf, and Olga Sorkine. 2012b. Smooth Shape-Aware Functions with Controlled Extrema. *Computer Graphics Forum* (proceedings of *EUROGRAPHICS/ACM SIGGRAPH Symposium on Geometry Processing*) 31, 5 (2012), 1577–1586.

* Pushkar Joshi, Mark Meyer, Tony DeRose, Brian Green, and Tom Sanocki. 2007. Harmonic coordinates for character articulation. *ACM TOG* 26, 3 (2007), 71–es.

* Tao Ju, Scott Schaefer, and Joe Warren. 2005. Mean value coordinates for closed triangular meshes. In *ACM Siggraph 2005 Papers*. 561–566.

* Torsten Langer, Alexander Belyaev, and Hans-Peter Seidel. 2006. Spherical barycentric coordinates. In *Symposium on Geometry Processing*. 81–88.

* Yaron Lipman, Johannes Kopf, Daniel Cohen-Or, and David Levin. 2007. GPU-assisted positive mean value coordinates for mesh deformations. In *SGP*. 117–123.

* Yaron Lipman, David Levin, and Daniel Cohen-Or. 2008. Green coordinates. *ACM ToG* 27, 3 (2008), 1–10.

* Élie Michel and Jean-Marc Thiery. 2023. Polynomial 2D Green Coordinates for Polygonal Cages. In *ACM SIGGRAPH 2023 Conference Proceedings*. 1–9.

* August Ferdinand Möbius. 1827. *Der barycentrische Calcul; ein neues Hülfsmittel zur analytischen Behandlung der Geometrie... Mit vier Kupfertafeln*. Barth.

* Oded Stein, Eitan Grinspun, Alec Jacobson, and Max Wardetzky. 2019. A mixed finite element method with piecewise linear elements for the biharmonic equation on surfaces. *arXiv preprint arXiv:1911.08029* (2019).

* Jean-Marc Thiery and Tamy Boubekeur. 2022. Green Coordinates for Triquad Cages in 3D. In *SIGGRAPH Asia 2022 Conference Papers*. 1–8.

* Jean-Marc Thiery, Pooran Memari, and Tamy Boubekeur. 2018. Mean value coordinates for quad cages in 3D. *ACM ToG* 37, 6 (2018), 1–14.

* Philipp von Radziewsky, Elmar Eisemann, Hans-Peter Seidel, and Klaus Hildebrandt. 2016. Optimized subspaces for deformation-based modeling and shape interpolation. *Computers & Graphics* 58 (2016), 128–138.

* Christoph Von-Tycowicz, Christian Schulz, Hans-Peter Seidel, and Klaus Hildebrandt. 2015. Real-time nonlinear shape interpolation. *ACM Transactions on Graphics* (TOG) 34, 3 (2015), 1–10.

* Yu Wang, Alec Jacobson, Jernej Barbič, and Ladislav Kavan. 2015. Linear subspace design for real-time shape deformation. *ACM Transactions on Graphics* (TOG) 34, 4 (2015), 1–11.

* Ofir Weber, Mirela Ben-Chen, Craig Gotsman, et al. 2009. Complex barycentric coordinates with applications to planar shape deformation. In *Computer Graphics Forum*, Vol. 28. Citeseer, 587.

* Ofir Weber, Roi Poranne, and Craig Gotsman. 2012. Biharmonic coordinates. In *Computer Graphics Forum*, Vol. 31. Wiley Online Library, 2409–2422.

* Juyong Zhang, Bailin Deng, Zishun Liu, Giuseppe Patanè, Sofien Bouaziz, Kai Hormann, and Ligang Liu. 2014. Local barycentric coordinates. *ACM Transactions on Graphics* (TOG) 33, 6 (2014), 1–12.
<figure><img src="image_13.png" /><figcaption>Fig. 10. *t* can be decomposed into 3 triangles (*T*₀, *T*₁, *T*₂) with (here) orientations *s*(*T*₀) = 1, *s*(*T*₁) = 1, *s*(*T*₂) = -1: *t* = *T*₀ + *T*₁ - *T*₂.</figcaption></figure>

A DERIVATION OF $\bar{\psi}_t(\eta)$

We detail here the derivation of Eq. (10), that gives the result of
$$ \int_{\xi \in t} \|\xi - \eta\| d\xi. $$

Following the lines of Lipman and colleagues [2008], we decom-
pose the triangle $t = (t^0, t^1, t^2)$ into three triangles $T_i := (p, t^i, t^{(i+1)\%3})$,
$p$ being the projection of $\eta$ onto the support plane of $t$: $p := \eta +$
$((t^0 - \eta) \cdot n_t)n_t$, and express the integral as

$$
\int_{\xi \in t} \|\xi - \eta\| d\xi = \sum_{i=0}^{2} s(T_i) \int_{\xi \in T_i} \|\xi - \eta\| d\xi, \quad (50)
$$

where "the sign of $T_i$" $s(T_i) \in \{-1; 1\}$ describes whether $T_i$ has
similar orientation than $t$ or not (see Fig. 10).

Derivation of $\int_{\xi \in T_i} \|\xi - \eta\| d\xi$ (Eq. (10)). To simplify notations, we focus on $T_0 = (p, t^0, t^1)$, but the derivation for the other two triangles $T_1$ and $T_2$ is obviously similarly obtained.

focus on $T_0 = (p, t^0, t^1)$, but the derivation for the other two triangles $T_1$ and $T_2$ is obviously similarly obtained.

Using radial coordinates $\xi \to (r, \theta)$ as $\xi = p + r(\cos(\theta)u_1 + \sin(\theta)u_2)$ (see Fig. 11-left), and noting *d* the unsigned distance from $\eta$ to the tangent plane supporting *t*, we obtain:

$$
\begin{align*}
\int_{\xi \in T_0} \|\xi - \eta\| d\xi &= \int_{\theta=0}^{\beta} \int_{r=0}^{R(\theta)} \sqrt{d^2 + r^2} \overbrace{rdrd\theta}^{d\xi(r,\theta)} \\
&= -\frac{d^3\beta}{3} + \frac{1}{3} \int_{\theta=0}^{\beta} \left(d^2 + R(\theta)^2\right)^{3/2} d\theta \\
&= -\frac{d^3\beta}{3} + \frac{1}{3} \int_{\theta=0}^{\beta} \left(d^2 + \frac{a^2}{\sin(\pi - \gamma_0 - \theta)^2}\right)^{3/2} d\theta \\
&= -\frac{d^3\beta}{3} + \frac{1}{3} \int_{\theta=\gamma-\beta}^{\gamma} \left(d^2 + \frac{a^2}{\sin(\theta)^2}\right)^{3/2} d\theta,
\end{align*}
$$

where (see Fig. 11-middle) $R(\theta) = \delta_0 \sin(\gamma_0)/\sin(\pi - \gamma_0 - \theta)$ is
obtained using the law of sines, $a := \delta_0 \sin(\gamma_0) = \delta_1 \sin(\gamma_1)$ (see

<figure><img src="image_14.png" /><figcaption>Fig. 11. Notations for computations within $T_0$. All lengths and angles are unsigned.</figcaption></figure>

Fig. 11-right) is simply the (unsigned) distance from *p* to the line (*t*⁰*t*¹), and the last integral corresponds to a simple change of variables *θ* → *π* − *γ*₀ − *θ* (further noting *γ* the angle *γ* := *π* − *γ*₀).

We focus on the remaining integral for now, and using a change of variables $x = \tan(\theta)$, we obtain:

$$
\int_{\theta=\beta-\gamma}^{\gamma} \left( d^2 + \frac{a^2}{\sin(\theta)^2} \right)^{3/2} d\theta = \qquad (51)
$$

$$
\int_{x=\tan(\gamma-\beta)}^{\tan(\gamma)} \left( (d^2 + a^2) + \frac{a^2}{x^2} \right)^{3/2} \frac{1}{1+x^2} dx \approx \int_{x=\tan(\gamma-\beta)}^{\tan(\gamma)} f(x)dx. \quad (52)
$$

A hidden difficulty in integrating *f* is that the sign of *x* (denoted by *s*(*x*)) appears in its expression, as

$$
f(x) := \left( (d^2 + a^2) + \frac{a^2}{x^2} \right)^{3/2} \frac{1}{1 + x^2} \quad (53)
$$

$$
= s(x) \left( (d^2 + a^2)x^2 + a^2 \right)^{3/2} \frac{1}{x^3(1+x^2)}. \quad (54)
$$

It is worth noting that integrating *f* on the angular domain [*γ* – *β*, *γ*] might result in a discontinuity in the range domain [tan(*γ* – *β*), tan(*γ*)] (if *π*/2 is contained within it). If the integration constant is not chosen appropriately, integrating *f* might therefore require to separate the domain [*γ* – *β*, *γ*] in two, and deal with the generated discontinuity, leading to unstable computations.

To avoid this pitfall, considering that *f* is even (i.e., *f*(-*x*) = *f*(*x*)
and accounting for possible change of *s*(*x*) within the integration
domain, we fix the constant values of its antiderivative *F* to make
it odd (in particular: *F*(+∞) = -*F*(-∞) = 0, corresponding to
*tan*(*θ*) = ±∞ ↔ *θ* = *π*/2 ± 0), leading to:

$$
F(x) := \int^{x} f(u) du = \tag{55}
$$

$$
s(x) \left[ d^3 \left( \acos \left( \frac{d}{V} \right) - \frac{\pi}{2} \right) - \frac{3ad^2 + a^3}{4} \log \left( \frac{(a+U)^2}{(d^2+a^2)x^2} \right) - \frac{a^2U}{2x^2} \right]
$$

with $U := \sqrt{d^2 x^2 + a^2 (1 + x^2)}$ and $V := \sqrt{(d^2 + a^2)(1 + x^2)}$.

A first expression can be obtained for our integral as

$$
\int_{\xi \in T_0} \|\xi - \eta\| d\xi = -\frac{d^3\beta}{3} + \frac{1}{3}(F(\tan(\gamma)) - F(\tan(\gamma - \beta))). \quad (56)
$$

We note that $\gamma = \pi - \gamma_0$ implies (as all angles are positive and belong to $[0, \pi]$ in our construction):

$$
\begin{align*}
\gamma - \beta &= \gamma_1, \\
\cos(\gamma) &= -\cos(\gamma_0), & \sin(\gamma) &= \sin(\gamma_0), & \tan(\gamma) &= -\tan(\gamma_0), \\
\cos(\gamma - \beta) &= \cos(\gamma_1), & \sin(\gamma - \beta) &= \sin(\gamma_1), & \tan(\gamma - \beta) &= \tan(\gamma_1).
\end{align*}
$$

Since *F* is odd, Eq. 56 can be transformed into

$$
\int_{\xi \in T_0} \|\xi - \eta\| d\xi = - \frac{d^3 \beta + F(\tan(\gamma_0)) + F(\tan(\gamma_1))}{3}. \quad (57)
$$

We further simplify the expression of *F* by expressing it w.r.t.
cosines and sines instead of tangents. For *x* = tan(*γ*<sub>*i*</sub>), and noting
$(C_i, S_i) = (\cos(\gamma_i), \sin(\gamma_i))$, we observe that

$$s(x) = s(C_i) \qquad (58)$$

$$a = \delta_i S_i, \quad \forall i = \{0; 1\} \qquad (59)$$

$$\sqrt{d^2 S_i^2 + a^2} = \sqrt{d^2 S_i^2 + \delta_i^2 S_i^2} = l_i S_i, \quad \forall i = \{0; 1\} \qquad (60)$$

$$\delta_i C_i = \epsilon_i \zeta_i, \text{ with } (\epsilon_0, \epsilon_1) = (1, -1) \qquad (61)$$

for $\zeta_i := (\eta - t_i) \cdot u_e$.

Using this, the three terms of $F$ (Eq. (55)) can be simplified into:

$$\begin{aligned}
F_3(x) &:= -s(x) \frac{3ad^2+a^3}{4} \log \left( \frac{\left(a + \sqrt{d^2x^2 + a^2(1+x^2)}\right)^2}{(d^2+a^2)x^2} \right) \\
&= -s(C_i) \frac{3ad^2+a^3}{4} \log \left( \frac{\left(a + \sqrt{d^2S_i^2/C_i^2 + a^2(1+S_i^2/C_i^2)}\right)^2 C_i^2}{(d^2+a^2)S_i^2} \right) \\
&= -s(C_i) \frac{3ad^2+a^3}{4} \log \left( \frac{\left(a|C_i| + \sqrt{d^2S_i^2 + a^2}\right)^2}{(d^2S_i^2+a^2) - a^2C_i^2} \right) \\
&= -s(C_i) \frac{3ad^2+a^3}{4} \log \left( \frac{\sqrt{d^2S_i^2 + a^2 + a|C_i|}}{\sqrt{d^2S_i^2 + a^2 - a|C_i|}} \right) \\
&= \frac{3ad^2+a^3}{4} \log \left( \frac{\sqrt{d^2S_i^2 + a^2 - aC_i}}{\sqrt{d^2S_i^2 + a^2 + aC_i}} \right) \\
&= \frac{3ad^2+a^3}{4} \log \left( \frac{\left(\sqrt{d^2S_i^2 + a^2} - aC_i\right)^2}{(d^2+a^2)S_i^2} \right) \\
&= \frac{3ad^2+a^3}{4} \log \left( \frac{(l_iS_i - \delta_iS_iC_i)^2}{(d^2+a^2)S_i^2} \right) \\
&= \epsilon_i \frac{3ad^2+a^3}{2} \log \left( \frac{l_i - \zeta_i}{\sqrt{d^2+a^2}} \right).
\end{aligned}$$

Note that, while the last step is presented extremely quickly, it is obtained by repeating the reasoning than the one detailed from line 2 to line 6 above.

$$\begin{aligned}
F_1(x) &:= s(x) d^3 \left( \operatorname{acos} \left( \frac{d}{\sqrt{(d^2 + a^2)(1 + x^2)}} \right) - \frac{\pi}{2} \right) \\
&= s(C_i) d^3 \left( \operatorname{acos} \left( \frac{d}{\sqrt{(d^2 + a^2)(1 + S_i^2/C_i^2)}} \right) - \frac{\pi}{2} \right) \\
&= s(C_i) d^3 \left( \operatorname{acos} \left( \frac{d|C_i|}{\sqrt{d^2 + a^2}} \right) - \frac{\pi}{2} \right) \\
&= d^3 \left( \operatorname{acos} \left( \frac{dC_i}{\sqrt{d^2 + a^2}} \right) - \frac{\pi}{2} \right),
\end{aligned}$$

$$\begin{aligned}
F_2(x) &:= -s(x) \frac{a^2 \sqrt{d^2 x^2 + a^2 (1 + x^2)}}{2x^2} \\
&= -s(C_i) \frac{a^2 \sqrt{d^2 S_i^2 / C_i^2 + a^2 (1 + S_i^2 / C_i^2) C_i^2}}{2S_i^2} \\
&= -s(C_i) \frac{a^2 \sqrt{d^2 S_i^2 + a^2 |C_i|}}{2S_i^2} \\
&= -\frac{a^2 \sqrt{d^2 S_i^2 + a^2 C_i}}{2S_i^2} \\
&= -\frac{a l_i \delta_i C_i}{2} \\
&= -\frac{a \epsilon_i l_i \zeta_i}{2}.
\end{aligned}$$

Focusing on $F_1$, summing up those for $\gamma_0$ and $\gamma_1$ results in:

$$\begin{aligned}
& -\frac{d^3}{3}\beta - F_1(\tan(\gamma_0)) - F_1(\tan(\gamma_1)) \\
&= -\frac{d^3}{3} \left( \beta + \operatorname{acos} \left( \frac{dC_0}{\sqrt{d^2 + a^2}} \right) - \frac{\pi}{2} + \operatorname{acos} \left( \frac{dC_1}{\sqrt{d^2 + a^2}} \right) - \frac{\pi}{2} \right) \\
&= \frac{d^3}{3} \underbrace{\left( \beta + \operatorname{acos} \left( \frac{dC_0}{\sqrt{d^2 + a^2}} \right) - \operatorname{acos} \left( \frac{dC_1}{\sqrt{d^2 + a^2}} \right) \right)}_{(\gamma_0 + \gamma_1) - \pi - \beta}, \\
&= \underbrace{\frac{d^3}{3} \left( \beta + \operatorname{acos} \left( \frac{dC_0}{\sqrt{d^2 + a^2}} \right) - \operatorname{acos} \left( \frac{dC_1}{\sqrt{d^2 + a^2}} \right) \right)}_{(\pi - \beta) - (\gamma_0 + \gamma_1)}, \\
&= \therefore \bar{\omega}_{T_0}(\eta)
\end{aligned}$$

making the unsigned, positive solid angle $\bar{\omega}_{T_0}(\eta)$ appear. Multiplying this unsigned solid angle by $s(T_0)$ and summing this expression for the three subtriangles $\{T_i\}$ results in the $\frac{d^3 \omega_t(\eta)}{3}$ term of Eq. (10) (which is true for all possible configurations of $\{s(T_i)\}_{i}$. Note that, while $d_t$ and $\omega_t(\eta)$ are signed in Eq. (10), their product is always positive as both terms have the same sign for all $\eta$.

Noting that $a$ in the first two terms above is unsigned as well, multiplying $a$ by $s(T_i)$ results in the signed $a_e^t$ term of Eq. 10, finishing the derivation of Eq. 10 by considering the contributions of the two remaining terms $F_2$ and $F_3$ for all three triangles $\{T_i\}$.
## B DERIVATION OF EQ. (18)

We first rewrite $∇\bar{ψ}$, and focus only on its components along $r_{e'}^t$, which is the rotated edge facing $v_t$ in $t$. Noting that $(u_{e_1}×u_{e_0}) = (\sin(\alpha)d_t/D_e)r_e^t + \lambda n_t$ (the component along $n_t$ being ignored when taking the dot product with $r_{e'}^t$), and using the law of sines to derive $(l_{e_0} + l_{e_1})\sin(\alpha) = \|e\|(\sin(\beta_0) + \sin(\beta_1)) = \|e\|(D_e/l_{e_0} + D_e/l_{e_1})$:

<figure><img src="image_15.png" /></figure>

$$ 
\begin{aligned}
\nabla_{\eta} \left( \int_{\xi \in t} \|\xi - \eta\| d\xi \right) (\eta) \cdot r_{e'}^t = & \\
& \frac{d_t^4}{3} \sum_{e \in t} \frac{-2\|e\|(l_{e_0} + l_{e_1})r_e^t \cdot r_{e'}^t}{((l_{e_0} + l_{e_1})^2 - \|e\|^2)(l_{e_0}l_{e_1})} \\
& + \sum_{e \in t} \frac{r_e^t \cdot r_{e'}^t}{6} \left( (3D_e^2) \log\left(\frac{l_{e_1} - \zeta_{e_1}}{l_{e_0} - \zeta_{e_0}}\right) - l_{e_1}\zeta_{e_1} + l_{e_0}\zeta_{e_0} \right) \\
& + \sum_{e \in t} \frac{a_e^t}{6} \left( (2d_t^2 + D_e^2) \left[ \frac{ue_1 - ue}{l_{e_1} - \zeta_{e_1}} - \frac{ue_0 - ue}{l_{e_0} - \zeta_{e_0}} \right] \right) \cdot r_{e'}^t \\
& + \sum_{e \in t} \frac{a_e^t}{6} \left( -u_{e_1}\zeta_{e_1} + u_{e_0}\zeta_{e_0} + (l_{e_0} - l_{e_1})ue \right) \cdot r_{e'}^t.
\end{aligned}
$$

Differentiating twice this expression gives the formula for

$$ \nabla_{\eta}^{T} \left( H_{\eta} \left( \int_{\xi \in t} \|\xi - \eta\| d\xi \right) r_{e'}^{t} \right) \text{ (required in Sec. 3.2), since, for a fixed vector } r \text{ (independent of } \eta \text{), the following holds:} $$

$$ \nabla_{\eta}^{T} \left( H_{\eta} \left( \int_{\xi \in t} \|\xi - \eta\| d\xi \right) r \right) = H_{\eta} \left( \nabla_{\eta}^{T} \left( \int_{\xi \in t} \|\xi - \eta\| d\xi \right) (\eta) \cdot r \right). $$

While the derivation is tedious (differentiating all terms twice, and refactoring them in order to find "as-small-as-possible expressions", requires several pages of text), it does actually not constitute a real challenge, and it can be even obtained using automatic differentiation tools from this point on. We omit it for conciseness.

## C SPECIAL CASES COMPUTATIONS

We detail here a list of special cases that are required for the evaluation of the harmonic/biharmonic coordinates on the cage. Those are essentially obtained by considering the general case and evaluating those formulas by limit.

**Strictly inside triangles:** For a point $\eta$ strictly inside the triangle $t$ (noting $\{e\}$ its 3D edges oriented CCW), **all formulas are well-behaved once the terms multiplied by $d_t = 0$ are removed.** We give the coordinates and their normal derivatives, as those exhibit

an interesting simplified structure:

$$ 
\begin{aligned}
\psi_t(\eta) &= 2|t| \sum_i C_{ei} \Gamma_{ti}^t(\eta) \\
\phi_{ti}^t(\eta) &= \frac{\Gamma_{ti}^t(\eta)}{2} \\
\bar{\psi}_t(\eta) &= \sum_{e \in t} \frac{a_e^t}{6} \left( \log\left(\frac{l_{e_1} - \zeta_{e_1}}{l_{e_0} - \zeta_{e_0}}\right) - l_{e_1}\zeta_{e_1} + l_{e_0}\zeta_{e_0} \right) \\
\bar{\phi}_{ti}^t(\eta) &= 0 \quad \forall i \\
\nabla \psi_t(\eta) \cdot n_t &= 1/2 \\
\nabla \phi_{ti}^t(\eta) \cdot n_t &= \frac{1}{2|t|} \det \left( n_t \times \sum_{e \in t} eC_e; e^i; n_t \right) \\
\nabla \bar{\psi}_t(\eta) \cdot n_t &= 0 \\
\nabla \bar{\phi}_{ti}^t(\eta) \cdot n_t &= \frac{1}{2} \Gamma_{ti}^t(\eta) \psi_t(\eta) + \lambda_2 \sum_{e \in t} \frac{a_e^t}{6} D_e^2 \left( [m_{e_1} - m_{e_0}] \cdot \mathbf{g}_i^t \right) \\
&\quad + \lambda_2 \sum_{e \in t} \frac{r_e^t \cdot \mathbf{g}_i^t}{6} \left( D_e^2 \log\left(\frac{l_{e_1} - \zeta_{e_1}}{l_{e_0} - \zeta_{e_0}}\right) - l_{e_1}\zeta_{e_1} + l_{e_0}\zeta_{e_0} \right) \\
&\quad + \lambda_2 \sum_{e \in t} \frac{a_e^t}{3} (\bar{D}_e \cdot \mathbf{g}_i^t) \log\left(\frac{l_{e_1} - \zeta_{e_1}}{l_{e_0} - \zeta_{e_0}}\right) \\
&\quad + \lambda_2 \sum_{e \in t} \frac{a_e^t}{6} \left( -u_{e_1}\zeta_{e_1} + u_{e_0}\zeta_{e_0} + (l_{e_0} - l_{e_1})ue \right) \cdot \mathbf{g}_i^t
\end{aligned}
$$

$C_e := \frac{1}{4\pi\|e\|} \log\left(\frac{l_{e_0}+l_{e_0}+\|e\|}{l_{e_0}+l_{e_0}-\|e\|}\right)$, and $e^i$ denoting the edge opposite to corner $i$ in triangle $t$.

**Strictly on cage vertices:** For a point $\eta = t^k$ being the $k^{th}$ corner of a triangle $t$, the coordinates are given by:

$$ 
\begin{aligned}
\psi_t(t^k) &= 2|t|C_{ek} \\
\phi_{ti}^t(t^k) &= \omega_t \delta_i^k \\
\bar{\psi}_t(t^k) &= \frac{a_e^t}{6} \left( \log\left(\frac{l_{e^k_1} - \zeta_{ek_1}}{l_{e^k_0} - \zeta_{ek_0}}\right) - l_{e^k_1}\zeta_{ek_1} + l_{e^k_0}\zeta_{ek_0} \right) \\
\bar{\phi}_{ti}^t(t^k) &= 0 \quad \forall i
\end{aligned}
$$

Note that, to compute $\phi_{ti}^t(t^k)$, one has to make sure to output the solid angle of $t$ seen from $\eta$, as $\eta$ tends to $t^k$. It is incorrect to set this value to $1/2$ (wrongly considering the point $t^k$ on the triangle), and this value should be set in order to obtain $\phi_{tk}(t^k) = \lim_{\epsilon \to 0^+} \sum_{t \in t^k} \omega_t^\epsilon(t^k)$, $\omega_t^\epsilon$ being the solid angle of $t$, whose corner has been displaced by $\epsilon$ along its normal: in the end, what is important is for $\phi_{tk}$ to be associated with the solid angle of the triangles adjacent to the considered corner, from the point of view of a 3D point inside the cage as it tends to $\eta = t^k$. In practice, we handle this case by limit in our implementation, by displacing $t^k$ slightly in order to make sure that $\eta$ does not lie on the considered triangles, as the solid angle does not depend on its exact placement (but only on its one-ring edges).
# D ALTERNATE CONSTRUCTIONS FOR CONSTRAINED BIHARMONIC DEFORMATIONS

## D.1 (1,2)-constrained biharmonic deformations

Since

$$f(\eta) = \sum_{i \in \mathcal{V}} \phi_i(\eta) v'_i + \sum_{j \in \mathcal{T}} \psi_j(\eta) n'_j + \sum_{i \in \mathcal{V}} \bar{\phi}_i(\eta) L'_i + \sum_{j \in \mathcal{T}} \bar{\psi}_j(\eta) D'_j$$

$$\frac{\partial f}{\partial n}(\eta) = \sum_{i \in \mathcal{V}} (\nabla \phi_i(\eta) \cdot n) v'_i + \sum_{j \in \mathcal{T}} (\nabla \psi_j(\eta) \cdot n) n'_j + \sum_{i \in \mathcal{V}} (\nabla \bar{\phi}_i(\eta) \cdot n) L'_i + \sum_{j \in \mathcal{T}} (\nabla \bar{\psi}_j(\eta) \cdot n) D'_j,$$

enforcing $f(p_k) = \sum_i \Gamma_i(p_k) v'_i$ (interpolation of position at sample $p_k$) and $\frac{\partial f}{\partial n_t}(q_k) = n'_t$ at sample $q_k$ of triangle $t$ (interpolation of the normal derivative on triangle $t$) leads to the following constraints:

$$\sum_{i} \Gamma_i(p_k) v'_i = \sum_{i \in \mathcal{V}} \phi_i(p_k) v'_i + \sum_{j \in \mathcal{T}} \psi_j(p_k) n'_j + \sum_{i \in \mathcal{V}} \bar{\phi}_i(p_k) L'_i + \sum_{j \in \mathcal{T}} \bar{\psi}_j(p_k) D'_j \\ n'_t = \sum_{i \in \mathcal{V}} (\nabla \phi_i(q_k) \cdot n_t) v'_i + \sum_{j \in \mathcal{T}} (\nabla \psi_j(q_k) \cdot n_t) n'_j + \sum_{i \in \mathcal{V}} (\nabla \bar{\phi}_i(q_k) \cdot n_t) L'_i + \sum_{j \in \mathcal{T}} (\nabla \bar{\psi}_j(q_k) \cdot n_t) D'_j.$$

Setting those constraints on all vertices and triangles leads to the following constraint in matrix form:

$$\left\{ \begin{array}{l} M_V V' = \Phi_V V' + \Psi_V N' + \bar{\Phi}_V L' + \bar{\Psi}_V D' \\ M_N N' = \partial_n \Phi_N V' + \partial_n \Psi_N N' + \partial_n \bar{\Phi}_N L' + \partial_n \bar{\Psi}_N D' \\ \left( \begin{array}{cc} \partial_n \bar{\Phi}_N & \partial_n \bar{\Psi}_N \\ \bar{\Phi}_V & \bar{\Psi}_V \end{array} \right) \begin{pmatrix} L' \\ D' \end{pmatrix} = \left( \begin{array}{cc} -\partial_n \Phi_N & M_N - \partial_n \Psi_N \\ M_V - \Phi_V & -\Psi_V \end{array} \right) \begin{pmatrix} V' \\ N' \end{pmatrix} \Rightarrow \\ \left( \begin{array}{c} L' \\ D' \end{array} \right) = A^{\dagger} B \left( \begin{array}{c} V' \\ N' \end{array} \right) =: \left( \begin{array}{cc} C_{11} & C_{12} \\ C_{21} & C_{22} \end{array} \right) \left( \begin{array}{c} V' \\ N' \end{array} \right). \end{array} \right.$$

Given deformation parameters ($V'$, $N'$), the matrix constraints give auxiliary parameters ($L'$, $D'$) for our biharmonic deformation function to match (in the least-squares sense) the Dirichlet and Neumann constraints at sampled locations ($\{p_k\}$, $\{q_k\}$).

Considering an input mesh $\mathcal{M}$ encoded using our unconstrained biharmonic coordinates (see Eq. (19)) as (in matrix form):

$$\mathcal{M}' = \Phi_{\mathcal{M}} V' + \Psi_{\mathcal{M}} N' + \bar{\Phi}_{\mathcal{M}} L' + \bar{\Psi}_{\mathcal{M}} D',$$

plugging our matrix constraints leads to the following expression for the (1,2)-regularized biharmonic coordinates (see Eq. (20)):

$$\mathcal{M}' = \underbrace{[\Phi_{\mathcal{M}} + \bar{\Phi}_{\mathcal{M}} C_{11} + \bar{\Psi}_{\mathcal{M}} C_{21}]}_{=:\alpha_{\mathcal{M}}^{(1,2)}} V' + \underbrace{[\Psi_{\mathcal{M}} + \bar{\Phi}_{\mathcal{M}} C_{12} + \bar{\Psi}_{\mathcal{M}} C_{22}]}_{=:\beta_{\mathcal{M}}^{(1,2)}} N'$$

## D.2 (1,3)- constrained biharmonic deformations, as in Weber et al. [2012]

Since

$$\Delta f(\eta) = \sum_{i \in \mathcal{V}} \phi_i(\eta) L'_i + \sum_{j \in \mathcal{T}} \psi_j(\eta) D'_j, \quad (62)$$

one can evaluate this expression at points $p_k \in \partial\Omega$, and derive a first constraint by enforcing linear interpolation of $L'$ at $p_k$ (averaging $L'$ at its corners):

$$\sum_{i} \Gamma_i(p_k) L'_i = \Delta f(p_k) = \sum_{i \in \mathcal{V}} \phi_i(p_k) L'_i + \sum_{j \in \mathcal{T}} \psi_j(p_k) D'_j.$$

Considering this constraint for a collection of samples $\{p_k\}$ leads to the following matrix constraint:

$$M_L L' = \Phi_L L' + \Psi_L D' \Rightarrow \\ D' = \Psi_L^\dagger (M_L - \Phi_L) L' =: AL', \quad (63)$$

$M_L$ denoting a mass matrix averaging, for each line $k$ corresponding to the $k^{th}$ sample, the quantities defined on the vertices (i.e., each row of $M_L$ has three non-zeros entries, that sum to 1, at the columns of the vertices indexed by the triangle on which $p_k$ lies).

Considering similarly the interpolation of the Dirichlet constraints at sparse locations $q_k \in \partial\Omega$ leads similarly to

$$\sum_{i} \Gamma_i(q_k) v'_i = \sum_{i \in \mathcal{V}} \phi_i(q_k) v'_i + \sum_{j \in \mathcal{T}} \psi_j(q_k) n'_j + \sum_{i \in \mathcal{V}} \bar{\phi}_i(q_k) L'_i + \sum_{j \in \mathcal{T}} \bar{\psi}_j(q_k) D'_j.$$

Stacking all constraints for a collection of sampled constraints $\{q_k\}$ leads to the following matrix constraint:

$$\begin{aligned} M_V V' &= \Phi_V V' + \Psi_V N' + \bar{\Phi}_V L' + \bar{\Psi}_V D' \\ &= \Phi_V V' + \Psi_V N' + [\bar{\Phi}_V + \bar{\Psi}_V A] L' \Rightarrow \\ L' &= [\bar{\Phi}_V + \bar{\Psi}_V A]^\dagger [(I - \Phi_V) V' - \Psi_V N'] =: BV' + CN'. \end{aligned} \quad (64)$$

Considering an input mesh $\mathcal{M}$ encoded using our unconstrained biharmonic coordinates (see Eq. (19)) as (in matrix form):

$$\mathcal{M}' = \Phi_{\mathcal{M}} V' + \Psi_{\mathcal{M}} N' + \bar{\Phi}_{\mathcal{M}} L' + \bar{\Psi}_{\mathcal{M}} D',$$

plugging our matrix constraints (Eqs. (63),(64)) leads to the following expression for Weber's (1,3)-regularized biharmonic coordinates (see Eq. (20)):

$$\mathcal{M}' = \underbrace{[\Phi_{\mathcal{M}} + \bar{\Phi}_{\mathcal{M}} B + \bar{\Psi}_{\mathcal{M}} AB]}_{=:\alpha_{\mathcal{M}}^{(1,3)}-Weber et al.} V' + \underbrace{[\Psi_{\mathcal{M}} + \bar{\Phi}_{\mathcal{M}} C + \bar{\Psi}_{\mathcal{M}} AC]}_{=:\beta_{\mathcal{M}}^{(1,3)}-Weber et al.} N'. \quad (65)$$

The main differences between our formulation (see Section 4.1.1) and Weber's formulation introduced in 2D in [Weber et al. 2012], is that Weber and colleagues dedicate the setting of $D'$ to enforce the Laplacian regularization constraints (Eq. (63)) before setting $L'$ to enforce the position regularization constraints (Eq. (64)), while we jointly optimize $L'$ and $D'$ to enforce both constraints sets. This minor difference in our construction results however in much better-behaved deformations, and offers means of balancing positional versus Laplacian regularizations (see Section 4.1.2).