Ellipsoidal Density-Equalizing Map for Genus-0 Closed Surfaces

Abstract

Surface parameterization is a fundamental task in geometry processing and plays an important role in many science and engineering applications. In recent years, the density-equalizing map, a shape deformation technique based on the physical principle of density diffusion, has been utilized for the parameterization of simply connected and multiply connected open surfaces. More recently, a spherical density-equalizing mapping method has been developed for the parameterization of genus-0 closed surfaces. However, for genus-0 closed surfaces with extreme geometry, using a spherical domain for the parameterization may induce large geometric distortion. In this work, we develop a novel method for computing density-equalizing maps of genus-0 closed surfaces onto an ellipsoidal domain. This allows us to achieve ellipsoidal area-preserving parameterizations and ellipsoidal parameterizations with controlled area change. We further propose an energy minimization approach that combines density-equalizing maps and quasi-conformal maps, which allows us to produce ellipsoidal density-equalizing quasi-conformal maps for achieving a balance between density-equalization and quasi-conformality. Using our proposed methods, we can significantly improve the performance of surface remeshing for genus-0 closed surfaces. Experimental results on a large variety of genus-0 closed surfaces are presented to demonstrate the effectiveness of our proposed methods.

1 Introduction

Surface parameterization is the process of mapping a complicated surface onto a simpler domain. It is closely related to many tasks in geometry processing including surface registration, texture mapping, surface remeshing, and shape analysis. In recent decades, numerous efforts have been devoted to the development of efficient surface parameterization algorithms with different desired properties [1, 2, 3, 4] for various science and engineering applications.

As genus-0 closed surfaces are topologically equivalent to a sphere, most prior works have focused on the problem of parameterizing genus-0 closed surfaces onto the unit sphere. For conformal parameterization, a vast number of spherical conformal parameterization methods have been developed over the past few decades based on linearization [5], harmonic energy minimization [6], Ricci flow [7], conformal curvature flow [8], quasi-conformal theory [9], partial welding [10] and Dirichlet energy minimization [11]. For area-preserving parameterization, the theory and computation of optimal mass transportation (OMT) have been extensively studied [12, 13]. Based on OMT, some area-preserving parameterization methods for genus-0 closed surfaces have been developed [14, 15, 16, 17]. Besides the above-mentioned conformal and area-preserving methods, some other prior works have also developed spherical parameterization algorithms with different mapping criteria [18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28]. While there is a topological equivalence between genus-0 closed surfaces and the sphere, the geometric difference between them may be large. For instance, mapping a genus-0 closed surface with an elongated shape to a sphere may induce large geometric distortion. Therefore, a few recent studies have focused on ellipsoidal parameterizations [29, 30, 31], in which an ellipsoid is used as the parameter domain for parameterizing genus-0 closed surfaces.

Density-equalizing maps [32] are a class of mappings with shape deformation produced based on prescribed density distribution. More specifically, given a planar domain with a density distribution prescribed in it, a density-equalizing map will deform the domain with the high-density regions enlarged and the low-density regions shrunk. Over the past few decades, density-equalizing maps have been widely used in cartogram creation and data visualization [33]. Also, different variants and improved algorithms have been developed [34, 35, 36]. Choi and Rycroft [37] developed a novel method for computing surface density-equalizing maps and demonstrated its effectiveness in surface parameterization for simply connected open surfaces. Later, the method was extended for other mapping problems for simply connected open surfaces [38, 39, 40] and volumetric domains [41] with applications to medical visualization and shape modeling. Lyu et al. [42] developed a surface parameterization method for multiply connected open surfaces by combining density-equalizing maps and quasi-conformal maps. More recently, they have developed spherical density-equalizing mapping methods for genus-0 closed surfaces [43]. Using their methods, spherical area-preserving parameterizations and spherical parameterizations with controlled area change can be achieved. However, as mentioned above, mapping genus-0 closed surfaces with extreme geometry onto a spherical domain may induce large geometric distortion. For instance, the spherical area-preserving parameterization of an elongated surface may possess low area distortion but extremely high angular distortion, which is undesirable for many practical applications.

Motivated by the above works, here we develop a novel method for computing ellipsoidal density-equalizing maps for genus-0 closed surfaces in $\mathbb{R}^{3}$. Specifically, given a genus-0 closed surface, we compute density-equalizing maps onto ellipsoidal domains with different prescribed radii based on a density distribution encoding the desired mapping effect (Fig. 1(a)). This allows us to easily achieve ellipsoidal area-preserving parameterizations and ellipsoidal parameterizations with controlled area change. We then further propose an algorithm for computing ellipsoidal density-equalizing quasi-conformal maps, which combines density-equalizing maps and quasi-conformal maps to achieve a balance between area and angle distortions. Moreover, throughout the algorithm, we can optimize the shape of the target ellipsoidal domain to further reduce overall geometric distortion (Fig. 1(b)). This provides us with an effective and automatic way of representing any genus-0 closed surface by an optimal ellipsoid. We apply our proposed methods for surface remeshing of genus-0 closed surfaces and demonstrate the improvement over prior approaches via various examples.

The organization of this paper is as follows. In Section 2, we introduce the mathematical background of our work. In Section 3, we describe our proposed methods for computing ellipsoidal density-equalizing maps (EDEM) and ellipsoidal density-equalizing quasi-conformal maps (EDEQ). In Section 4, we present experimental results of our proposed methods on various genus-0 closed surfaces. In Section 5, we describe the application of our proposed methods to genus-0 surface remeshing. Finally, in Section 6, we conclude this paper and discuss possible future works.

Figure 1: An illustration of the proposed ellipsoidal density-equalizing mapping method (EDEM) and ellipsoidal density-equalizing quasi-conformal mapping method (EDEQ). (a) Given any input genus-0 closed surface, we can apply the proposed EDEM method to compute ellipsoidal density-equalizing maps to ellipsoidal domains with different prescribed elliptic radii. (b) Given any input genus-0 closed surface (top left), we can start with an initial spherical parameterization (top right) and then apply the proposed EDEQ method to simultaneously optimize the elliptic radii of the domain and the mapping onto it, thereby achieving an ellipsoidal density-equalizing quasi-conformal map with minimal geometric distortions (bottom).

2 Mathematical background

Our proposed methods are primarily based on the theory of density-equalizing maps and quasi-conformal geometry. The relevant concepts are introduced in this section.

2.1 Density-equalizing maps

Gastner and Newmann [32] proposed a method for computing the density-equalizing maps for 2D planar domains based on the principle of density diffusion. Given a positive density function $\rho$ defined on the planar domain, the method produces shape deformation following the density gradient. More specifically, note that the advection equation is given by

$$\frac{\partial\rho}{\partial t}=-\nabla\cdot\mathbf{j},$$ (1)

where $\mathbf{j}=-\nabla\rho$ is the density flux by Fick’s law. This gives the diffusion equation

$$\frac{\partial\rho}{\partial t}=\Delta\rho.$$ (2)

Now, since the flux $\mathbf{j}$ can be expressed as $\mathbf{j}=\rho\mathbf{v}$, where $\mathbf{v}$ is the velocity field, we have

$$\mathbf{v}=\frac{\mathbf{j}}{\rho}=-\frac{\nabla\rho}{\rho}.$$ (3)

Therefore, the position of any tracer particle $\mathbf{r}$ at time $t$ can be traced by:

$$\mathbf{r}(t)=\mathbf{r}(0)+\int^{t}_{0}\mathbf{v}(\mathbf{r},\tau)\mathrm{d\tau}.$$ (4)

where $\mathbf{r}(0)$ is the initial position. Taking the time $t\to\infty$, the density $\rho$ will be fully equalized on the entire domain, and the resulting shape deformation is a density-equalizing map. In particular, as the shape deformation is induced by the density diffusion process, it is easy to see that high-density regions will be enlarged throughout the process and low-density regions will be shrunk (see Fig. 2 for an illustration).

Figure 2: An illustration of density-equalizing maps. The density flow enlarges the regions with high density and shrinks the regions with low density.

The above density-equalizing mapping method has been widely applied to cartogram creation and sociological data visualization. Specifically, the input density $\rho$ is commonly defined as some prescribed quantity (known as the “population”) per unit area, where the “population” can either be the actual population of a certain region on the geographical map, the average income of a region, or any other sociological data to be visualized.

In recent years, the density-equalizing mapping method has been further introduced to the field of surface parameterization. In particular, Choi et al. [37, 38] proposed methods for computing the density-equalizing maps for simply connected open surfaces in $\mathbb{R}^{3}$. Lyu et al. [42] developed parameterization methods for multiply connected open surfaces based on density-equalizing maps. More recently, Lyu et al. [43] developed a method for computing spherical density-equalizing maps for genus-0 closed surfaces.

2.2 Quasi-conformal theory

It is well-known that conformal maps [44] are angle-preserving. Intuitively, they map infinitesimal circles to infinitesimal circles. Quasi-conformal maps [45, 46] are a generalization of conformal maps taking infinitesimal circles to infinitesimal ellipses with bounded eccentricity. Mathematically, a quasi-conformal map $f:\overline{\mathbb{C}}\rightarrow\overline{\mathbb{C}}$ satisfies the Beltrami equation:

$$\frac{\partial f}{\partial\bar{z}}=\mu(z)\frac{\partial f}{\partial z}$$ (5)

for some complex-valued function $\mu$ satisfying $\|\mu\|_{\infty}<1$. Here, $\mu$ is called the Beltrami coefficient, which encodes important information about the conformal distortion of the map $f$. Specifically, $f$ is a conformal map if and only if $\mu=0$, as Eq. (5) becomes the Cauchy–Riemann equation. Intuitively, around a point $z_{0}\in\mathbb{C}$, the first order approximation of $f$ can be expressed as:

$$f(z)\approx f(z_{0})+f_{z}(z_{0})(z-z_{0})+f_{\bar{z}}(z_{0})\overline{(z-z_{0})}=f(z_{0})+f_{z}(z_{0})(z-z_{0}+\mu(z_{0})\overline{(z-z_{0})}).$$ (6)

The above approximation suggests that $f$ maps an infinitesimal circle centered at $z_{0}$ to an infinitesimal ellipse centered at $f(z_{0})$. Additionally, we can determine the angles and scaling factors of the maximal magnification and maximal shrinkage using $\mu$ (see Fig. 3). More specifically, the angle of maximal magnification is given by $\operatorname{arg}(\mu(z_{0}))/2$, with the magnifying factor $|f_{z}(z_{0})|(1+|\mu(z_{0})|)$. Also, the angle of maximal shrinkage is given by $(\operatorname{arg}(\mu(z_{0}))+\pi)/2$, with the shrinking factor $|f_{z}(z_{0})(|1-|\mu(z_{0})|)$. The maximal dilation of $f$ is $K(f)=\frac{1+\|\mu\|_{\infty}}{1-\|\mu\|_{\infty}}$.

