Smooth digital terrain modeling in irregular domains using finite element thin plate splines and adaptive refinement

1.   Introduction

This paper proposes a new method for computing spherical area-preserving mappings of topological spheres using a Riemannian optimization framework.

Area-preserving mappings can serve as parameterization of surfaces and have been applied to the field of computer graphics and medical imaging ^[6]. State-of-the-art methods for the computation of area-preserving parameterizations include diffusion-based methods ^[8], optimal mass transportation-based methods ^[18], and nonlinear optimization methods ^[20]. It is worth noting that most previous works consider area-preserving parameterization from simply connected open surfaces to planar regions. For the spherical area-preserving mapping of genus-zero closed surfaces, recently developed methods include the adaptive area-preserving parameterization method ^[7], the spherical optimal transportation mapping ^[10], and the stretch energy minimization method ^[21].

Riemannian optimization generalizes unconstrained, continuous optimization (defined on Euclidean space) following the same principle that Riemannian geometry is a generalization of Euclidean space. In traditional optimization methods, nonlinear matrix manifolds rely on the Euclidean vector space structure. In contrast, in the Riemannian optimization framework, algorithms and convergence analysis are formulated based on the language of differential geometry, and numerical linear algebra techniques are then used for implementation. Riemannian optimization focuses on matrix manifolds, which are manifolds in the sense of classical Riemannian geometry that admit a natural representation of their elements in the form of matrices. Retraction mappings play a vital role in this framework, as they are used to turn objective functions defined in Euclidean space into functions defined on an abstract manifold so that constraints are taken into account explicitly. A vast body of literature on the theory and implementation of Riemannian optimization now exists. Some of the earliest references for this field, especially for line-search methods, go back to Luenberger ^[15], Gabay ^[14], and Udrişte ^[19]. More recent literature that encompasses the findings and developments of the last three decades is ^[2,5,12].

Our Riemannian gradient descent (RGD) method involves classical components from computational geometry (simplicial surfaces and mappings) and Riemannian optimization (line search, retractions, and Riemannian gradients). To the best of the authors' knowledge, this is the first time this approach has been used in computational geometry. In this paper, we minimize the normalized stretch energy ^[20] for simplicial mappings subject to the constraint that the image of the mapping belongs to a power manifold of $n$ unit spheres embedded in $\mathbb{R}^{3}$.

1.1.   Contributions

In this paper, we propose an RGD method to compute spherical area-preserving mappings on a power manifold of unit spheres. In particular, the main contributions are as follows:

(i) We combine the tools from the Riemannian optimization framework and the components from computational geometry to propose an RGD method for computing spherical area-preserving mappings of topological spheres.

(ii) We explore two different line-search strategies: one using MATLAB's $\texttt{fminbnd}$, and the other using the quadratic/cubic interpolant approximation from [11,§ 6.3.2].

(iii) We conduct extensive numerical experiments on several mesh models to demonstrate the accuracy and efficiency of the algorithm.

(iv) We demonstrate the competitiveness and efficiency of our algorithm over three state-of-the-art methods for computing area-preserving mappings.

(v) We show that our algorithm is stable with respect to small perturbations in the initial mesh model.

(vi) We apply the algorithm to the concrete application of landmark-aligned surface registration between two human brain models.

1.2.   Outline of the paper

The remaining part of this paper is organized as follows. Section 2 introduces the main concepts on simplicial surfaces and mappings and presents the formulation of the objective function. Section 3 provides some preliminaries on the Riemannian optimization framework, briefly recalls the geometry of the unit sphere, and then describes the tools needed to perform optimization on the power manifold. Section 4 is about the RGD method and reports known convergence results. Section 5 discusses the extensive numerical experiments to evaluate our algorithm in terms of accuracy and efficiency, and provides a concrete application. Finally, we wrap our paper with conclusions and a future outlook in Section 6. The appendices contain other details about the calculations of the area of the simplicial mapping, the Hessian of the objective function, and the line-search procedure used in the RGD method.

1.3.   Notation

In this section, we list the notations and symbols adopted in order of appearance in the paper. Symbols specific to a particular section are usually not included in this list.

2.   Simplicial surfaces, mappings, and objective function

We introduce the simplicial surfaces and mappings in subsection 2.1 and the objective function in subsection 2.2.

2.1.   Simplicial surfaces and mappings

A simplicial surface ${\mathcal M}$ is the underlying set of a simplicial $2$-complex $\mathcal{K}({\mathcal M}) = {\mathcal F}({\mathcal M})\cup {\mathcal E}({\mathcal M})\cup\mathcal{V}({\mathcal M})$ composed of vertices

oriented triangular faces

and undirected edges

A simplicial mapping $f\colon {\mathcal M}\to \mathbb{R}^3$ is a particular type of piecewise affine mapping with the restriction mapping $f|_\tau$ being affine, for every $\tau\in {\mathcal F}({\mathcal M})$. We denote

The mapping $f$ can be represented as a matrix

or a vector

2.2.   The objective function

The authalic energy for simplicial mappings $f\colon {\mathcal M}\to \mathbb{R}^{3}$ is defined as ^[20]

where $E_{S}$ is the stretch energy defined as

where $I_{3}$ is the identity matrix of size 3-by-3, $\otimes$ denotes the Kronecker product, and $L_S(f)$ is the weighted Laplacian matrix $L_S(f)$. The latter is defined by

in which $\omega_S(f)$ is the modified cotangent weight defined as

with $\theta_{i, j}^k(f)$ being the angle opposite to the edge $f([v_i, v_j])$ at the point $f(v_k)$ on the image $f({\mathcal M})$, as illustrated in Figure 1.

It is proved that $E_A(f)\geq 0$ and the equality holds if and only if $f$ preserves the area [20,Corollary 3.4]. Due to the optimization process, the image area ${\mathcal A}(f)$ is not constant, hence we introduce a prefactor $| {\mathcal M}|/ {\mathcal A}(f)$ and define the normalized stretch energy as

