# A Robust Hessian-based Trust Region Algorithm for Spherical Conformal
Parameterizations

Zhong-Heng Tan School of Mathematics and Shing-Tung Yau Center, Southeast
University, Nanjing 211189; Nanjing Center for Applied Mathematics, Nanjing
211135, People’s Republic of China. Tiexiang Li Corresponding author. School
of Mathematics and Shing-Tung Yau Center, Southeast University, Nanjing
211189; Nanjing Center for Applied Mathematics, Nanjing 211135, People’s
Republic of China. Wen-Wei Lin Nanjing Center for Applied Mathematics, Nanjing
211135, People’s Republic of China; Department of Applied Mathematics, Yang
Ming Chiao Tung University, Hsinchu 300, Taiwan. Shing-Tung Yau Yau
Mathematical Sciences Center, Tsinghua University, Beijing 100084, China.

###### Abstract

Surface parameterizations are widely applied in computer graphics, medical
imaging and transformation optics. In this paper, we rigorously derive the
gradient vector and Hessian matrix of the discrete conformal energy for
spherical conformal parameterizations of simply connected closed surfaces of
genus-$0$. In addition, we give the sparsity structure of the Hessian matrix,
which leads to a robust Hessian-based trust region algorithm for the
computation of spherical conformal maps. Numerical experiments demonstrate the
local quadratic convergence of the proposed algorithm with low conformal
distortions. We subsequently propose an application of our method to surface
registrations that still maintains local quadratic convergence.

Key words— spherical conformal parameterization, conformal energy
minimization, Riemann surfaces of genus-0, Hessian matrix, local quadratic
convergence

MSC(2020) 49Q10, 52C26, 65D18, 65F05, 68U05

##  1 Introduction

The Poincaré-Klein-Koebe theorem is a fundamental theorem in Riemann geometry
that states that a simply connected Riemann surface $\mathcal{M}$ is
conformally equivalent to either a unit sphere $\mathbb{S}^{2}$, a complex
plane $\mathbb{C}$ or a unit disk $\mathbb{D}$. From a numerical point of
view, for a given discrete triangular mesh of $\mathcal{M}$, we require that
the conformal parameterization between $\mathcal{M}$ and the canonical shape
is a bijective map from $\mathcal{M}$ to $\mathbb{S}^{2}$ (Riemann surface of
genus-zero) or $\mathbb{D}$ (Riemann surface with a single boundary) while
minimizing the total angle distortion induced by the Dirichlet energy [12, 16,
49]. Conformal parameterizations, also known as angle-preserving
parameterizations, preserve the intersection angle of two arbitrary
intersecting curves on the surface up to the map. In other words, conformal
parameterizations preserve the local shapes of the surfaces. The development
of digital 3D object technologies has enabled the representation of smooth
surfaces in the real world through high-resolution meshes on computers. High-
resolution meshes have the ability to characterize the intricate geometrical
structure of surfaces. Nevertheless, high resolution also makes computations
and processes more challenging on complicated meshes. Recently, surface
parameterizations have received increasing attention due to their ability to
transform intricate geometries into simply shaped regions. It is clear that
the parameterization of a surface is not necessarily unique, even when the
target region is given. Conformal parameterization is one of the most commonly
used parameterization methods. As a result, it is used in several fields, such
as medical imaging [27, 39, 31, 42], texture mapping [32, 12], and
transformation optics [45, 22].

A detailed overview of surface conformal parameterizations can be found in
classical surveys [25, 9, 33, 10, 30]. Numerous methods for conformal
parameterizations have been proposed, including boundary first flattening
[41], discrete conformal equivalence [26], harmonic energy minimization [38]
and solving the Laplacian-Beltrami equation with a particular Dirac delta
function as the right-hand term [32]. In earlier years, based on the time-flow
technique, Jin et al. [36] and Yang et al. [46] proposed the discrete Ricci
flow with conformal circle packing metric, which is a negative gradient flow
of some convex energy and can be accelerated by the Newton method. However,
the evaluations of the related coefficients at each iterative step are
somewhat complicated due to the circle packing metric. Huang et al. [34]
proposed the quasi-implicit Euler method (QIEM) in view of the nonlinear heat
diffusion process with normalization on $\mathbb{S}^{2}$. However, the
convergence can be very slow, and there is thus far no theory to support its
convergence. With the great development of GPU-based computation in recent
years, parallelizable algorithms were subsequently proposed [14].

In 2015, Choi et al. [15] applied the quasi-conformal approach to spherical
parameterization and proposed the FLASH algorithm, which produces the
composition of the quasi-conformal maps with the same Beltrami coefficient to
obtain the conformal map. In 2019, Yueh et al. [48] proposed a north-south
hemisphere alternating iteration, called spherical conformal energy
minimization (SCEM), which alternatively maps the north and south poles to
infinity and fixes the corresponding hemisphere while updating the other
hemisphere. From numerical experiments, both FLASH and SCEM are satisfactory
with high accuracy and effectiveness compared to those of the other previously
mentioned methods. FLASH is a direct method that solves two double-sized
linear systems compared to SCEM, which utilizes stereographic projection to
transform the unit sphere $\mathbb{S}^{2}$ to the extended complex plane
$\overline{\mathbb{C}}$. The resulting map, which is composed of the inverse
stereographic projection and the computed conformal map, is the ideal
conformal parameterization from the given Riemann surface of genus zero to
$\mathbb{S}^{2}$. In addition, it simplifies the spherical constraint problem
in $\mathbb{R}^{3}$ into an unconstrained problem in $\overline{\mathbb{C}}$,
making it highly efficient. However, the stereographic projection maps the
north pole of $\mathbb{S}^{2}$ to infinity and others to $\mathbb{C}$.
Therefore, the computational error near the north pole becomes relatively
large in practical applications. To mitigate this error, the north-south
hemisphere alternatingly iterative SCEM reduces the conformal distortion near
the poles and has numerically sublinear convergence and asymptotically
$R$-linear convergence. However, larger computational errors in conformal
quantities are transferred to the junction of the equator. From a theoretical
point of view, both FLASH and SCEM consider minimizing the Dirichlet energy on
the extended complex plane $\overline{\mathbb{C}}$. From a numerical point of
view, both FLASH and SCEM may get larger angle distortions near the junction
region of the computational conformal map because they are connected by two
submaps via stereographic projection.

Different from FLASH and SCEM, in this paper, we solve the spherical conformal
energy optimization problem to obtain the conformal map by directly employing
the spherical coordinates for the representation of the conformal energy.
Unlike the stereographic projection, the spherical coordinate representation
is not a conformal map. We derive the related gradient vector and Hessian
matrix and develop a robust Hessian-based trust region (HBTR) algorithm that
has local quadratic convergence. The main contributions of this paper can be
divided into three folds.

  * •

We use the spherical coordinates to represent the conformal energy on
$\mathbb{S}^{2}$ and develop the HBTR algorithm to directly minimize the
discrete conformal energy, which is described as the difference between the
discrete Dirichlet energy on $\mathbb{S}^{2}$ and the image area of the
conformal map. In Section 5.3, numerical experiments demonstrate that the
discrete conformal energy, mean and standard deviation (SD) are reduced to
$1/4$ and $1/2$, respectively, when the mesh size is refined by a half.

  * •

We give the explicit derivation of the gradient vector and Hessian matrix and
present the significant sparsity of the Hessian matrix. This benefits the
development of fast computations for calculating the Newton iterations.
Combined with the trust region technique, we thus propose the robust HBTR
algorithm.

  * •

Numerical experiments and comparisons with existing algorithms demonstrate the
advantages of our method in terms of conformality and robustness. The
computational cost of the HBTR algorithm is of the same order of magnitude as
that of FLASH and SCEM but slightly more expensive for most models due to
spending a high percentage of time selecting the convergence region.
Furthermore, the bijectivity can almost be guaranteed because the conformal
energy is expressed by the unity of spherical coordinates. For few meshes
resulting in nonbijective maps, folding can be easily removed by performing a
postprocessing algorithm [24].

This paper is organized as follows. The Section 2 provides a brief review of
the conformal map and conformal energy. Then, in Section 3 we present the
theoretical derivation of the gradient vector and Hessian matrix. In Section
4, we describe the proposed algorithm for spherical conformal parameterization
that uses the sparsity property of the Hessian matrix. The numerical
performance and comparison with other methods are presented in Section 5. In
Section 6, we demonstrate an application of the algorithm to surface
registration. A concluding mark is given in Section 7.

The frequently used notations in this paper are listed here. Bold letters,
e.g., $\mathbf{a},\mathbf{s}$ denote vectors. $\mathbf{a}_{i}$ denotes the
$i$-th entry of $\mathbf{a}$. $\mathbf{1}_{m\times n}$ denotes the $m\times n$
matrix of all ones. The notation without subscript $\mathbf{1}$ denotes the
vector of all ones with proper dimension. $\mathbf{e}_{i}$ denotes the $i$-th
column of the identity matrix with the proper dimension.
$\big{[}A\big{]}_{ij}$ denotes the $(i,j)$-th entry of matrix $A$.
$\mathrm{diag}(\mathbf{a})$ denotes the diagonal matrix with the $(i,i)$-th
entry being $\mathbf{a}_{i}$. $[v_{i},v_{j}]$ denotes the edge formed by
$v_{i}$ and $v_{j}$. $[v_{i},v_{j},v_{k}]$ denotes the triangle formed by
$v_{i}$, $v_{j}$ and $v_{k}$, and $\big{|}[v_{i},v_{j},v_{k}]\big{|}$ denotes
the area of that triangle. The other notations are defined wherever they
appear.

##  2 Conformal map and Conformal Energy

In this section, we briefly review the continuous conformal map and the
conformal energy. Readers can refer to [28, 29, 35, 40, 32] for more details.

Let $\mathcal{M}=\mathbf{r}(u,v)$ be a surface and $f$ be a continuous and
bijective vector-valued map on $\mathcal{M}$, which maps $\mathcal{M}$ to
another surface $\widetilde{\mathcal{M}}=\tilde{\mathbf{r}}(u,v)$ with
$(u,v)\in\Omega\subset\mathbb{R}^{2}$. Let
$\gamma_{1},\gamma_{2}\subset\mathcal{M}$ be two arbitrary curves intersecting
at a point. Then, we call $f$ conformal if the intersecting angle of
$\gamma_{1}$ and $\gamma_{2}$ equals that of $f(\gamma_{1})$ and
$f(\gamma_{2})$ in $\widetilde{\mathcal{M}}$. Equivalently, $f$ is conformal
if and only if the first fundamental forms $I(u,v)$ and $\tilde{I}(u,v)$ of
$\mathcal{M}$ and $\widetilde{\mathcal{M}}$ with respect to $(u,v)$, i.e.,

$\displaystyle
I(u,v)=\begin{bmatrix}\langle\mathbf{r}_{u},\mathbf{r}_{u}\rangle&\langle\mathbf{r}_{u},\mathbf{r}_{v}\rangle\\\
\langle\mathbf{r}_{v},\mathbf{r}_{u}\rangle&\langle\mathbf{r}_{v},\mathbf{r}_{v}\rangle\\\
\end{bmatrix},\quad\tilde{I}(u,v)=\begin{bmatrix}\langle\tilde{\mathbf{r}}_{u},\tilde{\mathbf{r}}_{u}\rangle&\langle\tilde{\mathbf{r}}_{u},\tilde{\mathbf{r}}_{v}\rangle\\\
\langle\tilde{\mathbf{r}}_{v},\tilde{\mathbf{r}}_{u}\rangle&\langle\tilde{\mathbf{r}}_{v},\tilde{\mathbf{r}}_{v}\rangle\\\
\end{bmatrix},$

satisfy

$\displaystyle\tilde{I}(u,v)=\eta(u,v)I(u,v),$

where $\mathbf{r}_{u}:=\frac{\partial\mathbf{r}}{\partial u}$,
$\mathbf{r}_{v}:=\frac{\partial\mathbf{r}}{\partial v}$ and $\eta(u,v)$ is a
positive scalar function on $\Omega$. Then, the conformal energy functional
[35] of $f$ is defined as

$\displaystyle
E_{C}(f)=\frac{1}{2}\int_{\mathcal{M}}\|\nabla_{\mathcal{M}}f\|_{F}^{2}\mathrm{d}\sigma-\mathcal{A}(f),$
(2.1)

where $\nabla_{\mathcal{M}}$ is the tangential gradient, $\mathrm{d}\sigma$ is
the area element on $\mathcal{M}$ and $\mathcal{A}(f)$ is the area of the
image surface $\widetilde{\mathcal{M}}=f(\mathcal{M})$. Let us note that the
first term in (2.1) is the Dirichlet energy functional of $f$. It has been
proven that [35, 40]

  * •

$E_{C}(f)\geq 0$,

  * •

$E_{C}(f)=0$ if and only if $f$ is conformal.

From these properties, it is reasonable to adopt conformal energy minimization
(CEM) to obtain the conformal map. In the continuous scheme, if the target
surface $\widetilde{\mathcal{M}}$ is given, its area remains constant. Thus, a
method to solve (2.1) is to minimize the Dirichlet energy, which is a
quadratic functional. In the discrete scheme, the area is actually not a
constant. Introducing the discrete area term typically yields a better
conformal parameterization. It is easy to verify that the optimization problem
(2.1) has a trivial solution $f=constant$. In this case, all vertices shrink
into a point, and the area of the formed region becomes $0$, which violates
the requirement of $f$ being bijective. We refer to the trivial solution as
’degeneration’. Directly minimizing the Dirichlet energy without any
additional constraints often leads to degeneration in practical computations.
The area term is a natural penalty for the parameterization and hence can
greatly weaken the degeneration during iteration [35]. On the other hand,
introducing the area term can further decrease the conformal distortion. For
disk parameterizations, Yueh et al.[49] proposed a disk conformal
parameterization algorithm through Dirichlet energy minimization, while a
novel algorithm minimizing the disk conformal energy was later proposed by Kuo
et al. [37], in which they derived a particular and simple area expression of
a discrete disk represented by polar coordinates to design a fast algorithm.
As a consequence, conformal energy minimization [37] gives significantly less
conformal energy and angle distortion compared to those of the Dirichlet
energy minimization [49].

For the spherical parameterization, it is clearly feasible to introduce the
area term. However, it is a significant challenge that the area of the
discrete sphere, which is the sum of areas of all triangles formed by
vertices, cannot be simply expressed as in the disk case, even if both
vertices are of unit length. Simultaneously, the gradient vector and Hessian
matrix of the conformal energy also have complicated expressions. In the next
section, we carefully derive the expression of the conformal energy for the
spherical parameterization with spherical coordinates, along with its gradient
vector and Hessian matrix. The Hessian matrix has a special sparsity structure
as the Laplacian matrix, which guides the fast construction and computation
associated with it.

##  3 Discrete Conformal Energy on Closed Surfaces

Let $M$ be a discrete closed surface of genus-$0$ composed of triangles,
called a triangulation. Given that $M$ has $n$ vertices, we denote the vertex
set, edge set and triangle face set of $M$ by
$\mathcal{V}(M)=\\{v_{i},i=1,2,\cdots,n\\}$,
$\mathcal{E}(M)=\\{[v_{i},v_{j}]\\}$ and
$\mathcal{F}(M)=\\{T_{ijk}:=[v_{i},v_{j},v_{k}]\\},$ respectively.

In this paper, we aim to find a piecewise linear map $f$ that maps surface $M$
to a discrete unit sphere $\mathbb{S}^{2}$ conformally, in which all of the
vertices are unit vectors. Every point in a triangle can be represented as the
convex combination of its vertices via barycentric coordinates,

$\displaystyle
f(\hat{v})=\lambda_{i}(\hat{v})f(v_{i})+\lambda_{j}(\hat{v})f(v_{j})+\lambda_{k}(\hat{v})f(v_{k}),\text{
for }\hat{v}\in T_{ijk},$ (3.1)

where

