1 Introduction

A parameterization of a surface is a diffeomorphism that maps the surface onto a planar domain. A canonical coordinate system on the surface can be induced by parameterization. Surface parameterizations have been widely studied and applied in various tasks of computer graphics, such as surface registration, remeshing, morphing, alignment, and texture mapping.

An ideal parameterization usually preserves as much as possible of the geometric information such as length, angle, and area. A length-preserving map is also called an isometry. Classical parameterization methods that minimize the length distortion include the as-rigid-as-possible surface parameterization [33, 47], the most isometric parametrization [11, 23], and the isometric distortion minimization [37]. In addition, an area-preserving map is also called an authalic map or equiareal map. Classical computational methods for equiareal parameterizations include the stretch-minimizing method [39, 49], the Lie advection method [58], the discrete optimal mass transportation [12, 16, 45, 46, 57], the density-equalizing mapping [9], and the stretch energy minimization [50, 52, 53, 55].

In contrast, an angle-preserving map is called a conformal map. Classical computational methods for conformal parameterizations include the circle packing method [27, 44], the linear Laplace–Beltrami equation [6, 22], the angle-based flattening method [41, 43], the least squares conformal map [32], the holomorphic one-form method [18, 19, 29], the spectral conformal parameterization [26, 35], the discrete Ricci flow [28, 56], and the boundary first flattening [40]. In particular, for a simply connected open Riemann surface, the uniformization theorem guarantees that the shape of the domain can be a unit disk. The computation of the disk conformal mapping is a fundamental issue that has been widely studied. Several efficient algorithms have been developed in recent years, including the heat diffusion [17, 25], the fast disk conformal parameterization (FDCP) [10], the linear disk conformal parameterization (LDCP) [8], the conformal energy minimization (CEM) [30, 54], and the parallelizable global conformal parameterization [7].

Furthermore, some parameterization methods aim to reach a trade-off between minimizing the angle and area distortions, such as the polar decomposition method [36] and the balanced energy minimization [51].

The development of the algorithms and applications of surface parameterization can be found in several classic survey papers [5, 14, 15, 21, 24, 42].

Contribution In this paper, we develop a constructive disk conformal parameterization (CDCP) algorithm via the composition of conformal operations, namely, the spherical harmonic mapping of the doubly covered surface, a conformal mapping by the geodesic algorithm [34], and a Möbius transformation. The contributions and advantages of the proposed CDCP algorithm are threefold:

  1. (i)

    The bijectivity of the mapping computed by our CDCP algorithm is guaranteed under the Delaunay assumption of the triangular mesh.

  2. (ii)

    We prove that the number of unknowns of the associated linear system for the spherical harmonic mapping of a doubly covered surface can be reduced by roughly half so that the computational time cost is significantly reduced.

  3. (iii)

    According to numerical experiments, both the efficiency and accuracy of our CDCP algorithm are competitive or even improved compared to the other two state-of-the-art algorithms.

Notations In this paper, we use the following notations.

  • Bold letters, for instance, \({{\textbf {f}}}\), denote (complex) vectors.

  • Capital letters, for instance, L, denote matrices.

  • Typewriter letters, for instance, \({\mathtt {I}}\), \({\mathtt {B}}\), denote ordered sets of indices.

  • \(n_{{\mathtt {I}}}\) denotes the number of elements of the set \({\mathtt {I}}\).

  • \({{\textbf {f}}}_{i}\) denotes the ith entry of the vector \({{\textbf {f}}}\).

  • denotes the subvector of \({{\textbf {f}}}\) composed of \({{\textbf {f}}}_i\), for \(i\in {\mathtt {I}}\).

  • \(A_{i,j}\) denotes the (i, j)th entry of the matrix A.

  • \(A_{{\mathtt {I}},{\mathtt {J}}}\) denotes the submatrix of A composed of \(A_{i,j}\), for \(i\in {\mathtt {I}}\) and \(j\in {\mathtt {J}}\).

  • \({\mathbb {R}}\) denotes the set of real numbers.

  • \({\mathbb {S}}^n:=\{ {{\textbf {x}}}\in {\mathbb {R}}^{n+1} \,|\, \Vert {{\textbf {x}}}\Vert =1 \}\) denotes the n-sphere in \({\mathbb {R}}^{n+1}\).

  • \({\mathbb {C}}:=\{ u+\text {i}v \mid u,v\in {\mathbb {R}}\}\) denotes the complex plane.

  • \({\mathbb {H}}:=\{ u+\text {i}v \mid u\in {\mathbb {R}}, v>0 \}\) denotes the upper-half complex plane.

  • \({\text {Re}}(z)\) denotes the real part of the complex number z.

  • \({\text {Im}}(z)\) denotes the imaginary part of the complex number z.

  • \(\left[ {v}_0, \ldots , {v}_m\right] \) denotes the m-simplex of the vertices \({v}_0, \ldots , {v}_m\).

  • \(|[{v}_0, \ldots , {v}_m]|\) denotes the volume of the m-simplex \(\left[ v_0, \ldots , v_m\right] \).

  • \(\text {i}\) denotes the imaginary unit \(\sqrt{-1}\).

  • \(I_n\) denotes the identity matrix of size \(n\times n\).

  • \({{\textbf {0}}}\) denotes zero vectors and matrices of appropriate sizes.

Outline The remaining part of this paper is organized as follows. In Sect. 2, we briefly review the most related works to our CDCP algorithm. In Sect. 3, we introduce the discretization of surfaces and mappings by triangular meshes and piecewise affine mappings, respectively. In Sect. 4, we describe the proposed CDCP algorithm in detail. In Sect. 5, we prove the number of unknowns of the linear system involved in the spherical harmonic mapping of a doubly covered surface can be reduced by roughly half. In Sect. 6, we present numerical results of the proposed CDCP algorithm and compare them to the other two state-of-the-art algorithms in terms of efficiency and accuracy. Concluding remarks are given in Sect. 7.

2 Previous Works