To perform numerical optimization via the RGD method, we need to compute the Euclidean gradient. By applying the formula $\nabla E_S(f) = 2\, (I_3 \otimes L_S(f)) {\rm{vec}}(\mathbf{f})$ from [20,(3.6)], the gradient of $E(f)$ can be formulated as

The following proposition gives an explicit formula for calculating $\nabla {\mathcal A}(f)$.

Proposition 2.1 (Formula for $\nabla {\mathcal A}$). The gradient of ${\mathcal A}$ can be explicitly formulated as

Proof. By applying the explicit formulas $E_S(f|_\tau) = \frac{ {\mathcal A}(f|_\tau)^2}{|\tau|}$ and $\nabla E_S(f|_\tau) = 2 \, {\rm{vec}}(L_S(f|_\tau) \, \mathbf{f}_\tau)$ from ^[20], the chain rule yields

which is equivalent to (2.6). □

Other details about the calculation of ${\mathcal A}(f)$ are reported in Appendix 7.

3.   Riemannian optimization framework and geometry

The Riemannian optimization framework ^[2,5,12] solves constrained optimization problems where the constraints have a geometric nature. This approach utilizes the underlying geometric structure of the problems, which allows the constraints to be taken explicitly into account. The optimization variables are constrained to a smooth manifold, and the optimization is performed on that manifold. Typically, the manifolds considered are matrix manifolds, meaning there is a natural representation of their elements in matrix form. In particular, in this paper, the problem is formulated on a power manifold of $n$ unit spheres embedded in $\mathbb{R}^{3}$, and we use the RGD method for minimizing the cost function (2.4) on this power manifold.

Traditional optimization methods rely on the Euclidean vector space structure. For instance, the steepest descent method for minimizing a function $g\colon \mathbb{R}^{n} \to \mathbb{R}$ updates the current iterate $\mathbf{x}_{k}$ by moving in the direction $\mathbf{d}_{k}$ of the anti-gradient of $g$, by a step size $\alpha_{k}$ chosen according to an appropriate line-search rule. However, on a nonlinear manifold, the vector addition $\mathbf{x}_{k} + \alpha_{k} \mathbf{d}_{k}$ does not make sense in general due to the manifold curvature. We need the notion of tangent vectors and their length to generalize the steepest descent direction to a Riemannian manifold. The retraction mapping is also a vital tool in the Riemannian optimization framework. It is a mapping from the tangent space to the manifold, used to transform objective functions defined in Euclidean space into functions defined on a manifold while explicitly taking the constraints into account.

Similarly to the line-search method in Euclidean space, a line-search method in the Riemannian framework determines at a current iterate $\mathbf{x}_{k}$ on a manifold $M$ a search direction $\boldsymbol{\xi}$ on the tangent space $\mathrm{T}_{\mathbf{x}}M$. The next iterate $\mathbf{x}_{k+1}$ is then determined by a line search along a curve $\alpha \mapsto \text{R}_{ \mathbf{x}}(\alpha \boldsymbol{\xi})$ where $\text{R}_{ \mathbf{x}} \colon \mathrm{T}_{\mathbf{x}}M \to M$ is the retraction mapping. The procedure is then repeated for $\mathbf{x}_{k+1}$ taking the role of $\mathbf{x}_{k}$. Similarly to optimization methods in Euclidean space, search directions can be the negative of the Riemannian gradient, leading to the Riemannian steepest descent method. Other choices of search directions lead to other methods, e.g., Riemannian versions of the trust-region method ^[1] or BFGS ^[17].

In what follows, we introduce some fundamental geometry concepts used in Riemannian optimization, which are necessary to formulate the RGD method for our problems. We first briefly recall the geometry of the unit sphere $S^{2}$ embedded in $\mathbb{R}^{3}$, and then we switch to the power manifold of $n$ unit spheres, denoted by $\big(S^{2}\big)^{n}$.

3.1.   Geometry of the unit sphere $S^{2}$

The unit sphere $S^{2}$ is a Riemannian submanifold of $\mathbb{R}^{3}$ defined as

The Riemannian metric (inner product) on the sphere is inherited from the embedding space $\mathbb{R}^{3}$, i.e.,

where $\mathrm{T}_{\mathbf{x}}S^{2}$ is the tangent space to $S^{2}$ at $\mathbf{x} \in S^{2}$, defined as the set of all vectors orthogonal to $\mathbf{x}$ in $\mathbb{R}^{3}$, i.e.,

The projector $\text{P}_{ \mathrm{T}_{\mathbf{x}}S^{2}}$ from $\mathbb{R}^{3}$ onto the tangent space $\mathrm{T}_{\mathbf{x}}S^{2}$ is a mapping $\text{P}_{ \mathrm{T}_{\mathbf{x}}S^{2}} \colon \mathbb{R}^{3} \to \mathrm{T}_{\mathbf{x}}S^{2}$, defined by

In the following, points on the unit sphere are denoted by $\mathbf{f}_{\ell}$ (the vertices of the simplicial mapping $f$), and tangent vectors are represented by $\boldsymbol{\xi}_{\ell}$.

3.2.   Geometry of the power manifold $\big(S^{2}\big)^{n}$

We aim to minimize the function $E(f) = E(\mathbf{f}_{1}, \dots, \mathbf{f}_{n})$ (2.4), where each $\mathbf{f}_{\ell}$, $\ell = 1, \dots, n$, lives on the same manifold $S^{2}$. This leads us to consider the power manifold of $n$ unit spheres