$\displaystyle\lambda_{i}(\hat{v})=\frac{|[\hat{v},v_{j},v_{k}]|}{|T_{ijk}|},\quad\lambda_{j}(\hat{v})=\frac{|[v_{i},\hat{v},v_{k}]|}{|T_{ijk}|},\quad\lambda_{k}(\hat{v})=\frac{|[v_{i},v_{j},\hat{v}]|}{|T_{ijk}|}.$

After $f(v_{i}),i=1,2,\cdots,n$ is obtained, the discrete sphere is also
formed. Hence, we compute
$\mathbf{f}=[\mathbf{f}_{1}^{\top},\mathbf{f}_{2}^{\top},\cdots,\mathbf{f}_{n}^{\top}]^{\top}\in\mathbb{R}^{n\times
3}$ with $\mathbf{f}_{i}=f(v_{i})\in\mathbb{R}^{1\times 3}$ for the given
triangulation of $M$.

From the perspective of the conformal energy (2.1), the conformal map $f$ is
obtained by solving the optimization problem

$\displaystyle\mathbf{f}^{*}=\operatornamewithlimits{argmin}_{\mathbf{f}:\|\mathbf{f}_{i}\|=1,i=1,2,\cdots,n}E_{C}(f),$
(3.2)

where the discrete conformal energy $E_{C}(f)$ of $f$ on $M$ is given by

$\displaystyle E_{C}(f)=\frac{1}{2}\langle
L\mathbf{f},\mathbf{f}\rangle-A(f),$ (3.3)

in which $L$ is the Laplacian matrix defined as

$\displaystyle\big{[}L\big{]}_{ij}=\begin{cases}-w_{ij}&\text{if}~{}i\neq
j,~{}[v_{i},v_{j}]\in\mathcal{E}(M),\\\
\sum_{k\in\mathcal{N}(i)}w_{ik}&\text{if}~{}i=j,\\\
0&\text{if}~{}[v_{i},v_{j}]\notin\mathcal{E}(M),\end{cases}$ (3.4)

with cotangent weights $w_{ij}=\frac{1}{2}(\cot\alpha_{ij}+\cot\alpha_{ji})$.
Here, $\alpha_{ij}$ and $\alpha_{ji}$ are the angles opposite to the edge
$[v_{i},v_{j}]$ 111A closed surface has no boundary. Therefore, an edge must
correspond to $2$ opposite angles. and $\mathcal{N}(i)$ is the index set of
the adjacent vertices $v_{i}$, as shown in fig. 1. Readers can refer to [29,
Chap. 24] for a detailed derivation of the discrete Dirichlet energy on
triangulation. $A(f)$ in (3.3) is the area of the triangulation $f(M)$, which
can be written as the sum of areas of all triangles in $f(M)$, that is,
$A(f)=\sum_{T_{ijk}\in\mathcal{F}(M)}|f(T_{ijk})|.$

Let $v_{ij}=v_{i}-v_{j}$ and $\mathbf{f}_{ij}=\mathbf{f}_{i}-\mathbf{f}_{j}$
for simplicity, $\alpha_{ij}(f)$ and $\alpha_{ki}(f)$ be the angle,
respectively, opposite to the edges $[\mathbf{f}_{i},\mathbf{f}_{j}]$ and
$[\mathbf{f}_{k},\mathbf{f}_{i}]\in\mathcal{E}\big{(}f(M)\big{)}$, And let
$L(f)$ be a Laplacian matrix with respect to $f$ as $L$, by replacing the
cotangent weights as
$w_{ij}(f)=\frac{1}{2}(\cot\alpha_{ij}(f)+\cot\alpha_{ji}(f))$ in the target
sphere $\mathbb{S}^{2}$. $A(f)$ and its gradient with respect to $\mathbf{f}$
can be represented as in the following lemma.

######  Lemma 1.

The area of the image of $f$ and its gradient can be represented as

$\displaystyle A(f)=\frac{1}{2}\langle L(f)\mathbf{f},\mathbf{f}\rangle,$
$\displaystyle\nabla_{\mathbf{f}}A(f)=L(f)\mathbf{f}.$

###### Proof.

The area of a triangle $f(T_{ijk})\in\mathcal{F}(f(M))$ can be calculated by

$\displaystyle|f(T_{ijk})|$
$\displaystyle=\frac{1}{2}\|\mathbf{f}_{ij}\times\mathbf{f}_{jk}\|=\frac{\|\mathbf{f}_{ij}\times\mathbf{f}_{jk}\|^{2}}{4|f(T_{ijk})|}$
$\displaystyle=\frac{\langle\mathbf{f}_{jk}\times(\mathbf{f}_{ij}\times\mathbf{f}_{jk}),\mathbf{f}_{ij}\rangle}{4|f(T_{ijk})|}=\frac{\langle\mathbf{f}_{jk}\times(\mathbf{f}_{ki}\times\mathbf{f}_{ij}),\mathbf{f}_{ij}\rangle}{4|f(T_{ijk})|}$
(3.5)
$\displaystyle=\frac{\big{\langle}\langle\mathbf{f}_{ij},\mathbf{f}_{jk}\rangle\mathbf{f}_{ki}-\langle\mathbf{f}_{jk},\mathbf{f}_{ki}\rangle\mathbf{f}_{ij},\mathbf{f}_{ij}\big{\rangle}}{4|f(T_{ijk})|}$
(3.6)
$\displaystyle=-\frac{\langle\mathbf{f}_{ki},\mathbf{f}_{ij}\rangle\langle\mathbf{f}_{jk},\mathbf{f}_{jk}\rangle+\langle\mathbf{f}_{ij},\mathbf{f}_{jk}\rangle\langle\mathbf{f}_{ki},\mathbf{f}_{ki}\rangle+\langle\mathbf{f}_{jk},\mathbf{f}_{ki}\rangle\langle\mathbf{f}_{ij},\mathbf{f}_{ij}\rangle}{8|f(T_{ijk})|}$
$\displaystyle=\frac{1}{4}\big{(}\cot\alpha_{jk}(f)\|\mathbf{f}_{jk}\|^{2}+\cot\alpha_{ki}(f)\|\mathbf{f}_{ki}\|^{2}+\cot\alpha_{ij}(f)\|\mathbf{f}_{ij}\|^{2}\big{)},$
(3.7)

In the last equation (3.7), we use the formula
$\cot\alpha_{ij}(f)=-\frac{\langle\mathbf{f}_{jk},\mathbf{f}_{ki}\rangle}{2|T_{ijk}|}$.
Hence, the area $A(f)$ can be represented as

$\displaystyle A(f)$
$\displaystyle=\sum_{T_{ijk}\in\mathcal{F}(M)}|f(T_{ijk})|$
$\displaystyle=\frac{1}{4}\sum_{T_{ijk}\in\mathcal{F}(M)}\big{(}\cot\alpha_{jk}(f)\|\mathbf{f}_{jk}\|^{2}+\cot\alpha_{ki}(f)\|\mathbf{f}_{ki}\|^{2}+\cot\alpha_{ij}(f)\|\mathbf{f}_{ij}\|^{2}\big{)}$
$\displaystyle=\frac{1}{2}\sum_{ij:[v_{i},v_{j}]\in\mathcal{E}(M)}w_{ij}(f)\|\mathbf{f}_{ij}\|^{2}=\frac{1}{2}\langle
L(f)\mathbf{f},\mathbf{f}\rangle.$

Then, by using the derivation of (3.5)-(3.6), we have

$\displaystyle\frac{\partial|f(T_{ijk})|}{\partial\mathbf{f}_{i}}=$
$\displaystyle\frac{2\mathbf{f}_{jk}\times(\mathbf{f}_{ij}\times\mathbf{f}_{jk})}{8|f(T_{ijk})|}=\frac{\langle\mathbf{f}_{ij},\mathbf{f}_{jk}\rangle\mathbf{f}_{ki}-\langle\mathbf{f}_{jk},\mathbf{f}_{ki}\rangle\mathbf{f}_{ij}}{4|f(T_{ijk})|}$
$\displaystyle=$
$\displaystyle\frac{1}{2}\big{(}\cot\alpha_{ij}(f)\mathbf{f}_{ij}-\cot\alpha_{ki}(f)\mathbf{f}_{ki}\big{)}\in\mathbb{R}^{1\times
3},$ (3.8)

Hence, the derivative of $A(f)$ with respect to $\mathbf{f}_{i}$ is

$\displaystyle\frac{\partial
A(f)}{\partial\mathbf{f}_{i}}=\sum_{T_{ijk}\in\mathcal{F}(M)}\frac{\partial|f(T_{ijk})|}{\partial\mathbf{f}_{i}}=\frac{1}{2}\sum_{T_{ijk}\in\mathcal{F}(M)}\big{(}\cot\alpha_{ij}(f)\mathbf{f}_{ij}-\cot\alpha_{ki}(f)\mathbf{f}_{ki}\big{)}=\mathbf{e}_{i}^{\top}L(f)\mathbf{f},$

for $i=1,2,\cdots,n$. This completes the proof. ∎

(a) Opposite angles

(b) Adjacent vertices indices

Figure 1: Illustrations of opposite angles and adjacent vertices indices set.
(a) $\alpha_{ij}$ and $\alpha_{ji}$ are the pair of opposite angles of edge
$[v_{i},v_{j}]$; $\alpha_{jk}$ and $\alpha_{k^{\prime}j}$ are the opposite
angles of edge $[v_{j},v_{k}]$ in $T_{ijk}$ and $[v_{k^{\prime}},v_{j}]$ in
$T_{k^{\prime}ji}$, respectively. (b) $v_{j_{1}},\cdots,v_{j_{6}}$ are the
adjacent vertices of $v_{i}$. Hence,
$\mathcal{N}(i)=\\{j_{1},j_{2},\cdots,j_{6}\\}$.

From (3.3) and Lemma 1, we have

$\displaystyle E_{C}(f)=\frac{1}{2}\langle
D(f)\mathbf{f},\mathbf{f}\rangle,\quad\text{ with }D(f)=L-L(f),$ (3.9)
$\displaystyle\nabla_{\mathbf{f}}E_{C}(f)=D(f)\mathbf{f}.$ (3.10)

Clearly, $D(f)$ has the same sparsity structure as $L$ and $L(f)$, and the
corresponding cotangent weight is
$\big{[}D(f)\big{]}_{ij}:=\tilde{w}_{ij}=w_{ij}-w_{ij}(f)$.

######  Remark 1.

If $f^{*}$ is a conformal map, the inner angles of the triangles in $M$ are
identical to those in $f^{*}(M)$. We have $w_{ij}(f^{*})=w_{ij}$ and
$D(f)=\mathbf{0}$, which leads to $E_{C}(f^{*})=0$.

Let $\mathbf{f}=[\mathbf{x},\mathbf{y},\mathbf{z}]$ and
${\text{vec}}(\mathbf{f})=[\mathbf{x}^{\top},\mathbf{y}^{\top},\mathbf{z}^{\top}]^{\top}$
be its vectorization. Since the vertices of $\mathbf{f}$ are on the sphere, we
adopt the spherical coordinates to represent all vertices by

$\mathbf{x}=\cos\bm{\theta}\odot\sin{\bm{\phi}},\quad\mathbf{y}=\sin\bm{\theta}\odot\sin{\bm{\phi}},\quad\mathbf{z}=\cos{\bm{\phi}},$
(3.11a) and let
$\mathbf{u}=\cos\bm{\theta}\odot\cos{{\bm{\phi}}},\quad\mathbf{v}=\sin\bm{\theta}\odot\cos{{\bm{\phi}}},\quad\mathbf{w}=\sin{{\bm{\phi}}},$
(3.11b)

where $\bm{\theta},{{\bm{\phi}}}\in\mathbb{R}^{n}$ are the azimuth and
elevation angle vectors of the corresponding vertices of $\mathbf{f}$,
respectively, and $\odot$ is the Hadamard product. Similar to the
stereographic projection in [48, 15], the spherical coordinate representation
eliminates the unit-length constraint and transforms the Cartesian coordinate
$[\mathbf{x},\mathbf{y},\mathbf{z}]$ into the spherical coordinate
$[\bm{\theta},{{\bm{\phi}}}]$, reducing the input variable scale by $1/3$.
Furthermore, the spherical coordinate projects the whole sphere into a bounded
region $[0,2\pi]\times[0,\pi]$, which avoids the computational error near the
north pole [15] or the equator [48] by using the stereographic projection. In
the further discussion, we denote that
$\nabla:=[\nabla_{\bm{\theta}}^{\top},\nabla_{{\bm{\phi}}}^{\top}]^{\top}$
with respect to $(\bm{\theta},{{\bm{\phi}}})$ unless a special illustration in
the rest of the paper. The gradient of any scalar $a$ and Jacobi matrix of any
vector $\mathbf{a}$ with respect to $(\bm{\theta},{{\bm{\phi}}})$ are
represented as

$\nabla a=\left[\frac{\partial
a}{\partial\bm{\theta}_{1}},\cdots,\frac{\partial
a}{\partial\bm{\theta}_{n}},\frac{\partial
a}{\partial{{\bm{\phi}}}_{1}},\cdots,\frac{\partial
a}{\partial{{\bm{\phi}}}_{n}}\right]^{\top},\quad\nabla\mathbf{a}=\left[\frac{\partial\mathbf{a}}{\partial\bm{\theta}_{1}},\cdots,\frac{\partial\mathbf{a}}{\partial\bm{\theta}_{n}},\frac{\partial\mathbf{a}}{\partial{\bm{\phi}}_{1}},\cdots,\frac{\partial\mathbf{a}}{\partial{\bm{\phi}}_{n}}\right],$

respectively. Furthermore, the gradient of
$\mathbf{x},\mathbf{y},\mathbf{z},\mathbf{u},\mathbf{v},\mathbf{w}$ can be
represented by themselves,

$\displaystyle\begin{array}[]{l@{}l@{}l@{}l}&\nabla\mathbf{x}=[-\mathrm{diag}(\mathbf{y}),\mathrm{diag}(\mathbf{u})],&\nabla\mathbf{y}=[\mathrm{diag}(\mathbf{x}),\mathrm{diag}(\mathbf{v})],&\nabla\mathbf{z}=[\mathbf{0}_{n\times
n},-\mathrm{diag}(\mathbf{w})],\\\
&\nabla\mathbf{u}=[-\mathrm{diag}(\mathbf{v}),-\mathrm{diag}(\mathbf{x})],&\nabla\mathbf{v}=[\mathrm{diag}(\mathbf{u}),-\mathrm{diag}(\mathbf{y})],&\nabla\mathbf{w}=[\mathbf{0}_{n\times
n},\mathrm{diag}(\mathbf{z})].\end{array}$ (3.14)

We now vectorize the gradient (3.10) and let

$\displaystyle\mathbf{p}=D(f)\mathbf{x},\quad\mathbf{q}=D(f)\mathbf{y},\quad\mathbf{r}=D(f)\mathbf{z}.$
(3.15)

From (3.9), we have
$E_{C}(f)=\frac{1}{2}(\mathbf{x}^{\top}\mathbf{p}+\mathbf{y}^{\top}\mathbf{q}+\mathbf{z}^{\top}\mathbf{r})$.
Then, from (3.10) and (3.14), the gradient of $E_{C}(f)$ with respect to
$(\bm{\theta},{\bm{\phi}})$ is

$\displaystyle\mathbf{g}:=\nabla E_{C}$
$\displaystyle=\nabla{\text{vec}}(\mathbf{f})^{\top}\nabla_{{\text{vec}}(\mathbf{f})}E_{C}$
(3.16)
$\displaystyle=\nabla\mathbf{x}^{\top}\mathbf{p}+\nabla\mathbf{y}^{\top}\mathbf{q}+\nabla\mathbf{z}^{\top}\mathbf{r}=\begin{bmatrix}-\mathrm{diag}(\mathbf{y})\mathbf{p}+\mathrm{diag}(\mathbf{x})\mathbf{q}\\\
\mathrm{diag}(\mathbf{u})\mathbf{p}+\mathrm{diag}(\mathbf{v})\mathbf{q}-\mathrm{diag}(\mathbf{w})\mathbf{r}\end{bmatrix}\in\mathbb{R}^{2n\times
1}.$ (3.17)