Besides encoding important geometric information about the associated quasi-conformal map, the Beltrami coefficient is also related to the bijectivity of the mapping [45, 46]:

$\mathbf{Theorem}$ 2.1.

If $f$ is a $C^{1}$ mapping satisfying $\|\mu_{f}\|_{\infty}<1$, then $f$ is bijective.

Figure 3: An illustration of quasi-conformal maps. The Beltrami coefficient determines the conformality distortions

Moreover, the Beltrami coefficient of a composition of two quasi-conformal maps can be expressed in the following way. Let $f,g:\overline{\mathbb{C}}\rightarrow\overline{\mathbb{C}}$ be two quasi-conformal maps with Beltrami coefficient $\mu_{f}$ and $\mu_{g}$, respectively. The Beltrami coefficient of the composition map $g\circ f$ is given by

$$\mu_{g\circ f}=\dfrac{\mu_{f}+(\mu_{g}\circ f)\tau}{1+\overline{\mu_{f}}(\mu_{g}\circ f)\tau},$$ (7)

where $\tau=\overline{f_{z}}/f_{z}$. In particular, if $g$ is a conformal map, then $\mu_{g\circ f}=\mu_{f}$. In other words, composing a conformal map with a given quasi-conformal map will not change its Beltrami coefficient.

Lui et al. [47] proposed an algorithm called Linear Beltrami Solver (LBS) for efficiently reconstructing a quasi-conformal mapping based on a given Beltrami coefficient $\mu$. Specifically, for any given $\mu=\rho+i\tau$, the corresponding quasi-conformal map $f=u+iv$ can be computed by solving the following equations:

$$\left\{\begin{array}[]{l}\nabla\cdot(A\nabla u)=0\\ \nabla\cdot(A\nabla v)=0\end{array}\right.$$ (8)

where $A=\left(\begin{array}[]{cc}\alpha_{1}&\alpha_{2}\\ \alpha_{2}&\alpha_{3}\end{array}\right)$ with

$$\alpha_{1}=\frac{(\rho-1)^{2}+\tau^{2}}{1-\rho^{2}-\tau^{2}}\ \ \alpha_{2}=-\frac{2\tau}{1-\rho^{2}-\tau^{2}},\ \ \alpha_{3}=\frac{1+2\rho+\rho^{2}+\tau^{2}}{1-\rho^{2}-\tau^{2}}.$$ (9)

In the discrete case, Eq. (8) becomes two sparse symmetric positive definite linear systems and can be easily solved numerically. In the following, we denote the above method for reconstructing a quasi-conformal map $f$ from a Beltrami coefficient $\mu$ by $f=\textbf{LBS}(\mu)$.

3 Proposed methods

In this section, we first propose a novel algorithm for computing ellipsoidal density-equalizing maps (abbreviated as EDEM) of genus-0 closed surfaces onto a prescribed ellipsoid. Next, we propose another algorithm for computing ellipsoidal density-equalizing quasi-conformal maps (abbreviated as EDEQ), which allows us to simultaneously optimize the shape of the ellipsoidal domain and the mapping onto it to achieve minimal geometric distortions.

3.1 Ellipsoidal density-equalizing map (EDEM)

Consider a given genus-0 closed surface $\mathcal{M}$ in $\mathbb{R}^{3}$ discretized as a triangle mesh $\left(\mathcal{V},\mathcal{E},\mathcal{F}\right)$, where $\mathcal{V}$ is the set of vertices, $\mathcal{E}$ is the set of edges, and $\mathcal{F}$ is the set of triangle faces. Our goal is to compute a bijective ellipsoidal density-equalizing map $f:\mathcal{M}\rightarrow\mathcal{E}_{a,b,c}$, where $\mathcal{E}_{a,b,c}$ is a prescribed ellipsoid with elliptic radii $a,b,c$:

$$\mathcal{E}_{a,b,c}=\left\{(x,y,z)\in\mathbb{R}^{3}:\ \frac{x^{2}}{a^{2}}+\frac{y^{2}}{b^{2}}+\frac{z^{2}}{c^{2}}=1\right\}.$$ (10)

Here, the density $\rho$ involved in the density-equalizing mapping method will be defined by some prescribed positive quantity (denoted as the “population”) per unit area. By changing the “population”, we can easily control the shape deformation and achieve different desired parameterization results.

3.1.1 Initial ellipsoidal parameterization

First, we compute an initial ellipsoidal conformal parameterization $f_{0}:\mathcal{M}\rightarrow\mathcal{E}_{a,b,c}$.

To accomplish this, we apply the Fast Ellipsoidal Conformal Mapping (FECM) method [30], which is an efficient algorithm for computing an ellipsoidal conformal parameterization by leveraging quasi-conformal theory. Initially, the input genus-0 closed surface $\mathcal{M}$ is conformally mapped to the unit sphere $\mathbb{S}^{2}$. Then, the FECM method maps the sphere onto the complex plane using the stereographic projection, followed by a Möbius transformation step that maps two desired polar points of $\mathcal{M}$ to $0$ and $\infty$, respectively. However, it is important to note that the stereographic projection may lead to an uneven distribution of points. To address this issue, the FECM method incorporates an additional rescaling step to enhance the distribution of points.

Next, the FECM method looks for an inverse projection that maps the planar domain to the target ellipsoid $\mathcal{E}_{a,b,c}$ conformally. It is worth noting that while one can extend the ordinally stereographic projection and its inverse to the ellipsoidal case, the inverse ellipsoidal stereographic projection is not inherently conformal. Therefore, by utilizing the idea of quasi-conformal composition, the FECM method identifies a suitable additional quasi-conformal map on the plane and composes it with the inverse ellipsoidal stereographic projection to form a conformal map from the plane to the ellipsoid $\mathcal{E}_{a,b,c}$.

Finally, by combining all of the above-mentioned mappings, the FECM method produces an ellipsoidal conformal parameterization $f_{0}:\mathcal{M}\rightarrow\mathcal{E}_{a,b,c}$. Readers are referred to [30] for more details.

3.1.2 Density-equalizing map on ellipsoid

After obtaining an initial ellipsoidal parameterization, we handle the subsequent shape deformation problem on the ellipsoidal domain $\mathcal{E}_{a,b,c}$.

We define the initial density $\rho$ as the “population” per unit area, where the “population” is some prescribed positive quantity defined everywhere on the ellipsoid. In the following, we propose an iterative scheme for computing ellipsoidal density-equalizing maps on $\mathcal{E}_{a,b,c}$, with the bijectivity of the mappings enforced. As the density diffusion process is time-dependent, here we denote the position of a vertex $\mathbf{r}$ on $\mathcal{E}_{a,b,c}$ at time $t$ as $\mathbf{r}(t)$. The initial position of $\mathbf{r}$ is represented as $\mathbf{r}(0)$.

Analogous to the original density-equalizing mapping formulation [32], the diffusion equation of $\rho$ on the ellipsoid is given by:

$$\frac{\partial\rho(\mathbf{r}(t),t)}{\partial t}=\Delta\rho(\mathbf{r}(t),t),$$ (11)

where $\Delta$ is the Laplace–Beltrami operator. The velocity field induced by the density gradient on the ellipsoid is given by

$$\mathrm{v}(\mathbf{r}(t),t)=-\frac{\nabla\rho(\mathbf{r}(t),t)}{\rho(\mathbf{r}(t),t)}.$$ (12)

Then, the position of any vertex on the ellipsoid can be computed using the following equation:

$$\mathbf{r}(t)=\mathbf{r}(0)+\int^{t}_{0}v(\mathbf{r}(\tau),\tau)d\tau.$$ (13)

As time progresses towards infinity ($t\rightarrow\infty$), the density is equalized and we obtain the desired ellipsoidal density-equalizing map.

Note that in the discrete case, the initial density $\rho$ needs to be discretized on the ellipsoidal triangle mesh. Here, we first define the initial density $\rho^{0}_{\mathcal{F}}(T)=\frac{\text{Population}}{\text{Area}(T)}$ on every triangular face $T$ on the ellipsoid. The subsequent iterative process involves updating the positions of all vertices $\left\{\mathrm{r}_{i}(t_{n})\right\}_{i}$, where $i=1,2,\dots,|\mathcal{V}|$ and at different time points $t_{n}$, where $n=1,2,\dots$ with the time step size $\delta t=t_{n+1}-t_{n}$. In addition to the face density $\rho^{n}_{\mathcal{F}}$, we also denote the vertex density as $\rho^{n}_{\mathcal{V}}$. Note that $\rho^{n}_{\mathcal{V}}$ can be obtained from $\rho^{n}_{\mathcal{F}}$ by a simple formula [37]:

$$\rho^{n}_{\mathcal{V}}=M\rho^{n}_{\mathcal{F}},$$ (14)

where $M$ is a $|\mathcal{V}|\times|\mathcal{F}|$ face-to-vertex conversion matrix such that:

$$M_{ij}=\left\{\begin{array}[]{ll}\frac{\text{Area}(T_{j})}{\sum_{T\in N^{\mathcal{F}}(\mathbf{r}_{i})}\text{Area}(T)}&\text{ if $T_{j}$ is incident to $\mathbf{r}_{i}$},\\ 0&\text{ otherwise.}\end{array}\right.$$ (15)

Here, $N^{\mathcal{F}}(\mathbf{r}_{i})$ is the 1-ring neighborhood of the vertex $\mathbf{r}_{i}$.

Now, to solve Eq. (11) in the discrete case, note that the Laplace-Beltrami operator $\Delta_{n}$ at the $n$-th iteration can be discretized as:

$$\Delta_{n}=-A^{-1}_{n}L_{n}.$$ (16)

Here, $A_{n}$ is a diagonal matrix (known as the lumped mass matrix) of size $|\mathcal{V}|\times|\mathcal{V}|$ given by:

$$(A_{n})_{ii}=\frac{1}{3}\sum_{[\mathbf{r}_{i},\mathbf{r}_{j},\mathbf{r}_{k}]\in N^{\mathcal{F}}(\mathbf{r}_{i})}\text{Area}\left([\mathbf{r}_{i}(t_{n}),\mathbf{r}_{j}(t_{n}),\mathbf{r}_{k}(t_{n})]\right),$$ (17)

where $[\mathbf{r}_{i},\mathbf{r}_{j},\mathbf{r}_{k}]$ is a triangle in $N^{\mathcal{F}}(\mathbf{r}_{i})$. It is worth noting that every element of the lumped mass matrix $A_{n}$ is the area sum of the 1-ring neighborhood of the corresponding vertex at $t=t_{n}$. The matrix $L_{n}$ is a $|\mathcal{V}|\times|\mathcal{V}|$ symmetric sparse matrix with cotangent weights [48]:

$$(L_{n})_{ij}=\left\{\begin{array}[]{ll}-\frac{1}{2}(\cot\alpha_{ij}+\cot\beta_{ij})&\text{ if }[\mathbf{r}_{i},\mathbf{r}_{j}]\in\mathcal{E},\\ -\sum_{k\neq i}(L_{n})_{ik}&\text{ if }j=i,\\ 0&\text{ otherwise,}\end{array}\right.$$ (18)

where $\alpha_{ij}$ and $\beta_{ij}$ are the two angles opposite to the edge $[\mathbf{r}_{i},\mathbf{r}_{j}]$. Now, Eq. (11) can be solved using the semi-discrete backward Euler method:

$$\frac{\rho^{n+1}_{\mathcal{V}}-\rho^{n}_{\mathcal{V}}}{\delta t}=\Delta_{n}\rho^{n+1}_{\mathcal{V}},$$ (19)

from which we have

$$\rho^{n+1}_{\mathcal{V}}=(A_{n}+\delta tL_{n})^{-1}(A_{n}\rho^{n}_{\mathcal{V}}).$$ (20)

After solving Eq. (11), the velocity field can be updated using Eq. (12). More specifically, we first compute the discretized density gradient on the triangle mesh. For any triangle element $T=[\mathbf{r}_{i},\mathbf{r}_{j},\mathbf{r}_{k}]\in\mathcal{F}$, using the Whitney 0-forms, $\rho^{n+1}_{\mathcal{V}}$ can be interpolated at any point $x$ on the triangle $T$ by:

$$\rho^{n+1}_{\mathcal{V}}(x)=\rho^{n+1}_{\mathcal{V}}(\mathbf{r}_{i})\Psi^{W}_{i}(x)+\rho^{n+1}_{\mathcal{V}}(\mathbf{r}_{j})\Psi^{W}_{j}(x)+\rho^{n+1}_{\mathcal{V}}(\mathbf{r}_{k})\Psi^{W}_{k}(x),$$ (21)

where $\Psi^{W}_{i}(x)$, $\Psi^{W}_{j}(x)$, $\Psi^{W}_{k}(x)$ are hat functions and $\rho^{n+1}_{\mathcal{V}}(\mathbf{r}_{i})$, $\rho^{n+1}_{\mathcal{V}}(\mathbf{r}_{j})$, $\rho^{n+1}_{\mathcal{V}}(\mathbf{r}_{k})$ are the vertex densities at $\mathbf{r}_{i},\mathbf{r}_{j},\mathbf{r}_{k}$. By the property $\nabla\Psi^{W}_{i}=\frac{N\times e_{ij}}{2\text{Area}(T)}$ (see Ref. [37]), the density gradient $\nabla\rho^{n+1}_{\mathcal{F}}(T)$ on the triangle $T$ is given by:

$$\nabla\rho^{n+1}_{\mathcal{F}}(T)=\frac{\mathbf{n}_{\mathcal{F}}(T)\times\left(\rho^{n+1}_{\mathcal{V}}(\mathbf{r}_{i})e_{jk}+\rho^{n+1}_{\mathcal{V}}(\mathbf{r}_{j})e_{ki}+\rho^{n+1}_{\mathcal{V}}(\mathbf{r}_{k})e_{ij}\right)}{2\text{Area}(T)},$$ (22)

where $e_{ij}$, $e_{jk}$, $e_{ki}$ are the three directed edges $[\mathbf{r}_{i}(t_{n}),\mathbf{r}_{j}(t_{n})]$, $[\mathbf{r}_{j}(t_{n}),\mathbf{r}_{k}(t_{n})]$, $[\mathbf{r}_{k}(t_{n}),\mathbf{r}_{i}(t_{n})]$, and $\mathbf{n}_{\mathcal{F}}(T)$ is the outward unit normal vector of the triangle face $T$. Then, using Eq. (23), we can obtain $\nabla\rho^{n+1}_{\mathcal{V}}$ on every vertex by:

$$\nabla\rho^{n}_{\mathcal{V}}=M\nabla\rho^{n}_{\mathcal{F}}.$$ (23)

Therefore, the velocity field is given by:

$$\mathrm{v}^{n+1}_{\mathcal{V}}(\mathbf{r}_{i})=-\frac{\nabla\rho^{n+1}_{\mathcal{V}}(\mathbf{r}_{i})}{\rho^{n+1}_{\mathcal{V}}(\mathbf{r}_{i})}.$$ (24)

However, it is noteworthy that in the discrete case, the velocity in Eq. (24) may not be located in the tangent space. To resolve this issue, we add a projection step to project the velocity onto the admissible space. Since the equation of the ellipsoid $\mathcal{E}_{a,b,c}$ is

$$\frac{x^{2}}{a^{2}}+\frac{y^{2}}{b^{2}}+\frac{z^{2}}{c^{2}}=1,$$ (25)

it is easy to see that the outward unit normal vector of $\mathcal{E}_{a,b,c}$ is

$$\mathbf{n}_{\mathcal{V}}\left(\mathbf{r}(t)\right)=\frac{\left(\frac{2x}{a^{2}},\frac{2y}{b^{2}},\frac{2z}{c^{2}}\right)}{\left(\left(\frac{2x}{a^{2}}\right)^{2}+\left(\frac{2y}{b^{2}}\right)^{2}+\left(\frac{2z}{c^{2}}\right)^{2}\right)^{\frac{1}{2}}}=\frac{\left(\frac{x}{a^{2}},\frac{y}{b^{2}},\frac{z}{c^{2}}\right)}{\left(\left(\frac{x}{a^{2}}\right)^{2}+\left(\frac{y}{b^{2}}\right)^{2}+\left(\frac{z}{c^{2}}\right)^{2}\right)^{\frac{1}{2}}}.$$ (26)

Therefore, the projected velocity field $\widetilde{\mathrm{v}}^{n+1}_{\mathcal{V}}$ can be obtained by

$$\widetilde{\mathrm{v}}^{n+1}_{\mathcal{V}}(\mathbf{r}_{i})=\mathrm{v}^{n+1}_{\mathcal{V}}(\mathbf{r}_{i})-\left(\mathrm{v}^{n+1}_{\mathcal{V}}(\mathbf{r}_{i})\cdot\mathbf{n}_{\mathcal{V}}(\mathbf{r}_{i})\right)\mathbf{n}_{\mathcal{V}}(\mathbf{r}_{i}),$$ (27)

where $\mathbf{n}_{\mathcal{V}}(\mathbf{r}_{i})$ is the outward unit normal vector at the vertex $\mathbf{r}_{i}$. We can then update the vertex positions based on Eq. (13):

$$\mathbf{r}_{i}(t_{n+1})=\mathbf{r}_{i}(t_{n})-\delta t\ \widetilde{\mathrm{v}}^{n+1}_{\mathcal{V}}(\mathbf{r}_{i}).$$ (28)

Note that there may still be small numerical errors that cause the vertices to move outside the ellipsoidal surface. To ensure that all vertices remain on the ellipsoid exactly, we further divide the coordinates of each vertex $\mathbf{r}_{i}(t_{n+1})$ by $\sqrt{\frac{x^{2}}{a^{2}}+\frac{y^{2}}{b^{2}}+\frac{z^{2}}{c^{2}}}$, where $\mathbf{r}_{i}(t_{n+1})=(x,y,z)$. This ensures that all updated vertices satisfy the parametric equation in Eq. (10) and hence lie on the ellipsoid $\mathcal{E}_{a,b,c}$.

3.1.3 Ensuring the bijectivity of the ellipsoidal mapping throughout the iterative process

As mentioned in [42, 43], conventional density-equalizing mapping methods do not provide any guarantee of bijectivity throughout the diffusion process. Specifically, mesh overlaps may occur under the mapping update in Eq. (28) if the density gradient or the time step size is too large. To resolve this issue, we follow the idea in [43] and propose a correction scheme that ensures the bijectivity of the ellipsoidal mapping at each iteration.

We first introduce a rescaling transformation $h:\mathcal{E}_{a,b,c}\to\mathbb{S}^{2}$ defined by:

$$h\left(x,y,z\right)=\left(\frac{x}{a},\frac{y}{b},\frac{z}{c}\right),$$ (29)

where $a$, $b$, and $c$ are the radii of the ellipsoid. Using the rescaling transformation $h$, the initial ellipsoidal map $\mathbf{r}(0)$ and the ellipsoidal map $\mathbf{r}(t_{n})$ after the $n$-th iterative update can be mapped onto two unit spheres $\mathcal{S}^{0}:=h(\mathbf{r}(0))$ and $\mathcal{S}^{n}:=h(\mathbf{r}(t_{n}))$, respectively. It is easy to see that $h$ is a bijection and hence will not have any effect on the presence or absence of local mesh fold-overs in the ellipsoids.

Now, we can apply the spherical overlap correction scheme proposed by the SDEM work [43], which is designed for enforcing the bijectivity of a spherical mapping. Specifically, the scheme involves projecting the two spheres $\mathcal{S}^{0}$ and $\mathcal{S}^{n}$ onto the extended complex plane via the north-pole stereographic projection and the south-pole stereographic projection and correcting the local mesh fold-overs on the plane using the LBS method [47]. Readers are referred to [43] for the detailed description of the spherical overlap correction scheme. We denote the process of the spherical overlap correction scheme as a map $\tilde{g}_{n}:\mathcal{S}^{n}\to\mathbb{S}^{2}$, where the mapping result $\tilde{g}_{n}(\mathcal{S}^{n})$ is folding-free.

Finally, we apply the inverse rescaling transformation $h^{-1}$ to rescale the spherical mapping result to the ellipsoid. Again, it is easy to see that $h^{-1}$ is a bijection and will not have any effect on the presence or absence of mesh fold-overs. Therefore, the composition $h^{-1}\circ\tilde{g}_{n}\circ h$ is a mapping from $\mathcal{E}_{a,b,c}$ to $\mathcal{E}_{a,b,c}$ that can correct any local mesh fold-overs in $\mathbf{r}(t_{n})$. This completes the ellipsoidal overlap correction scheme.

3.1.4 Re-coupling the deformation and density

Recall that the proposed ellipsoidal density-equalizing mapping scheme creates shape deformations on a prescribed ellipsoid based on density diffusion. Specifically, the velocity field at each iteration is determined by the density and its gradient at the mesh vertices. As the computational procedure involves multiple discretization schemes, numerical errors may accumulate throughout the iterations. Also, the above-mentioned overlap correction scheme may alter the vertex positions when correcting the local mesh fold-overs, thereby leading to a discrepancy between the vertex positions and the actual density flow. To address this issue, we follow the idea in [43] and introduce an extra step to re-couple the density and deformation at the end of every iteration.

More specifically, at the end of the $n$-th iteration, we do not directly use the density obtained by solving the diffusion equation in Eq. (20) to continue the next iteration. Instead, we re-define the density $\rho^{n+1}_{\mathcal{F}}(T)$ on the triangle element $T=[\mathbf{r}_{i},\mathbf{r}_{j},\mathbf{r}_{k}]$ using the updated vertex positions $\mathbf{r}(t_{n+1})$ on the ellipsoid:

$$\rho^{n+1}_{\mathcal{F}}(T)=\frac{\text{Population}}{\text{Area}[\mathbf{r}_{i}(t_{n+1}),\mathbf{r}_{j}(t_{n+1}),\mathbf{r}_{k}(t_{n+1})]}.$$ (30)

Once the updated density $\rho^{n+1}_{\mathcal{F}}$ is obtained, we calculate $\rho^{n+1}_{\mathcal{V}}$ using Eq. (14) and proceed with the next iteration. Altogether, this additional re-coupling step allows us to reduce the accumulation of numerical errors as well as the discrepancy between the shape deformation and the density.

3.1.5 Summary

By integrating the initial ellipsoidal conformal parameterization (Section 3.1.1) and the diffusion iteration scheme (Section 3.1.2) with the overlap correction scheme (Section 3.1.3) and the re-coupling scheme (Section 3.1.4), we have the proposed EDEM algorithm for computing bijective ellipsoidal density-equalizing maps for genus-0 closed surfaces. The proposed algorithm is summarized in Algorithm 1. In practice, we set the step size $\delta t=0.1$, the stopping parameter $\epsilon=10^{-3}$, and the maximum number of iterations allowed $n_{\max}=300$.

Input: A genus-0 closed surface $\mathcal{M}$, a prescribed population, the elliptic radii $a,b,c$, the step size $\delta t$, the stopping parameter $\epsilon$, and the maximum number of iterations allowed $n_{\max}$.

Output: An ellipsoidal density-equalizing map $f:\mathcal{M}\to\mathcal{E}_{a,b,c}$.

1

2Compute an initial ellipsoidal conformal parameterization $f_{0}:\mathcal{M}\to\mathcal{E}_{a,b,c}$ [30];

3

4Compute the initial density $\rho^{0}_{\mathcal{F}}$ on $f_{0}(\mathcal{M})$ based on the prescribed population;

5

6Set $n=0$;

7

8repeat

9      

10      Compute $A_{n}$ and $L_{n}$ on $\mathcal{E}_{a,b,c}$ by Eq. (17) and Eq. (18);

11      

12      Obtain $\rho^{n+1}_{\mathcal{V}}$ by solving the diffusion equation Eq. (11);

13      

14      Compute the velocity field ${\mathbf{v}}^{n+1}_{\mathcal{V}}$ at every vertex by Eq. (24);

15      

16      Compute the projected velocity field $\widetilde{\mathbf{v}}^{n+1}_{\mathcal{V}}$ onto the ellipsoid using Eq. (27);

17      

18      Update the position of all vertices using Eq. (28) and enforce the vertices to remain on $\mathcal{E}_{a,b,c}$;

19      

20      Apply the overlap correction scheme in Section 3.1.3;

21      

22      Apply the re-coupling scheme in Section 3.1.4 to update $\rho^{n+1}_{\mathcal{F}}$;

23      

24      Update $n=n+1$;

25      

26until $\frac{\text{sd}\left(\rho^{n}_{\mathcal{V}}\right)}{\text{mean}\left(\rho^{n}_{\mathcal{V}}\right)}<\epsilon$  or  $n\geq n_{\max}$;

27The resulting bijective ellipsoidal density-equalizing map is $f:\mathcal{M}\rightarrow\mathcal{E}_{a,b,c}$ ;

Algorithm 1 Ellipsoidal density-equalizing map (EDEM)

3.2 Ellipsoidal density-equalizing quasi-conformal map (EDEQ)

While our EDEM algorithm can produce shape deformations and achieve prescribed area changes on the given ellipsoid $\mathcal{E}_{a,b,c}$, other geometric distortions such as conformal distortions may be uncontrolled. Also, note that the EDEM method can be applied to ellipsoids with any prescribed elliptic radii. It is therefore natural to ask whether one can further reduce the overall geometric distortions by optimizing both the shape of the ellipsoid and the associated ellipsoidal mappings.

As described in the DEQ method [42], the density-equalizing process can be considered as an energy minimization problem involving the prescribed density. Analogously, as our EDEM method aims to create a shape deformation following the density diffusion process on the ellipsoid, the density-equalizing effect of a map $f$ can be assessed using the following energy:

$$E_{\text{EDEM}}(f)=\int\left\|\frac{\nabla\rho}{\rho}\right\|^{2},$$ (31)

where $\rho$ is the density associated with the map $f$. Also, note that by quasi-conformal theory, the conformal distortion of a map $f$ can be represented using the norm of its Beltrami coefficient $\|\mu\|$, where a smaller $\|\mu\|$ implies a smaller eccentricity of the local ellipses and hence a smaller angle distortion. Therefore, we can assess the conformal distortion of a mapping by considering the following energy:

$$E_{\text{BC}}(f)=\int|\mu|^{2}.$$ (32)

In this section, our goal is to develop an ellipsoidal density-equalizing quasi-conformal mapping algorithm, which we abbreviate as EDEQ, to reduce both of the above energies $E_{\text{EDEM}}$ and $E_{\text{BC}}$. Now, it is noteworthy that the shape of the ellipsoid $\mathcal{E}_{a,b,c}$, characterized by the elliptic radii $a,b,c$, can effectively provide us with extra degrees of freedom in reducing the energies. Therefore, we consider minimizing the following combined energy $E$ with the radii $a,b,c$ and the map $f:\mathcal{M}\to\mathcal{E}_{a,b,c}$ being the variables:

$$E(a,b,c,f)=E_{\text{EDEM}}+\alpha E_{\text{BC}}=\int\left\|\frac{\nabla\rho}{\rho}\right\|^{2}+\alpha\int|\mu|^{2},$$ (33)

where $\alpha$ is a nonnegative weighting parameter for balancing the density-equalizing error and the conformal distortion. The existence of the minimizer of $E$ can be proved by following the arguments in [42, 19].

3.2.1 Decoupling the combined energy

To simplify the optimization process, we first decouple the minimization problem of the combined energy $E$ into two subproblems of minimizing $E_{1}$ and $E_{2}$ below:

$$E_{1}(f)=\int_{\mathcal{E}_{a,b,c}}\left\|\frac{\nabla\rho}{\rho}\right\|^{2},$$ (34)

$$E_{2}(a,b,c)=\int_{\mathcal{E}_{a,b,c}}\left\|\frac{\nabla\rho}{\rho}\right\|^{2}+\alpha\int_{\mathcal{E}_{a,b,c}}|\mu|^{2}.$$ (35)

More specifically, the energy $E_{1}$ in Eq. (34) aims to equalize the density distribution across a fixed ellipsoidal domain with radii $a,b,c$, promoting a more uniform density throughout the iterations. The energy $E_{2}$ in Eq. (35) focuses on optimizing the shape of the ellipsoid to minimize the combined geometric distortion. We can then utilize our developed computational scheme in the EDEM method for the subproblem $E_{1}$ and focus on updating the ellipsoidal geometry in the subproblem $E_{2}$.

3.2.2 Subproblem $E_{1}$

We first consider the descent direction of $E_{1}=\int_{\mathcal{E}_{a,b,c}}\left\|\frac{\nabla\rho}{\rho}\right\|^{2}$, with the elliptic radii $a,b,c$ fixed. As in the formulation of the EDEM method in Section 3.1, the velocity field can be computed by $\mathbf{v}=-\frac{\nabla\rho}{\rho}$. Also, to preserve the given ellipsoidal shape, we remove the normal component $\mathbf{v}^{\perp}=\left(\mathbf{v}\cdot\mathbf{n}\right)\mathbf{n}$, where $\mathbf{n}$ is the outward normal unit vector. The descent direction of $E_{1}$ is then given by:

$$dE_{1}=\mathbf{v}-\mathbf{v}^{\perp}.$$ (36)

Hence, we have the following iterative scheme for the energy $E_{1}$:

$$f_{n+1}=f_{n}+\delta tdE^{n}_{1},$$ (37)

where $\delta t$ is the time step size, $f_{n}$ is the map at the $n$-th iteration, and $dE^{n}_{1}$ is the descent direction of $E_{1}$ at the $n$-th iteration.

3.2.3 Subproblem $E_{2}$

To minimize the energy $E_{2}$, we first consider the Beltrami coefficient term $\int_{\mathcal{E}_{a,b,c}}|\mu|^{2}$. Recall that the Beltrami coefficient is a complex-valued function defined on the complex plane. To extend it for assessing the conformal distortion of an ellipsoidal map, we use a triangle-based approach as follows. Let $T^{0}=[\mathbf{r}_{i}(t_{0}),\mathbf{r}_{j}(t_{0}),\mathbf{r}_{k}(t_{0})]$ be a triangle element on the initial ellipsoidal conformal parameterization and $T^{n}=[\mathbf{r}_{i}(t_{n}),\mathbf{r}_{j}(t_{n}),\mathbf{r}_{k}(t_{n})]$ be the corresponding triangle at the $n$-th iteration. We rigidly map both $T^{0}$ and $T^{n}$ onto the complex plane. Specifically, this can be done by fixing one of the vertices of each triangle at the origin, one of the triangle edges containing this vertex onto the real axis, and then computing the position of the remaining vertex based on the edge lengths and angles of the triangle. We denote the resulting triangles on the complex plane as $\widetilde{T}^{0}$ and $\widetilde{T}^{n}$. Then, we can consider the mapping between $\widetilde{T}^{0}$ and $\widetilde{T}^{n}$ as a quasi-conformal map $f_{\widetilde{T}^{n}}=u+iv$. The Beltrami coefficient of $f_{\widetilde{T}^{n}}$ can be obtained by:

$$\mu_{\widetilde{T}^{n}}=\frac{\left(u_{x}-v_{y}\right)+i\left(v_{x}+u_{y}\right)}{\left(u_{x}+v_{y}\right)+i\left(v_{x}-u_{y}\right)}.$$ (38)

From the above formula, we can get the norm of the Beltrami coefficient $\|\mu_{\widetilde{T}^{n}}\|$ and use it to represent the conformal distortion between $T^{0}$ and $T^{n}$. We can then repeat the above procedure for all triangular faces of the surface to get the term $\int_{\mathcal{E}_{a,b,c}}|\mu|^{2}$.

Next, we consider optimizing the shape of the ellipsoid (i.e. the radii $a,b,c$) based on the energy $E_{2}$. To simplify the formulation, we keep $a$ unchanged and focus on optimizing the two other elliptic radii $b$ and $c$ using an iterative procedure. We first define the radius step sizes along $b$ and $c$ as $\delta b$ and $\delta c$, respectively. Now, there are nine possible combinations of the radii in the form of $(a,b+k_{b}\delta b,c+k_{c}\delta c)$ with $k_{b},k_{c}=0,1,-1$. We can then consider the following energies associated with the nine combinations, including the original one:

$$\begin{array}[]{l}E_{2}(a,b,c)=\displaystyle\int_{\mathcal{E}_{a,b,c}}\left(\left\|\frac{\nabla\rho}{\rho}(a,b,c)\right\|^{2}+\alpha\left|\mu(a,b,c)\right|^{2}\right),\end{array}$$ (39)

the ones with only one among $b$ and $c$ changed:

$$\begin{array}[]{l}\displaystyle E_{2}(a,b+\delta b,c)=\int_{\mathcal{E}_{a,b+\delta b,c}}\left(\left\|\frac{\nabla\rho}{\rho}(a,b+\delta b,c)\right\|^{2}+\alpha|\mu(a,b+\delta b,c)|^{2}\right),\\ \displaystyle E_{2}(a,b,c+\delta c)=\int_{\mathcal{E}_{a,b,c+\delta c}}\left(\left\|\frac{\nabla\rho}{\rho}(a,b,c+\delta c)\right\|^{2}+\alpha|\mu(a,b,c+\delta c)|^{2}\right),\\ \displaystyle E_{2}(a,b-\delta b,c)=\int_{\mathcal{E}_{a,b-\delta b,c}}\left(\left\|\frac{\nabla\rho}{\rho}(a,b-\delta b,c)\right\|^{2}+\alpha|\mu(a,b-\delta b,c)|^{2}\right),\\ \displaystyle E_{2}(a,b,c-\delta c)=\int_{\mathcal{E}_{a,b,c-\delta b}}\left(\left\|\frac{\nabla\rho}{\rho}(a,b,c-\delta c)\right\|^{2}+\alpha|\mu(a,b,c-\delta c)|^{2}\right),\end{array}$$ (40)

and the ones with both $b$ and $c$ changed:

$$\begin{array}[]{l}\displaystyle E_{2}(a,b+\delta b,c+\delta c)=\int_{\mathcal{E}_{a,b+\delta b,c+\delta c}}\left(\left\|\frac{\nabla\rho}{\rho}(a,b+\delta b,c+\delta c)\right\|^{2}+\alpha|\mu(a,b+\delta b,c+\delta c)|^{2}\right),\\ \displaystyle E_{2}(a,b+\delta b,c-\delta c)=\int_{\mathcal{E}_{a,b+\delta b,c-\delta c}}\left(\left\|\frac{\nabla\rho}{\rho}(a,b+\delta b,c-\delta c)\right\|^{2}+\alpha|\mu(a,b+\delta b,c-\delta c)|^{2}\right),\\ \displaystyle E_{2}(a,b-\delta b,c+\delta c)=\int_{\mathcal{E}_{a,b-\delta b,c+\delta c}}\left(\left\|\frac{\nabla\rho}{\rho}(a,b-\delta b,c+\delta c)\right\|^{2}+\alpha|\mu(a,b-\delta b,c+\delta c)|^{2}\right),\\ \displaystyle E_{2}(a,b-\delta b,c-\delta c)=\int_{\mathcal{E}_{a,b-\delta b,c-\delta c}}\left(\left\|\frac{\nabla\rho}{\rho}(a,b-\delta b,c-\delta c)\right\|^{2}+\alpha|\mu(a,b-\delta b,c-\delta c)|^{2}\right).\end{array}$$ (41)

Here, the terms $\frac{\nabla\rho}{\rho}(\tilde{a},\tilde{b},\tilde{c})$ and $\mu(\tilde{a},\tilde{b},\tilde{c})$ are computed based on the updated ellipsoid $\mathcal{E}_{\tilde{a},\tilde{b},\tilde{c}}$. By selecting the combination with the smallest energy value, we can update the shape of the ellipsoid and reduce $E_{2}$. The shape update procedure is summarized in Algorithm 2.

Input: A ellipsoid $\mathcal{E}_{a,b,c}$, a prescribed population, and the step sizes $\delta b$ and $\delta c$.

Output: An ellipsoidal triangle mesh $\mathcal{E}_{\tilde{a},\tilde{b},\tilde{c}}$.

1

2Compute the energies $E_{2}(a,b,c)$, $E_{2}(a,b+\delta b,c)$, $E_{2}(a,b,c+\delta c)$, $E_{2}(a,b-\delta b,c)$, $E_{2}(a,b,c-\delta c)$, $E_{2}(a,b+\delta b,c+\delta c)$, $E_{2}(a,b+\delta b,c-\delta c)$, $E_{2}(a,b-\delta b,c+\delta c)$, $E_{2}(a,b-\delta b,c-\delta c)$;

3

4Select $(\tilde{a},\tilde{b},\tilde{c})$ that gives the smallest energy;

5

6Update the ellipsoid as $\mathcal{E}_{\tilde{a},\tilde{b},\tilde{c}}$ ;

7

Algorithm 2 Shape update of the ellipsoidal domain

3.2.4 Summary

In Sections 3.2.2 and 3.2.3, we discussed how to solve the subproblems $E_{1}$ and $E_{2}$. Combining these approaches, we have our proposed ellipsoidal density-equalizing quasi-conformal map (EDEQ) algorithm.

Specifically, analogous to the EDEM algorithm, here we start by computing an initial ellipsoidal conformal parameterization (Section 3.1.1) with radii $(a,b,c)$. The initial density and the initial combined energy can be obtained using the initial ellipsoidal parameterization. Then, to solve the subproblem $E_{1}$, we deform the initial ellipsoid based on the descent direction (Section 3.2.2) iteratively. The overlap correction scheme (Section 3.1.3) and the re-coupling scheme (Section 3.1.4) are also applied at each iteration. We repeat the above process for a certain number of iterations (set by a prescribed parameter $K$) so that the energy $E_{1}$ defined based on the initial elliptic radii is sufficiently reduced. We then include an additional shape update step (Section 3.2.3) after every $K$ iteration to modify the ellipsoidal radii by solving the subproblem $E_{2}$. With the elliptic radii updated, we solve the subproblem $E_{1}$ as described above again and repeat the above process. Also, note that a finer adjustment of the elliptic radii will be needed as the iterations continue. Therefore, instead of using fixed radius step sizes $\delta b,\delta c$ in the shape update algorithm, we add a scaling factor $0.9^{m}$ to the prescribed step sizes $\delta b,\delta c$ at the $m$-th shape update step. We repeat the iterations until $|\frac{E^{n+1}-E^{n}}{E^{n}}|<\epsilon$ for some stopping parameter $\epsilon$. The final EDEQ map is obtained by $f=f_{n}$, from which we also get the final optimal elliptic radii $(a,b,c)$.

The proposed EDEQ algorithm is summarized in Algorithm 3. In practice, the initial radii $(a_{0},b_{0},c_{0})$ can be chosen arbitrarily, the time step size is set to be $\delta t=0.1$, the initial radius step sizes are set to be $\delta b=0.1$ and $\delta c=0.1$, the number of iterations for each set of fixed radii is set to be $K=5$, the error threshold is set to be $\epsilon=10^{-5}$, and the maximum number of iterations is set to be $n_{\text{max}}=300$.

Input: A genus-0 closed surface $\mathcal{M}$, a prescribed population, the initial elliptic radii $(a_{0},b_{0},c_{0})$, the step size $\delta t$, the radius step sizes $\delta b$ and $\delta c$, the number of iterations for each set of fixed radii $K$, the stopping parameter $\epsilon$, and the maximum number of iterations allowed $n_{\max}$.

Output: An ellipsoidal density-equalizing quasi-conformal map $f:\mathcal{M}\to\mathcal{E}_{a,b,c}$ with the optimal ellipsoid $\mathcal{E}_{a,b,c}$.

1

2Set $(a,b,c)=(a_{0},b_{0},c_{0})$;

3

4Compute an initial ellipsoidal conformal parameterization $f_{0}:\mathcal{M}\to\mathcal{E}_{a,b,c}$ [30];

5

6Compute the initial density $\rho^{0}_{\mathcal{F}}$ on $f_{0}(\mathcal{M})$ based on the prescribed population;

7

8Set $n=0$;

9

10repeat

11      

12      Compute the descent direction $dE^{n}_{1}$ by Eq. (36);

13      

14      Update the mapping $f_{n+1}=f_{n}+\delta tdE^{n}_{1}$;

15      

16      Apply the overlap correction scheme in Section 3.1.3;

17      

18      Apply the re-coupling scheme in Section 3.1.4 to update $\rho^{n+1}_{\mathcal{F}}$;

19      

20      if $n=km$ for some positive integer $m$ then

21             Apply the shape update method (Algorithm 2) with the radius step sizes $0.9^{m}\delta b$ and $0.9^{m}\delta c$ and obtain the updated elliptic radii $(\tilde{a},\tilde{b},\tilde{c})$;

22            

23            Set $(a,b,c)=(\tilde{a},\tilde{b},\tilde{c})$;

24            

25      

26      Update $n=n+1$;

27      

28until $|\frac{E^{n+1}-E^{n}}{E^{n}}|<\epsilon$  or  $n\geq n_{\max}$;

29The resulting bijective ellipsoidal density-equalizing quasi-conformal map is $f:\mathcal{M}\rightarrow\mathcal{E}_{a,b,c}$, where $(a,b,c)$ are the optimal elliptic radii;

Algorithm 3 Ellipsoidal density-equalizing quasi-conformal map (EDEQ)

4 Experiments

In this section, we present experimental results to demonstrate the effectiveness of our proposed EDEM and EDEQ algorithms. The algorithms are implemented using MATLAB R2021a on the Windows platform. All experiments are conducted on a computer with an Intel(R) Core(TM) i9-12900 2.40 GHz processor and 32GB memory. The surface meshes are from online mesh repositories [49]. All surfaces are discretized in the form of triangular meshes.

4.1 Ellipsoidal density-equalizing map

To test our proposed EDEM algorithm, we first consider mapping an ellipsoidal surface with two different density distributions, including a discontinuous density distribution (Fig. 4(a)) and a continuous density distribution (Fig. 4(b)). In both examples, it can be observed that the input densities are highly non-uniform. By applying the proposed EDEM method, we obtain the corresponding ellipsoidal density-equalizing maps. It can be observed that the domains with a high initial density are enlarged while the domains with a low initial density are shrunk, which shows that our method can accurately produce shape deformations on the ellipsoid based on the prescribed density. Also, from the histograms of the initial and final densities, it can be observed that the density is effectively equalized.

Figure 4: Ellipsoidal density-equalizing maps of ellipsoidal surfaces. Each row shows one example. (a) An example with discontinuous input density. (b) An example with continuous input density. Left to right: The initial ellipsoidal surface color-coded with the initial density, the final EDEM result color-coded with the initial density, the histogram of the initial density, and the histogram of the final density.

Next, we consider some more complicated ellipsoidal mapping examples. Fig. 5(a) shows an ellipsoidal mesh with non-uniform distributed triangle elements. Specifically, the mesh is denser at the top part of the ellipsoid and much coarser at the bottom part. We define different populations on the ellipsoid such that the initial density of all the triangular regions is three times the density of all the pentagonal regions. From our EDEM mapping result, it can be observed that different regions are enlarged or shrunk accordingly. This shows that our method is capable of computing ellipsoidal density-equalizing maps on a non-uniform mesh. One may also wonder whether our method can be applied to more extreme ellipsoidal geometries and more extreme density distributions. In Fig. 5(b), we consider a more elongated ellipsoid with a more extreme density distribution prescribed on it, where the maximum and minimum density values differ by 10 times. It can be observed from our EDEM mapping result that different regions are effectively enlarged or shrunk based on the input density. This suggests that our method is capable of handling the extreme ellipsoidal geometry and extreme input density. Fig. 5(c) shows another ellipsoid with non-uniform mesh elements and a complicated continuous density distribution. Again, we can see that our proposed method produces an ellipsoidal density-equalizing map successfully. Lastly, in Fig. 5(d) we test our EDEM algorithm on a spherical surface (i.e. an ellipsoid with all three elliptic radii being identical). It can be observed from the mapping result that our method also works well for the spherical case.

Figure 5: Additional examples of ellipsoidal density-equalizing maps of ellipsoidal surfaces. For each example, the left figure shows the input surface color-coded with the initial density, and the right figure shows the EDEM result color-coded with the initial density. (a) An ellipsoid with non-uniformly distributed mesh elements and a prescribed discontinuous density. (b) An elongated ellipsoid with a prescribed continuous density. (c) An ellipsoid with non-uniformly distributed mesh elements and a complex prescribed density. (d) A sphere with a prescribed continuous density.

For a more quantitative analysis, we record the computational time, radii of the ellipsoid, variance of the initial density, variance of the final density, and number of overlaps for all the above ellipsoidal mapping examples in Table 1. It can be observed that our EDEM method can be effectively applied to a wide range of ellipsoidal shapes with different elliptic radii. Specifically, for ellipsoids with different radii and input densities, the variances of the final density in the EDEM result are all close to 0. This shows that all EDEM results are highly density-equalizing. Also, it can be observed that all mapping results are folding-free. Overall, the experiments have demonstrated the versatility of the EDEM method.

Table 1: The performance of our EDEM algorithm. For each surface, we record the number of triangle elements, the computational time, the elliptic radii $(a,b,c)$ of the ellipsoidal domain, the variance of the normalized initial density $\widetilde{\rho}_{1}=\frac{\rho_{1}}{\text{Mean}(\rho_{1})}$ (where $\rho_{1}$ is the initial vertex density), the variance of the normalized final density $\widetilde{\rho}_{2}=\frac{\rho_{2}}{\text{Mean}(\rho_{2})}$ (where $\rho_{2}$ is the final vertex density), and the number of overlaps.

Surface	# Faces	Time (s)	$\mathbf{(a,b,c)}$	$\text{Var}(\widetilde{\rho}_{1})$	$\text{Var}(\widetilde{\rho}_{2})$	# Overlaps
Ellipsoid 1 (Fig. 4(a))	7808	1.8171	(1,2,4)	0.1176	0.0032	0
Ellipsoid 2 (Fig. 4(b))	7808	4.6973	(1,2,4)	0.2084	0.0028	0
Ellipsoid 3 (Fig. 5(a))	20480	1.4869	(1,1,1.5)	0.3333	0.0137	0
Ellipsoid 4 (Fig. 5(b))	18904	0.3811	(1,0.6,2)	0.2135	0.0055	0
Ellipsoid 5 (Fig. 5(c))	20480	2.6060	(1,1,1.5)	0.2580	0.0010	0
Ellipsoid 6 (Fig. 5(d))	18904	0.6203	(1,1,1)	0.2940	0.0068	0

Next, we consider computing the ellipsoidal density-equalizing maps for more general genus-0 closed surfaces. In Fig. 6, we set the population as the face area of the original mesh and apply our EDEM algorithm to achieve ellipsoidal area-preserving parameterizations. For the initial parameterization step, we map each surface conformally onto the ellipsoid using the FECM method [30]. Then, we apply the EDEM method to compute the ellipsoidal density-equalizing parameterizations. To assess the area-preserving property of the ellipsoidal parameterizations, we consider the area distortion for each triangular face $T$ as follows:

$$d_{\text{area}}(f)(T)=\log\left(\frac{\text{Area}(f(T))/\sum_{T^{\prime}\in\mathcal{F}}\text{Area}(f(T^{\prime}))}{\text{Area}(T)/\sum_{T^{\prime}\in\mathcal{F}}\text{Area}(T^{\prime})}\right).$$ (42)

Here $f$ denotes the ellipsoidal parameterization and $\mathcal{F}$ is the set of all triangular faces. The normalization factors in the numerator and denominator ensure that a perfectly area-preserving mapping would result in $d_{\text{area}}\equiv 0$. In Fig. 6, we present various genus-0 surface models, the initial ellipsoidal parameterizations, the ellipsoidal parameterizations obtained by our EDEM method, and the initial and final area distortions. It can be observed that the EDEM method performs very well for genus-0 closed surfaces with complex structures. Specifically, by visually comparing the triangle element distributions in the initial ellipsoidal parameterizations and the EDEM parameterization results, it is easy to see that the EDEM method achieves a more uniform distribution on the prescribed ellipsoid. By comparing the initial and final area distortion histograms, we can also see that the area distortion is significantly reduced by our EDEM method. For a more quantitative analysis, Table 2 provides detailed statistics on the performance of our EDEM algorithm for ellipsoidal area-preserving parameterization. We can see that the EDEM method can efficiently reduce the area distortion by over 95% on average when compared to the initial parameterization, and the bijectivity of the parameterization is well-preserved in all examples. Altogether, the experiments show that our EDEM method is capable of computing bijective ellipsoidal area-preserving parameterizations for a large variety of genus-0 closed surfaces.

Figure 6: Ellipsoidal area-preserving parameterization of genus-0 closed surfaces obtained by the EDEM method. Each row shows one example. (a) The David model. (b) The Buddha model. (c) The hippocampus model. (d) The twisted ball model. Left to right: The input surface mesh, the initial ellipsoidal conformal parameterization obtained by the FECM method [30], the EDEM result, the histogram of the logged area ratio $d_{\text{area}}$ of the initial ellipsoidal parameterization, and the histogram of the logged area ratio $d_{\text{area}}$ of the final ellipsoidal parameterization.

Table 2: The performance of our EDEM algorithm for ellipsoidal area-preserving parameterization. For each surface, we record the number of triangle elements, the computational time, the elliptic radii $(a,b,c)$ of the target ellipsoidal domain, mean and standard deviation of the initial area distortion $|d_{\text{area}}(f_{0})|$ and the final area distortion $|d_{\text{area}}(f)|$, and the number of overlaps.

Surface	# Faces	Time (s)	$\mathbf{(a,b,c)}$	$|d_{\text{area}}(f_{0})|$		$|d_{\text{area}}(f)|$		# Overlaps
				Mean	SD	Mean	SD
David (Fig. 6(a))	21338	2.6982	(1,1.2,1.4)	0.7078	0.5248	0.0132	0.0136	0
Buddha (Fig. 6(b))	50002	10.5964	(1,1,1.4)	0.4365	0.3907	0.0366	0.0648	0
Hippocampus (Fig. 6(c))	12000	3.1002	(1,1,2.2)	1.3635	1.4740	0.0515	0.0973	0
Twisted Ball (Fig. 6(d))	38620	7.5060	(1,0.8,1.4)	0.8998	0.8195	0.0239	0.0289	0
Chinese Lion 1 (Fig. 1(a))	10000	10.1311	(1,1,0.8)	1.8306	1.8851	0.1407	0.2081	0
Chinese Lion 2 (Fig. 1(a))	10000	9.4919	(1,1,1.1)	1.6529	1.7923	0.1180	0.1759	0
Chinese Lion 3 (Fig. 1(a))	10000	7.9418	(1,1,1.4)	1.5394	1.8125	0.1137	0.1892	0

Table 3: The results obtained by our EDEM algorithm with different elliptic radii. Here, for different choices of the elliptic radii $(a,b,c)$ of the target ellipsoid, we record the mean and standard deviation of the final area distortion $|d_{\text{area}}(f)|$, the mean value of the norm of the Beltrami coefficient $|\mu|$, and the number of overlaps.

$\mathbf{a}$	$\mathbf{b}$	$\mathbf{c}$	$\text{Mean}(|d_{\text{area}}(f)|)$	$\text{SD}(|d_{\text{area}}(f)|)$	Mean$(|\mu|)$	# Overlaps
1	1	1	0.0895	0.1684	0.4443	0
	0.7	1	0.1349	0.3481	0.4620	0
	1.2	1	0.1042	0.2063	0.4446	0
	1	0.8	0.0996	0.1838	0.4709	0
	1	1.3	0.0816	0.1518	0.4017	0
	1	1.5	0.0730	0.1322	0.3730	0
	1	2.2	0.0515	0.0973	0.2862	0
	1.3	1.6	0.1484	0.3474	0.4026	0
	1.2	1.8	0.1340	0.3179	0.3942	0
	1.1	2	0.0975	0.2178	0.3718	0

It is natural to ask whether the performance of our EDEM method is affected by the shape of the target ellipsoid. Here, we map the Hippocampus model to various target ellipsoids with different elliptic radii using our EDEM method and analyze the area and angle distortions of the resulting mappings. In Table 3, we report the mean and variance of the area distortion, as well as the mean value of the norm of the Beltrami coefficient $|\mu|$. It can be observed that for all combinations of the elliptic radii, all mapping results are folding-free and possess low area distortion. This shows that our EDEM method is capable of computing area-preserving parameterizations of genus-0 closed surfaces onto different prescribed ellipsoids. However, if we also take the angle distortion into consideration, then some differences can be observed. For instance, if we increase the value of $c$, we can see that the angle distortion is further reduced. This experiment suggests that the overall ellipsoidal shape plays an important role and motivates the need for the proposed EDEQ method, which optimizes the ellipsoidal geometry and achieves ellipsoidal density-equalizing quasi-conformal maps.

Figure 7: Ellipsoidal density-equalizing quasi-conformal maps of ellipsoidal surfaces. Each row shows one example. (a) An example with discontinuous input density. (b) An example with continuous input density. Left to right: The initial ellipsoidal surface color-coded with the initial density, the final EDEQ result color-coded with the initial density, the histogram of the initial density, and the histogram of the final density.

4.2 Ellipsoidal density-equalizing quasi-conformal map

Next, we test our EDEQ algorithm for computing the ellipsoidal density-equalizing quasi-conformal mappings. In Fig. 7(a), we define different populations at different regions on the ellipsoid to obtain a discontinuous initial density. Note that the maximum density is three times larger than the minimum one. It can be observed in the EDEQ mapping result that the shape of the ellipsoid is changed under the mapping process, while the higher-density domain expands and the lower-density domain contracts. Comparing the initial and final density histograms, it can be observed that the density is highly equalized. Fig. 7(b) shows another example with a continuous initial density defined on an ellipsoid. Under the EDEQ algorithm, the ellipsoidal shape is stretched along the z-direction and becomes more elongated, and different regions are expanded or shrunk according to the given density. From the density histograms, it can again be observed that the mapping result is highly density-equalizing.

Figure 8: Ellipsoidal density-equalizing quasi-conformal maps of ellipsoidal surfaces. For each example, the left figure shows the input surface color-coded with the initial density, and the right figure shows the EDEQ result color-coded with the initial density. (a) An ellipsoid with non-uniformly distributed mesh elements and a prescribed discontinuous density. (b) An elongated ellipsoid with a prescribed continuous density. (c) An ellipsoid with non-uniformly distributed mesh elements and a complex prescribed density. (d) A sphere with a prescribed continuous density.

In Fig. 8, we further apply the EDEQ method to the same set of ellipsoidal examples in Fig. 5. It can be observed that our EDEQ method produces ellipsoidal density-equalizing mapping results with both the ellipsoidal geometry and the vertex positions optimized. For a more quantitative analysis, we record the initial and final elliptic radii, the variance of the initial and final densities, the mean value of the Beltrami coefficient $\mu$, and the number of overlaps for all the above-mentioned examples in Table 4. From the table, we can see that the ellipsoidal shapes are effectively optimized under the EDEQ method. Furthermore, in addition to possessing a highly density-equalizing effect as reflected in the variance of the final density, the mapping results also achieve a low conformal distortion. Specifically, note that the norm of the Beltrami coefficient $|\mu|$ normally ranges from 0 to 1, where $|\mu|=0$ indicates that the mapping is conformal. As the value of $|\mu|$ is small for all examples, we can see that the EDEQ method not only achieves density equalization but also reduces the conformal distortion of the ellipsoidal mappings.

Table 4: The performance of our proposed EDEQ algorithm. For each surface, we record the number of triangle elements, the initial elliptic radii $(a_{0},b_{0},c_{0})$ of the ellipsoid, the final elliptic radii $(a,b,c)$ of the ellipsoidal shape obtained by EDEQ, the variance of the normalized initial density $\widetilde{\rho}_{1}=\frac{\rho_{1}}{\text{Mean}(\rho_{1})}$ (where $\rho_{1}$ is the initial vertex density), the variance of the normalized final density $\widetilde{\rho}_{2}=\frac{\rho_{2}}{\text{Mean}(\rho_{2})}$ (where $\rho_{2}$ is the final vertex density), the mean value of norm of the Beltrami coefficient $|\mu|$, and the number of overlaps.

Surface	# Faces	$\mathbf{(a_{0},b_{0},c_{0})}$	$\mathbf{(a,b,c)}$	$\text{Var}(\widetilde{\rho}_{1})$	$\text{Var}(\widetilde{\rho}_{2})$	mean($|\mu|$)	# Overlaps
Ellipsoid 1 (Fig. 7(a))	7808	(1,1,2)	(1,1.8973,2.8163)	0.2503	0.0049	0.1434	0
Ellipsoid 2 (Fig. 7(b))	7808	(1,2,4)	(1,1.9888,4.3110)	0.2084	0.0029	0.1953	0
Ellipsoid 3 (Fig. 8(a))	20480	(1,1,1.5)	(1,1,1.4250)	0.333	0.0153	0.1852	0
Ellipsoid 4 (Fig. 8(b))	18904	(1,0.6,2)	(1,0.9821,2.5162)	0.2135	0.0082	0.2285	0
Ellipsoid 5 (Fig. 8(c))	20480	(1,1,1.5)	(1,1.0310,1.0783)	0.2580	0.0011	0.1679	0
Ellipsoid 6 (Fig. 8(d))	18904	(1,1,1)	(1,1.3095,0.9100)	0.2940	0.0072	0.2267	0

Next, we consider using the EDEQ method to compute the ellipsoidal area-preserving parameterizations for genus-0 closed surfaces. For each surface in Fig. 9, we first conformally map it onto an ellipsoid with some arbitrary elliptic radii using the FECM method [30]. Then, we apply the EDEQ algorithm with the population set as the original face area to obtain the ellipsoidal area-preserving parameterizations. It can be observed that the ellipsoidal shapes are changed significantly under the EDEQ method. Comparing the initial and final area distortion histograms, we can see that the mapping results are highly area-preserving.

Figure 9: Ellipsoidal area-preserving quasi-conformal parameterization of genus-0 closed surfaces. Each row shows one example. (a) The Duck model. (b) The Hippocampus model. (c) The Lion-Vase model. Left to right: The input surface mesh, the initial ellipsoidal conformal parameterization, the final EDEQ result, the histogram of the logged area ratio $d_{\text{area}}$ of the initial ellipsoidal parameterization, and the histogram of the logged area ratio $d_{\text{area}}$ of the final ellipsoidal parameterization.

Figure 10: The shape change of the ellipsoidal parameterization under the proposed EDEQ method. Each row shows one example. (a) The Duck model. (b) The Hippocampus model. (c) The Lion-Vase model. Left to right: The input surface mesh, the initial ellipsoidal conformal parameterization, the results after 10, 50, and 100 iterations, and the final EDEQ result.

As the choice of the initial elliptic radii is arbitrary, we can even start with a unit sphere and run the EDEQ method to eventually obtain an optimal ellipsoidal parameterization. In Fig. 10, we consider several genus-0 closed surface models and compute an initial spherical conformal parameterization. Then, we apply the proposed EDEQ algorithm to compute an optimal ellipsoidal area-preserving parameterization. It can be observed that throughout the EDEQ iterative process, both the ellipsoidal shapes and the vertex positions on the ellipsoids are changed gradually. For a more quantitative analysis, Table 5 summarizes the performance of our EDEQ algorithm in terms of the final radii, the initial and final area distortions, the conformal distortion, and the number of overlaps. Comparing the initial and final area distortions, it can be observed that the area distortion is significantly reduced in all examples. The mean of the norm of the Beltrami coefficient $|\mu|$ is also small, indicating that the conformal distortion is low. Besides, the number of overlaps is 0 for all examples, which shows that our mapping results are all folding-free. In Table 6, We further compare the performance of the EDEQ method and the EDEM method in terms of the area and conformal distortions, where both methods use the same initial elliptic radii. It can be observed that the EDEQ method achieves a significantly lower conformal distortion when compared to the EDEM method, while the area distortion remains comparable. Altogether, the above experiments demonstrate the effectiveness of our proposed EDEQ method for ellipsoidal density-equalizing quasi-conformal maps.

Table 5: The performance of our EDEQ algorithm for ellipsoidal area-preserving parameterization. For each surface, we start with an initial spherical conformal parameterization and then apply the EDEQ algorithm. The number of triangle elements, the final elliptic radii $(a,b,c)$, the mean and standard deviation of the initial area distortion $|d_{\text{area}}(f_{0})|$ and the final area distortion $|d_{\text{area}}(f)|$, the mean value of the norm of the Beltrami coefficient $|\mu|$, and the number of overlaps are recorded.

Surface	# Faces	$\mathbf{(a,b,c)}$	$|d_{\text{area}}(f_{0})|$		$|d_{\text{area}}(f)|$		mean($|\mu|$)	# Overlaps
			Mean	SD	Mean	SD
Duck (Fig. 10(a))	20000	(1,1.8797,1.8521)	1.0561	1.2190	0.0241	0.0379	0.2216	0
Hippocampus (Fig. 10(b))	12000	(1,0.1886,1.8818)	2.5651	1.7272	0.1028	0.1990	0.2602	0
Lion-Vase (Fig. 10(c))	38620	(1,0.3833,1.6136)	0.8264	0.7195	0.1032	0.1749	0.2306	0
David (Fig. 1(b))	21338	(1,1.0062,1.5126)	0.5652	0.4349	0.0125	0.0133	0.1279	0

Table 6: Comparison between the EDEM method and the EDEQ method. We record the mean and standard deviation of the final area distortion $|d_{\text{area}}(f)|$, the mean of the norm of the Beltrami coefficient $|\mu|$, and the number of overlaps for each example and each method.

Surface	EDEM				EDEQ
	Mean of $|d_{\text{area}}(f)|$	SD of $|d_{\text{area}}(f)|$	Mean$(|\mu|)$	# Overlaps	Mean of $|d_{\text{area}}(f)|$	SD of $|d_{\text{area}}(f)|$	Mean$(|\mu|)$	# Overlaps
Duck	0.0268	0.0507	0.2467	0	0.0241	0.0379	0.2216	0
Hippocampus	0.0895	0.1684	0.4443	0	0.1028	0.1990	0.2602	0
Lion-Vase	0.0909	0.1604	0.3137	0	0.1032	0.1749	0.2306	0

5 Application to genus-0 surface remeshing

Note that one important application of surface parameterization is surface remeshing [50]. Specifically, after computing the parameterization of a given surface, one can generate a high-quality mesh on the parameter domain and use the inverse of the parameterization mapping to map the mesh onto the input surface, thereby yielding a remeshed surface. For instance, Sheffer and de Sturler [51] utilized the angle-based flattening method to compute parameterizations for surface remeshing. Praun and Hoppe [52] computed spherical parameterization for surface remeshing. Later, different point cloud parameterization methods [53, 54] have been developed and applied to surface meshing. Analogously, our proposed EDEM and EDEQ methods can be easily applied to genus-0 surface remeshing.

Figure 11: The surface remeshing results of the Pig model obtained by SCM [9], SDEM [43], ECM [30], the proposed EDEM method, and the proposed EDEQ method. For each surface, two different views are provided.

Figure 12: The surface remeshing results of the Bear model obtained by SCM [9], SDEM [43], ECM [30], the proposed EDEM method, and the proposed EDEQ method. For each surface, two different views are provided.

Figure 13: The surface remeshing results of the Hippocampus model obtained by SCM [9], SDEM [43], ECM [30], the proposed EDEM method, and the proposed EDEQ method. For each surface, two different views are provided.

Figure 14: The surface remeshing results of the Lion-vase model obtained by SCM [9], SDEM [43], ECM [30], the proposed EDEM method, and the proposed EDEQ method. For each surface, two different views are provided.

In Figs. 11, 12, 13, and 14, we consider various genus-0 surface models and compare the remeshing results obtained by different parameterization methods for genus-0 closed surfaces, including spherical conformal mapping (SCM) [9], spherical density-equalizing mapping (SDEM) [43], ellipsoidal conformal mapping (ECM) [30], our ellipsoidal density-equalizing mapping (EDEM) method, and our ellipsoidal density-equalizing quasi-conformal mapping (EDEQ) method. For each genus-0 surface model $\mathcal{M}$ in $\mathbb{R}^{3}$ and each parameterization method, denote the parameterization mapping by $f:\mathcal{M}\rightarrow\mathcal{T}$, where the parameter domain $\mathcal{T}$ is either the unit sphere (for SCM and SDEM) or an ellipsoid (for ECM, EDEM and EDEQ). Then we construct a regular triangular mesh $\mathcal{T}_{\Delta}$ on $\mathcal{T}$ using the DistMesh method [55]. Finally, using the inverse mappings $f^{-1}:\mathcal{T}\to\mathcal{M}$, the regular mesh $\mathcal{T}_{\Delta}$ can be mapped back onto $\mathcal{M}$, yielding a remeshed surface $f^{-1}(\mathcal{T}_{\Delta})$. For a fair comparison, we enforce the number of vertices of the triangular mesh generated by DistMesh (and hence the number of vertices of the final remeshing result) for each method to be approximately 8500 (with variation $<\pm 1\%$).

From the remeshing results in Figs. 11, 12, 13, and 14, it can be observed that the two conformal approaches (SCM and ECM) lead to significantly non-uniform remeshing results. Moreover, the shape features near the sharp regions are lost in the SCM and ECM results. This can be explained by the fact that the conformal mappings preserve angles but may yield large area distortions. By contrast, the remeshing results obtained by the density-equalizing approaches (SDEM, EDEM, and EDEQ) are much more uniform and have a higher mesh quality. The features located near the sharp regions are well-preserved.

Figure 15: Comparison for the surface remeshing results between the SDEM method [43], the proposed EDEM method, and the proposed EDEQ method. Each row shows one example. For each method, we show a zoom-in of the remeshed surface obtained by the method to visualize the triangle quality.

To further compare the remeshing results obtained using the density-equalizing approaches, Fig. 15 shows the zoom-in images of the SDEM, EDEM and EDEQ remeshing results for the Pig, Bear, Hippocampus, and Lion-Vase models. From the zoom-in images, we can see that the EDEQ remeshing results are more uniform than the ones by the other two methods, especially at the regions corresponding to certain sharp features in the surface models. Moreover, comparing the mesh qualities in terms of the angles of the triangles, it can be observed that the SDEM remeshing results contain more skinny and irregular triangle elements when compared to the other two methods. This can be explained by the fact that the SDEM method maps all surfaces onto the unit sphere regardless of their geometries, which may lead to large geometric distortions. Also, while the EDEM method allows us to map the surfaces onto a prescribed ellipsoid, we may need to carefully choose a suitable ellipsoidal shape in order to reduce both the area and angle distortions of the mappings. By contrast, the EDEQ method optimizes both the ellipsoidal shape and the mapping onto it, giving the highest flexibility among the three methods. Therefore, it gives the most uniformly distributed and regular triangle elements in the remeshing results.

For a more quantitative analysis, we evaluate the surface remeshing performance of different approaches by considering the shape and size variations of the triangle elements in the remeshed surfaces. Ideally, the triangle elements on the remeshed surface should be as uniform as possible in terms of their size (i.e. the area of the triangles) and shape (i.e. the regularity of the triangles). To assess the size variation of the remeshed surface, we define the size variation measure $\delta_{\text{size}}$ as follows:

$$\delta_{\text{size}}=\log\left(\frac{A_{\text{max}}}{A_{\text{min}}}\right),$$ (43)

where $A_{\text{max}}$ is the maximum triangle area and $A_{\text{min}}$ is the minimum triangle area in the remeshed surface. It is easy to see that $\delta_{\text{size}}\geq 0$ for any remeshed surface, and the equality holds if and only if $A_{\text{max}}=A_{\text{min}}$, i.e. all triangle elements of the remeshed surface are equal in size. This shows that $\delta_{\text{size}}$ can effectively capture the size variation in the remeshed surface.

To assess the shape variation, we first define the face regularity $R_{i}$ for each face $T_{i}$ of the remeshed surface as follows:

$$R_{i}=\sum^{3}_{j=1}\left|\frac{e^{i}_{j}}{e^{i}_{1}+e^{i}_{2}+e^{i}_{3}}-\frac{1}{3}\right|,$$ (44)

where $e^{i}_{1},e^{i}_{2},e^{i}_{3}$ represent the length of the three edges of $T_{i}$. Note that $R_{i}=0$ if and only if $T_{i}$ is an equilateral triangle. Then, we define the shape variation measure $\delta_{\text{shape}}$ of the entire remeshed surface as

$$\delta_{\text{shape}}=\underset{i}{\text{mean}}(R_{i}).$$ (45)

It is easy to see that $\delta_{\text{shape}}\geq 0$ for any remeshed surface, and the equality holds if and only if all triangles in the remeshed surface are equilateral. Therefore, $\delta_{\text{shape}}$ can effectively quantify the overall shape variation of the remeshed surface.

Table 7 records the values of $\delta_{\text{size}}$ and $\delta_{\text{shape}}$ for all remeshing results shown in Figs. 11, 12, 13, and 14. For the two conformal approaches (SCM and ECM), it can be observed that the values of $\delta_{\text{shape}}$ are small but the values of $\delta_{\text{size}}$ are extremely large, which means the remeshing results obtained by those two methods are highly non-uniform in size. By contrast, the values of $\delta_{\text{size}}$ for the three density-equalizing approaches (SDEM, EDEM, and EDEQ) are small. This demonstrates the advantage of utilizing density-equalizing maps for surface remeshing. Furthermore, by considering both $\delta_{\text{size}}$ and $\delta_{\text{shape}}$ for the three density-equalizing methods, it can be observed that the EDEQ method gives the best overall remeshing performance among the three methods, which also matches our observations in Fig. 15.

Table 7: Comparison between the surface remeshing results obtained by SCM [9], ECM [30], SDEM [43], EDEM, and EDEQ. For each method, we record the shape variation measure $\delta_{\text{shape}}$ and the size variation measure $\delta_{\text{size}}$ of the remeshed surface.

Surface	Method	Shape variation $\delta_{\text{shape}}$	Size variation $\delta_{\text{size}}$
Pig (Fig. 11)	SCM	0.0499	4.5846
	ECM	0.0505	4.7123
	SDEM	0.1562	3.3251
	EDEM	0.1493	2.4706
	EDEQ	0.1478	1.4218
Bear (Fig. 12)	SCM	0.0535	3.7168
	ECM	0.0526	3.9181
	SDEM	0.1296	1.9349
	EDEM	0.1282	1.4401
	EDEQ	0.1263	1.3977
Hippocampus (Fig. 13)	SCM	0.0546	7.7204
	ECM	0.0549	7.0796
	SDEM	0.2593	4.1060
	EDEM	0.2428	1.3482
	EDEQ	0.1568	1.2508
Lion-Vase (Fig. 14)	SCM	0.0587	5.5245
	ECM	0.0576	5.4088
	SDEM	0.1744	2.9125
	EDEM	0.1583	1.7410
	EDEQ	0.1249	1.7081

6 Discussion

In this work, we have proposed a novel method for computing bijective ellipsoidal density-equalizing maps (EDEM) for genus-0 closed surfaces. Using this approach, we can efficiently obtain parameterizations with different desired mapping effects, including ellipsoidal area-preserving parameterizations and ellipsoidal parameterizations with controlled area change. Also, by considering a combined energy involving both a density-equalizing term and a quasi-conformal term, we can achieve ellipsoidal density-equalizing quasi-conformal maps (EDEQ) balancing the angle and area distortions. Experimental results on various ellipsoidal surfaces and genus-0 closed surface models have demonstrated the effectiveness of our proposed methods. Our proposed methods can be easily applied to surface remeshing of genus-0 closed surfaces, with the remeshing quality significantly improved when compared to other prior approaches.

We remark that while our proposed methods have considered both the density-equalizing property and the quasi-conformality, landmark-matching constraints have not been incorporated into our formulations. In the future, we plan to further extend our methods to facilitate landmark-matching surface parameterization and surface registration. Our models may also be combined with curvature-based methods to achieve curvature-based surface registrations with controlled area change for surfaces with extreme geometries.