with the metric of $S^{2}$ extended elementwise. In the remaining part of this section, we present the tools from Riemannian geometry needed to generalize gradient descent to this manifold. The projector onto the tangent space is used to compute the Riemannian gradient, as it will be explained in Section 4. The projection onto the power manifold turns points of $\mathbb{R}^{n\times 3}$ into points of $\big(S^{2}\big)^{n}$. Finally, the retraction turns an objective function defined on $\mathbb{R}^{n \times 3}$ into an objective function defined on the manifold $\big(S^{2}\big)^{n}$.

3.2.1.   Projection onto the tangent space

As seen above, the projector onto the tangent space to the unit sphere at the point $\mathbf{x}$ is given by (3.1). Here, the points are denoted by $\mathbf{f}_{\ell} \in \mathbb{R}^{3}$, $\ell = 1, \dots, n$, so we write

It clearly changes for every vertex $\mathbf{f}_{\ell}$. The projector from $\mathbb{R}^{n\times 3}$ onto the tangent space at $\mathbf{f}$ to the power manifold $\big(S^{2}\big)^{n}$ is a mapping

and can be represented by a block diagonal matrix of size $3n \times 3n$, i.e.,

In practice, we never actually create this matrix. Instead, we implement an efficient version using vectorized operations (MATLAB's $\texttt{bsxfun}$).

3.2.2.   Projection onto $\big(S^{2}\big)^{n}$

The projection of a single vertex $\mathbf{f}_{\ell}$ from $\mathbb{R}^{3}$ to the unit sphere $S^{2}$ is given by the normalization

Hence, the projection of the whole of $\mathbf{f}$ onto the power manifold $\big(S^{2}\big)^{n}$ is given by

defined by

Again, this representative matrix is only shown for illustrative purposes; in the actual implementation, we use row-wise normalization of $\mathbf{f}$.

3.2.3.   Retraction onto $\big(S^{2}\big)^{n}$

A retraction is a mapping from the tangent space to the manifold used to turn tangent vectors into points on the manifold and functions defined on the manifold into functions defined on the tangent space; see ^[2,3] for more details.

The retraction of a single tangent vector $\boldsymbol{\xi}_{\ell}$ from $\mathrm{T}_{ \mathbf{f}_{\ell}} S^{2}$ to $S^{2}$ is a mapping $\text{R}_{ \mathbf{f}_{\ell}} \colon \mathrm{T}_{ \mathbf{f}_{\ell}} S^{2} \to S^{2}$, defined by [2,Example 4.1.1]

Figure 2 provides an illustration of the retraction from $\mathrm{T}_{ \mathbf{f}_{\ell}} S^{2}$ to $S^{2}$.

For the power manifold $\big(S^{2}\big)^{n}$, the retraction of all the tangent vectors $\boldsymbol{\xi}_{\ell}$, $\ell = 1, \dots, n$, is a mapping

defined by

Again, this retraction is implemented by row-wise normalization of $\mathbf{f} + \boldsymbol{\xi}$.

4.   Riemannian gradient descent method

We are now in the position of introducing the RGD method. In this section, we will first provide the formula for the Riemannian gradient on the power manifold $\big(S^{2}\big)^{n}$. Then, we will explain the RGD method by providing its pseudocode, and finally, we will recall the known theoretical results that ensure the convergence of RGD.

The Riemannian gradient of the objective function $E$ in (2.4) is given by the projection onto $\mathrm{T}_{ \mathbf{f}} \big(S^{2}\big)^{n}$ of the Euclidean gradient of $E$, namely,

where $\nabla E$ is explicitly formulated in (2.5). This is always the case for embedded submanifolds; see [2,§ 3.6.1]. Figure 3 illustrates the difference between the Euclidean and the Riemannian gradient for one point on the unit sphere.

Given an initial iterate $f^{(0)} \in \big(S^{2}\big)^{n}$, the RGD algorithm generates a sequence of iterates $\lbrace f^{(k)} \rbrace$ as follows. At each iteration $k = 0, 1, 2, \ldots$, it chooses a search direction $\mathbf{d}^{(k)} = -\text{grad } E(f^{(k)})$ in the tangent space $\mathrm{T}_{ \mathbf{f}^{(k)}} \big(S^{2}\big)^{n}$ such that the sequence $\lbrace \mathbf{d}^{(k)} \rbrace$ is gradient related. Then, the new point $f^{(k+1)}$ is chosen such that it satisfies the sufficient decrease condition

where $c > 0$ and $\alpha_{k}$ is the step size for the given $\mathbf{d}^{(k)}$.

The RGD method on the power manifold $\big(S^{2}\big)^{n}$ is summarized in Algorithm 1. It has been adapted from [2,p. 63]. In practice, in our numerical experiments of Section 5, the initial mapping $\mathbf{f}^{(0)}\in \big(S^{2}\big)^{n}$ is computed by applying the fixed-point iteration (FPI) method of ^[21] until the first increase in energy is detected; see subsection 5.3, Algorithm 2. The line-search procedure used in line 7 is described in detail in Appendix 9.

Known convergence results

The RGD method has theoretically guaranteed convergence. For the reader's convenience, we report the two main results on the convergence of RGD. The first result is about the convergence of Algorithm 1 to critical points of the objective function. The second result is about the convergence of RGD to a local minimizer with a line-search technique.

Theorem 4.1 ([2,Theorem 4.3.1]). Let $\lbrace f^{(k)} \rbrace$ be an infinite sequence of iterates generated by Algorithm 1. Then, every accumulation point $f^{(\ast)}$ of $\lbrace f^{(k)} \rbrace$ is a critical point of the cost function $E$.

Remark 4.1. We are implicitly saying that a sequence can have more than one accumulation point; for example, from a sequence $\lbrace f^{(k)} \rbrace$, we may extract two subsequences such that they have two distinct accumulation points.

The proof of Theorem 4.1 can be done by contradiction, but it remains pretty technical, so we refer the interested reader to [2,p. 65]. It should be pointed out that Theorem 4.1 only guarantees the convergence to critical points. It does not tell us anything about their nature, i.e., whether the critical points are local minimizers, local maximizers, or saddle points. However, if the smallest eigenvalue of the Hessian of $E$ at $f^{(\ast)}$, $\lambda_{H, \min} > 0$, then the critical point $f^{(\ast)}$ is a local minimizer of $E$. Under this assumption (i.e., that $\lambda_{H, \min} > 0$), [2,§ 4.5.2] gives an asymptotic convergence bound for Algorithm 1. Indeed, the following result uses the smallest and largest eigenvalues of the Hessian of the objective function at a critical point $f^{(\ast)}$.

Theorem 4.2 ([2,Theorem 4.5.6]). Let $\lbrace f^{(k)} \rbrace$ be an infinite sequence of iterates generated by Algorithm 1, converging to a point $f^{(\ast)}$. (By Theorem 4.1, $f^{(\ast)}$ is a critical point of $E$.) Let $\lambda_{H, \min}$ and $\lambda_{H, \max}$ be the smallest and largest eigenvalues of the Hessian of $E$ at $f^{(\ast)}$. Assume that $\lambda_{H, \min} > 0$ (hence $f^{(\ast)}$ is a local minimizer of $E$). Then, given $r$ in the interval $(r_{\ast}, 1)$ with $r_{\ast} = 1-\min\left(2\sigma \bar{\alpha}\lambda_{H, \min}, 4\sigma(1-\sigma)\beta \frac{\lambda_{H, \min}}{\lambda_{H, \max}} \right)$, there exists an integer $K \geq 0$ such that

for all $k\geq K$, where $c$ is the parameter in (4.2) of Algorithm 1. Note that $0 < r_{\ast} < 1$ since $\beta, \sigma \in (0, 1)$.

As noted in [2,p. 71], in typical cases of Algorithm 1, the constant $c$ in the descent condition (4.2) equals 1, hence (4.3) reduces to $E(f^{(k+1)}) - E(f^{(\ast)}) \leq r\left(E(f^{(k)}) - E(f^{(\ast)}) \right)$.

5.   Numerical experiments

In this section, we demonstrate the convergence behavior of the RGD method using twelve mesh models and two line-search techniques. We present numerical results in tables and provide convergence plots. The computational time cost of the RGD in every table does not include the time cost for the initial mapping by the FPI method. In subsection 5.2, we introduce a correction for bijectivity, which helps to unfold the folding triangles. We then compare RGD to the FPI method of ^[21] in subsection 5.3 and the adaptive area-preserving parameterization method of ^[7] in subsection 5.4. In subsection 5.5, we compare our RGD to the spherical optimal transportation mapping proposed by Cui et al. ^[10]. Then, in subsection 5.6, we show that the algorithm is robust to noise by starting the algorithm from an initial guess with small perturbations. Finally, in subsection 5.7, we apply our algorithm to the concrete application of surface registration of two brain models.

We conducted all our experiments on a laptop Lenovo ThinkPad T460s, with Windows 10 Pro and MATLAB R2021a installed, with Intel Core i7-6600 CPU, 20GB RAM, and Mesa Intel HD Graphics 520. The benchmark triangular mesh models used in our numerical experiments are shown in Figure 4, arranged per increasing number of vertices and faces, from the top left to bottom right.

5.1.   Convergence behavior of RGD

To provide the RGD method with a good initial mapping, we first apply the FPI method of ^[21], described in subsection 5.3. Specifically, we apply the FPI until the first increase in energy occurs. We adopted two different line-search strategies: one that uses MATLAB's $\texttt{fminbnd}$ and another that uses the quadratic/cubic interpolant of [11,§ 6.3.2], described in Appendix 9.

In all the experiments, we monitor the authalic energy $E_A(f) := E_S(f) - {\mathcal A}(f)$ defined in subsection 2.2 instead of the normalized stretch energy (2.4) because when $E_{A} = 0$, we know from the theory that $f$ is an area-preserving mapping. Strictly speaking, in practice, we never obtain a mapping that exactly preserves the area, but we obtain a mapping that is area-distortion minimizing, since $E_{A}$ is never identically zero. We also monitor the ratio between the standard deviation and the mean of the area ratio. This quantity has been considered in ^[8]. Finally, the computational time is always reported in seconds, and #Fs denotes the number of folding triangles.

Tables 1 and 2 report the numerical results for RGD for minimizing the normalized stretch energy $E$, when run for a fixed maximum number of iterations (10 on the left and 100 on the right). Table 1 is for the RGD that uses the $\texttt{fminbnd}$ line-search strategy, while Table 2 is for the RGD that uses the quadratic/cubic interpolant from [11,§ 6.3.2]. Similarly to the tables, Figures 5 and 6 illustrate the convergence behavior of RGD in the same setting when run for 100 iterations.

A comparison of the computational times between the two line-search strategies shows that the line-search strategy with a quadratic/cubic interpolant (Table 2) is much more efficient than the line-search strategy that uses MATLAB's $\texttt{fminbnd}$ (Table 1). In many cases, the former even yields more accurate results than the latter. This is particularly evident for the Art Statuette (1 111.67 s versus 654.57 s) and the Bimba Statue (2 320.19 s versus 775.36 s) mesh models. The numerical results show that, in general, RGD can still decrease energy as the number of iterations increases.

Table 3 presents the results using a different stopping criterion instead of a fixed maximum number of iterations. Specifically, we keep track of the following quantity

which represents the relative improvement in the authalic energy $E_{A}(f)$ between two successive iterates, and we stop the RGD method if the value of $\Delta E_{A}^{(k)}$ is smaller than $10^{-3}$.

It is apparent from Table 3 that RGD stops at different numbers of iterations for different mesh models. Comparing with Table 2 above, it suggests that the improvement after doing more iterations is not significant for a few mesh models, especially Lion Statue, Art Statuette, and Bimba Statue. Table 3 also shows the values of the monitored quantities SD/Mean, $E_{A}(f)$, and #Fs before and after the bijectivity correction described in subsection 5.2. In all cases where the bijectivity correction is applied, the number of folding triangles #Fs decreases or reduces to zero. Moreover, the post-processing for bijectivity correction improves the SD/Mean value in all the cases. For the Right Hand mesh model, the SD/Mean value even improves significantly from 0.2063 to 0.1071 before and after the bijectivity correction. In one case (Gargoyle), the value of authalic energy $E_{A}(f)$ also improves.

Figure 7 shows the convergence behavior for the first three smallest mesh models considered, namely Right Hand, David Head, and Cortical Surface, when the algorithm is run for many more iterations; here, 10 000 iterations. It shows that RGD keeps decreasing the authalic energy $E_{A}$, albeit very slowly. The values of SD/Mean are also improved for all three mesh models. The corresponding numerical results are reported in Table 4; results for 100 iterations are also reported for easier comparison.

We compute with MATLAB the eigenvalues of the Hessian matrix at the minimizer. Table 5 reports on the smallest eigenvalue of the Hessian of the initial mapping (produced by the FPI method) and eigenvalues of the Hessian of the mapping when the prescribed maximum number of RGD iterations is achieved. The eigenvalues of the Hessian are computed by the MATLAB built-in function $\texttt{eigs}$ with the option $\texttt{smallestabs}$ and the number of eigenvalues to compute being 2 000. We observe that the Hessian eigenvalues are significantly closer to zero after running the RGD method, as we expected.

Figure 8 displays the spherical mappings resulting from applying our RGD algorithm to the benchmark mesh models considered, offering a qualitative insight into the goodness of the mappings obtained.

5.2.   Correction for bijectivity

The proposed RGD method does not guarantee the produced mapping is bijective. To remedy this drawback, we introduce a post-processing method to ensure bijectivity in the mapping. This is achieved by employing a modified version of the mean value coordinates, as described in ^[13].

Suppose we have a spherical mapping $f\colon {\mathcal M}\to S^{2}$ that may not be bijective. First, we map the spherical image to the extended complex plane $\overline{\mathbb{C}} = \mathbb{C}\cup\{\infty\}$ by the stereographic projection

Denote the mapping $h = \varPi_{ S^{2}}\circ f$ and the complex-valued vector $\mathbf{h}_i = h(v_i)$, for $v_i\in\mathcal{V}({\mathcal M})$. The folding triangular faces in the southern hemisphere are now mapped in $\mathbb{D}\subset\overline{\mathbb{C}}$, which can be unfolded by solving the linear systems

where $\mathtt{I} = \{ i \mid | h(v_i) | < r \}$ denotes the vertex index set with $r$ being a value slightly larger than $1$, e.g., $r = 1.2$, $\mathtt{B} = \{1, \ldots, n\} \backslash \mathtt{I}$, and $L_{M}$ is the Laplacian matrix defined as

with $\omega_M(h)$ being a variant of the mean value weight ^[13] defined as

in which $\varphi_{i, j}^k(h)$ is the angle opposite to the edge $[h(v_j), h(v_k)]$ at the point $h(v_i)$ on $h({\mathcal M})$, as illustrated in Figure 9. Then, the mapping is updated by replacing $\mathbf{h}_{\mathtt{I}}$ with $\widetilde{\mathbf{h}_{\mathtt{I}}}$. Next, an inversion

is performed to reverse the positions of the southern and northern hemispheres. Then, the linear system (5.2) is solved again to unfold the triangular faces originally located in the northern hemisphere. We denote the updated mapping as $\widetilde{h}$. Ultimately, the corrected mapping is given by $\varPi_{ S^2}^{-1}\circ \widetilde{h}$, where $\varPi_{ S^2}^{-1}$ denotes the inverse stereographic projection

In our numerical experiments, we perform the bijectivity correction both after the FPI method and after the RGD method, if needed.

5.3.   Comparison with the FPI method

To validate the algorithm, we compare it with the numerical results on the same models obtained using the FPI method. This method was proposed by Yueh et al. ^[21] to calculate area-preserving boundary mapping while parameterizing 3-manifolds.

The FPI method of stretch energy minimization (SEM) is performed as follows. First, a spherical harmonic mapping $f^{(0)}\colon {\mathcal M}\to S^2$ is computed by the Laplacian linearization method ^[4]. Suppose a spherical mapping $f^{(k)}$ is computed. Then, we map the spherical image to the extended complex plane by the stereographic projection (5.1) together with an inversion (5.4). Then, we define the mapping $h^{(k)}\colon {\mathcal M}\to\mathbb{C}\cup\{\infty\}$ as

The mapping $h^{(k)}$ is represented as the complex-valued vector by $\mathbf{h}^{(k)}_i = h^{(k)}(v_i)$. Next, we define the interior vertex index set as

which collects the indices of vertices in the circle centered at the origin of radius $r$. Other vertices are defined as boundary vertices and the associated index set is defined as

Then, the interior mapping is updated by solving the linear system

while the boundary mapping remains the same, i.e., $\mathbf{h}^{(k)}_{\mathtt{B}^{(k+1)}} = \mathbf{h}^{(k)}_{\mathtt{B}^{(k)}}$. The updated spherical mapping $f^{(k+1)}\colon {\mathcal M}\to S^2$ is computed by

where $\varPi_{ S^2}^{-1}$ is the inverse stereographic projection (5.5). This procedure is summarized by Algorithm 2.

Table 6 reports the numerical results for the FPI method applied to the twelve benchmark mesh models considered. Two different stopping criteria are considered: increase in authalic energy (columns 2 to 6) and maximum number of iterations (columns 7 to 10). From Table 6, it appears that performing more iterations of the FPI method does not necessarily improve the mapping $f$. In most cases, except for David Head and Cortical Surface, the authalic energy and the values of SD/Mean increase. This motivates us to use the FPI method only to calculate a good initial mapping and then switch to RGD.

5.4.   Comparison with the adaptive area-preserving parameterization

In this section, we compare the numerical results of our RGD method with the adaptive area-preserving parameterization for genus-zero closed surfaces proposed by Choi et al. ^[7]. The computational procedure is summarized as follows. First, the mesh is punctured by removing two triangular faces $\tau_1$ and $\tau_2$ that share a common edge. Then, the FPI of the SEM ^[22] is applied to compute an area-preserving initial mapping ${g}_0\colon {\mathcal M}\backslash\{\tau_1, \tau_2\}\to\mathbb{D} := \{\mathbf{x}\in\mathbb{R}^2 \mid \|\mathbf{x}\|_2\leq 1\}$. The Beltrami coefficient of the mapping ${g}_0$ is denoted by $\mu_{g_0}$. Next, a quasi-conformal map $g\colon {\mathcal M}\backslash\{\tau_1, \tau_2\}\to\mathbb{D}$ with $\|\mu_g\|_\infty = \|\lambda\mu_{g_0}\|_\infty < 1$ is constructed ^[9]. The scaling factor $\lambda\in[0, 1]$ is chosen to be $0.2$ in practice. Then, the optimal mass transport mapping $h_r\colon {\mathcal M}\backslash\{\tau_1, \tau_2\}\to r\mathbb{D}$ is computed by the method proposed by Zhao et al. ^[23], with the optimal radius $r$ satisfying

The final spherical area-preserving parameterization is obtained by the composition mapping $f = \varPi_{ S^{2}}^{-1}\circ h_r \circ rg$, where $\varPi_{ S^{2}}^{-1}$ denotes the inverse stereographic projection (5.5).

Table 7 reports on the numerical results. The algorithm of Choi et al. ^[7] was run with the default parameters found in the code package^*. The only satisfactory results are those for the David Head mesh model. However, even in this case, the values of the monitored quantities (SD/Mean, $E_{A}(f)$, and computation time) do not compete with those obtained by running our RGD method for only ten iterations; compare with the fourth column of Table 2. Our method is much more efficient and accurate and provides mappings with better bijectivity properties.

^*Available at https://www.math.cuhk.edu.hk/ptchoi/files/adaptiveareapreservingmap.zip.

5.5.   Comparison with the spherical optimal transportation mapping

In this section, we compare the numerical results of our RGD method with the spherical optimal transportation mapping proposed by Cui et al. ^[10]. The computational procedure can be summarized as follows. First, the algorithm calculates a conformal map from a genus-zero closed surface ${\mathcal M}$ to the unit sphere and then normalizes the area of ${\mathcal M}$. Iteratively, a power diagram is computed from the sphere and the power radii. The gradient of the energy is calculated from the vertices, followed by the Hessian matrix, which depends on the source measure. A linear equation is solved at every step to update the power radii, and a line search strategy is used. The step length parameter is chosen so that all cells of the power diagram are non-degenerated. If some cells are degenerated, the step length is reduced by half, and the power diagram is computed again until all cells are non-degenerated. After the loop, the centroid of each power cell in the power diagram is computed, and the mapping from each vertex to each centroid of the power cell is returned. The algorithmic details can be found in [10,Algorithm 1].

Here, we run the algorithm of ^[10] with the default parameters found in the code package^†. The target mesh for the mesh model was adapted to be its spherical harmonic mapping as suggested in ^[10].

^†Available at https://www3.cs.stonybrook.edu/gu/software/SOT/index.html

From Table 8, it appears that the method of ^[10] is capable, at least when working, of computing bijective spherical mappings within a relatively moderate number of iterations. However, the area-preserving properties, as measured by SD/Mean and $E_{A}(f)$, are inferior to those obtained with our RGD method for all four mesh models considered.

5.6.   Numerical stability

In this section, we investigate the numerical stability of our scheme. To this aim, we introduce Gaussian random noise to the vertex normal of each vertex in every mesh model according to a given value of noise variance $\sigma_{\mathrm{noise}}$. We then re-run the entire algorithm, i.e., we first perform a few iterations of the FPI method (Algorithm 2) to obtain a mapping that is used as an initial mapping for the RGD method (Algorithm 1). We then calculate a relative error on the authalic energy, as defined below in (5.6). We repeat this procedure for different values of noise variance $\sigma_{\mathrm{noise}}$.

Figure 10 shows the original mesh model of the Lion Statue (panel (a)) and two noisy versions (panels (b) and (c)).

We compute the following relative error on the authalic energy

where ${E}_{A}\big(\tilde{f}\big)$ and $E_{A}(f)$ are the authalic energies after 100 iterations of RGD for the mesh model with and without noise, respectively, and $\widetilde{v}$ and $v$ denote the coordinates of the vertices of the mesh model with and without noise, respectively.

Tables 9 and 10 report the numerical results for all the noisy mesh models, with $\sigma_{\mathrm{noise}} = 1 \times 10^{-3}$ and $5 \times 10^{-3}$, respectively. The values of authalic energy for the non-noisy mesh models from Table 2 are reported in the second last column for easier comparison. We observe that the values for SD/Mean and $E_{A}\big(\tilde{f}\big)$ remain bounded and reasonable with respect to the original mesh models, demonstrating that our method is stable to noise.

5.7.   Registration problem between two brain surfaces

A registration mapping between surfaces ${\mathcal M}_0$ and ${\mathcal M}_1$ refers to a bijective mapping $g \colon {\mathcal M}_0\to {\mathcal M}_1$. An ideal registration mapping keeps important landmarks aligned while preserving specified geometry properties. In this section, we demonstrate a framework for the computation of landmark-aligned area-preserving parameterizations of genus-zero closed surfaces.

Suppose a set of landmark pairs $\{(p_i, q_i) \mid p_i\in {\mathcal M}_0, \ q_i\in {\mathcal M}_1\}_{i = 1}^m$ is given. The goal is to compute an area-preserving simplicial mapping $g\colon {\mathcal M}_0\to {\mathcal M}_1$ that satisfies $g(p_i) \approx q_i$, for $i = 1, \ldots, m$. First, we compute area-preserving parameterizations $f_0\colon {\mathcal M}_0\to S^{2}$ and $f_1\colon {\mathcal M}_1\to S^{2}$ of surfaces ${\mathcal M}_0$ and ${\mathcal M}_1$, respectively. The simplicial registration mapping $h \colon S^{2} \to S^{2}$ that satisfies $h\circ f_0(p_i) = f_1(q_i)$, for $i = 1, \ldots, m$, can be carried out by minimizing the registration energy

Let

be the matrix representation of $h$. The gradient of $E_R$ with respect to $\mathbf{h}$ can be formulated as

where $\mathbf{r}$ is the matrix of the same size as $\mathbf{h}$ given by

In practice, we define the midpoints $c_i$ of each landmark pairs on $S^{2}$ as

for $i = 1, \ldots, m$, and compute $h_0$ and $h_1$ on $S^{2}$ that satisfy $h_0\circ f_0(p_i) = c_i$ and $h_1\circ f_1(q_i) = c_i$, respectively. Ultimately, the registration mapping $g\colon {\mathcal M}_0\to {\mathcal M}_1$ is obtained by the composition mapping $g = f_1^{-1} \circ h_1^{-1} \circ h_0 \circ f_0$. Figure 11 schematizes this composition of functions for the landmark-aligned morphing process from one brain to another.

A landmark-aligned morphing process from ${\mathcal M}_0$ to ${\mathcal M}_1$ can be constructed by the linear homotopy $H\colon {\mathcal M}_0\times[0, 1]\to\mathbb{R}^3$ defined as

In Figure 12, we demonstrate the morphing process from one brain to another brain by four snapshots at four different values of $t$. The brain surfaces are obtained from the source code package in ^[9].

6.   Conclusions and outlook

In this paper, we introduced an RGD method for computing spherical area-preserving mappings of genus-zero closed surfaces. Our approach combines the tools of Riemannian optimization and computational geometry to develop a method based on the minimization of the normalized stretch energy. The proposed algorithm has theoretically guaranteed convergence and is accurate, efficient, and robust. We tested two different line-search strategies and conducted extensive numerical experiments on various mesh models to demonstrate the algorithm's stability and effectiveness. By comparing with three existing methods for computing area-preserving mappings, we demonstrated that our algorithm is more efficient than these state-of-the-art methods. Moreover, we show that our approach is stable and robust even when the mesh model undergoes small perturbations. Finally, we applied our algorithm to the practical problem of landmark-aligned surface registration between two human brain models.

There are some directions in which we could conduct further research. Specifically, we would like to enhance the speed of convergence of the algorithm we have proposed while keeping the computational cost low. One potential way to improve this would be to use appropriate Riemannian generalizations of the conjugate gradient method or the limited memory BFGS (L-BFGS) method, as suggested in ^[17]. Another research direction may target genus-one or higher genus closed surfaces.

Appendix A.   Calculation of the gradient of the image area

The image area of the simplicial mapping $f$ can be calculated as follows. Let $\tau = [v_i, v_j, v_k]$ and denote $f_{ij}^\ell = f_{i}^\ell - f_{j}^\ell$, $f_{ik}^\ell = f_{i}^\ell - f_{k}^\ell$, and $f_{jk}^\ell = f_{ij}^\ell - f_{k}^\ell$, for $\ell = 1, 2, 3$. The image area of a simplicial mapping $f$ can be formulated as

where

The functionals ${\mathcal A}_{12}$, ${\mathcal A}_{13}$ and ${\mathcal A}_{23}$ measure the image area of mappings $\varPi_{P_{xy}}\circ f$, $\varPi_{P_{xz}}\circ f$ and $\varPi_{P_{yz}}\circ f$, respectively. The partial derivatives of ${\mathcal A}(f|_{\tau})$ can be formulated as

Appendix B.   Calculation of the Hessian of the stretch energy functional

In this appendix, we describe the calculation of the Hessian of the stretch energy functional. Let $\tau = [v_i, v_j, v_k]$. From [20,Theorem 3.5 and (3.1)], we know that

where

and $\mathbf{f}|_\tau^\ell = (\mathbf{f}_i^\ell, \mathbf{f}_j^\ell, \mathbf{f}_k^\ell)^\top$, $\ell = 1, 2, 3$. A direct calculation yields that the Hessian matrix

can be formulated as

where $\mathbf{h}_{ijk}^\ell = (f_j^\ell-f_k^\ell, f_k^\ell-f_i^\ell, f_i^\ell-f_j^\ell) ^{\top}$, $\ell = 1, 2, 3$.

Appendix C.   Line-search procedure of the RGD method

In this appendix, we describe the line-search procedure used in our RGD method. We start by detailing how to compute the derivative appearing in the sufficient decrease condition, and then we describe the interpolant line-search strategy from [11,§ 6.3.2].

In Euclidean space $\mathbb{R}^{n}$, the steepest descent method updates a current iterate by moving in the direction of the anti-gradient, by a step size chosen according to an appropriate line-search rule. The step size is related to the sufficient decrease condition ^[16]

In the Euclidean setting, the univariate function $\phi(\alpha)$ in the line-search procedure is

where $E$ is the objective function considered, and $\mathbf{f}^{(k)}$ is the current iterate. However, this function changes in the Riemannian optimization framework because we cannot directly perform the vector addition $\mathbf{f}^{(k)} + \alpha \, \mathbf{d}^{(k)}$ without leaving the manifold.

Let $\mathbf{f}^{(k)}$ be a point of $\big(S^{2}\big)^{n}$ at the $k$th iteration of our RGD algorithm, and $\text{R}_{ \mathbf{f}^{(k)}}$ the retraction at $\mathbf{f}^{(k)}$, as defined in subsection 3.2.3. Let the objective function be $E \colon \big(S^{2}\big)^{n} \to \mathbb{R}$. For a fixed point $\mathbf{f}^{(k)}$ and a fixed tangent vector $\mathbf{d}^{(k)}$, we introduce the vector-valued function $\boldsymbol{\psi} \colon \mathbb{R} \to \big(S^{2}\big)^{n}$, defined by $\alpha \mapsto \text{R}_{ \mathbf{f}^{(k)}}(\alpha\, \mathbf{d}^{(k)})$. Hence, by composition of functions, we have $\phi \colon \mathbb{R} \to \mathbb{R}$, defined by

Note that $\phi(0) = E(\boldsymbol{\psi}(0)) = E(\mathbf{f}^{(k)}) =: E^{(k)}$, due to the definition of retraction.

To evaluate the sufficient decrease condition (9.1), we need to calculate $\phi'(0)$, and to compute $\phi'(0)$ we need to calculate the derivative of the retraction $\text{R}_{ \mathbf{f}^{(k)}} (\alpha \, \mathbf{d}^{(k)})$ with respect to $\alpha$. The derivative of $\phi(\alpha)$ is given by the chain rule

and then we evaluate it at $\alpha = 0$, i.e.,

In general, the retraction and the derivative $\boldsymbol{\psi}'(\alpha)$ depend on the choice of the manifold. The differentials of a retraction provide the so-called vector transports, and the derivative of the retraction on the unit sphere appears in [2,§ 8.1.2].

In the following formulas, we omit the superscript $^{(k)}$ referring to the current iteration of RGD and assume that the line-search direction $\mathbf{d}$ is also partitioned as the matrix $\mathbf{f}$; see (2.1). Recalling the retraction on the power manifold $\big(S^{2}\big)^{n}$ from (3.3), we can write $\boldsymbol{\psi}(\alpha)$ as

The derivative $\boldsymbol{\psi}'(\alpha)$ can be computed via the formula (cf. [2,Example 8.1.4])

At $\alpha = 0$, this simplifies into

This value is needed to evaluate the sufficient decrease condition (9.1).

Quadratic/cubic approximation

At every step of our RGD algorithm, we want to satisfy the sufficient decrease condition (9.1). Here, we adopt the safeguarded quadratic/cubic approximation strategy described in [11,§ 6.3.2].

In the following, we let $\alpha_{k}$ and $\alpha_{k-1}$ denote the step lengths used at iterations $k$ and $k-1$ of the optimization algorithm, respectively. We denote the initial guess using $\alpha_{0}$. We suppose that the initial guess is given; alternatively, one can use [16,(3.60)] as initial guess, i.e.,

If $\alpha_{0}$ satisfies the sufficient decrease condition (9.1), i.e.,

then $\alpha_{0}$ is accepted as step length, and we terminate the search. Otherwise, we build a quadratic approximation $\phi_{\mathrm{q}}(\alpha)$ of $\phi(\alpha)$ using the information we have, that is, $\phi(0)$, $\phi'(0)$, and $\phi(\alpha_{0})$. The quadratic model is

The new trial value $\alpha_{1}$ is defined as the minimizer of this quadratic, i.e.,

We terminate the search if the sufficient decrease condition is satisfied at $\alpha_{1}$, i.e.,

Otherwise, we need to backtrack again. We now have four pieces of information about $\phi(\alpha)$, so it is desirable to use all of them. Hence, at this and any subsequent backtrack steps during the current iteration of RGD, we use a cubic model $\phi_{\mathrm{c}}(\alpha)$ that interpolates the four pieces of information $\phi(0)$, $\phi'(0)$, and the last two values of $\phi(\alpha)$, and set $\alpha_{k}$ to the value of $\alpha$ at which $\phi_{\mathrm{c}}(\alpha)$ has its local minimizer.

Let $\alpha_{\mathrm{prev}}$ and $\alpha_{\mathrm{2prev}}$ be the last two previous values of $\alpha_{k}$ tried in the backtrack procedure. The cubic that fits $\phi(0)$, $\phi'(0)$, $\phi(\alpha_{\mathrm{prev}})$, and $\phi(\alpha_{\mathrm{2prev}})$ is

where

The local minimizer of this cubic is given by [11,(6.3.18)]

and we set $\alpha_{k}$ equal to this value. If necessary, this process is repeated, using a cubic interpolant of $\phi(0)$, $\phi'(0)$, and the two most recent values of $\phi$, namely, $\phi(\alpha_{\mathrm{prev}})$ and $\phi(\alpha_{\mathrm{2prev}})$, until a step size that satisfies the sufficient decrease condition is located. Numerical experiments in Section 5 demonstrate the usage and effectiveness of this line-search technique.

Use of AI tools declaration

The authors declare they have not used Artificial Intelligence (AI) tools in the creation of this article.

Author contributions

Marco Sutti and Mei-Heng Yueh: conceptualization, data curation, formal analysis, funding acquisition, investigation, methodology development, project administration, resource provision, software development, supervision, validation, visualization, writing of the original draft, review, and editing. All authors of this article have been contributed equally. All authors have read and approved the final version of the manuscript for publication.

Acknowledgments

The work of the first author was supported by the National Center for Theoretical Sciences in Taiwan (R.O.C.) under the NSTC grant 112-2124-M-002-009-. The work of the second author was supported by the National Science and Technology Council and the National Center for Theoretical Sciences in Taiwan.

Conflict of interest

The authors declare that they have no conflict of interest.