Moreover, the Hessian matrix is also obtained

$\displaystyle H:=\nabla\nabla E_{C}$
$\displaystyle=(\nabla\mathbf{x}^{\top}\nabla\mathbf{p}+\nabla\mathbf{y}^{\top}\nabla\mathbf{q}+\nabla\mathbf{z}^{\top}\nabla\mathbf{r})$
(3.18)
$\displaystyle+\begin{bmatrix}-\mathrm{diag}(\mathbf{p})\nabla\mathbf{y}+\mathrm{diag}(\mathbf{q})\nabla\mathbf{x}\\\
\mathrm{diag}(\mathbf{p})\nabla\mathbf{u}+\mathrm{diag}(\mathbf{q})\nabla\mathbf{v}-\mathrm{diag}(\mathbf{r})\nabla\mathbf{w}\end{bmatrix}\in\mathbb{R}^{2n\times
2n},$

where

$\displaystyle\nabla\mathbf{p}=D(f)\nabla\mathbf{x}-\sum_{\ell=1}^{n}\mathbf{x}_{\ell}\nabla(L(f)\mathbf{e}_{\ell}),\
\nabla\mathbf{q}=D(f)\nabla\mathbf{y}-\sum_{\ell=1}^{n}\mathbf{y}_{\ell}\nabla(L(f)\mathbf{e}_{\ell}),\
\nabla\mathbf{r}=D(f)\nabla\mathbf{z}-\sum_{\ell=1}^{n}\mathbf{z}_{\ell}\nabla(L(f)\mathbf{e}_{\ell}).$

Via the chain rule

$\displaystyle\nabla(L(f)\mathbf{e}_{\ell})=$
$\displaystyle\nabla_{\mathbf{x}}(L(f)\mathbf{e}_{\ell})\nabla\mathbf{x}+\nabla_{\mathbf{y}}(L(f)\mathbf{e}_{\ell})\nabla\mathbf{y}+\nabla_{\mathbf{z}}(L(f)\mathbf{e}_{\ell})\nabla\mathbf{z}\in\mathbb{R}^{n\times
2n},$

we have

$\begin{bmatrix}\sum_{\ell=1}^{n}\mathbf{x}_{\ell}\nabla(L(f)\mathbf{e}_{\ell})\\\
\sum_{\ell=1}^{n}\mathbf{y}_{\ell}\nabla(L(f)\mathbf{e}_{\ell})\\\
\sum_{\ell=1}^{n}\mathbf{z}_{\ell}\nabla(L(f)\mathbf{e}_{\ell})\end{bmatrix}=\begin{bmatrix}\sum_{\ell=1}^{n}\mathbf{x}_{\ell}\nabla_{\mathbf{x}}(L(f)\mathbf{e}_{\ell})&\sum_{\ell=1}^{n}\mathbf{x}_{\ell}\nabla_{\mathbf{y}}(L(f)\mathbf{e}_{\ell})&\sum_{\ell=1}^{n}\mathbf{x}_{\ell}\nabla_{\mathbf{z}}(L(f)\mathbf{e}_{\ell})\\\
\sum_{\ell=1}^{n}\mathbf{y}_{\ell}\nabla_{\mathbf{x}}(L(f)\mathbf{e}_{\ell})&\sum_{\ell=1}^{n}\mathbf{y}_{\ell}\nabla_{\mathbf{y}}(L(f)\mathbf{e}_{\ell})&\sum_{\ell=1}^{n}\mathbf{y}_{\ell}\nabla_{\mathbf{z}}(L(f)\mathbf{e}_{\ell})\\\
\sum_{\ell=1}^{n}\mathbf{z}_{\ell}\nabla_{\mathbf{x}}(L(f)\mathbf{e}_{\ell})&\sum_{\ell=1}^{n}\mathbf{z}_{\ell}\nabla_{\mathbf{y}}(L(f)\mathbf{e}_{\ell})&\sum_{\ell=1}^{n}\mathbf{z}_{\ell}\nabla_{\mathbf{z}}(L(f)\mathbf{e}_{\ell})\end{bmatrix}\begin{bmatrix}\nabla\mathbf{x}\\\
\nabla\mathbf{y}\\\ \nabla\mathbf{z}\end{bmatrix}.$

(3.19)

Plugging (3.19) into (3.18), we can rewrite the Hessian matrix as

$\displaystyle H$ $\displaystyle=\begin{bmatrix}\nabla\mathbf{x}\\\
\nabla\mathbf{y}\\\
\nabla\mathbf{z}\end{bmatrix}^{\top}(L_{1}-L_{2})\begin{bmatrix}\nabla\mathbf{x}\\\
\nabla\mathbf{y}\\\
\nabla\mathbf{z}\end{bmatrix}-\begin{bmatrix}\mathrm{diag}(\mathbf{p}\odot\mathbf{x}+\mathbf{q}\odot\mathbf{y})&\mathrm{diag}(\mathbf{p}\odot\mathbf{v}-\mathbf{q}\odot\mathbf{u})\\\
\mathrm{diag}(\mathbf{p}\odot\mathbf{v}-\mathbf{q}\odot\mathbf{u})&\mathrm{diag}(\mathbf{p}\odot\mathbf{x}+\mathbf{q}\odot\mathbf{y}+\mathbf{r}\odot\mathbf{z})\end{bmatrix}$
(3.20) $\displaystyle\equiv H_{1}-H_{2}-K,$ (3.21)

where

$\displaystyle H_{s}$
$\displaystyle=\left[\nabla\mathbf{x}^{\top},\nabla\mathbf{y}^{\top},\nabla\mathbf{z}^{\top}\right]L_{s}\begin{bmatrix}\nabla\mathbf{x}\\\
\nabla\mathbf{y}\\\ \nabla\mathbf{z}\end{bmatrix},s=1,2,$ $\displaystyle K$
$\displaystyle=\begin{bmatrix}\mathrm{diag}(\mathbf{p}\odot\mathbf{x}+\mathbf{q}\odot\mathbf{y})&\mathrm{diag}(\mathbf{p}\odot\mathbf{v}-\mathbf{q}\odot\mathbf{u})\\\
\mathrm{diag}(\mathbf{p}\odot\mathbf{v}-\mathbf{q}\odot\mathbf{u})&\mathrm{diag}(\mathbf{p}\odot\mathbf{x}+\mathbf{q}\odot\mathbf{y}+\mathbf{r}\odot\mathbf{z})\end{bmatrix}$

with

(3.22)

In the matrix $L_{2}$, the Jacobian matrix of each column of $L(f)$ should be
considered. For the $\ell$-th column entries $L(f)\mathbf{e}_{\ell}$, it is
clear that $\left[L(f)\right]_{i\ell}\neq 0$ if and only if
$i\in\mathcal{N}(\ell)$ or $i=\ell$ by the definition of the Laplacian matrix
in (3.4). The nondiagonal entries are negative cotangent weights, and the
diagonal entries are the sums of the cotangent weights. Hence, the Jacobian
matrix of $L(f)\mathbf{e}_{\ell}$ is formed by the gradient of cotangent
weights $w_{ij}(f)$, especially the gradient of cotangent functions
$c_{ij}(f):=\cot\alpha_{ij}(f)$. Before discussing entries of $\nabla
L\mathbf{e}_{\ell}$, we first give a lemma for the sparsity structure of the
Jacobian block
$\sum_{\ell=1}^{n}\mathbf{a}_{\ell}\nabla_{\mathbf{b}}(L(f)\mathbf{e}_{\ell})$,
$\mathbf{a},\mathbf{b}=\mathbf{x},\mathbf{y},\mathbf{z}$.

######  Lemma 2.

Each subblock
$\sum_{\ell=1}^{n}\mathbf{a}_{\ell}\nabla_{\mathbf{b}}(L(f)\mathbf{e}_{\ell})$
of $L_{2}$ in (3.22) with
$\mathbf{a},\mathbf{b}=\mathbf{x},\mathbf{y},\mathbf{z}$ is of identical
sparsity structure to $L$.

###### Proof.

By the cotangent formula in (3.8), one can observe that the gradients with
respect to $\mathbf{x},\mathbf{y},\mathbf{z}$ are the $1$st, $2$nd and $3$rd
entries of that with respect to $\mathbf{f}$, respectively. Hence,
$\sum_{\ell=1}^{n}\mathbf{a}_{\ell}\nabla_{\mathbf{b}}(L(f)\mathbf{e}_{\ell})$
has identical sparsity for all
$\mathbf{a},\mathbf{b}=\mathbf{x},\mathbf{y},\mathbf{z}$. Here, we denote the
nonzero indices set of
$\mathbf{e}_{i}^{\top}\nabla_{\mathbf{b}}(L(f)\mathbf{e}_{j})$ as

$S(i,j)=\\{k~{}|~{}\big{[}\nabla_{\mathbf{b}}(L(f)\mathbf{e}_{j})\big{]}_{ik}\neq
0\\}.$

Clearly, the nonzero indices set of
$\mathbf{e}_{i}^{\top}\sum_{\ell=1}^{n}\mathbf{a}_{\ell}\nabla_{\mathbf{b}}(L(f)\mathbf{e}_{\ell})$
is $\bigcup_{\ell=1}^{n}S(i,\ell)$.

  * •

For $i\neq j$ and $[v_{i},v_{j}]\notin\mathcal{E}(M)$, the entries are $0$,
and therefore, the gradients are also $0$. Hence,

$\displaystyle S(i,j)=\varnothing.$ (3.23)

  * •

For $i\neq j$ and $[v_{i},v_{j}]\in\mathcal{E}(M)$, by the definition of
cotangent weight $w_{ij}(f)=\frac{1}{2}\big{(}c_{ij}(f)+c_{ji}(f)\big{)}$, we
have
$\mathbf{e}_{i}^{\top}\nabla_{\mathbf{b}}(L(f)\mathbf{e}_{j})=-\nabla_{\mathbf{b}}w_{ij}(f)=-\frac{1}{2}\big{(}\nabla_{\mathbf{b}}c_{ij}(f)+\nabla_{\mathbf{b}}c_{ji}(f)\big{)}$.
As shown in Figure 1a, one can see that $w_{ij}(f)$ relates to only $4$
vertices
$\mathbf{f}_{i},\mathbf{f}_{j},\mathbf{f}_{k},\mathbf{f}_{k^{\prime}}$.
Therefore, the entries of $\nabla_{\mathbf{b}}w_{ij}(f)$ are $0$ except for
the $i,j,k,k^{\prime}$-th entries; that is,

$\displaystyle S(i,j)=\\{i,j,k,k^{\prime}\\}\subset\mathcal{N}(i)\cup\\{i\\}.$
(3.24)

  * •

For $i=j$, we have
$\mathbf{e}_{i}^{\top}\nabla_{\mathbf{b}}(L(f)\mathbf{e}_{i})=\nabla_{\mathbf{b}}\big{[}L(f)\big{]}_{ii}=\sum_{j\in\mathcal{N}(i)}\nabla_{\mathbf{b}}w_{ij}(f)$.
Hence, we can easily verify that
$\mathbf{e}_{i}^{\top}\nabla_{\mathbf{b}}(L_{f}\mathbf{e}_{i})$ relates to the
whole adjacent vertices of $\mathbf{f}_{i}$ from Figure 1b. It immediately
follows that

$\displaystyle S(i,i)=\mathcal{N}(i)\cup\\{i\\}.$ (3.25)

Combining (3.23), (3.24) and (3.25), we have
$S(i,\ell)\subset\mathcal{N}(i)\cup\\{i\\}=S(i,i),\text{ for }\ell\neq i.$ It
follows that $\bigcup_{\ell=1}^{n}S(i,\ell)=S(i,i)=\mathcal{N}(i)\cup\\{i\\}.$
Since the index set of nonzero entries of $\mathbf{e}_{i}^{\top}L$ is
$\mathcal{N}(i)\cup\\{i\\}$, the lemma is obtained immediately. ∎

The Lemma 2 demonstrates that
$\big{[}\sum_{\ell=1}^{n}\mathbf{a}_{\ell}\nabla_{\mathbf{b}}(L(f)\mathbf{e}_{\ell})\big{]}_{ij}\neq
0$ if and only if $[v_{i},v_{j}]\in\mathcal{E}(M)$ or $i=j$ for every $i,j$.
As a result, $L_{2}$ in (3.22) is stacked by $9$ matrices with the same
sparsity as $L$ in $3\times 3$ form. Additionally, one can further observe
that $H$ in (3.20) is also stacked by $4$ matrices with the same sparsity as
$L$ in $2\times 2$ form. We summarize this idea as the following theorem.

######  Theorem 1.

$H$ in (3.20) is of identical sparsity to $\mathbf{1}_{2\times 2}\otimes L.$

###### Proof.

By the Hessian matrix representation of $H$ in (3.20) and the Lemma 2, the
proof is obtained. ∎

Theorem 1 demonstrates the high sparsity of the Hessian matrix. This property
confirms the feasibility of practically solving the large-scale linear system

$\displaystyle H\mathbf{s}=-\mathbf{g},$ (3.26)

which inspires us to apply a Newton-type algorithm for solving the
optimization problem (3.2). Furthermore, in the representation of $H$ in
(3.20),
$L_{1}-L_{2}=\nabla_{{\text{vec}}(\mathbf{f})}\nabla_{{\text{vec}}(\mathbf{f})}E_{C}$
is the Hessian matrix of the conformal energy with respect to
$\mathbf{x},\mathbf{y},\mathbf{z}$. The term
$\big{[}\nabla\mathbf{x}^{\top},\nabla\mathbf{y}^{\top},\nabla\mathbf{z}^{\top}\big{]}^{\top}$
in (3.20) is the Jacobian matrix of $\mathbf{x},\mathbf{y},\mathbf{z}$ with
respect to $\bm{\theta}$ and $\bm{\phi}$, which is stacked by $3\times 2$
diagonal matrices. The matrix $K$ in (3.20) is also stacked by $2\times 2$
diagonal matrices. These structures are invariant for arbitrary
parameterization, and therefore, the sparsity structure of Hessian matrix $H$
is also invariant. Hence, this property is also available on other
parameterizations of closed surfaces, including other expressions of spheres
and other target regions.

We proceed further to analyze the entries of Hessian matrix $H$. The entries
of $L_{1}$ are derived in (3.15) and (3.9), and the entries of $K$ in (3.20)
are derived by (3.11a)-(3.11b) and (3.15). We now give the specific derivation
for entries of $L_{2}$ in (3.20). First, we give the gradients of

$c_{ij}(f)=-\frac{\langle\mathbf{f}_{ki},\mathbf{f}_{jk}\rangle}{2\left|f(T_{ijk})\right|},\quad
c_{jk}(f)=-\frac{\langle\mathbf{f}_{ij},\mathbf{f}_{ki}\rangle}{2\left|f(T_{ijk})\right|},\quad
c_{ki}(f)=-\frac{\langle\mathbf{f}_{jk},\mathbf{f}_{ij}\rangle}{2\left|f(T_{ijk})\right|}$

in the triangle $[\mathbf{f}_{i},\mathbf{f}_{j},\mathbf{f}_{k}]$, where
$c_{ij}(f):=\cot\alpha_{ij}(f)$ as before. Since
$c_{ij}(f),c_{jk}(f),c_{ki}(f)$ only relates to vertices
$\mathbf{f}_{i},\mathbf{f}_{j},\mathbf{f}_{k}$, their gradients with respect
to other vertices are $0$. The gradients with respect to $\mathbf{f}_{j}$ are
calculated by