Now, we review the most related works to the proposed CDCP algorithm, namely, the double covering technique for computing disk conformal parameterization in Sect. 2.1, and the zipper algorithm for computing conformal mappings between Jordan regions in Sect. 2.2.

2.1 Double Covering Technique

The double covering technique for computing disk conformal parameterization was first proposed by Gu and Yau [19]. Given a simply connected open surface \({\mathcal {M}}\) and let \({\mathcal {M}}^-\) be an identical copy of \({\mathcal {M}}\) with the reverse orientation. The double covering of \({\mathcal {M}}\) is a closed surface \(\widetilde{{\mathcal {M}}}\) composed of \({\mathcal {M}}\) and \({\mathcal {M}}^-\) with their boundaries \(\partial {\mathcal {M}}\) and \(\partial {\mathcal {M}}^-\) identified. Based on the fact that harmonic mappings between closed surfaces are conformal, the heat diffusion can be applied to conformally map the closed surface \(\widetilde{{\mathcal {M}}}\) to \({\mathbb {S}}^2\). Finally, the stereographic projection and the Möbius transformation are used to obtain the desired disk conformal parameterization, as illustrated in Fig. 1.

The double covering technique nicely avoids the boundary decision problem on the computation of the disk conformal parameterization; however, the main drawback of this approach is that the number of triangular faces and interior vertices of the mesh becomes doubled so that the computational time cost would significantly increase. To remedy this issue, we applied the structure of the associated discrete Laplace–Beltrami operator of the doubly covered triangular mesh and reduce roughly half the number of unknowns in the linear system, which is introduced in Sect. 5. In addition, to ensure the boundary vertices are precisely mapped on \({\mathbb {S}}^1\), we apply the geodesic algorithm [34] to do so, which is introduced in Sect. 2.2.

Fig. 1
figure 1

An illustration of the algorithm [19] for the computation of the disk conformal mapping using the heat diffusion with the double covering technique

Full size image

2.2 Geodesic Algorithm

Fig. 2
figure 2

An illustration of the geodesic algorithm, where \(\zeta _j=\varphi _{j-2}(z_j)\), for \(j=3, \ldots , n\)

Full size image

Suppose \(z_1, \ldots , z_n \in {\mathbb {C}}\) counterclockwise forms a piecewise affine curve that enclose a region \(\varOmega \subset {\mathbb {C}}\). The fundamental issue is to find a conformal map that maps \(\varOmega \) to \({\mathbb {H}}\). In the 1980s, an elementary algorithm for computing such conformal mappings was proposed independently by Kühnau and Marshall, but the convergence was not proved yet. Later, the convergence of the geodesic algorithm, a variant of the zipper algorithm, is proven by Marshall and Rohde [34]. They use a sequence of conformal mappings to map a specified region \(\varOmega \) to the upper-half plane \({\mathbb {H}}\). In particular, each boundary point of \(\varOmega \) is mapped on the real axis. Now, we introduce in detail the geodesic algorithm [34].

First, we map the interval \([z_1, z_2]\) onto the imaginary axis by