$\displaystyle\frac{\partial}{\partial\mathbf{f}_{j}}c_{ij}(f)$
$\displaystyle=\frac{1}{|f(T_{ijk})|}\left[-\frac{\partial|f(T_{ijk})|}{\partial\mathbf{f}_{j}}c_{ij}(f)-\frac{1}{2}\mathbf{f}_{ki}\right],$
(3.27a) $\displaystyle\frac{\partial}{\partial\mathbf{f}_{j}}c_{jk}(f)$
$\displaystyle=\frac{1}{|f(T_{ijk})|}\left[-\frac{\partial|f(T_{ijk})|}{\partial\mathbf{f}_{j}}c_{jk}(f)+\frac{1}{2}\mathbf{f}_{ki}\right],$
(3.27b) $\displaystyle\frac{\partial}{\partial\mathbf{f}_{j}}c_{ki}(f)$
$\displaystyle=\frac{1}{|f(T_{ijk})|}\left[-\frac{\partial|f(T_{ijk})|}{\partial\mathbf{f}_{j}}c_{ki}(f)+\frac{1}{2}(\mathbf{f}_{jk}-\mathbf{f}_{ij})\right],$
(3.27c)

with

$\frac{\partial|f(T_{ijk})|}{\partial\mathbf{f}_{i}}=\frac{1}{2}\big{(}c_{ij}(f)\mathbf{f}_{ij}-c_{ki}(f)\mathbf{f}_{ki}\big{)},$
(3.28a)
$\frac{\partial|f(T_{ijk})|}{\partial\mathbf{f}_{j}}=\frac{1}{2}\big{(}c_{jk}(f)\mathbf{f}_{jk}-c_{ij}(f)\mathbf{f}_{ij}\big{)},$
(3.28b)
$\frac{\partial|f(T_{ijk})|}{\partial\mathbf{f}_{k}}=\frac{1}{2}\big{(}c_{ki}(f)\mathbf{f}_{ki}-c_{jk}(f)\mathbf{f}_{jk}\big{)}.$
(3.28c)

The others can be obtained by rotating the subscript $i,j,k$ in turn.

Now, we derive the entries of matrix $L_{2}$ in (3.20). Let us consider the
block
$\sum_{\ell=1}^{n}\mathbf{a}_{\ell}\nabla_{\mathbf{b}}(L(f)\mathbf{e}_{\ell})$
for $\mathbf{a},\mathbf{b}=\mathbf{x},\mathbf{y},\mathbf{z}$, the $(i,j)$-th
entry of which is

$\displaystyle\left[\sum_{\ell=1}^{n}\mathbf{a}_{\ell}\nabla_{\mathbf{b}}(L(f)\mathbf{e}_{\ell})\right]_{ij}$
$\displaystyle=\sum_{\ell=1}^{n}\mathbf{a}_{\ell}\frac{\partial}{\partial\mathbf{b}_{j}}\big{[}L(f)\big{]}_{i\ell}=\mathbf{a}_{i}\frac{\partial}{\partial\mathbf{b}_{j}}\big{[}L(f)\big{]}_{ii}+\sum_{\ell\in\mathcal{N}(i)}\mathbf{a}_{\ell}\frac{\partial}{\partial\mathbf{b}_{j}}\big{[}L(f)\big{]}_{i\ell}$
$\displaystyle=\mathbf{a}_{i}\frac{\partial}{\partial\mathbf{b}_{j}}\sum_{\ell\in\mathcal{N}(i)}w_{i\ell}(f)-\sum_{\ell\in\mathcal{N}(i)}\mathbf{a}_{\ell}\frac{\partial}{\partial\mathbf{b}_{j}}w_{i\ell}(f)$
$\displaystyle=\sum_{\ell\in\mathcal{N}(i)}\mathbf{a}_{i\ell}\frac{\partial}{\partial\mathbf{b}_{j}}w_{i\ell}(f).$

(i) For $j\neq i$, there are only three vertices $v_{j},v_{k},v_{k^{\prime}}$
related to $v_{j}$ among all adjacent vertices of $v_{i}$, as shown in Figure
1. Therefore,

$\displaystyle\left[\sum_{\ell=1}^{n}\mathbf{a}_{\ell}\nabla_{\mathbf{b}}(L(f)\mathbf{e}_{\ell})\right]_{ij}=\sum_{\ell\in\\{j,k,k^{\prime}\\}}\mathbf{a}_{i\ell}\frac{\partial}{\partial\mathbf{b}_{j}}w_{i\ell}(f)$
$\displaystyle=$
$\displaystyle\mathbf{a}_{ik}\frac{\partial}{\partial\mathbf{b}_{j}}w_{ik}(f)+\mathbf{a}_{ik^{\prime}}\frac{\partial}{\partial\mathbf{b}_{j}}w_{ik^{\prime}}(f)+\mathbf{a}_{ij}\frac{\partial}{\partial\mathbf{b}_{j}}w_{ij}(f)$
$\displaystyle=$
$\displaystyle\frac{1}{2}\left[\left(\mathbf{a}_{ij}\frac{\partial}{\partial\mathbf{b}_{j}}c_{ij}(f)-\mathbf{a}_{ki}\frac{\partial}{\partial\mathbf{b}_{j}}c_{ki}(f)\right)+\left(\mathbf{a}_{ik^{\prime}}\frac{\partial}{\partial\mathbf{b}_{j}}c_{ik^{\prime}}(f)-\mathbf{a}_{ji}\frac{\partial}{\partial\mathbf{b}_{j}}c_{ji}(f)\right)\right].$
(3.29)

One can see that the terms in the first bracket relate only to triangle
$T_{ijk}$, while the terms in the second bracket relate to triangle
$T_{k^{\prime}ji}$, which are of the same form. This characteristic is also
similar to $L$, whose cotangent weight is $w_{ij}=\frac{1}{2}(c_{ij}+c_{ji})$
with $c_{ij},c_{ji}$ also relating to $T_{ijk}$ and $T_{k^{\prime}ji}$,
respectively. Hence, we only give the specific representation of the first
bracket of (3.29). The second bracket is obtained similarly. By
(3.27a)-(3.27c), we have

$\displaystyle\mathbf{a}_{ij}\frac{\partial}{\partial\mathbf{b}_{j}}c_{ij}(f)+\mathbf{a}_{ik}\frac{\partial}{\partial\mathbf{b}_{j}}c_{ki}(f)$
$\displaystyle=$
$\displaystyle-\frac{1}{2|f(T_{ijk})|}\left[\frac{2\partial|f(T_{ijk})|}{\partial\mathbf{a}_{i}}\frac{2\partial|f(T_{ijk})|}{\partial\mathbf{b}_{j}}+(2\mathbf{a}_{ki}\mathbf{b}_{jk}-\mathbf{a}_{jk}\mathbf{b}_{ki})\right].$
(3.30)

(ii) For $j=i$, the related vertices are the whole adjacent vertices of
$v_{i}$ and itself, i.e., $\mathcal{N}(i)\cup\\{i\\}$. Hence,

$\displaystyle\left[\sum_{\ell=1}^{n}\mathbf{b}_{\ell}\nabla_{\mathbf{a}}(L(f)\mathbf{e}_{\ell})\right]_{ii}=\sum_{\ell\in\mathcal{N}(i)}\mathbf{a}_{i\ell}\frac{\partial}{\partial\mathbf{b}_{i}}w_{i\ell}(f)$
$\displaystyle=$ $\displaystyle\frac{1}{2}\sum_{\\{j,k\\}\in
S_{\mathcal{E}}(i)}\left(\mathbf{a}_{ij}\frac{\partial}{\partial\mathbf{b}_{i}}c_{ij}(f)-\mathbf{a}_{ki}\frac{\partial}{\partial\mathbf{b}_{i}}c_{ki}(f)\right),$

where
$S_{\mathcal{E}}(i)=\big{\\{}\\{j,k\\}~{}|~{}j,k\in\mathcal{N}(i),[v_{j},v_{k}]\in\mathcal{E}(M)\big{\\}}$.
Similarly, via the gradient formulas (3.27a)-(3.27c), we have

$\displaystyle\mathbf{a}_{ij}\frac{\partial}{\partial\mathbf{b}_{i}}c_{ij}(f)-\mathbf{a}_{ki}\frac{\partial}{\partial\mathbf{b}_{i}}c_{ki}(f)$
$\displaystyle=$
$\displaystyle-\frac{1}{2|f(T_{ijk})|}\left[\frac{2\partial|f(T_{ijk})|}{\partial\mathbf{a}_{i}}\frac{2\partial|f(T_{ijk})|}{\partial\mathbf{b}_{i}}+\mathbf{a}_{jk}\mathbf{b}_{jk}\right].$
(3.31)

We now have derived the entries of matrix $L_{1}-L_{2}$, which is the Hessian
matrix of the conformal energy with respect to Cartesian coordinates
$[\mathbf{x},\mathbf{y},\mathbf{z}]$. The following theorem shows that the
nullity of $L_{1}-L_{2}$ is $3$.

######  Theorem 2.

$L_{1}-L_{2}$ is a symmetric Laplacian matrix with its null space having
orthogonal basis
$[\mathbf{1}^{\top},\mathbf{0}^{\top},\mathbf{0}^{\top}]^{\top}$,
$[\mathbf{0}^{\top},\mathbf{1}^{\top},\mathbf{0}^{\top}]^{\top}$ and
$[\mathbf{0}^{\top},\mathbf{0}^{\top},\mathbf{1}^{\top}]^{\top}$.

###### Proof.

Since $L_{1}-L_{2}$ is the Hessian matrix of the conformal energy with respect
to Cartesian coordinates $[\mathbf{x},\mathbf{y},\mathbf{z}]$, it is obviously
symmetric. By the gradient formulas (3.27a)-(3.27c), it is easily seen that

$\displaystyle\left(\frac{\partial}{\partial\mathbf{f}_{i}}+\frac{\partial}{\partial\mathbf{f}_{j}}+\frac{\partial}{\partial\mathbf{f}_{k}}\right)c_{ij}(f)=0.$

Therefore, by using the representation according to the adjacent vertices, we
have

$\displaystyle\sum_{j=1}^{n}\left[\sum_{\ell=1}^{n}\mathbf{a}_{\ell}\nabla_{\mathbf{b}}(L(f)\mathbf{e}_{\ell})\right]_{ij}$
$\displaystyle=\frac{1}{2}\sum_{j\in\mathcal{N}(i)}\mathbf{a}_{ij}\left(\frac{\partial}{\partial\mathbf{b}_{i}}+\frac{\partial}{\partial\mathbf{b}_{j}}+\frac{\partial}{\partial\mathbf{b}_{k}}\right)c_{ij}(f)$
$\displaystyle-\frac{1}{2}\sum_{j\in\mathcal{N}(i)}\mathbf{a}_{ji}\left(\frac{\partial}{\partial\mathbf{b}_{i}}+\frac{\partial}{\partial\mathbf{b}_{j}}+\frac{\partial}{\partial\mathbf{b}_{k^{\prime}}}\right)c_{ji}(f)=0,~{}i=1,2,\cdots,n.$

Additionally, by the equality
$\mathbf{1}^{\top}(\nabla_{\mathbf{b}}L(f)\mathbf{e}_{\ell})=\nabla_{\mathbf{b}}(\mathbf{1}^{\top}L(f)\mathbf{e}_{\ell})=0$,
we also have

$\sum_{i=1}^{n}\left[\sum_{\ell=1}^{n}\mathbf{a}_{\ell}\nabla_{\mathbf{b}}(L(f)\mathbf{e}_{\ell})\right]_{ij}=\left[\sum_{\ell=1}^{n}\mathbf{a}_{\ell}\nabla_{\mathbf{b}}(\mathbf{1}^{\top}L(f)\mathbf{e}_{\ell})\right]\mathbf{e}_{j}=0,~{}j=1,2,\cdots,n.$

We conclude that the sum of each row and column of
$\sum_{\ell=1}^{n}\mathbf{a}_{\ell}\nabla_{\mathbf{b}}(L(f)\mathbf{e}_{\ell})$
are zero, which guarantees the $3$ orthogonal basis of the null space. ∎

Finally, we focus on

$\displaystyle
H_{2}:=\left[\nabla\mathbf{x}^{\top},\nabla\mathbf{y}^{\top},\nabla\mathbf{z}^{\top}\right]L_{2}\begin{bmatrix}\nabla\mathbf{x}\\\
\nabla\mathbf{y}\\\
\nabla\mathbf{z}\end{bmatrix}=\begin{bmatrix}H_{2,\bm{\theta}\bm{\theta}}&H_{2,\bm{\theta}\bm{\phi}}\\\
H_{2,\bm{\phi}\bm{\theta}}&H_{2,\bm{\phi}\bm{\phi}}\end{bmatrix}$

defined in (3.20) and (3.21). The formulas (3.30) and (3.31) show that the
entries of $L_{2}$ can be rewritten as the inner products of two terms related
to $\mathbf{a}$ and $\mathbf{b}$, respectively. Taking (3.31) as an example,
we have

$\displaystyle\mathbf{a}_{ij}\frac{\partial}{\partial\mathbf{b}_{i}}c_{ij}(f)+\mathbf{a}_{ik}\frac{\partial}{\partial\mathbf{b}_{i}}c_{ki}(f)=-\frac{1}{2|f(T_{ijk})|}\left\langle\left[\frac{2\partial|f(T_{ijk})|}{\partial\mathbf{a}_{i}},\mathbf{a}_{jk}\right],\left[\frac{2\partial|f(T_{ijk})|}{\partial\mathbf{b}_{i}},\mathbf{b}_{jk}\right]\right\rangle,$

From gradient formulas (3.28a)-(3.28c), we can find that the first and second
terms are only associated with $\mathbf{a}$ and $\mathbf{b}$, respectively.
Additionally, (3.30) can also be written as a similar inner product form.
Then, since $\mathbf{f}_{i}$ depends only on
$(\bm{\theta}_{i},{\bm{\phi}}_{i})$ for $i=1,2,\cdots,n$, the Jacobian
matrices of $\mathbf{x},\mathbf{y},\mathbf{z}$ with respect to
$\bm{\theta},{\bm{\phi}}$ are diagonal. It is easy to express the products of
Jacobian matrices and $L_{2}$. Without loss of generality, we discuss only
entry $\big{[}H_{2,\bm{\theta}\bm{\phi}}\big{]}_{ii}$,

$\displaystyle\big{[}H_{2,\bm{\theta}\bm{\phi}}\big{]}_{ii}=$
$\displaystyle\sum_{\\{j,k\\}\in
S_{\mathcal{E}}(i)}\sum_{\mathbf{a},\mathbf{b}=\mathbf{x},\mathbf{y},\mathbf{z}}\left(\mathbf{a}_{ij}\frac{\partial}{\partial\mathbf{b}_{i}}c_{ij}(f)-\mathbf{a}_{ki}\frac{\partial}{\partial\mathbf{b}_{i}}c_{ki}(f)\right)\frac{\partial\mathbf{a}_{i}}{\partial\bm{\theta}_{i}}\frac{\partial\mathbf{b}_{i}}{\partial{\bm{\phi}}_{i}}$
$\displaystyle=$ $\displaystyle-\sum_{\\{j,k\\}\in
S_{\mathcal{E}}(i)}\frac{1}{2|f(T_{ijk})|}\sum_{\mathbf{a},\mathbf{b}=\mathbf{x},\mathbf{y},\mathbf{z}}\left\langle\frac{\partial\mathbf{a}_{i}}{\partial\bm{\theta}_{i}}\left[\frac{2\partial|f(T_{ijk})|}{\partial\mathbf{a}_{i}},\mathbf{a}_{jk}\right],\frac{\partial\mathbf{b}_{i}}{\partial{\bm{\phi}}_{i}}\left[\frac{2\partial|f(T_{ijk})|}{\partial\mathbf{b}_{i}},\mathbf{b}_{jk}\right]\right\rangle$
$\displaystyle=$ $\displaystyle-\sum_{\\{j,k\\}\in
S_{\mathcal{E}}(i)}\frac{1}{2|f(T_{ijk})|}\left\langle\sum_{\mathbf{a}=\mathbf{x},\mathbf{y},\mathbf{z}}\frac{\partial\mathbf{a}_{i}}{\partial\bm{\theta}_{i}}\left[\frac{2\partial|f(T_{ijk})|}{\partial\mathbf{a}_{i}},\mathbf{a}_{jk}\right],\sum_{\mathbf{b}=\mathbf{x},\mathbf{y},\mathbf{z}}\frac{\partial\mathbf{b}_{i}}{\partial{\bm{\phi}}_{i}}\left[\frac{2\partial|f(T_{ijk})|}{\partial\mathbf{b}_{i}},\mathbf{b}_{jk}\right]\right\rangle$
$\displaystyle=$ $\displaystyle-\sum_{\\{j,k\\}\in
S_{\mathcal{E}}(i)}\frac{1}{2|f(T_{ijk})|}\left\langle\frac{\partial\mathbf{f}_{i}}{\partial\bm{\theta}_{i}}^{\top}\left[\frac{2\partial|f(T_{ijk})|}{\partial\mathbf{f}_{i}},\mathbf{f}_{jk}\right],\frac{\partial\mathbf{f}_{i}}{\partial{\bm{\phi}}_{i}}^{\top}\left[\frac{2\partial|f(T_{ijk})|}{\partial\mathbf{f}_{i}},\mathbf{f}_{jk}\right]\right\rangle.$

The partial differential terms in (3.28a)-(3.28c) show that they are the
linear combination of $\mathbf{f}$. Hence, only the sum with respect to
$\mathbf{x},\mathbf{y},\mathbf{z}$ is necessary, such as

$\begin{bmatrix}\mathbf{a}^{i}_{jk}\\\
\mathbf{b}^{i}_{jk}\end{bmatrix}:=\left(\frac{\partial\mathbf{f}_{i}}{\partial(\bm{\theta}_{i},{\bm{\phi}}_{i})}\right)^{\top}\mathbf{f}_{jk}=\left[\begin{aligned}
-&\mathbf{y}_{i}\mathbf{x}_{jk}+\mathbf{x}_{i}\mathbf{y}_{jk}\\\
&\mathbf{u}_{i}\mathbf{x}_{jk}+\mathbf{v}_{i}\mathbf{y}_{jk}-\mathbf{w}_{i}\mathbf{z}_{jk}\end{aligned}\right],$
(3.32a) $\begin{bmatrix}\mathbf{a}^{i}_{ki}\\\
\mathbf{b}^{i}_{ki}\end{bmatrix}:=\left(\frac{\partial\mathbf{f}_{i}}{\partial(\bm{\theta}_{i},{\bm{\phi}}_{i})}\right)^{\top}\mathbf{f}_{ki}=\left[\begin{aligned}
-&\mathbf{y}_{i}\mathbf{x}_{ki}+\mathbf{x}_{i}\mathbf{y}_{ki}\\\
&\mathbf{u}_{i}\mathbf{x}_{ki}+\mathbf{v}_{i}\mathbf{y}_{ki}-\mathbf{w}_{i}\mathbf{z}_{ki}\end{aligned}\right],$
(3.32b) $\begin{bmatrix}\mathbf{a}^{i}_{ij}\\\
\mathbf{b}^{i}_{ij}\end{bmatrix}:=\left(\frac{\partial\mathbf{f}_{i}}{\partial(\bm{\theta}_{i},{\bm{\phi}}_{i})}\right)^{\top}\mathbf{f}_{ij}=\left[\begin{aligned}
-&\mathbf{y}_{i}\mathbf{x}_{ij}+\mathbf{x}_{i}\mathbf{y}_{ij}\\\
&\mathbf{u}_{i}\mathbf{x}_{ij}+\mathbf{v}_{i}\mathbf{y}_{ij}-\mathbf{w}_{i}\mathbf{z}_{ij}\end{aligned}\right].$
(3.32c)

The others are defined similarly by modifying the superscripts and subscripts.
Here, $\mathbf{a}^{i}_{\cdot\cdot}$ denotes the $\partial\bm{\theta}_{i}$
term, and $\mathbf{b}^{i}_{\cdot\cdot}$ denotes the $\partial{\bm{\phi}}_{i}$
term. Thus, we have

$\displaystyle\frac{\partial\mathbf{f}_{i}}{\partial\bm{\theta}_{i}}^{\top}\left[\frac{2\partial|f(T_{ijk})|}{\partial\mathbf{f}_{i}},\mathbf{f}_{jk}\right]=\frac{\partial\mathbf{f}_{i}}{\partial\bm{\theta}_{i}}^{\top}$
$\displaystyle\big{[}c_{ij}(f)\mathbf{f}_{ij}-c_{ki}(f)\mathbf{f}_{ki},\mathbf{f}_{jk}\big{]}$
$\displaystyle=$
$\displaystyle\big{[}c_{ij}(f)\mathbf{a}_{ij}^{i}-c_{ki}(f)\mathbf{a}_{ki}^{i},\mathbf{a}_{jk}^{i}\big{]}.$

We can see that they differ only in $\mathbf{f}_{\cdot\cdot}$ and
$\mathbf{a}_{\cdot\cdot}^{i},\mathbf{b}_{\cdot\cdot}^{i}$. Hence, the entries
of $H_{2}$ differ from only those of $L_{2}$ with the notations
$\mathbf{x}_{\cdot\cdot},\mathbf{y}_{\cdot\cdot},\mathbf{z}_{\cdot\cdot}$ in
(3.30) and (3.31) replaced by $\mathbf{a}_{\cdot\cdot}^{i}$ and
$\mathbf{b}_{\cdot\cdot}^{i}$. Taking $H_{2,\bm{\theta}\bm{\phi}}$ as an
example, the diagonal and nondiagonal entries are

$\displaystyle\big{[}H_{2,\bm{\theta}\bm{\phi}}\big{]}_{ii}=$
$\displaystyle-\sum_{\\{j,k\\}\in
S_{\mathcal{E}}(i)}\frac{1}{4|f(T_{ijk})|}\left[\frac{2\partial|f(T_{ijk})|}{\partial\bm{\theta}_{i}}\frac{2\partial|f(T_{ijk})|}{\partial{\bm{\phi}}_{i}}+\mathbf{a}^{i}_{jk}\mathbf{b}^{i}_{jk}\right],$
$\displaystyle\big{[}H_{2,\bm{\theta}\bm{\phi}}\big{]}_{ij}=$
$\displaystyle-\frac{1}{4|f(T_{ijk})|}\left[\frac{2\partial|f(T_{ijk})|}{\partial\bm{\theta}_{i}}\frac{2\partial|f(T_{ijk})|}{\partial{\bm{\phi}}_{j}}+(2\mathbf{a}_{ki}^{i}\mathbf{b}_{jk}^{j}-\mathbf{a}_{jk}^{i}\mathbf{b}_{ki}^{j})\right]$
$\displaystyle-\frac{1}{4|f(T_{k^{\prime}ji})|}\left[\frac{2\partial|f(T_{k^{\prime}ji})|}{\partial\bm{\theta}_{i}}\frac{2\partial|f(T_{k^{\prime}ji})|}{\partial{\bm{\phi}}_{j}}+(2\mathbf{a}_{ik^{\prime}}^{i}\mathbf{b}_{k^{\prime}j}^{j}-\mathbf{a}_{k^{\prime}j}^{i}\mathbf{b}_{ik^{\prime}}^{j})\right],$

respectively, where

$\frac{\partial|f(T_{ijk})|}{\partial\bm{\theta}_{i}}=\frac{1}{2}\big{(}c_{ij}(f)\mathbf{a}^{i}_{ij}-c_{ki}(f)\mathbf{a}^{i}_{ki}\big{)},$
(3.33a)
$\frac{\partial|f(T_{ijk})|}{\partial{\bm{\phi}}_{i}}=\frac{1}{2}\big{(}c_{ij}(f)\mathbf{b}^{i}_{ij}-c_{ki}(f)\mathbf{b}^{i}_{ki}\big{)}.$
(3.33b)

The other entries are obtained similarly by replacing
$\partial\mathbf{x}_{i}$, $\partial\mathbf{y}_{i}$, $\partial\mathbf{z}_{i}$
and $\mathbf{x}_{\cdot\cdot}$, $\mathbf{y}_{\cdot\cdot}$,
$\mathbf{z}_{\cdot\cdot}$ by $\partial\bm{\theta}_{i},\partial{\bm{\phi}}_{i}$
and $\mathbf{a}_{\cdot\cdot}^{i},\mathbf{b}_{\cdot\cdot}^{i}$ as in
(3.32a)-(3.32c), respectively. Now, we have obtained the entries of Hessian
matrix $H$ in (3.20).

######  Theorem 3.

Let $H=\begin{bmatrix}H_{\bm{\theta}\bm{\theta}}&H_{\bm{\theta}{\bm{\phi}}}\\\
H_{{\bm{\phi}}\bm{\theta}}&H_{{\bm{\phi}}{\bm{\phi}}}\end{bmatrix}$ with
$H_{\bm{\theta}{\bm{\phi}}}=H_{{\bm{\phi}}\bm{\theta}}^{\top}$. The diagonal
and nondiagonal entries of the blocks are