$$\begin{aligned} \varphi _1(z) = {\left\{ \begin{array}{ll} 1 &{} \text {if }z=\infty ,\\ \infty &{} \text {if }z=z_1,\\ \sqrt{\frac{z-z_2}{z-z_1}} &{} \text {otherwise}. \end{array}\right. } \end{aligned}$$
(1)

The mapping \(\varphi _1\) is composed of \((\varphi _1)^2: z \mapsto \frac{z-z_2}{z-z_1}\) and the square root mapping with the branch \(\sqrt{-1}=\text {i}\). The mapping \((\varphi _1)^2\) is a Möbius transformation that maps the interval \([z_1,z_2]\) onto the nonpositive real axis, and the square root mapping is holomorphic with a nonvanishing first derivative. As a result, \(\varphi _1\) is conformal. For convenience, we denote \(\zeta _3 = \varphi _1(z_3)\).

Remark 1

For \(z\in \varOmega \), we write \(z-z_1=r_1e^{\text {i}\theta _1}\) and \(z-z_2=r_2e^{\text {i}\theta _2}\) as the polar forms, where \(\theta _1, \theta _2\in [-\pi ,\pi ]\). As a result, \(\varphi _1(z)\) can be written as

$$\begin{aligned} \varphi _1(z)=\sqrt{\frac{z-z_2}{z-z_1}}=\sqrt{\frac{r_2}{r_1}}e^{\text {i}\frac{\theta _2-\theta _1}{2}}. \end{aligned}$$

Since \(z_1, z_2\) are on the same line, the argument \(\frac{\theta _1-\theta _2}{2}\in [0, \frac{\pi }{2}]\). Namely, \(\varphi _1(\varOmega )\) is in the first quadrant, as illustrated in Fig. 2. In particular, \({\text {Re}}\zeta _3>0\).

On the image of \(\varphi _1\), the interval \(\varphi _1([z_2, z_3]) = [0,\zeta _3]\). Next, we map the interval \([0,\zeta _3]\) onto the imaginary axis by

$$\begin{aligned} \varphi _2(z) = {\left\{ \begin{array}{ll} \text {i} &{} \text {if }z=\varphi _1(z_2),\\ 0 &{} \text {if }z=\varphi _1(z_3),\\ \sqrt{{\mathscr {M}}_{\zeta _3}(z)^2-1} &{} \text {otherwise}, \end{array}\right. } \end{aligned}$$
(2)

where \({\mathscr {M}}_c:\overline{{\mathbb {C}}}\rightarrow \overline{{\mathbb {C}}}\) is the Möbius transformation with respect to \(c\in {\mathbb {C}}\) as

$$\begin{aligned} {\mathscr {M}}_c(z)= {\left\{ \begin{array}{ll} \frac{{\text {Re}}(c)}{\text {i}{\text {Im}}(c)} &{} \text {if }z=\infty ,\\ \infty &{} \text {if }z=\frac{-|c|^2}{\text {i}{\text {Im}}(c)},\\ \frac{{\text {Re}}(c)z}{|c|^2+\text {i}{\text {Im}}(c)z} &{} \text {otherwise}. \end{array}\right. } \end{aligned}$$

Remark 2

The intuition of \(\varphi _2\) is stated as follows. Note that \({\mathscr {M}}_{\zeta _3}\) maps 0, \(\zeta _3\), \(\frac{-|\zeta _3|^2}{\text {i}{\text {Im}}(\zeta _3)}\) to 0, 1, \(\infty \), respectively. From a geometrical perspective, a semicircular arc from 0 to \(\frac{-|\zeta _3|^2}{\text {i}{\text {Im}}(\zeta _3)}\) is mapped to the non-negative real line \({\mathbb {R}}_{\ge 0}\) from 0 to \(\infty \), as illustrated in Fig. 2. Then, the interval \({\mathscr {M}}_{\zeta _3}([0, \zeta _3])\) is mapped to the imaginary axis by \(z\mapsto \sqrt{z^2-1}\).

Next, suppose we have mapped the point \(z_{j}\) to 0 and would like to map the interval \([0, \zeta _{j+1}]\) onto the imaginary axis, where

$$\begin{aligned} \zeta _{j+1} := \varPhi _{j-1}(z_{j+1}):= \varphi _{j-1} \circ \cdots \circ \varphi _{1}(z_{j+1}). \end{aligned}$$

The conformal map

$$\begin{aligned} \varphi _j(z) = {\left\{ \begin{array}{ll} \text {i} &{} \text {if }z=\varPhi _{j-1}(z_j),\\ 0 &{} \text {if }z=\varPhi _{j-1}(z_{j+1}),\\ \sqrt{{\mathscr {M}}_{\zeta _{j+1}}(z)^2-1} &{} \text {otherwise}, \end{array}\right. } \end{aligned}$$
(3)

\(j=3, \ldots , n-1\), can be applied to do so.

Finally, we define

$$\begin{aligned} \varphi _{n}(z) = {\left\{ \begin{array}{ll} \varPhi _{n-1}(z_1)^2&{} \text {if }z=\infty , \\ \infty &{} \text {if }z=\varPhi _{n-1}(z_1),\\ \left( \frac{\varPhi _{n-1}(z_1)z}{\varPhi _{n-1}(z_1)-z}\right) ^2 &{} \text {otherwise}. \end{array}\right. } \end{aligned}$$
(4)

Then, the composition mapping

$$\begin{aligned} \varPhi _n:=\varphi _{n} \circ \cdots \circ \varphi _{1} \end{aligned}$$
(5)

is the desired conformal map that maps \(\varOmega \) onto \({\mathbb {H}}\), where \(\varphi _1, \ldots , \varphi _n\) are defined in (1), (2), (3) and (4), respectively. These mappings have been applied to the conformal welding procedure for the computation of conformal parameterizations [7].

Ultimately, the disk conformal parameterization can be achieved by the Möbius transformation \(\psi _{\mathbb {D}}:{\mathbb {H}}\rightarrow {\mathbb {D}}\) defined as

$$\begin{aligned} \psi _{\mathbb {D}}(z) = {\left\{ \begin{array}{ll} 1 &{} \text {if }z=\infty ,\\ \infty &{} \text {if }z=\overline{\varPhi _n(0)},\\ \frac{z-\varPhi _n(0)}{z-\overline{\varPhi _n(0)}} &{} \text {otherwise}. \end{array}\right. } \end{aligned}$$
(6)

As a result, \(\psi _{\mathbb {D}}\circ \varPhi _n:\varOmega \rightarrow {\mathbb {D}}\) is a conformal mapping that maps \(\varOmega \) onto \({\mathbb {D}}\) and preserves the origin.

3 Discrete Surfaces, Mappings, and Operators

In this section, we introduce the discretization of surfaces, mappings, and the Laplace–Beltrami operator in Sects. 3.1, 3.2 and 3.3, respectively.

3.1 Discrete Surfaces

The discrete surface considered in this paper is the triangular mesh of genus-zero with a boundary. A triangular mesh \({\mathcal {M}}\) is consisted of a set of n vertices

$$\begin{aligned} {\mathcal {V}}({\mathcal {M}}) = \left\{ {v}_s \equiv \left( {v}_s^1, {v}_s^2, {v}_s^3 \right) ^\top \in {\mathbb {R}}^3 \right\} _{s=1}^n \end{aligned}$$

with coordinates in \({\mathbb {R}}^3\), and a set of oriented triangular faces

$$\begin{aligned} {\mathcal {F}}({\mathcal {M}}) = \left\{ \left[ {v}_i, {v}_j, {v}_k \right] \subset {\mathbb {R}}^3 \text { for some vertices }\left\{ {v}_i, {v}_j, {v}_k \right\} \subset {\mathcal {V}}({\mathcal {M}})\right\} , \end{aligned}$$

where the bracket \(\left[ {v}_i, {v}_j, {v}_k \right] \) denotes the convex hull or the 2-simplex of affinely independent set of points \(\{ {v}_i, {v}_j, {v}_k \}\). In addition, the set of non-oriented edges of \({\mathcal {M}}\) is denoted as

$$\begin{aligned} {\mathcal {E}}({\mathcal {M}}) = \left\{ \left[ {v}_i, {v}_j\right] \,|\, \left[ {v}_i, {v}_j, {v}_k \right] \in {\mathcal {F}}({\mathcal {M}}) \text { for some }{v}_k\in {\mathcal {V}}({\mathcal {M}}) \right\} . \end{aligned}$$

The union \({\mathcal {V}}({\mathcal {M}})\cup {\mathcal {E}}({\mathcal {M}})\cup {\mathcal {F}}({\mathcal {M}})\) forms a homogeneous simplicial 2-complex.

3.2 Discrete Mappings

A discrete mapping \(f:{\mathcal {M}}\rightarrow {\mathbb {C}}\) on \({\mathcal {M}}\) is a piecewise affine map defined as

$$\begin{aligned} f(v) = {\left\{ \begin{array}{ll} {{\textbf {f}}}_i &{} \text{ if } v \text{ is } \text{ the } i\text{ th } \text{ vertex } v_i\in {\mathcal {V}}({\mathcal {M}}),\\ \lambda _i \, {{\textbf {f}}}_i + \lambda _j \, {{\textbf {f}}}_j + \lambda _k \, {{\textbf {f}}}_k &{} \text{ if } v\in [v_i, v_j, v_k]\in {\mathcal {F}}({\mathcal {M}}), \end{array}\right. } \end{aligned}$$

with coefficients \(\lambda _i = \frac{|[{v}, {v}_j, {v}_k]|}{|[{v}_i, {v}_j, {v}_k]|}\), \(\lambda _j = \frac{|[{v}_i, {v}, {v}_k]|}{|[{v}_i, {v}_j, {v}_k]|}\) and \(\lambda _k = \frac{|[{v}_i, {v}_j, {v}]|}{|[{v}_i, {v}_j, {v}_k]|}\) being the barycentric coordinates of v in the triangle \([{v}_i, {v}_j, {v}_k]\). Here the notation \(|[{v}_0, \ldots , {v}_m]|\) denotes the volume of the m-simplex \(\left[ v_0, \ldots , v_m\right] \). In particular, \(|[v_i, v_j, v_k]|\) and \(|[v_i, v_j]|\) denote the area and length of the triangle \([v_i, v_j, v_k]\) and interval \([v_i, v_j]\), respectively. It is common to denote the space of piecewise affine maps on \({\mathcal {M}}\) as \(C^{PL}({\mathcal {M}})\) [17]. A map \(f\in C^{PL}({\mathcal {M}})\) can be represented as a complex-valued vector

$$\begin{aligned} {{\textbf {f}}} \equiv ({{\textbf {f}}}_1, \ldots , {{\textbf {f}}}_n)^\top = (f(v_1), \ldots , f(v_n))^\top \in {\mathbb {C}}^{n}. \end{aligned}$$

A discrete conformal mapping on a triangular mesh \({\mathcal {M}}\) is a piecewise affine mapping \(f:{\mathcal {M}}\rightarrow {\mathbb {D}}\) that preserves the angles, that is,

$$\begin{aligned} {\left\{ \begin{array}{ll} \angle (u, v, w) = \angle (f(u), f(v), f(w)),\\ \angle (w, u, v) = \angle (f(w), f(u), f(v)),\\ \angle (v, w, u) = \angle (f(v), f(w), f(u)), \end{array}\right. } \end{aligned}$$

for every \([u,v,w]\in {\mathcal {F}}({\mathcal {M}})\), where

$$\begin{aligned} \angle (u, v, w) = \arccos \left( \frac{|[u,v]|^2 + |[v,w]|^2 - |[w,u]|^2}{2|[u,v]|\,|[v,w]|}\right) . \end{aligned}$$

However, in general, such mapping does not exist because the sum of neighboring angles of a vertex is not necessarily \(2\pi \). Due to this fact, a discrete conformal mapping usually refers to a piecewise affine mapping that minimizes the conformality distortion in some sense. One of the most commonly considered conformality distortion is the mean and standard deviation (SD) of the angular distortions defined as

$$\begin{aligned} {\mathcal {D}}_{\mathrm{angle}}(f, \tau , v) = |\angle (u, v, w) - \angle (f(u), f(v), f(w))|, \end{aligned}$$
(7)

where \(\tau =[u,v,w]\in {\mathcal {F}}({\mathcal {M}})\) and \(v\in {\mathcal {V}}({\mathcal {M}})\cap \tau \). A conformal map has both the mean and SD of the angular distortions being 0.

3.3 Discrete Laplace–Beltrami Operator

The discrete Laplace–Beltrami operator on \(C^{PL}({\mathcal {M}})\) is a matrix \({L}_{\mathcal {M}}\) defined as

$$\begin{aligned}{}[{L}_{\mathcal {M}}]_{i,j} = {\left\{ \begin{array}{ll} -\frac{1}{2}\left( \cot \theta _{i,j}+\cot \theta _{j,i}\right) &{}\text{ if } [{v}_i,{v}_j]\not \subset \partial {\mathcal {M}},\\ -\frac{1}{2}\cot \theta _{i,j} &{}\text{ if } [{v}_i,{v}_j]\subset \partial {\mathcal {M}},\\ -\sum _{k\ne i} [{L}_{\mathcal {M}}]_{i,k} &{}\text{ if } i = j, \\ 0 &{}\text{ otherwise } \end{array}\right. } \end{aligned}$$
(8)

in which \(\theta _{i,j}\) and \(\theta _{j,i}\) are the two angles opposite to the edge \([{v}_i,{v}_j]\) connecting vertices \({v}_i\) and \({v}_j\) on \({\mathcal {M}}\), as illustrated in Fig. 3.

Fig. 3
figure 3

An illustration for the cotangent weights

Full size image

The formula of \({L}_{\mathcal {M}}\) is derived from the discrete Dirichlet energy

$$\begin{aligned} E_{\mathcal {M}}({{\textbf {f}}}) = \frac{1}{2} \, {{\textbf {f}}}^* {L}_{\mathcal {M}} \, {{\textbf {f}}}, \end{aligned}$$
(9)

where \({{\textbf {f}}}^*\) denotes the conjugate transpose of \({{\textbf {f}}}\) and \({L}_{\mathcal {M}}\) is the matrix in (8). The discrete Dirichlet energy (9) is derived by substituting \(f\in C^{PL}({\mathcal {M}})\) into the Dirichlet energy \({\mathcal {E}}_{\mathcal {M}}(f) = \int _{{\mathcal {M}}} \Vert \nabla f\Vert ^2 \text {d} v_{\mathcal {M}}\) and omitting those singularities of measure zero at vertices and edges of the triangular mesh \({\mathcal {M}}\). A critical point of (9) is called a discrete harmonic mapping that satisfies \(\nabla E_{\mathcal {M}}({{\textbf {f}}}) = {{\textbf {0}}}\), i.e., \({L}_{\mathcal {M}} \, {{\textbf {f}}} = {{\textbf {0}}}\). The derivation in detail can be found in [20, Section 10.1]. Under a given boundary condition , the interior mapping is uniquely determined by the linear system

(10)

Remark 3

In the finite element method, the Laplace–Beltrami operator for triangular mesh is discretized as \(\triangle _{\mathcal {M}} \approx M^{-1} L_{\mathcal {M}}\), where \({L}_{\mathcal {M}}\) is the cotangent Laplacian defined as in (8), and M is a mass matrix for the area normalization usually being a diagonal matrix with entries being areas of Voronoi or barycentric cells of vertices [38]. As minimizing the discrete Dirichlet energy (9) is the main concern, \(\triangle _{\mathcal {M}}\) is discretized by \(L_{\mathcal {M}}\), i.e., \(M=I_n\). In particular, for the boundary value problem of the Laplace–Beltrami equation (10), the solution is independent of the choices of vertex mass.

4 Numerical Methods

In this section, we introduce our new algorithm for the computation of the disk conformal mapping. Suppose \({\mathcal {M}}\) is a simply connected open triangular mesh of n vertices with \({\mathtt {I}}\) and \({\mathtt {B}}\) being the sets of indices of interior and boundary vertices, respectively. Let \(\widetilde{{\mathcal {M}}}\) be the doubly covered surface of \({\mathcal {M}}\), which is a genus-zero closed surface. First, we find the central triangular face of the mesh \({\mathcal {M}}\) by applying the heat diffusion, which is introduced in Sect. 4.1. Then, we compute a spherical harmonic map \({\widetilde{h}}:\widetilde{{\mathcal {M}}}\rightarrow \overline{{\mathbb {C}}}\), which is introduced in Sect. 4.2. Then, the restriction map \(h:={\widetilde{h}}|_{{\mathcal {M}}}\) is a conformal map that maps \({\mathcal {M}}\) onto a region \(\varOmega =h({\mathcal {M}})\subset \overline{{\mathbb {C}}}\). Finally, we apply the geodesic algorithm [34] to map \(\varOmega \) onto \({\mathbb {D}}\subset {\mathbb {C}}\) and obtain the disk conformal map, which is introduced in Sect. 4.3.

4.1 Central Triangular Face

First, we find the surface center of \({\mathcal {M}}\) by the heat diffusion

$$\begin{aligned} \frac{\partial f(v,t)}{\partial t} = \triangle _{\mathcal {M}} f(v,t) \end{aligned}$$
(11)

with the initial function being

$$\begin{aligned} f(v,0) = {\left\{ \begin{array}{ll} 1 &{} \text{ if } v\in \partial {\mathcal {M}}, \\ 0 &{} \text{ otherwise, } \end{array}\right. } \end{aligned}$$

and a suitable terminal time \(t_1\). By applying the implicit Euler method, the heat diffusion (11) is discretized as a sequence of linear systems

$$\begin{aligned} \left( I_n + \delta t L_{\mathcal {M}} \right) {{\textbf {f}}}^{(k+1)} = {{\textbf {f}}}^{(k)}, \end{aligned}$$
(12)

where \(\delta t\) is an appropriate step size of the time step, \({{\textbf {f}}}^{(0)} = {{\textbf {1}}}_{\partial {\mathcal {M}}}\), and

$$\begin{aligned} {{\textbf {f}}}^{(k)} = \left( f(v_1, k\delta t), \ldots , f(v_n,k\delta t) \right) ^\top , ~ \text { for }k=0, 1, \ldots , \frac{t_1}{\delta t}. \end{aligned}$$

In particular, by choosing both the terminal time \(t_1\) and the time step size \(\delta t\) in (12) to be \(10^4\), the iteration (12) becomes the linear system

$$\begin{aligned} \left( I_n + 10^4 L_{\mathcal {M}} \right) {{\textbf {f}}} = {{\textbf {1}}}_{\partial {\mathcal {M}}}, \end{aligned}$$
(13)

where the matrix \(L_{\mathcal {M}}\) is defined in (8) and the right-hand-side vector \({{\textbf {1}}}_{\partial {\mathcal {M}}}\) is the indicator vector defined as

$$\begin{aligned} ({{\textbf {1}}}_{\partial {\mathcal {M}}})_i = {\left\{ \begin{array}{ll} 1 &{} \text{ if } v_i\in \partial {\mathcal {M}}, \\ 0 &{} \text{ otherwise. } \end{array}\right. } \end{aligned}$$

The central vertex is then defined as

$$\begin{aligned} v_c = \mathop {{\text {argmin}}}\limits _{v_i\in {\mathcal {V}}({\mathcal {M}})} {{\textbf {f}}}_i, \end{aligned}$$
(14)

and the central triangular face is defined as

$$\begin{aligned} \tau _c = \mathop {{\text {argmin}}}\limits _{[v_i, v_j, v_k]\in {\mathcal {F}}({\mathcal {M}})} ({{\textbf {f}}}_i+{{\textbf {f}}}_j+{{\textbf {f}}}_k). \end{aligned}$$
(15)

Remark 4

A suitable terminal time \(t_1\) should be sufficiently large so that each entry of \({{\textbf {f}}}\) is far from zero. In addition, \(t_1\) should also be far from infinity so that entries of \({{\textbf {f}}}\) would be far from constant. Based on these two facts, the value \(10^4\) is a nice option. The result would be the same or similar if we change the value of \(t_1\) into, e.g., \(10^3\) or \(10^5\), as shown in Fig. 4. On the other hand, the smaller time step size \(\delta t\) is used, the more linear systems need to be solved to reach the terminal time. In that case, the sparse Cholesky decomposition of the symmetric matrix \(I_n + \delta t L_{{\mathcal {M}}}\) should be first computed so that each linear system of the form (12) would cost merely two triangular systems, which can be computed efficiently by forward and backward substitutions.

Fig. 4
figure 4

The surface centers of the mesh models Ear, Foot, Chinese Lion and Stanford Bunny computed by (14) with \(t_1\) being \(10^3\), \(10^4\) and \(10^5\), respectively

Full size image

4.2 Spherical Harmonic Map

The spherical harmonic map can be computed by solving the Laplace–Beltrami equation [6, 22]

$$\begin{aligned} \triangle _{\widetilde{{\mathcal {M}}}} {\widetilde{h}} = \left( \frac{\partial }{\partial u} - \text {i} \frac{\partial }{\partial v} \right) \delta _p, \end{aligned}$$
(16)

where \(\delta _p\) is the Dirac delta function with p being the center of \({\mathcal {M}}\). In practice, we let \(\mathtt {P} = \{i,j,k\}\) be the index set given by the central face \(\tau _c = [v_i,v_j,v_k]\) in (15). Under the piecewise affine approximation, (16) is represented as a linear system

$$\begin{aligned} L_{\widetilde{{\mathcal {M}}}} \widetilde{{{\textbf {h}}}} = \widetilde{{{\textbf {r}}}}, \end{aligned}$$
(17)

where the discrete Laplace–Beltrami operator with respect to \(\widetilde{{\mathcal {M}}}\) and the right-hand-side vector are represented as

respectively, and \({{\textbf {r}}}\) is the vector with only 3 nonzero entries

Since harmonic maps between genus-zero closed surfaces are conformal [48], the restriction is a discrete conformal mapping that maps \({\mathcal {M}}\) to \(\varOmega =h({\mathcal {M}})\subset {\overline{{\mathbb {C}}}}\).

Remark 5

Under the Delaunay assumption of the triangular mesh, the bijectivity of the mapping computed by (17) is guaranteed by applying [13, Theorem 6.7], because the 3-point boundary condition would form a triangle—a convex shape.

Remark 6

The angular distortion of the map \(\widetilde{{{\textbf {h}}}}\) would be relatively large near the constrained points . Due to this observation, we choose the triangular face at the center of the surface \({\mathcal {M}}\) as the index set \(\mathtt {P}\) so that the angular distortion would be relatively small at . In our approach, the heat diffusion (11) with a sufficient large terminal time \(t_1\) is applied to find such a triangular face due to the fact that it would approximate \(\triangle _{\mathcal {M}} f = 0\), which is consistent with (16).

Noting that the matrix \(L_{\widetilde{{\mathcal {M}}}}\) in the linear system (17) is of size \(2n_{{\mathtt {I}}}+n_{{\mathtt {B}}}\) and almost half of blocks are replica of the other blocks, which can be reduced into a linear system of size \(n_{{\mathtt {I}}}+n_{{\mathtt {B}}}\) as

$$\begin{aligned} L_{{\mathcal {M}}} {{\textbf {h}}} = {{\textbf {r}}}. \end{aligned}$$
(18)

The solutions of the linear systems (17) and (18) are equivalent up to a translation and a scaling by \(\frac{1}{2}\). A proof for this statement is given in Corollary 1 in Sect. 5.

In addition, the matrix \({L}_{{\mathcal {M}}}\) in 18 is singular with the kernel being

$$\begin{aligned} \ker ({L}_{{\mathcal {M}}}) = \text {span}\{(1, \ldots , 1)^\top \}. \end{aligned}$$

In practice, we fix , then the system 18 becomes nonsingular with a unique solution.

4.3 Disk Conformal Map

Due to the symmetry of the doubly covered triangular mesh, the boundary map is close to a circle. To ensure the boundary vertices are mapped to \({\mathbb {S}}^1\), we apply the geodesic algorithm [34] introduced in Sect. 2.2 to do so. Namely, the boundary map is updated by

where \(\varPhi _n\) and \(\psi _{\mathbb {D}}\) are given in (5) and (6), respectively. Then, the interior map is computed by solving the linear system

To further reduce the angular distortions of the resulting mapping f, the parameterization can be improved by the Möbius transformation \(\psi _\omega :{\mathbb {D}}\rightarrow {\mathbb {D}}\) defined as

$$\begin{aligned} \psi _\omega (z) = \frac{z-\omega }{1-{\overline{\omega }}z} \end{aligned}$$
(19)

in which \(\omega \) is chosen to be

$$\begin{aligned} \omega = \mathop {{\text {argmin}}}\limits _{\omega \in {\mathbb {D}}} \mathop {{\text {mean}}}\limits _{\begin{array}{c} \tau \in {\mathcal {F}}({\mathcal {M}})\\ v\in {\mathcal {V}}({\mathcal {M}})\cap \tau \end{array}} {\mathcal {D}}_{\mathrm{angle}}(\psi _\omega \circ f, \tau , v), \end{aligned}$$
(20)

where \({\mathcal {D}}_{\mathrm{angle}}\) is the angular distortion defined in (7). The optimization problem (20) can be efficiently solved by the Nelder–Mead simplex method [31].

The algorithm for the computation of disk conformal mappings is summarized in Algorithm 1.

figure a

5 Equivalence of Boundary Maps

Now, we prove the equivalence of the boundary maps computed by (17) and (18) up to a translation and a scaling by \(\frac{1}{2}\).

Theorem 1

Given two linear systems of block forms

$$\begin{aligned} \begin{bmatrix} A &{} C \\ C^\top &{} 2 B &{} C^\top \\ &{} C &{} A \end{bmatrix} \begin{bmatrix} \widetilde{{{\textbf {f}}}}\\ \widetilde{{{\textbf {g}}}}\\ \widetilde{{{\textbf {h}}}} \end{bmatrix} = \begin{bmatrix} {{\textbf {b}}}\\ {{\textbf {0}}}\\ {{\textbf {0}}} \end{bmatrix} \end{aligned}$$
(21)

and

$$\begin{aligned} \begin{bmatrix} A &{} C \\ C^\top &{} B \end{bmatrix} \begin{bmatrix} {{\textbf {f}}}\\ {{\textbf {g}}} \end{bmatrix} = \begin{bmatrix} {{\textbf {b}}}\\ {{\textbf {0}}} \end{bmatrix}, \end{aligned}$$
(22)

where A is nonsingular. Suppose \((\widetilde{{{\textbf {f}}}}^\top , \widetilde{{{\textbf {g}}}}^\top , \widetilde{{{\textbf {h}}}}^\top )^\top \) is a solution to (21). Then, \(\widetilde{{{\textbf {g}}}}^\top = \frac{1}{2}{{\textbf {g}}}^\top \) for some solution \(({{\textbf {f}}}^\top , {{\textbf {g}}}^\top )^\top \) to (22).

Proof

Suppose \((\widetilde{{{\textbf {f}}}}^\top , \widetilde{{{\textbf {g}}}}^\top , \widetilde{{{\textbf {h}}}}^\top )^\top \) is a solution to (21). Then, \(\widetilde{{{\textbf {f}}}}\), \(\widetilde{{{\textbf {g}}}}\) and \(\widetilde{{{\textbf {h}}}}\) satisfies

figure b

From (23a) and the assumption that A is nonsingular,

$$\begin{aligned} \widetilde{{{\textbf {f}}}} = A^{-1} ({{\textbf {b}}} - C\widetilde{{{\textbf {g}}}}) = A^{-1} {{\textbf {b}}} - A^{-1}C\widetilde{{{\textbf {g}}}}. \end{aligned}$$
(24)

Similarly, from (23c), we have

$$\begin{aligned} \widetilde{{{\textbf {h}}}} = -A^{-1} C\widetilde{{{\textbf {g}}}}. \end{aligned}$$
(25)

Substituting (24) and (25) into (23b), we obtain

$$\begin{aligned} C^\top (A^{-1} {{\textbf {b}}} - 2A^{-1}C\widetilde{{{\textbf {g}}}}) + 2B\widetilde{{{\textbf {g}}}} = {{\textbf {0}}}. \end{aligned}$$

That is,

$$\begin{aligned} (B - C^\top A^{-1} C)(2\widetilde{{{\textbf {g}}}}) = -C^\top A^{-1}{{\textbf {b}}}. \end{aligned}$$
(26)

On the other hand, the linear system (22) can be written as

figure c

From (27a) and the assumption that A is nonsingular,

$$\begin{aligned} {{\textbf {f}}} = A^{-1} ({{\textbf {b}}} - C{{\textbf {g}}}) = A^{-1} {{\textbf {b}}} - A^{-1}C{{\textbf {g}}}. \end{aligned}$$
(28)

Substituting (28) into (27b), we obtain

$$\begin{aligned} B {{\textbf {g}}} = -C^\top A^{-1} {{\textbf {b}}}+C^\top A^{-1}C{{\textbf {g}}}. \end{aligned}$$

That is,

$$\begin{aligned} (B - C^\top A^{-1} C){{\textbf {g}}} = -C^\top A^{-1}{{\textbf {b}}}. \end{aligned}$$
(29)

Therefore, from (26) and (29), there is a solution \(({{\textbf {f}}}^\top , {{\textbf {g}}}^\top )^\top \) to (22) that satisfies \({{\textbf {g}}}=2\widetilde{{{\textbf {g}}}}\) as desired. \(\square \)

By substituting \(L_{{\mathtt {I,I}}}\), \(L_{{\mathtt {I,B}}}\), and \(L_{{\mathtt {B,B}}}\) into A, B, and C in Theorem 1, respectively, we have the desired consequence as follows.

Corollary 1

The boundary maps and computed by the linear systems (17) and (18), respectively, satisfy , where c is a constant.

Corollary 1 indicates that we can solve the reduced linear system (18) of size \(n_{{\mathtt {I}}}+n_{{\mathtt {B}}}-1\) instead of solving the linear system (17) of size \(2n_{{\mathtt {I}}}+n_{{\mathtt {B}}}-1\). As a result, the computational time cost for computing the boundary map is significantly reduced.

Remark 7

Corollary 1 could also explain the reason that the initial mapping of the CEM algorithm in [30, 54] is already very close to conformal so that the number of iteration steps for the convergence is usually less than 20.

6 Numerical Experiments

In this section, we demonstrate the numerical results of the proposed CDCP algorithm for the computation of the disk conformal parameterizations of simply connected open triangular meshes. In addition, we compare the accuracy and the efficiency of the CDCP algorithm with other state-of-the-art algorithms, FDCP and LDCP. Furthermore, in Sect. 6.3, we discuss the utility of the Möbius transformation performed in Lines 9, 10 of Algorithm 1.

Some of the triangular mesh models are obtained from TurboSquid [3], AIM@SHAPE shape repository [1], the Stanford 3D scanning repository [2], and Sketchfab [4]. All numerical experiments are performed in MATLAB on a personal laptop with a 2.30 GHz CPU and 16 GB RAM.

In Fig. 5, we present the mesh models Max Planck, Buddha and Nefertiti Statue and their disk conformal parameterizations computed by the proposed CDCP algorithm. The color represents the normal map of the surface. In Fig. 6, we present the histogram of the angular distortions defined in (7). We see that most of the angular distortions are less than 2 degrees, which is quite satisfactory.

Fig. 5
figure 5

The mesh models Max Planck, Buddha and Nefertiti Statue and their disk conformal parameterizations computed by the CDCP algorithm

Full size image
Fig. 6
figure 6

The histogram of angular distortions of the disk conformal parameterizations computed by the CDCP algorithm for a Face Lin, b Knit Cap Man, c Bimba Statue, d Max Planck, e Buddha, and f Nefertiti Statue

Full size image

In the following, we compare the accuracy and efficiency of the proposed CDCP algorithm with other state-of-the-art algorithms in Sect. 6.1. In addition, we discuss the choices of the central triangular face and the utility if Möbius transformations in Sects. 6.2 and 6.3, respectively.

6.1 Comparison with State-of-the-Art Algorithms

In Table 1, we present the mean and SD of the angular distortion, defined in (7), of the parameterizations of benchmark mesh models computed by the FDCP [10], LDCP [8], and the proposed CDCP algorithm. The data in Table 1 is also presented in Fig. 7 for a better visualization. From Fig. 7, we see that the mean and SD of the angular distortions of the proposed CDCP algorithm is similar or the best compared to the other two state-of-the-art algorithms.

Table 1 The mean and SD of angular distortions of the disk conformal maps computed by the FDCP, LDCP, and CDCP algorithms
Full size table
Fig. 7
figure 7

The a mean and b SD of angular distortions (in degree) of the conformal parameterization by FDCP, LDCP, and the CDCP algorithms

Full size image

In Table 2, we present the computational time cost of disk conformal parameterizations of benchmark mesh models by FDCP, LDCP, and the CDCP algorithm. In Fig. 8, we further present the relationship between the number of vertices and the computational time cost. From Fig. 8, we see that the proposed CDCP algorithm has the smallest computational time cost compared to the other two state-of-the-art algorithms. In particular, for a mesh model with roughly 1 million vertices, it would cost merely 33 seconds to compute its disk conformal parameterization, which is satisfactory.

It is worth noting that all the demonstrated mesh models do not satisfy the Delaunay assumption. Nevertheless, all the resulting parameterizations computed by the CDCP algorithm are bijective.

In summary, the proposed CDCP algorithm has similar accuracy and significantly improved efficiency compared to the other two state-of-the-art algorithms. The files of mesh models, mapping results, and MATLAB executables of the CDCP algorithm are available at http://tiny.cc/CDCP.

Table 2 The computational time costs (sec.) of the conformal parameterizations computed by the FDCP, LDCP, and CDCP algorithms
Full size table
Fig. 8
figure 8

a The computational time costs for the conformal parameterizations by FDCP, LDCP, and the CDCP algorithms. b The relationship between the numbers of vertices and the computational time costs for the conformal parameterizations by FDCP, LDCP, and the CDCP algorithms

Full size image

6.2 Choices of the Central Triangular Face

As in Remark 6, the central triangular face is better chosen to be as far from the boundary as possible. In our approach, the heat diffusion (11) is applied to find such a central triangular face. Another heuristic method by directly finding the triangular face closest to the centroid of the surface could also deal with the task. We now compare these two methods for finding a central triangular face.

From Fig. 9, we observe that the vertex closest to the centroid of the surface could be distant from the actual surface center. It is crucial to know whether the choices of the central triangular faces would affect the accuracy of the resulting mapping in terms of angular distortions.

Fig. 9
figure 9

The surface centers of the mesh models Ear, Foot, Chinese Lion and Stanford Bunny computed by the directly finding the vertex closest to the centroid and the heat diffusion with \(t_1=10^4\), respectively

Full size image

In Fig. 10, we present the mean and SD of the angular distortion of the mappings computed by the CDCP algorithm with \(\tau _c\) being the face closest to the surface centroid and the surface center (obtained by the heat diffusion), respectively, without the Möbius correction. We observe that the choices of the central faces do not significantly affect the angular distortions for most of the demonstrated mesh models, however, for the very coarse mesh model Ear, the CDCP with \(\tau _c\) being the face closest to the surface center is significantly better. Thus, we adopted the surface center obtained by the heat diffusion, although it would cost solving a linear system.

Fig. 10
figure 10

The a mean and b SD of the angular distortions (in degree) with \(\tau _c\) being the face closest to the surface centroid and the surface center (obtained by the heat diffusion), respectively

Full size image

6.3 Utility of Möbius Transformations

The mappings in (5) and (6) are conformal, however, the transformation performed in Lines 6 and 7 in Algorithm 1 are piecewise affine approximations of (5) and (6), respectively. As a result, a piecewise affine approximation of the Möbius transformation in (19) can be applied to further improve the parameterization in terms of angular distortions in (7).

It is crucial to know how much the Möbius transformation in (19) can improve the angular distortions and how much time the optimization problem 20 costs.

In Fig. 11, we present the mean and SD of the angular distortions of the conformal parameterization by the CDCP algorithm with and without Möbius transformation (19). We observe that the Möbius transformation (19) could merely slightly improve the angular distortions for most of the benchmark mesh models, however, for the models Face Ho and Face Lin, the Möbius transformation (19) significantly improves the angular distortions.

On the other hand, in Fig. 12, we present the computational time costs for the conformal parameterizations by the CDCP algorithm with and without Möbius transformation (19). We observe that solving the optimization problem 20 would cost even more time than that of the other procedures. In particular, without Möbius transformation (19), it would costs about 14 secs. to compute the conformal parameterization for the mesh model Nefertiti Statue by the CDCP algorithm, but solving the optimization problem 20 for the Möbius transformation (19) would costs about 30 secs.

Fig. 11
figure 11

The a mean and b SD of the angular distortions (in degree) of the conformal parameterization by the CDCP algorithm with and without Möbius transformation (19)

Full size image
Fig. 12
figure 12

a The computational time costs for the conformal parameterizations by the CDCP algorithm with and without Möbius transformation (19). b The relationship between the numbers of vertices and the computational time costs for the conformal parameterizations by the CDCP algorithm with and without Möbius transformation (19)

Full size image

In summary, when the angular distortion is the main concern, the Möbius transformation (19) can be applied to slightly improve the parameterization. Otherwise, Lines 9–10 of Algorithm 1 can be omitted so that the efficiency is highly improved.

7 Concluding Remarks

In this paper, we have developed an efficient constructive algorithm called CDCP, combined with the spherical harmonic mapping, geodesic algorithm, and Möbius transformation, for the computation of disk conformal parameterizations of simply connected open surfaces. Although the double covering is applied, we have proved that the number of unknowns can be kept the same by applying the structure of the associated linear system so that the computational time cost can be significantly reduced. Numerical results indicate that the proposed CDCP algorithm has similar accuracy and the best efficiency compared to the other two state-of-the-art algorithms. In particular, it would cost less than 15 sec. on a personal laptop for the computation of a disk conformal parameterization of a simply connected open triangular mesh of roughly 1 million vertices. Such promising results show the potential of the CDCP algorithm on real-time applications of disk conformal parameterizations.