$\left\\{\begin{array}[]{r@{}l}\big{[}H_{\bm{\theta}\bm{\theta}}\big{]}_{ii}&=\sum_{\\{j,k\\}\in
S_{\mathcal{E}}(i)}\frac{1}{4|f(T_{ijk})|}\left[\left(\frac{2\partial|f(T_{ijk})|}{\partial\bm{\theta}_{i}}\right)^{2}+\big{(}\mathbf{a}^{i}_{jk}\big{)}^{2}\right]\\\\[4.0pt]
&+\sum_{j\in\mathcal{N}(i)}\tilde{w}_{ij}(\mathbf{x}_{i}\mathbf{x}_{j}+\mathbf{y}_{i}\mathbf{y}_{j}),\\\\[2.0pt]
\big{[}H_{\bm{\theta}\bm{\theta}}\big{]}_{ij}&=\frac{1}{4|f(T_{ijk})|}\left[\frac{2\partial|f(T_{ijk})|}{\partial\bm{\theta}_{i}}\frac{2\partial|f(T_{ijk})|}{\partial{\bm{\theta}}_{j}}+(2\mathbf{a}_{ki}^{i}\mathbf{a}_{jk}^{j}-\mathbf{a}_{jk}^{i}\mathbf{a}_{ki}^{j})\right]\\\\[4.0pt]
&+\frac{1}{4|f(T_{k^{\prime}ji})|}\left[\frac{2\partial|f(T_{k^{\prime}ji})|}{\partial\bm{\theta}_{i}}\frac{2\partial|f(T_{k^{\prime}ji})|}{\partial{\bm{\theta}}_{j}}+(2\mathbf{a}_{ik^{\prime}}^{i}\mathbf{a}_{k^{\prime}j}^{j}-\mathbf{a}_{k^{\prime}j}^{i}\mathbf{a}_{ik^{\prime}}^{j})\right]\\\\[4.0pt]
&-\tilde{w}_{ij}(\mathbf{x}_{i}\mathbf{x}_{j}+\mathbf{y}_{i}\mathbf{y}_{j});\\\
\end{array}\right.$
$\left\\{\begin{array}[]{r@{}l}\big{[}H_{\bm{\phi}\bm{\phi}}\big{]}_{ii}&=\sum_{\\{j,k\\}\in
S_{\mathcal{E}}(i)}\frac{1}{4|f(T_{ijk})|}\left[\left(\frac{2\partial|f(T_{ijk})|}{\partial\bm{\phi}_{i}}\right)^{2}+\big{(}\mathbf{b}^{i}_{jk}\big{)}^{2}\right]\\\\[4.0pt]
&+\sum_{j\in\mathcal{N}(i)}\tilde{w}_{ij}(\mathbf{x}_{i}\mathbf{x}_{j}+\mathbf{y}_{i}\mathbf{y}_{j}+\mathbf{z}_{i}\mathbf{z}_{j}),\\\\[4.0pt]
\big{[}H_{\bm{\phi}\bm{\phi}}\big{]}_{ij}&=\frac{1}{4|f(T_{ijk})|}\left[\frac{2\partial|f(T_{ijk})|}{\partial\bm{\phi}_{i}}\frac{2\partial|f(T_{ijk})|}{\partial{\bm{\phi}}_{j}}+(2\mathbf{b}_{ki}^{i}\mathbf{b}_{jk}^{j}-\mathbf{b}_{jk}^{i}\mathbf{b}_{ki}^{j})\right]\\\\[4.0pt]
&+\frac{1}{4|f(T_{k^{\prime}ji})|}\left[\frac{2\partial|f(T_{k^{\prime}ji})|}{\partial\bm{\phi}_{i}}\frac{2\partial|f(T_{k^{\prime}ji})|}{\partial{\bm{\phi}}_{j}}+(2\mathbf{b}_{ik^{\prime}}^{i}\mathbf{b}_{k^{\prime}j}^{j}-\mathbf{b}_{k^{\prime}j}^{i}\mathbf{b}_{ik^{\prime}}^{j})\right]\\\\[4.0pt]
&-\tilde{w}_{ij}(\mathbf{u}_{i}\mathbf{u}_{j}+\mathbf{v}_{i}\mathbf{v}_{j}+\mathbf{w}_{i}\mathbf{w}_{j});\\\
\end{array}\right.$
$\left\\{\begin{array}[]{r@{}l}\big{[}H_{\bm{\phi}\bm{\theta}}\big{]}_{ii}&=\sum_{\\{j,k\\}\in
S_{\mathcal{E}}(i)}\frac{1}{4|f(T_{ijk})|}\left[\frac{2\partial|f(T_{ijk})|}{\partial\bm{\phi}_{i}}\frac{2\partial|f(T_{ijk})|}{\partial{\bm{\theta}}_{i}}+\mathbf{b}^{i}_{jk}\mathbf{a}^{i}_{jk}\right]\\\\[4.0pt]
&+\sum_{j\in\mathcal{N}(i)}\tilde{w}_{ij}(\mathbf{v}_{i}\mathbf{x}_{j}-\mathbf{u}_{i}\mathbf{y}_{j}),\\\\[2.0pt]
\big{[}H_{\bm{\phi}\bm{\theta}}\big{]}_{ij}&=\frac{1}{4|f(T_{ijk})|}\left[\frac{2\partial|f(T_{ijk})|}{\partial\bm{\phi}_{i}}\frac{2\partial|f(T_{ijk})|}{\partial{\bm{\theta}}_{j}}+(2\mathbf{b}_{ki}^{i}\mathbf{a}_{jk}^{j}-\mathbf{b}_{jk}^{i}\mathbf{a}_{ki}^{j})\right]\\\\[4.0pt]
&+\frac{1}{4|f(T_{k^{\prime}ji})|}\left[\frac{2\partial|f(T_{k^{\prime}ji})|}{\partial\bm{\phi}_{i}}\frac{2\partial|f(T_{k^{\prime}ji})|}{\partial{\bm{\theta}}_{j}}+(2\mathbf{b}_{ik^{\prime}}^{i}\mathbf{a}_{k^{\prime}j}^{j}-\mathbf{b}_{k^{\prime}j}^{i}\mathbf{a}_{ik^{\prime}}^{j})\right]\\\\[4.0pt]
&-\tilde{w}_{ij}(\mathbf{v}_{i}\mathbf{x}_{j}-\mathbf{u}_{i}\mathbf{y}_{j});\end{array}\right.$
(3.34)

respectively.

Finally, we present a theorem for the proposed algorithm in Section 4, which
demonstrates singularity and the corresponding eigenpair of $H$.

######  Theorem 4.

The null space of $H$ defined in (3.20) has a basis
$[\mathbf{1}^{\top},\mathbf{0}^{\top}]^{\top}$, where
$\text{dim}(\mathbf{1})=\text{dim}(\mathbf{0})=n$.

###### Proof.

To prove the assertion, we need to indicate only that

$\displaystyle\mathbf{e}_{i}^{\top}H_{\bm{\theta}\bm{\theta}}\mathbf{1}=0,\quad\text{
and
}\quad\mathbf{e}_{i}^{\top}H_{{\bm{\phi}}\bm{\theta}}\mathbf{1}=0,\quad\text{for
}i=1,2,\cdots,n.$

We first consider $H_{{\bm{\phi}}\bm{\theta}}$. By (3.34) we have

$\displaystyle\mathbf{e}_{i}^{\top}H_{{\bm{\phi}}\bm{\theta}}\mathbf{1}=\sum_{\\{j,k\\}\in
S_{\mathcal{E}}(i)}\frac{1}{4|f(T_{ijk})|}$
$\displaystyle\Bigg{[}\frac{2\partial|f(T_{ijk})|}{\partial{\bm{\phi}}_{i}}\cdot
2\left(\frac{\partial}{\partial\bm{\theta}_{i}}+\frac{\partial}{\partial\bm{\theta}_{j}}+\frac{\partial}{\partial\bm{\theta}_{k}}\right)|f(T_{ijk})|$
$\displaystyle+\mathbf{b}^{i}_{jk}\mathbf{a}^{i}_{jk}+2\mathbf{b}_{ki}^{i}\mathbf{a}_{jk}^{j}-\mathbf{b}_{jk}^{i}\mathbf{a}_{ki}^{j}+2\mathbf{b}_{ij}^{i}\mathbf{a}_{jk}^{k}-\mathbf{b}_{jk}^{i}\mathbf{a}_{ij}^{k}\Bigg{]}$
(3.35)

It is easy to verify that

$\displaystyle\mathbf{a}_{ij}^{i}=\mathbf{a}_{ij}^{j},\quad\mathbf{a}_{jk}^{j}=\mathbf{a}_{jk}^{k},\quad\mathbf{a}_{ki}^{k}=\mathbf{a}_{ki}^{i}.$
(3.36)

Combining (3.36) with (3.33a), we have

$\displaystyle\left(\frac{\partial}{\partial\bm{\theta}_{i}}+\frac{\partial}{\partial\bm{\theta}_{j}}+\frac{\partial}{\partial\bm{\theta}_{k}}\right)|f(T_{ijk})|=0.$

Hence, by using (3.36) again, (3.35) becomes

$\displaystyle\mathbf{e}_{i}^{\top}H_{{\bm{\phi}}\bm{\theta}}\mathbf{1}=$
$\displaystyle\sum_{\\{j,k\\}\in
S_{\mathcal{E}}(i)}\frac{1}{4|f(T_{ijk})|}\left[\mathbf{b}^{i}_{jk}\mathbf{a}^{i}_{jk}+2\mathbf{b}_{ki}^{i}\mathbf{a}_{jk}^{j}-\mathbf{b}_{jk}^{i}\mathbf{a}_{ki}^{j}+2\mathbf{b}_{ij}^{i}\mathbf{a}_{jk}^{k}-\mathbf{b}_{jk}^{i}\mathbf{a}_{ij}^{k}\right]$
$\displaystyle=$ $\displaystyle\sum_{\\{j,k\\}\in
S_{\mathcal{E}}(i)}\frac{1}{4|f(T_{ijk})|}\mathbf{b}^{i}_{jk}\left(\mathbf{a}_{ij}^{i}+\mathbf{a}^{i}_{jk}+\mathbf{a}_{ki}^{i}\right)=0.$

Similarly, we also have
$\mathbf{e}_{i}^{\top}H_{\bm{\theta}\bm{\theta}}\mathbf{1}=0$. Notably,
$\mathbf{e}_{i}^{\top}H_{{\bm{\phi}}{\bm{\phi}}}\mathbf{1}\neq 0$. Therefore,
the theorem is proved. ∎

######  Remark 2.

Geometrically, the one-dimensional null space of $H$ reveals that the
conformal energy is invariant up to a rotation along the latitude, while that
along the longitude is not characterized, since it is not linearly related to
${\bm{\phi}}$.

##  4 Hessian-Based Trust Region Algorithm

In this section, we develop an HBTR algorithm to minimize the conformal energy
for the computation of the conformal map from a closed surface of genus-$0$ to
a unit sphere. Here, we review the optimization problem,

$\displaystyle\min E_{C}(\bm{\theta},{\bm{\phi}}):=\frac{1}{2}\langle
D(f)\mathbf{f},\mathbf{f}\rangle,\qquad\mathbf{f}=[\cos\bm{\theta}\odot\sin{\bm{\phi}},\sin\bm{\theta}\odot\sin{\bm{\phi}},\cos{\bm{\phi}}].$
(4.1)

We loosen the box constraint
$(\bm{\theta},{\bm{\phi}})\in[0,2\pi]^{n}\times[0,\pi]^{n}$ and consider (4.1)
as an unconstrained problem. In Section 3, the gradient vector and Hessian
matrix of $E_{C}({\bm{\theta}},{\bm{\phi}})$ are derived in (3.16) and Theorem
3, respectively. The sparsity of the Hessian matrix shown in Theorem 1
guarantees the feasibility of fast computation associated with $H$.
Specifically, we solve the large-scale sparse linear system to obtain the
Newton direction: $H\mathbf{s}=-\mathbf{g}$ as in (3.26). As Theorem 4
demonstrated, $H$ is singular, and the general solution can be expressed as a
spherical solution with arbitrary rotation along latitude according to Remark
2. Therefore, we fix the first entry of $\mathbf{s}$ and let

$H=\begin{bmatrix}h_{11}&\mathbf{h}_{1}^{\top}\\\
\mathbf{h}_{1}&\hat{H}\end{bmatrix},\quad\mathbf{s}=\begin{bmatrix}0\\\
\hat{\mathbf{s}}\end{bmatrix},\quad\mathbf{g}=\begin{bmatrix}g_{1}\\\
\hat{\mathbf{g}}\end{bmatrix}.$

Then, we solve the linear system

$\displaystyle\hat{H}\hat{\mathbf{s}}=-\hat{\mathbf{g}}.$ (4.2)

It is easy to verify that $H\mathbf{s}=-\mathbf{g}$ holds for
$\mathbf{s}=[0,\hat{\mathbf{s}}^{\top}]^{\top}$. Geometrically, this approach
means that the longitude of the first vertex always remains invariant and that
meaningless rotation of the sphere is avoided during the iteration.

######  Remark 3.

$H$ can be rearranged from $L\otimes\mathbf{1}_{2\times 2}$ to
$\mathbf{1}_{2\times 2}\otimes L$. The rearranged matrix is constructed by the
$2\times 2$ block in $L$ form such that the block LU decomposition is
appropriate for the fast computation of the linear system (4.2).

Therefore, it is natural to use the Newton-type method to solve the
optimization problem (4.1). However, the pure Newton method is insufficient
for solving the optimization problem directly because the Newton method is
well known to have local quadratic convergence. Unfortunately, it is almost
impractical to directly seek an initial guess sufficiently close to the ideal
solution for (4.1), especially when facing surfaces with high-curvature
regions or complicated shapes. Additionally, the conformal energy of (4.1)
with respect to $(\bm{\theta},{\bm{\phi}})$ is nonconvex and nonlinear. Hence,
the quadratic convergence generally disappears at the beginning of iterations,
and it might take much time to reach the neighborhood with quadratic
convergence via the Newton step. Moreover, importantly, the Hessian matrix $H$
is not uniformly positive semidefinite; that is, $H$ may be indefinite at some
points $(\bm{\theta},{\bm{\phi}})$. Indeed, negative curvature occurs
frequently during iteration in practical experiments. Consequently, it is
possible that the Newton direction may not be a descent direction.

To overcome this drawback, we introduce the negative gradient direction, which
is a descent direction. More specifically, we search the trial step
$\mathbf{d}$ from the $2$D subspace spanned by the Newton direction
$\mathbf{s}$ and the gradient direction $\mathbf{g}$. If $\mathbf{s}$ is not a
descent direction, it is a negative curvature direction, which is still
beneficial information for the choice of descent direction. Moreover, there
must be a descent direction in this $2$D subspace since
$-\mathbf{g}\in\text{span}(\mathbf{s},\mathbf{g})$. When the iterative point
is in the convergence neighborhood, the Newton direction guarantees the
quadratic convergence of the algorithm. We utilize the trust region method to
search the trial step from $\text{span}(\mathbf{s},\mathbf{g})$. In other
words, during each iterative step, we consider the trust region subproblem
proposed by Shultz et al.[43]

$\displaystyle\begin{array}[]{cl}\min&\mathbf{d}^{\top}H\mathbf{d}+\mathbf{g}^{\top}\mathbf{d}\\\
\text{s.t.}&\|\mathbf{d}\|\leq\mit\Delta,\
\mathbf{d}\in\text{span}(\mathbf{s},\mathbf{g}),\end{array}$ (4.5)

where $\mit\Delta$ is the trust region radius. The $2$D optimization problem
(4.5) is easy to solve and costs less time. The global and local convergence
of the trust region method with subproblem (4.5) has been proven by [43], and
the practical experiment has been verified later in [13].

For the error measurement, it is worth noting that
$E_{C}(\bm{\theta},{\bm{\phi}})$ is a periodic function. Since we loosen the
box constraint $(\bm{\theta},{\bm{\phi}})\in[0,2\pi]^{n}\times[0,\pi]^{n}$,
$(\bm{\theta},{\bm{\phi}})$ might move larger than $2\pi$, while the vertices
on the sphere move much less. Hence, it is inappropriate to check the length
of step $\mathbf{d}$. Inspired by [47], we adopt the optimal rotation and
measure the error with

$\displaystyle\delta^{(k)}=\min_{R\in
SO(3)}\left\|\mathbf{f}^{(k+1)}-\mathbf{f}^{(k)}R\right\|^{2},$ (4.6)

where $SO(3)=\left\\{R\in\mathbb{R}^{3\times
3}~{}|~{}R^{\top}R=I,~{}\det(R)=1\right\\}$ and $\mathbf{f}^{(k)}$ are the
vertices on the sphere at the $k$-th iteration. The error measurement
eliminates the error from rotation. In addition, it avoids the miscalculation
caused by the periodicity of $E_{C}(\bm{\theta},{\bm{\phi}})$.

We summarize the proposed Hessian-based trust region (HBTR) algorithm for the
optimization problem (4.1).

Algorithm 1 HBTR for Spherical CEM Problem

0: Triangulation $M$ with vertices $\\{v_{i},i=1,2,\cdots,n\\}$, tolerance
$\varepsilon$.

0: $\mathbf{f}\in\mathbb{R}^{n\times 3}$ inducing the conformal map $f$ as in
(3.1).

1: Set $k=0$ and $\delta^{(0)}=+\infty$.

2: Compute the initial guess $\mathbf{f}$ and the corresponding conformal
energy $E_{C}^{(0)}$ in (3.3).

3: while $\delta>\varepsilon$ do

4: Compute the gradient vector $\mathbf{g}$ and the Hessian matrix $H$ by
(3.16) and Theorem 3.

5: Solve the linear system $\hat{H}\hat{\mathbf{s}}=-\hat{\mathbf{g}}$ via
block LU decomposition to get the Newton direction $\mathbf{s}$.

6: Solve the trust region subproblem (4.5) to get the trial step $\mathbf{d}$.

7: Let $E\leftarrow E_{C}\big{(}(\bm{\theta},{\bm{\phi}})+\mathbf{d}\big{)}$.
If $E_{C}^{(k)}>E$, set $k\leftarrow k+1$ and update
$\displaystyle(\bm{\theta},{\bm{\phi}})\leftarrow(\bm{\theta},{\bm{\phi}})+\mathbf{d},$
$\displaystyle E_{C}^{(k)}\leftarrow E.$

8: Compute the error $\delta^{(k)}$ by (4.6) and tune the trust region radius
$\mit\Delta$.

9: end while

######  Remark 4.

Similar to the disk parameterization in [47], the bijectivity of the resulting
map can almost be guaranteed under the spherical coordinate representation. If
not, the folding triangles can be removed by the mean value coordinate [24].

##  5 Numerical Experiments

In this section, we describe the numerical performance of our proposed HBTR
method for spherical conformal parameterization on several triangulation
models. All experimental programs are executed in MATLAB R2021a on a personal
computer with a 2.50 GHz CPU and 64 GB RAM. Most of the triangulation models
are taken from AIM@SHAPE shape repository [2], ALICE [1], Gu’s personal
website [7], the Stanford 3D scanning repository [5], Human Connectome Project
[4], and TurboSquid [6]. The triangulations of brain cortical surfaces are
generated from BraTS datasets [11] via library JIGSAW [20, 19, 18, 21, 17] and
toolbox Iso2Mesh [23, 44, 8]. From the abovementioned benchmarks, we take the
triangulation models for experiments as shown in Figure 2 and present their
basic information for numbers of vertices and faces in Table 1. Notably, no
folding occurs on the $8$ models by HBTR. We apply the SCEM algorithm proposed
in [48] for the initial guess in the experiments. Among the vast experiments,
the SCEM can stably provide a great initial guess in a very short time, which
is appropriate for the HBTR algorithm.

Mesh | $\\#V$ | $\\#F$ | Mesh | $\\#V$ | $\\#F$
---|---|---|---|---|---
Apple | 17839 | 35674 | Fandisk | 6475 | 12946
Arnold | 14530 | 29056 | Horse | 21013 | 42022
Brain | 32160 | 64316 | Planck | 51108 | 102212
Bunny | 55684 | 111364 | Venus | 14303 | 28602
Table 1: The features of triangulation models with $\\#V$ and $\\#F$ being the number of vertices and triangle faces, respectively. |  |  |
---|---|---|---
Apple | Arnold | Brain | Bunny
|  |  |
Fandisk | Horse |  Planck | Venus
Figure 2: The triangulation models for the experiments.

###  5.1 Convergence behavior and conformal distortion

We first present the convergence behavior of the proposed Algorithm 1. Figure
3 shows the relationship between the number of iterations $k$ and conformal
energy $E_{C}^{(k)}$, the infinity norm of gradient
$\|\mathbf{g}^{(k)}\|_{\infty}$ and the error $\delta^{(k)}$ by Algorithm 1
for models in Figure 2. As shown, the conformal energy $E_{C}^{(k)}$ with the
scale on the right decreases linearly first. Meanwhile,
$\|\mathbf{g}^{(k)}\|_{\infty}$ and $\delta^{(k)}$ with the scale on the left
remain stable within a range. At this moment, the gradient and Newton
directions are utilized for the trial step. Then, $E_{C}^{(k)}$ tends to level
off, while $\|\mathbf{g}^{(k)}\|_{\infty}$ and $\delta^{(k)}$ descend sharply.
More specifically, $\|\mathbf{g}^{(k)}\|_{\infty}$ becomes $10^{-10}$ on most
of the models when the iteration loops terminate, implying that the iteration
stops at a critical point and the algorithm converges. Moreover,
$\delta^{(k)}$ descends in quadratic order. Taking Brain as an example, we
have

$\delta^{(25)}=7.4\times 10^{-4},\quad\delta^{(26)}=9.3\times
10^{-6},\quad\delta^{(27)}=5.9\times 10^{-10}.$

Obviously, it shows the quadratic convergence of the HBTR algorithm.

|
---|---
|

Figure 3: The relationship between the number of iterations $k$ and the
conformal energy $E_{C}^{(k)}$, the infinity norm of gradient
$\|\mathbf{g}^{(k)}\|_{\infty}$ and the error $\delta^{(k)}$ by Algorithm 1.
The $x$-axes represent the number of iterations. The left $y$-axes represent
the conformal energies, and the right $y$-axes represent infinity norms of
gradient and errors, respectively.

Figures 4 and 5 present the histograms of absolute angle distortion (degree)
of each angle, denoted as $|\alpha_{jk}-\alpha_{jk}(f)|$ for $v_{i}$ in
$T_{ijk}$, and Beltrami coefficients $\mu$ [15] of triangular faces,
respectively. If $|\mu|=0$, the map $f$ is conformal. The subplots in the
upper right of each histogram are the front and back of the angle distortion
distributions/Beltrami coefficient distributions on the resulting unit
spheres, respectively. One can see that most of angle distortions are less
than $5$ degrees and most of Beltrami coefficients are less than $0.1$ as
well, guaranteeing the conformal performance of the HBTR algorithm. Angle
distortion has similar performance to the Beltrami coefficient. Furthermore,
Figure 6 shows the absolute value of the discrete Gauss curvature $|\kappa|$,
the average angle distortion $\epsilon_{\alpha}$ and the average of the norm
of the Beltrami coefficient $\overline{|\mu|}$ at each vertex for the models
Arnold, Brain, Fandisk and Horse, denoted as

$\displaystyle\kappa(v_{i})=$ $\displaystyle 2\pi-\sum_{\\{j,k\\}\in
S_{\mathcal{E}}(i)}\alpha_{jk},$ $\displaystyle\epsilon_{\alpha}(v_{i})=$
$\displaystyle\frac{1}{|\mathcal{N}(i)|}\sum_{\\{j,k\\}\in
S_{\mathcal{E}}(i)}|\alpha_{jk}-\alpha_{jk}(f)|,$
$\displaystyle\overline{|\mu|}(v_{i})=$
$\displaystyle\frac{1}{|\mathcal{N}(i)|}\sum_{T_{ijk}\ni i}|\mu(T_{ijk})|,$

where $|\mathcal{N}(i)|$ is the number of adjacent vertices of $v_{i}$ and
$\mu(T_{ijk})$ is the Beltrami coefficient on triangle $T_{ijk}$. We
approximately find
$\overline{|\mu|}\propto\epsilon_{\alpha}\propto\sqrt{|\kappa|}$ in these $4$
models except Brain, which approximately satisfies
$\overline{|\mu|}\propto\epsilon_{\alpha}\propto|\kappa|$. The figures are
plotted according to the relationships. The textures illustrate the high
similarity of high curvatures, large angle distortions and large Beltrami
coefficient distributions. In other words, the large angle distortion and the
large Beltrami coefficient regions are mainly at those with high curvatures,
such as the ears, eyes and nose of Arnold, corners of Fandisk, and ears and
legs of Horse. The small angle distortion and small Beltrami coefficient
regions are roughly at those with low curvature. Therefore, the HBTR performs
relatively poorly at vertices with high curvature, which is an issue in our
future work. Additionally, the angle distortion and the Beltrami coefficient
have diffusion trends. Taking Fandisk as an example, high-curvature regions
are on corners and edges, and the curvatures in other regions are mostly $0$.
Furthermore, the angle distortions and the Beltrami coefficients diffuse from
the corners to the adjacent regions gradually. These phenomena also occur in
the other models, which are not shown.

|
---|---

Figure 4: Histograms of angle distortions on triangulation models.

|
---|---

Figure 5: Histograms of the Beltrami coefficients in the triangulation models. |  |  |  |
---|---|---|---|---
(a) Gauss curvatures
|  |  |
(b) Angle distortions |
|  |  |  |
(c) Beltrami coefficients |
Figure 6: (a) Gauss curvature distributions (top), (b) angle distortion
distributions, and (c) Beltrami coefficient distributions (bottom) for models
Arnold, Brain, Fandisk and Horse.

###  5.2 Comparison with state-of-the-art algorithms

In this subsection, we compare the accuracy of the proposed HBTR algorithm
with that of two state-of-the-art algorithms for spherical conformal
parameterization, namely, FLASH [15] and the SCEM [48]. The algorithm FLASH,
simply put, via the stereographic projection, applies the composition of two
quasi-conformal maps to construct the ideal conformal map, which is not an
iterative algorithm. The MATLAB program of FLASH is obtained from Choi’s
website [3]. The SCEM algorithm, as mentioned in Section 2, adopts the north-
south hemisphere alternating iteration to compute the conformal map. The
maximum number of iterations of the HBTR and SCEM is $500$. The loop
termination condition of SCEM is that the difference between the conformal
energies of two consecutive iteration steps is less than $10^{-9}$. For the
HBTR algorithm, we set $\varepsilon=10^{-9}$, as in Figure 3, which guarantees
convergence.

Table 2 shows the comparison of conformal energies between FLASH, SCEM and
HBTR. We can see that the conformal energies by HBTR are smallest among all
models. In the view of angle distortions, it is observed that the HBTR
algorithm has well performance for $50$-th percentile and $75$-th percentile
as FLASH and SCEM. In addition, among the $8$ testing examples, FLASH, SCEM
and HBTR did not produce foldings for spherical conformal maps.

Mesh | Conformal Energy | Angle Distortion
---|---|---
$50$-th percentile |  $75$-th percentile
FLASH[15] | SCEM[48] | HBTR | FLASH[15] | SCEM[48] | HBTR | FLASH[15] | SCEM[48] | HBTR
Apple | 3.74e-04 | 3.49e-04 | 3.43e-04 | 0.162 | 0.156 | 0.153 | 0.280 | 0.268 | 0.264
Arnold | 7.82e-03 | 9.07e-03 | 4.32e-03 | 0.676 | 0.675 | 0.659 | 1.232 | 1.173 | 1.233
Brain | 2.69e-02 | 2.36e-02 | 2.33e-02 | 0.821 | 0.729 | 0.729 | 1.849 | 1.694 | 1.694
Bunny | 1.93e-03 | 1.90e-03 | 1.85e-03 | 0.343 | 0.339 | 0.363 | 0.759 | 0.746 | 0.735
Fandisk | 3.50e-02 | 1.39e-02 | 1.34e-02 | 1.071 | 0.936 | 0.936 | 1.974 | 1.646 | 1.640
Horse | 4.58e-03 | 4.30e-03 | 4.18e-03 | 1.027 | 1.020 | 1.010 | 2.292 | 2.283 | 2.283
Planck | 1.71e-03 | 1.19e-03 | 1.13e-03 | 0.311 | 0.245 | 0.231 | 0.578 | 0.457 | 0.446
Venus | 6.66e-03 | 3.54e-03 | 3.27e-03 | 0.666 | 0.453 | 0.442 | 1.186 | 0.831 | 0.821

Table 2: Comparison of conformal energies and angle distortions between FLASH,
SCEM and HBTR.

###  5.3 Convergence behavior of the discrete scheme

In this subsection, we check the numerical convergence of the discrete
conformal energy of the resulting map to the continuous energy in (2.1). We
consider two ellipsoids with semiaxis lengths of $(1.1,1,0.9)$ and
$(2.0,1,0.3)$ and generate triangular meshes with different resolutions. The
basic information of the meshes is in Table 3, where $h$ represents the
maximum diameter of triangles in the mesh. Then, we use FLASH, SCEM and HBTR
to compute conformal parameterizations. The conformal energies, means and SDs
of angle distortion are used to measure the conformal distortion of the
algorithms. Figure 7 shows the relationship between the measurements and $h$.
The $x$-axis represents $h$ and $y$-axis represents the conformal energy, mean
and SD of angle distortion, respectively. We can see that the conformal
energy, mean and SD do not stably decrease as $h$ decreases for FLASH, while
$E_{C}=O(h^{2})$ and angle distortions are linearly related to $h$ for SCEM
and HBTR. Specifically, as $h\to\frac{1}{2}h$, $E_{C}\to\frac{1}{4}E_{C}$, the
means and SDs are reduced by half. Therefore, SCEM and HBTR are robust in the
respective of convergence of the discrete scheme. The conformal energies,
means and SDs are lowest for HBTR compared with those of FLASH and SCEM,
demonstrating the robustness and accuracy of HBTR.

$\\#V$ | 642 | 2562 | 10242 | 40962 | 163842 | 655362
---|---|---|---|---|---|---
$\\#F$ | 1280 | 5120 | 20480 | 81920 | 327680 | 1310720
$h$ | $(1.1,1,0.9)$ | 0.1796 | 0.0901 | 0.0451 | 0.0226 | 0.0113 | 0.0056
$(2.0,1,0.3)$ | 0.3192 | 0.1607 | 0.0806 | 0.0403 | 0.0202 | 0.0101
Table 3: The ellipsoids meshes for checking the convergence of the discrete
scheme.

(a) Conformal energy $(1.1,1,0.9)$ (b) Mean $(1.1,1,0.9)$ (c) Standard
deviation $(1.1,1,0.9)$

(d) Conformal energy $(2.0,1,0.3)$ (e) Mean $(2.0,1,0.3)$ (f) Standard
deviation $(2.0,1,0.3)$

Figure 7: The relationship between conformal energy, mean and standard
deviation of angle distortion and $h$ of algorithms FLASH, SCEM and HBTR for
$2$ ellipsoids. The figures from the top row to the bottom row are for
ellipsoids $(1.1,1,0.9)$ and $(2.0,1,0.3)$, respectively.

###  5.4 Removement of folding triangles

HBTR does not necessarily guarantee the bijectivity of the resulting map; that
is, folding triangles may occur in the image region $\mathbb{S}^{2}$. In this
subsection, we apply a postprocessing method, named mean value coordinates
[24], to remove the folding triangles. Let $L_{MV}$ be a Laplacian matrix
defined as

$\displaystyle\big{[}L_{MV}\big{]}_{ij}=\begin{cases}-w_{MV,ij}&\text{if}~{}i\neq
j,~{}[v_{i},v_{j}]\in\mathcal{E}(M),\\\
\sum_{k\in\mathcal{N}(i)}w_{MV,ik}&\text{if}~{}i=j,\\\
0&\text{if}~{}[v_{i},v_{j}]\notin\mathcal{E}(M),\end{cases}$ (5.1)

with

$w_{MV,ij}=\frac{\tan(\alpha_{jk}/2)+\tan(\alpha_{k^{\prime}j}/2)}{\|v_{ij}\|},$

where $\alpha_{jk}$ and $\alpha_{k^{\prime}j}$ are the angles opposite to
vertex $v_{i}$ in triangles $T_{ijk}$ and $T_{k^{\prime}ji}$, respectively, as
shown in Figure 1a. The postprocessing method is concluded in Algorithm 2.

Algorithm 2 Mean value coordinates for removing folding triangles

0: Triangulation $M$ with vertices $\\{v_{i},i=1,2,\cdots,n\\}$, and
$\mathbf{f}\in\mathbb{R}^{n\times 3}$ obtained by Algorithm 1.

0: $\mathbf{f}\in\mathbb{R}^{n\times 3}$ inducing the conformal map $f$ as in
(3.1) guaranteeing the bijectivity.

1: Generate Laplacian matrix $L_{MV}$ defined in (5.1).

2: Search the folding triangles in $f(M)$. Let
$\mathtt{I}=\\{i~{}|~{}T_{ijk}\text{ is a folding triangle.}\\}$ be the index
set of vertices contained in the folding triangles and
$\mathtt{O}=\\{1,2,\cdots,n\\}\setminus\mathtt{I}$.

3: while $\mathtt{I}\neq\varnothing$ do

4: Select the face center of an unfolding triangle as north pole and perform
the stereographic projection $\Pi(\mathbf{f})\to\mathbf{h}$.

5: Update the vertices $\mathbf{h}_{\mathtt{I}}$ by solving the linear system,
$\displaystyle\big{[}L_{MV}\big{]}_{\mathtt{I}\mathtt{I}}\mathbf{h}_{\mathtt{I}}=-\big{[}L_{MV}\big{]}_{\mathtt{I}\mathtt{O}}\mathbf{h}_{\mathtt{O}}.$
(5.2)

6: Perform the inverse stereographic projection
$\Pi^{-1}(\mathbf{h})\to\mathbf{f}$ and update $\mathtt{I}$ and $\mathtt{O}$.

7: end while

Unlike $L$ in conformal energy in (3.3), $L_{MV}$ is not symmetric. However,
its weights $w_{MV}$ must be positive. Therefore, this approach can guarantee
the bijectivity of the modified map. We present $2$ examples, the resulting
maps of which by HBTR are not bijective. Then, we use Algorithm 2 to remove
the folding triangles. Table 4 shows the conformal energies, angle distortions
and the number of folding triangles before and after the removal. The
conformal energies decrease slightly, and the angle distortions are almost
unchanged, while the folding triangles disappear.

Mesh | $\\#V$ | $\\#F$ | Conformal energy | Mean of angle distortion |  $\\#$folding
---|---|---|---|---|---
Bimba | 502575 | 1005146 | 1.132e-3/ 1.113e-3 | 5.328e-3/ 5.328e-3 | 24/ 0
RightBrain | 163842 | 327680 | 2.806e-2/ 2.799e-2 | 2.252e-2/ 2.252e-2 | 34/ 0
Table 4: The result of postprocessing for removing the folding triangles. The
left and right of ’/’ are the values before and after the postprocessing,
respectively.

##  6 Application to surface registrations

Given a fixed surface ${M}_{0}$ and a series of moving surfaces
${{M}_{t},t=1,2,\cdots}$, surface registration aims to find bijective maps
from the moving surfaces to the fixed surface. It is broadly applied in
computer vision and medical imaging. The goal of surface registration is to
transform surfaces from different sources into one coordinate system.
Therefore, the registration should ensure that the predominant features in the
fixed surface correspond to those in the target surfaces, which are often
expressed as landmarks in practical applications. It is generally not easy to
manage the registration because of the complicated structure of surfaces. With
the help of parameterization, we can transform the closed fixed surface
${M}_{0}$ into a unit sphere $\mathbb{S}^{2}$ via conformal map $f_{0}$ and
then register the moving surfaces to the obtained unit sphere $\mathbb{S}^{2}$
via registration map $f_{reg,t}$. As a result, the map
$\tilde{f}_{reg,t}=f_{0}^{-1}\circ{f}_{reg,t}$ is the registration map from
${M}_{t}$ to ${M}_{0}$. To obtain $f_{reg,t}$, we consider the optimization
problem

$\displaystyle\min
E_{B}(\bm{\theta},{\bm{\phi}}):=E_{C}(\bm{\theta},{\bm{\phi}})+\lambda
E_{reg}(\bm{\theta},{\bm{\phi}}),$ (6.1)

where $E_{reg}$ is the registration loss and $\lambda$ is its parameter. The
representation of $E_{reg}$ depends on the expression of predominant features.
The most common representation is the landmark-based registration

$\displaystyle E_{reg}=\frac{1}{2|S_{L}|}\sum_{i\in
S_{L}}\|\mathbf{f}_{t,i}-\mathbf{f}_{0,i}\|_{F}^{2},$ (6.2)

where $\mathbf{f}_{t,i}$ are the feature vertices on surface $M_{t}$ and
$S_{L}$ and $|S_{L}|$ are the indices set and the number of the landmark
vertices, respectively. Landmark-based registration aims to align the landmark
vertices such that the features of the surfaces are also aligned.

The conformal energy term guarantees the conformality of the map, while the
registration loss term aligns the features of the surfaces. Therefore, the
combination of the conformal energy and the registration loss results in a
conformal (as possible) registration map. The conformal registration map is an
elastic registration and preserves the local shape of the surface, which is
widely used in the field of medical imaging. For the optimization problem
(6.1), it is easy to derive the gradient vector and Hessian matrix of
registration loss generally. Benefiting from their simple representations, we
can also utilize HBTR to solve it. Moreover, the conformal energy is invariant
up to arbitrary rotation on $\mathbb{S}^{2}$. Therefore, we introduce an
optimal rotation to further decrease the registration loss. Based on Algorithm
1, we present the following spherical conformal registration algorithm.

Algorithm 3 HBTR for spherical conformal registration

0: Fixed surface $M_{0}$ and moving surface $M_{t}$ with landmarks,
registration parameter $\lambda$, tolerance $\varepsilon$.

0: The conformal registration map $\tilde{f}_{reg,t}$.

1: Compute the spherical conformal map of $M_{0}$ by Algorithm 1, denoted as
$f_{0}$.

2: Set $k=0$ and $\delta^{(0)}=+\infty$.

3: Compute the initial guess $\mathbf{f}$ of $M_{t}$.

4: Compute the optimal rotation $R$ according to the landmarks, update
$\mathbf{f}\leftarrow\mathbf{f}R$ and compute the corresponding term
$E_{B}^{(0)}$ in (6.1).

5: while $\delta>\varepsilon$ do

6: Compute gradient vector $\mathbf{g}$ and Hessian matrix $H$ of (6.1).

7: Solve the linear system $H\mathbf{s}=-\mathbf{g}$ via block LU
decomposition.

8: Solve the trust region subproblem (4.5) to get the trial step $\mathbf{d}$.

9: Compute the optimal rotation $R$ according to the landmarks to update
$\mathbf{d}$.

10: Let $E\leftarrow E_{B}\big{(}(\bm{\theta},{\bm{\phi}})+\mathbf{d}\big{)}$.
If $E_{B}^{(k)}>E$, set $k\leftarrow k+1$ and update
$\displaystyle(\bm{\theta},{\bm{\phi}})\leftarrow(\bm{\theta},{\bm{\phi}})+\mathbf{d},$
$\displaystyle E_{B}^{(k)}\leftarrow E.$

11: Compute the error $\delta^{(k)}$ by (4.6) and tune the trust region radius
$\mit\Delta$.

12: end while

13: Let
$\mathbf{f}_{reg,t}=[\cos\bm{\theta}\odot\sin{\bm{\phi}},\sin\bm{\theta}\odot\sin{\bm{\phi}},\cos{\bm{\phi}}]$.
Compute $\tilde{f}_{reg,t}=f_{0}^{-1}\circ f_{reg,t}$, which is the
registration map of $M_{t}$.

To present the registration performance of our method, we take $5$ right brain
cortex meshes RBrain0 - RBrain4 from the Human Connectome Project [4] as an
example, which are shown in the top row of Figure 8, in which the regions in
different colors represent different parts of the brain. We select $3$
landmark curves for each brain, which are in red, green and cyan,
respectively, as shown in the top row of Figure 8. We select RBrain0 as the
fixed surface and register RBrain1 - RBrain4 to it. The parameter $\lambda$ in
(6.1) is chosen as $1,5,10$ successively. The middle row shows the resulting
spheres by solving the registration problem (6.1) with the landmark curves,
respectively, with $\lambda=5$. The bottom row shows the corresponding
registered brains. The high similarity of the landmark curves between the
fixed brain and registered brains illustrates the well performance of our
method. Notably, all registration maps are bijective. Table 5 presents their
conformal energies, angle distortions and registration losses. As the
parameter $\lambda$ increases, the conformal energies and angle distortions
remain low.

|  |  |  |
---|---|---|---|---
|  |  |  |
|  |  |  |
RBrain0 (fixed) | RBrain1 | RBrain2 | RBrain3 | RBrain4
Figure 8: The right brains for registration, the registered spheres and the registered brains from top to bottom with $\lambda=5$. The curves in red, green and cyan are the landmark curves. Mesh | $\lambda$ | Conformal | Angle distortion | Registration
---|---|---|---|---
energy | Mean | SD | loss
RBrain1 | 1 | 1.31e-03 | 3.56e-01 | 3.13e-01 | 4.95e-03
5 | 3.41e-03 | 5.53e-01 | 5.27e-01 | 3.55e-03
10 | 4.73e-03 | 6.28e-01 | 6.21e-01 | 1.39e-03
RBrain2 | 1 | 1.18e-03 | 3.50e-01 | 3.20e-01 | 4.35e-03
5 | 3.00e-03 | 5.31e-01 | 5.43e-01 | 3.90e-03
10 | 8.31e-03 | 7.55e-01 | 8.92e-01 | 3.00e-03
RBrain3 | 1 | 9.98e-04 | 3.28e-01 | 2.89e-01 | 3.87e-03
5 | 2.34e-03 | 4.69e-01 | 4.42e-01 | 2.66e-03
10 | 5.28e-03 | 5.86e-01 | 6.55e-01 | 1.83e-03
RBrain4 | 1 | 1.05e-03 | 3.35e-01 | 3.02e-01 | 2.67e-03
5 | 2.63e-03 | 4.93e-01 | 4.82e-01 | 3.04e-03
10 | 6.65e-03 | 7.38e-01 | 7.41e-01 | 2.40e-03
Table 5: The conformal energies, angle distortions and registration losses of
registered brain surfaces with $\lambda=1,5,10$, respectively.

##  7 Conclusions

In this paper, we employ spherical coordinates and directly solve the
spherical CEM problem for the computation of the surface conformal
parameterization. Then, we give the explicit derivations of the gradient
vector and the Hessian matrix of the discrete conformal energy, which
preserves the sparsity as the Laplacian matrix. Due to the sparsity of the
Hessian matrix, the robust algorithm, called HBTR, is developed to solve the
spherical CEM problem. HBTR sufficiently combines the local quadratic
convergence and continuing descent advantages of the gradient and the Newton
directions. The numerical experiments actually show the conformality, the
robustness and the local quadratic convergence of the HBTR. For the discrete
scheme, we also present the quadratic convergence of the discrete conformal
energy to the continuous scheme. Since the gradient vector and Hessian matrix
of the registration loss have simple representations, we utilize the HBTR to
propose a modified version of HBTR for the application to surface
registrations. Significantly, the modified algorithm 3 also has quadratic
convergence, suggesting its potential for extension to other applications of
our method.

## Acknowledgements

T. Li was supported in parts by the National Natural Science Foundation of
China (NSFC) 12371377. W.-W. Lin was partially supported by the Ministry of
Science and Technology (MoST 110-2115-M-A49-004), Taiwan. This work was
partially supported by National Centre of Theoretical Sciences (NCTS) in
Taiwan. We thank Tianhe-2 and the Big Data Computing Center in Southeast
University, China, for their support of our use of their computing resources.

## References

  * [1] ALICE.  http://alice.loria.fr/.  (2016).
  * [2] Digital Shape Workbench - Shape Repository.  http://visionair.ge.imati.cnr.it/ontologies/shapes/.  (2016).
  * [3] Gary P. T. Choi’s website.  https://www.math.cuhk.edu.hk/~ptchoi/.  (2023).
  * [4] Human Connectome Project.  http://www.humanconnectome.org.  (2023).
  * [5] The Stanford 3D Scanning Repository.  http://graphics.stanford.edu/data/3Dscanrep/.  (2023).
  * [6] TurboSquid.  https://www.turbosquid.com/.  (2023).
  * [7] David Xianfeng Gu’s website.  http://www3.cs.stonybrook.edu/~gu/, (2017).
  * [8] Iso2Mesh.  http://iso2mesh.sourceforge.net, (2018).
  * [9] Emil Praun Alla Sheffer and Kenneth Rose.  Mesh parameterization methods and their applications.  Foundations and Trends in Computer Graphics and Vision, 2(2):105–171, 2006.
  * [10] Pierre Alliez, Giuliana Ucelli, Craig Gotsman, and Marco Attene.  Recent advances in remeshing of surfaces.  In Shape Analysis and Structuring, pages 53–82. Springer Berlin Heidelberg, 2008.
  * [11] Spyridon Bakas, Ujjwal Baid, Keyvan Farahani, et al.  The International Brain Tumor Segmentation (BraTS) Cluster of Challenges, 2023.
  * [12] Lévy Bruno, Petitjean Sylvain, Ray Nicolas, and Maillot Jérome.  Least squares conformal maps for automatic texture atlas generation.  ACM Transactions on Graphics, 21(3):362–371, 2002.
  * [13] Richard H. Byrd, Robert B. Schnabel, and Gerald A. Shultz.  Approximate solution of the trust region problem by minimization over two-dimensional subspaces.  Mathematical Programming, 40(1-3):247–263, 1988.
  * [14] Gary P. T. Choi, Yusan Leung-Liu, Xianfeng Gu, and Lok Ming Lui.  Parallelizable global conformal parameterization of simply-connected surfaces via partial welding.  SIAM Journal on Imaging Sciences, 13(3):1049–1083, 2020.
  * [15] Pui Tung Choi, Ka Chun Lam, and Lok Ming Lui.  Flash: Fast landmark aligned spherical harmonic parameterization for genus-0 closed brain surfaces.  SIAM Journal on Imaging Sciences, 8(1):67–94, 2015.
  * [16] Mathieu Desbrun, Mark Meyer, and Pierre Alliez.  Intrinsic parameterizations of surface meshes.  Computer Graphics Forum, 21(3):209–218, 2002.
  * [17] Darren Engwirda.  Locally optimal Delaunay-refinement and optimisation-based mesh generation.  PhD thesis, School of Mathematics and Statistics, The University of Sydney, 2014.
  * [18] Darren Engwirda.  Voronoi-based point-placement for three-dimensional delaunay-refinement.  Procedia Engineering, 124:330–342, 2015.
  * [19] Darren Engwirda.  Conforming restricted delaunay mesh generation for piecewise smooth complexes.  Procedia Engineering, 163:84–96, 2016.
  * [20] Darren Engwirda.  Generalised primal-dual grids for unstructured co-volume schemes.  Journal of Computational Physics, 375:155–176, 2018.
  * [21] Darren Engwirda and David Ivers.  Off-centre steiner points for delaunay-refinement on curved surfaces.  Computer-Aided Design, 72:157–171, 2016.
  * [22] Hossein Eskandari.  Strictly conformal transformation optics for directivity enhancement and unidirectional cloaking of a cylindrical wire antenna.  Scientific Reports, 12(1), 2022.
  * [23] Qianqian Fang and David A. Boas.  Tetrahedral mesh generation from volumetric binary and grayscale images.  In 2009 IEEE International Symposium on Biomedical Imaging: From Nano to Macro. IEEE, 2009.
  * [24] Michael S. Floater.  Mean value coordinates.  Computer Aided Geometric Design, 20(1):19–27, 2003.
  * [25] Michael S. Floater and Kai Hormann.  Surface parameterization: a tutorial and survey.  In Advances in Multiresolution for Geometric Modelling, pages 157–186. Springer Berlin Heidelberg, 2005.
  * [26] Mark Gillespie, Boris Springborn, and Keenan Crane.  Discrete conformal equivalence of polyhedral surfaces.  ACM Transactions on Graphics, 40(4):1–20, 2021.
  * [27] Xianfeng Gu, Yalin Wang, Tony F. Chan, Paul M. Thompson, and Shing-Tung Yau.  Genus zero surface conformal mapping and its application to brain surface mapping.  IEEE Transactions on Medical Imaging, 23(8):949–958, 2004.
  * [28] Xianfeng Gu and Shing-Tung Yau.  Computational Conformal Geometry, volume 3 of Advanced Lectures in Mathematics.  International Press and Higher Education Press, 2007.
  * [29] Xianfeng Gu and Shing-Tung Yau.  Computational Conformal Geometry.  International Press and Higher Education Press, 2020.
  * [30] Xianfeng David Gu, Wei Zeng, Feng Luo, and Shing-Tung Yau.  Numerical computation of surface conformal mappings.  Computational Methods and Function Theory, 11(2):747–787, 2011\.
  * [31] Krishna Chaitanya Gurijala, Rui Shi, Wei Zeng, Xianfeng Gu, and Arie Kaufman.  Colon flattening using heat diffusion riemannian metric.  IEEE Transactions on Visualization and Computer Graphics, 19(12):2848–2857, 2013.
  * [32] Steven Haker, Sigurd Angenent, Allen Tannenbaum, Ron Kikinis, Guillermo Sapiro, and Michael Halle.  Conformal surface parameterization for texture mapping.  IEEE Transactions on Visualization and Computer Graphics, 6(2):181–189, 2000.
  * [33] Kai Hormann, Bruno Lévy, and Alla Sheffer.  Mesh parameterization: Theory and practice.  In ACM SIGGRAPH Course Notes, 2007.
  * [34] Wei-Qiang Huang, Xianfeng David Gu, Tsung-Ming Huang, Song-Sun Lin, Wen-Wei Lin, and Shing-Tung Yau.  High performance computing for spherical conformal and riemann mappings.  Geometry, Imaging and Computing, 1(2):223–258, 2014.
  * [35] John E. Hutchinson.  Computing conformal maps and minimal surfaces.  Proceedings of the Centre for Mathematics and its Applications, 26:140–161, 1991.
  * [36] Miao Jin, Junho Kim, Feng Luo, and Xianfeng Gu.  Discrete surface ricci flow.  IEEE Transactions on Visualization and Computer Graphics, 14(5):1030–1043, 2008.
  * [37] Yueh-Cheng Kuo, Wen-Wei Lin, Mei-Heng Yueh, and Shing-Tung Yau.  Convergent conformal energy minimization for the computation of disk parameterizations.  SIAM Journal on Imaging Sciences, 14(4):1790–1815, 2021.
  * [38] Rongjie Lai, Zaiwen Wen, Wotao Yin, Xianfeng Gu, and Lok Ming Lui.  Folding-free global conformal mapping for genus-0 surfaces by harmonic energy minimization.  Journal of Scientific Computing, 58(3):705–725, 2014.
  * [39] Jingxin Nie, Tianming Liu, Gang Li, Geoffrey Young, Ashley Tarokh, Lei Guo, and Stephen T.C. Wong.  Least-square conformal brain mapping with spring energy.  Computerized Medical Imaging and Graphics, 31(8):656–664, 2007\.
  * [40] Ulrich Pinkall and Konrad Polthier.  Computing discrete minimal surfaces and their conjugates.  Experimental Mathematics, 2(15):15–36, 1993.
  * [41] Rohan Sawhney and Keenan Crane.  Boundary first flattening.  ACM Transactions on Graphics, 37(1):5:1–5:14, 2017.
  * [42] Rui Shi, Wei Zeng, Zhengyu Su, Jian Jiang, Hanna Damasio, Zhonglin Lu, Yalin Wang, Shing-Tung Yau, and Xianfeng Gu.  Hyperbolic harmonic mapping for surface registration.  IEEE Transactions on Pattern Analysis and Machine Intelligence, 39(5):965–980, 2017.
  * [43] Gerald A. Shultz, Robert B. Schnabel, and Richard H. Byrd.  A family of trust-region-based algorithms for unconstrained minimization with strong global convergence properties.  SIAM Journal on Numerical Analysis, 22(1):47–67, 1985.
  * [44] Anh Phong Tran, Shijie Yan, and Qianqian Fang.  Improving model-based functional near-infrared spectroscopy analysis using mesh-based anatomical and light-transport models.  Neurophotonics, 7(01):1, 2020.
  * [45] Lin Xu and Huanyang Chen.  Conformal transformation optics.  Nature Photonics, 9(1):15–23, 2014.
  * [46] Yong-Liang Yang, Ren Guo, Feng Luo, Shi-Min Hu, and Xianfeng Gu.  Generalized discrete ricci flow.  Computer Graphics Forum, 28(7):2005–2014, 2009.
  * [47] Mei-Heng Yueh, Tsung-Ming Huang, Tiexiang Li, Wen-Wei Lin, and Shing-Tung Yau.  Projected gradient method combined with homotopy techniques for volume-measure-preserving optimal mass transportation problems.  Journal of Scientific Computing, 88(3), 2021.
  * [48] Mei-Heng Yueh, Tiexiang Li, Wen-Wei Lin, and Shing-Tung Yau.  A novel algorithm for volume-preserving parameterizations of 3-manifolds.  SIAM Journal on Imaging Sciences, 12(2):1071–1098, 2019.
  * [49] Mei-Heng Yueh, Wen-Wei Lin, Chin-Tien Wu, and Shing-Tung Yau.  An efficient energy minimization for conformal parameterizations.  Journal of Scientific Computing, 73(1):203–227, 2017.