Subdivision connectivity remeshing via Teichmüller extremal map

Curvilinear surfaces in 3D Euclidean spaces are commonly represented by triangular meshes. The structure of the triangulation is important, since it affects the accuracy and efficiency of the numerical computation on the mesh. Remeshing refers to the process of transforming an unstructured mesh to one with desirable structures, such as the subdivision connectivity. This is commonly achieved by parameterizing the surface onto a simple parameter domain, on which a structured mesh is built. The 2D structured mesh is then projected onto the surface via the parameterization. Two major tasks are involved. Firstly, an effective algorithm for parameterizing, usually conformally, surface meshes is necessary. However, for a highly irregular mesh with skinny triangles, computing a folding-free conformal parameterization is difficult. The second task is to build a structured mesh on the parameter domain that is adaptive to the area distortion of the parameterization while maintaining good shapes of triangles. This paper presents an algorithm to remesh a highly irregular mesh to a structured one with subdivision connectivity and good triangle quality. We propose an effective algorithm to obtain a conformal parameterization of a highly irregular mesh, using quasi-conformal Teichmuller theories. Conformality distortion of an initial parameterization is adjusted by a quasi-conformal map, resulting in a folding-free conformal parameterization. Next, we propose an algorithm to obtain a regular mesh with subdivision connectivity and good triangle quality on the conformal parameter domain, which is adaptive to the area distortion, through the landmark-matching Teichmuller map. A remeshed surface can then be obtained through the parameterization. Experiments have been carried out to remesh surface meshes representing real 3D geometric objects using the proposed algorithm. Results show the efficacy of the algorithm to optimize the regularity of an irregular triangulation.


(Communicated by Sung Ha Kang)
Abstract. Curvilinear surfaces in 3D Euclidean spaces are commonly represented by triangular meshes. The structure of the triangulation is important, since it affects the accuracy and efficiency of the numerical computation on the mesh. Remeshing refers to the process of transforming an unstructured mesh to one with desirable structures, such as the subdivision connectivity. This is commonly achieved by parameterizing the surface onto a simple parameter domain, on which a structured mesh is built. The 2D structured mesh is then projected onto the surface via the parameterization. Two major tasks are involved. Firstly, an effective algorithm for parameterizing, usually conformally, surface meshes is necessary. However, for a highly irregular mesh with skinny triangles, computing a folding-free conformal parameterization is difficult. The second task is to build a structured mesh on the parameter domain that is adaptive to the area distortion of the parameterization while maintaining good shapes of triangles. This paper presents an algorithm to remesh a highly irregular mesh to a structured one with subdivision connectivity and good triangle quality. We propose an effective algorithm to obtain a conformal parameterization of a highly irregular mesh, using quasi-conformal Teichmüller theories. Conformality distortion of an initial parameterization is adjusted by a quasi-conformal map, resulting in a folding-free conformal parameterization. Next, we propose an algorithm to obtain a regular mesh with subdivision connectivity and good triangle quality on the conformal parameter domain, which is adaptive to the area distortion, through the landmark-matching Teichmüller map. A remeshed surface can then be obtained through the parameterization. Experiments have been carried out to remesh surface meshes representing real 3D geometric objects using the proposed algorithm. Results show the efficacy of the algorithm to optimize the regularity of an irregular triangulation.
1. Introduction. With the advance of image acquisition technologies, 3D geometric objects from the real world can now be effectively captured. The captured geometric data are commonly represented by triangular meshes. With these digital geometric data, systematic shape analysis of the geometric objects and numerical computations on them can be carried out. Applications can be found in different areas, such as medical imaging, computer graphics and computer visions. In order to ensure the accuracy and efficiency of the shape analysis and numerical computations, a high quality triangular mesh that represents the geometric object is crucial. The quality of a mesh often depends on 1. the structure of the connectivity; 2. the shapes of triangular faces and 3. geometric error from the original surface. For example, a surface mesh with lots of skinny triangles causes numerical errors in scientific computing. It also affects the accuracies of computing geometric quantities, such as curvatures, which further affects the accuracies of shape analysis results. Furthermore, it is commonly desirable to have a mesh with a good structure of the connectivity, such as the subdivision connectivity. A subdivision connectivity mesh automatically provides a sequence of meshes with different refinement levels capturing different levels of details. This facilitates multi-resolution algorithms, such as wavelets [24,36], progressive transmission [21] and multiresolution editing [46].
In practical situations, triangulation from 3D raw geometric data are usually unstructured. It is therefore necessary to develop algorithms to improve the quality of the triangulation from an irregular mesh. Such a process is called remeshing. Remeshing is commonly done by surface parameterization. The surface is firstly parameterized onto a simple parameter domain, for instance, a 2D rectangle. A regular triangulation on the parameter domain is then built, which can be projected onto the original surface via the parameterization. To preserve the regularity of the triangulation after the projection, the parameterization has to preserve the local geometry as much as possible. A popular choice is the conformal parameterization. For surface meshes with good triangle quality, conformal parameterizations can be obtained using conventional parameterization algorithms. However, for highly irregular meshes with skinny triangles, conformal parameterizations are generally difficult to compute (see Table 1). For example, skinny triangles may cause the discrete Laplace-Beltrami operator, which is involved in most conformal parameterization algorithms, to become singular. As a result, most algorithms fail to obtain the parameterization of those meshes (such as the human face and Lion meshes as shown in Figure 1). Even if the parameterization can be obtained, a folding-free mapping cannot be guaranteed. Flips of triangles under the parameterization can often be observed. This causes troubles for remeshing, since the projected mesh onto the original surface through a non-bijective parameterization may contains flips of triangles. Developing an effective algorithm to conformally and bijectively parameterize highly irregular surface meshes is necessary. With the conformal parameterization, the second major task is to build a structured mesh on the conformal parameter domain. The 2D structured mesh can then be projected onto the surface via the parameterization to produce a remeshed surface. However, one important issue is that the conformal parameterization often introduces area distortions. When an arbitrary regular mesh on the conformal parameter domain is projected onto the surface, the resolution at some regions may get coarsen. This results in a serious loss in the geometric details (see Figure 6). An adaptive triangulation on the parameter domain according to the area distortion under the parameterization is thus necessary.
In this paper, we present an algorithm to remesh a highly irregular mesh to a structured one with subdivision connectivity and good triangle quality. We first develop an effective algorithm to obtain a conformal parameterization of a highly irregular mesh, using quasi-conformal Teichmüller theories. The irregular mesh is firstly embedded in R 2 using the Tutte's embedding. This initial parameterization is expected to introduce conformality distortions. To fix the conformality distortions, the initial parameterization is composed with a quasi-conformal map from the initial parameter domain to the unit disk. A folding-free conformal parameterization can be obtained. Next, we propose a method to build a regular mesh with subdivision connectivity and good triangle quality on the conformal parameter domain, which is adaptive to the area distortions under the conformal parameterization. The basic idea is to transform a regular mesh with a constrained Teichmüller map (T-Map) that matches landmarks consistently. The landmarks are chosen such that their distribution follows the area distortion of the conformal parameterization. A T-Map is the "most conformal" map subject to the landmark constraints. It transforms the regular mesh to another regular mesh, which follows the area density of the conformal parameter domain. With this triangulation, a remeshed surface can be constructed through the projection by the parameterization, which preserves the original geometry well. To test the effectiveness of the proposed algorithm, experiments have been carried out on surface meshes extracted from real 3D geometric objects. Results show the efficacy of our proposed algorithm to parameterize and optimize the regularity of an irregular triangular mesh.
In short, the contributions of this paper are two-folded. Firstly, we propose an effective algorithm to robustly compute conformal parameterizations of irregular meshes, even for meshes with skinny triangles, using quasi-conformal theories. Secondly, we propose an efficient algorithm to obtain an adaptive regular mesh on the conformal parameter domain based on the area distortion of the parameterization, using the landmark-matching Teichmüller maps (T-Maps). The algorithm allows us to produced surface remesh of the original surface with subdivision connectivity and good triangle quality.
The rest of the paper is organized as follows. In section 2, some related works are presented. In section 3, the basic mathematical background is described. The proposed algorithm for surface remeshing is explained in details in section 5. Experimental results are reported in section 6. The paper is concluded in section 7.
2. Previous work. Surface remeshing is an important pre-processing in computer graphics and scientific computing on surfaces. Many remeshing algorithms have been proposed by various research groups. Existing remeshing algorithms can mainly be divided into two categories, namely, 1. the parameterization approach and 2. the explicit approach. The parameterization approach parameterizes the surface mesh onto a simple parameter domain, on which a structured mesh is built. Surface remeshing is then achieved by projecting the structured mesh onto the surface [33,17,6,1,2] For example, Praun et al. [33] proposed to remesh genus-0 closed surfaces by mapping it to the spherical domain conformally. Hormann et al. [17] proposed to remesh triangulated topologically disk-like surfaces by parameterizing the surface onto a planar domain, using the most isometric parameterization strategy (MIPS). Eck et al. [6] presented a remeshing algorithm based on partitioning the meshes into several triangular regions followed by parameterizing each regions using harmonic maps. Gu et al. [11] proposed to remesh an arbitrary surface onto a completely regular structure, through cutting the mesh along a network of edge paths, and parametrize the resulting single chart onto a square. Alliez et al. [1] proposed an isotropic remeshing algorithm with locally uniform edges, through building a weighted centroidal Voronoi tessellation in a conformal parameter space, where the specified density function is used for weighting. Later, Alliez et al. [2] proposed an algorithm that uses curvature directions to drive the remeshing process. The algorithm can produce meshes ranging from isotropic to anisotropic, from coarse to dense, and from uniform to curvature adapted.
Another category of remeshing algorithms is the explicit mesh modification approach, in which vertices are progressively adjusted until it matches some specified properties [32,41,3,7,8,15,16,34,37]. For example, Peyrè et al. [32] proposed a fast algorithm for the remeshing of a surface with a uniform or adaptive distribution, which is based on iteratively choosing the farthest point according to a weighted distance on the surface. Yan et al. [41] proposed a fast isotropic remeshing method, based on an efficient algorithm for the Restricted Voronoi Diagram (RVD) for computing the centroidal Voronoi tessellation (CVT). Later, centroidal Voronoi tessellation in universal covering spaces of manifold surfaces was proposed in [35], which can be used for remeshing surfaces of arbitrary topologies. Chen et al. [3] proposed a parameterization-free remeshing algorithm by progressively optimizing an initial resampled mesh through alternatively recovering the Delaunay mesh and moving each vertex to the centroid of its 1-ring neighborhood.
In this work, surface remeshing is carried out based on the conformal parameterization. Different algorithms for conformal parameterization has been recently proposed. For instance, Haker et al. proposed a finite element approximation of conformal parameterization in [14]. They linearized the Laplace-Beltrami operator and solved the sparse linear system for conformally parameterizing brain surfaces. In [23], Lévy et al. proposed a parameterization method by approximating the Cauchy-Riemann equations using the least-squares method. Desbrun et al. [4] introduced the intrinsic parameterizations which minimize the distortion of different intrinsic measures of the surface patches. In [18], Hurdal and Stephenson proposed a discrete mapping approach for spherical conformal parameterization which uses circle packing to produce "flattened" images of cortical surfaces on the sphere, the Euclidean plane, and the hyperbolic plane. Gu and Yau [12,13] introduced a nonlinear algorithm for spherical conformal parameterization. They performed the optimization in the tangent space of the sphere by gradient descent. The computation is more stable and accurate. In [22], Lai et al. reported an approach to obtain a folding-free global conformal mapping. Curvature flow methods, which deform the Riemannian metric conformally to the uniformization metric, have also recently been studied to obtain conformal parameterizations of surface with arbitrary topologies [19,42,43].
Quasi-conformal theories will also be applied in our proposed remeshing algorithm. The computation of quasi-conformal mappings has been recently studied to obtain smooth 1-1 correspondences with bounded conformality distortion [25,26,44,28]. For example, Lui et al. [26] proposed to compute quasi-conformal registration between hippocampal surfaces based on the holomorphic Beltrami flow method, which matches geometric quantities (such as curvatures) and minimizes the conformality distortion [25]. Wei et al. [44] proposed the Quasi-Yamabe flow method to compute quasi-conformal mapping for high-genus surfaces. Quasi-conformal mapping that matches landmarks consistently has also been proposed [25,45]. In [25], the authors proposed to compute the brain landmark-matching registration, which minimizes L 2 norm of the Beltrami coefficients. Wei et al. [45] also proposed to compute quasi-conformal mappings for feature matching face registration. The Beltrami coefficient associated to a landmark points matching parameterization is approximated. However, either exact landmark matching or the bijectivity of the mapping cannot be guaranteed, especially when very large deformations occur. Later, the extremal QC mapping, which minimizes the conformality distortion has been proposed. Zorin et al. [40] proposes a least square algorithm to compute mapping between connected domains with given Dirichlet condition defined on the whole boundaries. The extremal mapping is obtained by minimizing a least square Beltrami energy, which is non-convex. The algorithm can obtain an extremal mapping when initialization is carefully chosen. Recently, Lui et al. [27] proposed to compute the unique T-Map between simply-connected Riemann surfaces of finite type. The convergence of the algorithm has also been proven in [29]. The proposed algorithm was applied for landmark-based surface registration. Later, Ng et al. [30] extended the algorithm to compute the quasi-conformal extremal map between multiply-connected domains.
3. Basic mathematical background. In this work, conformal and quasi-conformal geometry theories will be applied. We describe briefly the basic mathematical theories mostly related to our proposed models in this section. For details, we refer the readers to [9].
Let Ω 1 and Ω 2 be simply-connected domains in C. A map f : Ω 1 → Ω 2 is conformal if it satisfies the Cauchy-Riemann equation: where f = u + iv. In term of metric, a conformal map f preserves the metric up to a multiplicative factor. Mathematically, where z and w are the coordinates on Ω 1 and Ω 2 respectively. λ is called the conformal factor. An immediate consequence of this property is that a conformal map preserves angles. Let S be a Riemann surface. The conformal parameterization of S can also be defined.
A natural generalization of the conformal map is the quasi-conformal map. A quasi-conformal map is an orientation preserving homeomorphism with bounded conformality distortions, in the sense that their first order approximations takes small circles to small ellipses of bounded eccentricity [9]. Mathematically, f : Ω 1 → Ω 2 is quasi-conformal provided that it satisfies the Beltrami equation: for some complex-valued function µ satisfying ||µ|| ∞ < 1. µ is called the Beltrami coefficient, which is a measure of non-conformality. µ f measures how far the map is deviated from a conformal map. µ ≡ 0 if and only if f is conformal. Infinitesimally, around a point p, f may be expressed with respect to its local parameter as follows: Obviously, f is not conformal if and only if µ(p) = 0. A quasi-conformal map f maps a small circle to a small ellipse. From µ(p), we can determine the directions of maximal magnification and shrinking and the amount of their distortions as well. Specifically, the angle of maximal magnification is arg(µ(p))/2 with magnifying factor 1+|µ(p)|; The angle of maximal shrinking is the orthogonal angle (arg(µ(p))− π)/2 with shrinking factor 1 − |µ(p)|. Thus, the Beltrami coefficient µ gives us all the information about the properties of the map. The maximal dilation of f is given by: Given a smooth BC µ : Ω 1 → C with µ ∞ < 1. There is always a diffeomorphism of C that satisfies the equation (3) [9]. 4. Problem setting. We will now describe our problem mathematically. Suppose S is a Riemann surface embedded in R 3 . In this work, S is assumed to be simplyconnected, although our proposed algorithm can be readily extended to surfaces of arbitrary topologies.
Let K S be the irregular triangular mesh approximating S. By convention, a mesh contains vertices, edges and faces. The intersection of two faces may either be a vertex, an edge or an empty set. An edge belongs to one face if it is a boundary edge, or else it belongs to two faces. We write Our goal is to build a mesh K S = {V K S , T K S } with good connectivity structure and good triangle quality, which is a good approximation of K S . More precisely, we construct a subdivision connectivity mesh with good triangle quality to approximate K S . A subdivision connectivity mesh is obtained from a sequential subdivision of a base mesh. In other words, we propose to construct a sequence of meshes . This mesh with the subdivision connectivity is automatically endowed with a sequence of meshes { K j S } L j=1 with different refinement details capturing different levels of details.
In addition to having the subdivision connectivity, we require K S to have good shapes of triangular faces. The triangle quality can be quantitatively measured by the radius-ratio defined as follow [31]: is the radius of the inscribed circle and R cir i is the radius of the circumscribed of the i-th triangle. The radius-ratio ρ has the range in 0 ≤ ρ ≤ 1. An equilateral triangle has ρ = 1 and a degenerated triangle has ρ = 0.
We can now formulate our problem as follows. Given an irregular mesh K S , we look for a sequence of meshes { K j S } L j=1 such that K S := K L S closely approximate K S and the radius-ratio ρ i is close to 1 for all T i ∈ T K S , 5. Remeshing via Teichmüller extremal map. In this section, we will describe our proposed remeshing algorithm via Teichmüller extremal map in details.

5.1.
Overview of the algorithm. The algorithm can be divided into three main steps.
1. Develop an effective algorithm for a folding-free conformal parameterization of a highly irregular surface mesh: the surface mesh with a highly irregular triangulation is conformally mapped onto a unit disk. 2. Develop an efficient algorithm to build an adaptive mesh on the parameter domain: a regular triangulation is built on the conformal parameter domain, which is adaptive to the area distortions under the conformal parameterization. This is done via Teichmüller extremal map. 3. Project the adaptive mesh to the surface: the adaptive regular triangulation on the conformal parameter domain is projected to the original surface.

5.2.
Robust conformal parameterization. Let K S be a triangular mesh approximating a surface S embedded in R 3 . Suppose K S is highly irregular, in the sense that both the triangle quality and connectivity condition are not satisfactory. Our first task is to obtain a conformal parameterization of K S onto a unit disk. We denote the conformal parameterization by Φ : Parameterizing an irregular mesh conformally is challenging. Existing algorithms are sensitive to the quality of the triangulation. For example, for meshes with skinny triangles (that is, one of the three inner angles of the triangle is close to π), the cotangent formula that approximates the Laplace-Beltrami operator usually leads to a singular matrix. As a result, algorithms, which rely on the cotangent formula, fail to obtain the parameterization of highly irregular meshes. Even if the parameterization can be obtained, foldings or flips of triangular faces can often be found. This causes troubles for remeshing, since it induces flips of triangles in the projected mesh onto the original surface through the non-bijective parameterization. Thus, we usually get stuck in the dilemma of a "chicken and egg" problem. To remesh an irregular mesh, we propose to look for a conformal parameterization onto the 2D domain and carry out the remeshing procedure on 2D. However, with a highly irregular triangulation, obtain the conformal parameterization is challenging.
We propose a robust algorithm to compute the conformal parameterization, which can deal with bad meshes. Our strategy is to embed K S to a regular foldingfree mesh K Ω in R 2 . The projection is not necessarily conformal. The initial embedding is then adjusted by computing a quasi-conformal map from K Ω onto another mesh K Ω , such that its composition with the initial embedding is conformal.
Denote the collection of all vertices of . A conformal parameterization of K S depends on its angle structure, which captures the information of the angles of each triangular faces. The angle structure depends on the lengths of each edges, which is often called the discrete metric defined as follows.
Definition 5.1. The discrete metric on a triangular mesh is a positive real-valued function l : E K S → R + defined on the collection of all edges that satisfies the following triangle inequality: where [u, w] denotes the edge joining the vertices u ∈ V and w ∈ V , and v i , v j and v k forms a triangular face.
The discrete metric gives information about the shapes of triangular faces. The robustness of most existing parameterization algorithms depends on the shape of the triangular faces, namely, whether the radius-ratio of faces are close to 1. In our case, the triangle quality of K S is bad and there are many skinny triangles. As a result, the cotangent formula that approximates the Laplace-Beltrami operator is usually singular. Parameterization algorithms that depend on the cotangent formula fail to embed the irregular mesh onto the 2D domain.
In order to find an initial embedding, we first assign K S with a fake discrete metric by setting l ≡ 1. In other words, we assume all edges have lengths equal to 1. Hence, all triangular faces are presumed to be equilateral. With this fake metric, we conformally parameterize K S onto a domain K Ω in R 2 using any existing conformal parameterization algorithms. Note that under this fake metric, the triangle quality is good enough for most parameterization algorithms to work. In this paper, we choose to use the Ricci flow algorithm with the fake discrete metric to parameterize K S [19]. The idea of Ricci flow is to conformally deform the surface metric g(t) = (g ij (t)) according to its induced Gaussian curvature K(t) to its uniformization metric through the heat flow equation on the metric: where g(0) is the initial fake discrete metric and g(∞) is the desired uniformization metric. From the uniformization metric, the surface can be conformally embedded onto a unit disk D. For details, please refer to [19]. Since all faces of K S are equilateral under the fake discrete metric, the faces are mapped to (almost) equilateral triangles in R 2 under the obtained initial parameterization (which is conformal under the fake metric). Thus, K Ω is a triangular mesh in R 2 with good triangle quality. For example, Figure 3(a) shows the parameterization of the foot mesh with the fake discrete metric. Note that the parameter domain, as shown in (B), is a triangular mesh of D with good triangle quality (the radius-ratio is close to 1 for all faces). This initial parameterization φ : K Ω → K S under the fake discrete metric is guaranteed to be bijective. As a result, we obtain a folding-free and non-degenerated 1 mesh K Ω in R 2 , which is homeomorphic to K S . 2 A folding-free and non-degenerated mesh K Ω (equivalently a homeomorphic φ) is crucial in our algorithm, which will be made clear in the following.
However, under the original discrete metric of K S , φ is not conformal, since angles are not preserved under the mapping. This is obvious since the shape of faces of K S are highly irregular, whereas the shape of faces of K Ω are (almost) equilateral. Our goal is to obtain a conformal parameterization of K S onto a simple domain in R 2 .
Our strategy is to compose our initial parameterization φ with a quasi-conformal map g that fixes the conformality distortion. More precisely, suppose the Beltrami coefficient of φ is given by µ φ . Let g : K Ω → D be a quasi-conformal map with Beltrami coefficient µ g . Then, the Beltrami coefficient of φ • g −1 is given by: Note that the above formula is valid when φ and g are orientation homeomorphism. In this case, both ||µ g || ∞ < 1 and ||µ φ || ∞ < 1. This is guaranteed since φ is a folding-free and homeomorphic. Hence, a homeomorphic φ is crucial.
Note that both φ and φ • g −1 are mapping from a 2D domain to the surface K S in R 3 . Their Beltrami coefficients can be computed by locally parameterizing K S . More specifically, let D ⊂ K Ω and ψ : φ(D) → Σ ⊂ R 2 be a local chart of φ(D). Then the Beltrami coefficient of φ on D is defined as the Beltrami coefficient of ψ • φ. We can show that this definition is independent of the choice of ψ and hence it is well-defined. The Beltrami coefficient of φ • g −1 can be similarly defined.
Motivated by the above theorem, we first compute the Beltrami coefficient µ φ of the initial parameterization φ. We then compute a quasi-conformal map g from K Ω to D. The composition map φ • g −1 is then a conformal parameterization of K S .
The associated quasi-conformal map g : K Ω → D can be obtained by solving the Beltrami's equation: (12) ∂g ∂z = µ ∂g ∂z In this work, we apply the Beltrami holomorphic flow (BHF) algorithm to solve Equation (12). As we shall notice, BHF is an iterative scheme to obtain the quasiconformal map, which does not rely on solving the Laplace-Beltrami equation discretized by the cotangent formula. Hence, the robustness of the algorithm is not affected by the singularity of the discrete Laplace-Beltrami operator. On the other hand, the Beltrami coefficient is iteratively adjusted, whose supreme norm is strictly less than 1. Thus, a folding-free homeomorphic conformal paramterization can be obtained.
We will explain the BHF algorithm briefly. For details, we refer the readers to [30]. The basic idea of the BHF algorithm is to iteratively deform an initial quasiconfomral map g 0 : K Ω → D to our desired quasi-conformal map g ∞ whose Beltrami coefficient is µ. The initial quasi-conformal map is chosen as the identity map in this work, whose Beltrami coefficeint is given by µ 0 ≡ 0. We proceed to adjust µ 0 and find its associated quasi-conformal map, which is closer to our target map g ∞ . As the Beltrami coefficient varies, its associated quasi-conformal map varies. Assume µ 0 varies by ω, and assume its associated quasiconformal map g µ0+ω varies by V. In other words, g µ0+ω (z) = g 0 (z) + V(z). The variation V can be found by solving the Beltrami equation. Let ν = µ 0 + ω. Define the differential operator A by: It follows from the Beltrami equation that In other words, finding V is equivalent to solving the partial differential equation (14) subject to the boundary condition that (15) (g 0 + V)| ∂KΩ = ∂D Note that the point-wise correspondence between ∂K Ω and ∂D is not required in (15). Equation (15) simply means V has to be tangential to ∂D.
To obtain g ∞ with Beltrami coefficient µ, we iteratively deform g 0 to our desired map. More specifically, we first set the initial map g 0 as the harmonic map whose boundary condition is given by the arc-length parameterization. We then flow g 0 to g 1 whose Beltrami coefficient is close to ν 1 := (1 − )µ 0 + µ ( > 0). This can be done by solving Equation (14) by putting ν = ν 1 with the boundary constraint to obtain V 0 . Note that V 0 on the boundaries is restricted to be tangential to ∂D. We get a new quasi-conformal map g 1 := g 0 + V 0 whose Beltrami coefficient is denoted by µ 1 . Suppose at the n th iteration, we have the quasi-conformal map g n with Beltrami coefficient µ n . We then flow g n to g n+1 whose Beltrami coefficient is close to ν n := (1 − )µ n + µ ( > 0). This is again done by solving Equation (14) by putting ν = ν n and f µ = f n with the boundary constraint to obtain V n . Set g n+1 := g n + V n whose Beltrami coefficient is denoted by µ n+1 . Note that in each step, can be chosen so that ||µ n+1 − µ|| ∞ is minimized. In practice, we choose = 1. A sequence of quasi-conformal maps {g n } ∞ n=1 is obtained, whose Beltrami coefficients converge to µ.
Once the quasi-conformal map g is computed, our desired conformal parameterization of K S can be obtained from the composition map Φ := φ • g −1 : Intuitively, the BHF algorithm iteratively adjusts the vertex positions of K Ω to a new mesh K Ω while keeping the connectivity, such that the map between K Ω and K S becomes conformal. Note that if ||µ|| ∞ > 1, the deforming mesh can become degenerated during the BHF process. As a result, the BHF algorithm will fail. Thus, the condition of ||µ|| ∞ < 1 is crucial. Again, this is guaranteed since the initial parameterization φ is homeomorphic.
Our algorithm to compute a conformal parameterization of an irregular triangular mesh can now be summarized as follows.

Algorithm 1 : (Robust Conformal parameterization)
Input : Irregular triangular mesh: K S . Output : Conformal parameterization Φ : K Ω → K S 1. Using Ricci flow, obtain an initial parameterization ϕ : K Ω → K S of K S with the fake discrete metric l ≡ 1. Compute the Beltrami coefficient µ φ of φ. 2. Using BHF, obtain a quasi-conformal map g : Our experimental results show that the proposed method can conformally parameterize highly irregular surface meshes with skinny triangular faces effectively. Details will be shown in Section 6.

Adaptive regular triangulation on D.
Once the irregular triangular surface mesh is parameterized conformally, a regular triangular mesh with subdivision connectivity can be constructed on the parameter domain. The structured triangulation on the parameter domain can then be projected onto the original surface mesh. Since the parameterization is conformal, the regularity of the triangulation is well-preserved. Hence, the original surface can be remeshed to a regular triangular mesh.
However, an important issue is the area distortion under the conformal parameterization. Although a conformal map preserves angles, it does not preserve area. Thus, a large number of triangular faces of the original surface mesh may possibly be squeezed to a small region on the conformal parameter domain (see Figure 3 and 4). If a uniform mesh is built on the parameter domain, some regions would have insufficient faces to approximate the original surface. As a result, geometric features might be lost after remeshing. A regular triangulation of the parameter domain that is adaptive to the area distortion is therefore necessary.
In this work, we propose an effective and efficient algorithm to obtain the adaptive regular triangulation of the parameter domain based on the centroidal Voronoi tessellation and the constrained Teichmüller mapping. We will first briefly describe the idea of centroidal Voronoi tessellation (CVT) and the constrained Teichmüller mapping (T-Map). Our proposed Teichmüller adaptive remeshing algorithm will then be explained in details.

Centroidal Voronoi tessellation (CVT).
In this work, centroidal Voronoi tessellation (CVT) will be used. We explain the concept of CVT briefly. We refer the readers to [5] for details.
CVT is a special kind of Voronoi tessellation. Given a set S and k elements z i in S (i = 1, 2, ..., k). The idea of Voronoi tessellation is to divide S into k subsets V 1 , V 2 , ..., V k such that: z i 's are called the generators and V i 's are called the Voronoi region with respect to z i .
Given a region V in R n and a density function ρ(w) defined on V . The mass centroid z * of V is given by: A centroidal Voronoi tessellation (CVT) is a Voronoi tessellation whose generators z i 's are the same as the mass centroids z * i 's of each Voronoi regions V i 's. The CVT generators are local uniformly distributed according to the density function ρ(w). This property was conjectured by Gersho [10] and was later proven for the 2D cases [38].
CVT is used in this work to find sample points, called landmarks, which are distributed according to the area distortion of the conformal parameterization.

Landmark-matching Teichmüller mapping (T-Map). Teichmüller map (T-Map) is used to project a regular mesh onto the conformal parameter domain.
We briefly describe the idea of the T-Map. For details, we refer the readers to [9,27,30] Suppose Ω 1 and Ω 2 are two simply-connected domains in R 2 . Let {p i ∈ Ω 1 } n i=1 and {q i ∈ Ω 2 } n i=1 be the corresponding landmarks on Ω 1 and Ω 2 respectively. One might be interested in looking for a bijective map f : Ω 1 → Ω 2 satisfying f (p i ) = q i (i = 1, 2, ..., n), which minimizes the conformality distortion. Mathematically, we look for such a landmark-matching bijective map f such that: (17) ||µ(f )|| ∞ ≤ ||µ(g)|| ∞ over all landmark-matching bijective map g : Ω 1 → Ω 2 , where µ(f ) and µ(g) are the Beltrami coefficients of f and g respectively. f is called the extremal map.
Intuitively, an extremal map is the "most conformal" map subject to the landmark constraints. The extremal map is closely related to the Teichmüller map (T-Map). A T-Map is a quasi-conformal map whose Beltrami coefficient is of the form: where 0 ≤ k < 1 is a positive constant, varphi : Ω 1 → C is a holomorphic function on Ω 1 . In other words, a T-Map has a uniform conformality distortion (|µ| = k).

Given any prescribed landmark correspondences ({p
, there exists a unique T-Map f : Ω 1 → Ω 2 satisfying f (p i ) = q i for i = 1, 2, ..., n, which is also the extremal map. Therefore, the landmark-matching T-Map is the "most conformal" bijective map satisfying the prescribed landmark constraints. This property is particularly important in our case, since a regular triangulation is required to be mapped to the parameter domain that preserves the angle structure.

5.3.3.
Teichmüller adaptive remeshing. In Section 5.2, we introduce an algorithm to compute a conformal parameterization Φ := φ • g −1 : K Ω → K S of the irregular mesh K S approximating a surface S. To remesh K S using Φ, we need to build a regular triangulation on D and project it to S for remeshing. Since Φ is conformal, the angle structure of the regular triangulation on D is well-preserved on S. Although Φ preserves angles, it does not preserve area. Hence, the regular mesh on Ω must be adaptive to the area distortion under Φ. Otherwise, geometric loss of S may be resulted.
We develop a method to build an adaptive mesh on D, which is called the Teichmüller adaptive remeshing. Our goal is to build a regular mesh on D according to the area distortion under Φ. In other words, more vertices are places on regions with larger squeezing.
We start with choosing k landmark points are uniformly distributed on K S . These uniform sampled landmark points give information about the density of the triangular faces of the conformal parameter domain K Ω . To get the landmark points, our strategy is to apply CVT to get k generators distributed according to the density function, which is given by the area distortion under φ. The density function ρ : K Ω → R + is defined as: where N (v i ) is the 1-ring neighbourhood faces of v i and A(T j ) is the area of the triangular face T j . The idea of adjusting the vertex positions according to the area distortion using CVT was firstly proposed in [32]. In [32], the positions of all vertices on the parameter domain were adjusted using CVT. Hence, the method was computationally expensive. In our current work, only a small amount of sample points (k is much less than the total number of vertices) are chosen. The sample points are considered as landmarks to give us information about the density of the conformal parameter domain. T-Map will then be used to map a regular mesh to the conformal parameter domain that follows the density based on the landmarks.
We initially take k points on K Ω and choose them as the initial generators of the CVT. We then use the Lloyd's method to compute the CVT with respect to the density function ρ, which is a simple iterative scheme as follows: and set z i = z * i for i = 1, 2, ..., k; 4. Go back to Step 2 until convergence.
After the CVT associated to the k generators {p i } k i=1 is computed, the Delaunay triangulation K CV T of {p i } k i=1 can be constructed. K CV T has a good connectivity, whose vertex degrees are usually of degree 6. However, the geometry of it may not be optimized. In some cases, sharp triangular faces may occur.
To further improve the geometry of the triangulation while preserving the structure of the connectivity, we compute the Tutte's embedding ϕ : CV T ⊂ D with good triangle quality while having the same connectivity as K CV T is then obtained. K tutte CV T is a regular triangulation with good geometry and connectivity. We can further refine K tutte CV T by sub-division to get a sequence of refined meshes: Since K tutte CV T has good triangle quality and connectivity, so are K tutte,i CV T for i = 1, 2, ..., L.
K tutte,i CV T gives a good triangular mesh approximating the parameter domain D. However, the vertices of K tutte,i CV T does not follow the density function ρ. We propose to deform K tutte,i CV T to another mesh K i regular of D, whose vertices follow ρ while preserving the regular angle structure of K tutte,i CV T as much as possible. Ideally, it is desirable to get a conformal map to deform K tutte,i CV T according to the density function ρ. But practically, such a conformal map does not exist. Thus, we look for an extremal Teichmüller map, which minimizes the conformality distortion. Mathematically, we search for a diffeomorphism T : D → D that minimizes ||µ(T )|| ∞ while the density of T ( K tutte CV T ) resembles ρ as good as possible. Note that the k generators give us information about the density of Φ. The above problem can be solved by finding a landmark-matching T-Map. More specifically, we search for an extremal T-Map T : D → D satisfying the landmark constraints: (21) T (Φ(p i )) = p i , i = 1, 2, ..., k.
T is then the optimized conformal map to transform K tutte,i CV T to another mesh K i regular := T ( K tutte,i CV T ) that follows the required density, while preserving the triangle quality of K tutte,i CV T as much as possible. Note again that the extremal T-Map T is computed on K tutte,i CV T , which is regular. The problem of numerical inaccuracies due to irregular triangulation can be avoided.
To compute T , the Quasi-conformal (QC) iteration can be used [27]: 1. Start with an initial landmark-matching map T 0 : K tutte,i CV T → D; 2. Suppose T n is obtained at the n th iteration. Compute ν n = µ(f n ). Let and L is the Laplace smoothing operator.
Set T n+1 = LBS LM (µ n ) where LBS LM (µ n ) find the closest landmarkmatching map T n+1 satisfying the Beltrami equation corresponding to µ n (in the L 2 -sense). More precisely, suppose µ n := ρ n + iτ n and T n+1 = u n+1 + iv n+1 . The Beltrami equation can be reformulated as follows.
The convergence of the above iterative scheme is proven in [29].
The proposed Teichmüller adaptive remeshing of the parameter domain can now be summarized as follows. Algorithm 2 : (Teichmüller adaptive remeshing) Input : Conformal parameterization Φ : Compute the Tutte's embedding of K CV T into D to obtain K tutte CV T . 4. Refine K tutte CV T by subdivision to get a sequence of subdivision meshes { K tutte,j CV T } L j=1 . 5. Compute the extremal T-Map T j : K tutte,j CV T → D such that T (Φ(p i )) = p i for i = 1, 2, ..., k. Compute the Teichmüller adaptive mesh K j regular = T j ( K tutte,j CV T ) for j = 1, ..., L.
In some situations, it may be desirable to keep the mesh topology (connectivity and vertex number) of the original surface mesh. In this case, we can choose k vertices {v i } k i=1 on K Ω as the sample points, which are distributed according to the area distortion using CVT. K Ω is then used as the base mesh. The Teichmüller mapT is then computed to deform K Ω such that it is adaptive to the density function. More precisely, we computeT : D → D such that:T (f −1 (v i )) = v i , where f : K Ω → K Ω is the quasi-conformal map as described in section 5.2. The algorithm can be summarized as follows.
Algorithm 3 : (Teichmüller adaptive remeshing with fixed mesh topology) Input : Conformal parameterization φ : K Ω ⊂ D → K S . Output : Teichmüller adaptive mesh K regular on D with fixed mesh topology 1. Compute the conformal factor ρ of φ 2. Using CVT, obtain k vertices f : K Ω → K Ω is the quasi-conformal map as described in section 5.2. Compute the Teichmüller adaptive mesh K regular =T (K Ω ).
The regular mesh K regular can then be projected to the surface S using ϕ to get a remeshed surface K remesh . Again, since the parameterization ϕ is conformal, the projected triangulation K remesh is a regular triangular mesh approximating the surface S.

Numerical implementation.
In this subsection, we briefly describe the numerical implementation of the proposed algorithms.
The computation of conformal parameterization of the irregular surface mesh requires approximating the Beltrami coefficient of the initial parameterization ϕ. The Beltrami coefficient µ is related to the Riemannian metric on S. Suppose S has a metric (23) ds 2 = Edx 2 + 2F dy 2 + Gdy 2 .
under the parameterization. Let dz = dx + idy and dz = dx − idy. The metric can be written as: Assume that ϕ(u, v) = (x(u, v), y(u, v), z(u, v)). Then: In the discrete setting, ϕ : K Ω → K S is a piecewise linear homeomorphism. The first order derivatives of x, y and z can be easily computed, which are constants on each triangular faces. E, F, G and hence µ can then be computed. Therefore, in the discrete setting, the Beltrami coefficient is piecewise constant on each triangular face.
Next, when applying CVT for the Teichmüller adaptive remeshing, the Voronoi tessellation of the k generators has to be computed. To ensure the efficiency, we use a quantized version of Voronoi tessellation. Let {p i } k i=1 be the chosen k generators in D. We compute the quantized Voronoi cell V q i (1 ≤ i ≤ k) as follows: Besides, the mass centroid z * has to be computed. Using the integral formula (16) to compute z * is time-consuming. We simplify the approximation of z * as follows: where |V q i | denotes the number of elements in V q i . We leave the details of the numerical implementations for Ricci flow, BHF and QC iterations but refer readers to related references. The numerical implementation for Ricci flow, BHF and QC iterations can be found in [19], [30] and [27] respectively. 6. Experimental results. We test the proposed remeshing algorithm on different surface meshes embedded in R 3 . All our experiments are carried out on a machine with the following configuration: AMD CPU 3.2 GHz and 16GB RAM. The algorithms are implemented using MATLAB. In this section, we report the remeshing results on six surface meshes, namely, 1. foot, 2. hand, 3. human face, 4. Venus, 5. mask and 6. lion. The original surface meshes of 1, 2, 3 and 4 are shown in Figure  2. The original meshes of 5 and 6 are shown in Figure 1. The surface meshes are visualized using a software called MeshLab.
Furthermore, to quantitatively measure the quality of a triangular mesh, we consider a triangle quality measurement, the radius-ratio [31]: (29) τ where R ins i is the radius of the inscribed circle and R cir i is the radius of the circumscribed circle of the i-th triangle. The radius-ratio τ i satisfies 0 ≤ τ i ≤ 1. An equilateral triangle has τ i = 1 and a degenerated triangle has τ i = 0.  6.1. Surface parameterization. We first test Algorithm 1 to parameterize irregular surface meshes. Figure 3 shows the parameterization results of the irregular surface mesh of a foot, whose original surface mesh is shown in Figure 2(A). The Tutte's embedding ϕ : K Ω → K S is firstly computed by endowing a fake metric (with all edge lengths equal to 1), which is shown in Figure 3(A). Note that the triangulation of the Tutte's embedding K Ω is regular and has no flipping triangles. From K Ω , we compute a quasi-conformal map g : K Ω → K Ω whose Beltrami coefficient is given by the conformality distortion of ϕ. The composition map φ = ϕ • g −1 is our desired conformal parameterization of K S , which is shown in Figure 3(B). Figure 4 shows another parameterization results of the irregular surface mesh of the Venus surface, whose original surface mesh is shown in Figure 2(D). The Tutte's embedding is shown in Figure 4(A) and the conformal parameterization is shown in Figure 4(B).  Table 1 gives the quantitative comparison of the parameterization results of the six surface meshes with other state-of-the-art approaches. The mean, standard deviation and minimum of the radius-ratio τ i for each mesh are recorded in the table, which show the triangle quality of the triangular meshes. The quantitative measurement is taken as the mean of the Beltrami coefficient of the parameterization. A parameterization is conformal if the value is 0. As shown in the table, our proposed method succeeded in computing the conformal parameterization for all six irregular surface meshes with good conformality. Yamabe flow [20] and inverse distance Ricci flow [42] often fail on the irregular meshes. The conventional Ricci flow [19] is robust in computing the conformal parameterization, however, the conformality is not satisfactory. The double-covering approach, which converts the surface to a genus-0 closed surface and parameterize it using spherical harmonic map (with cotangent formula) [13], is also sensitive to the mesh quality. It fails on the hand and Venus meshes. Also, the conformalities for the successful cases are worse than ours. These experimental results demonstrate that our proposed parameterization method is robust and accurate for computing conformal parameterizations of highly irregular surface meshes.   6.2. Teichmüller adaptive remeshing. Once the conformal parameterization is obtained, the irregular surface mesh can be remeshed by projecting a regular triangulation on the parameter domain onto the surface using the obtained parameterization. However, the regular mesh on the parameter domain should be adaptive to the area distortion under the conformal parameterization. Otherwise, geometric losses may occur. Figure 6(A) and (B) show the surface remeshing results of the foot and hand respectively by projecting a uniform mesh on D onto the surfaces. Note that since the uniform mesh on D is not adaptive to the area distortions of the conformal parameterizations, serious geometric losses near the fingers and toes can be observed. Figure 7 shows the Teichmüller adaptive remeshing results on the parameter domain of the foot. The conformal parameterization of the foot surface is shown in Figure 7(A). The parameterization is conformal but introduces area distortion. Using CVT, we pick sparse sample points according to the density function on the parameter domain and construct the Delaunay triangulation of the sample points to obtain the base mesh, which is shown in 7(B). In 7(C), we build a denser mesh by subdividing. Note that the geometry of the mesh is not optimized. Some triangular faces have sharp angles. Also, the mesh is not smooth. Jumps across the edges of the base mesh can be obviously seen. Using the T-Map, we compute the Teichmüller adaptive mesh on the conformal parameter domain, which is shown in 7(D). The obtained mesh is comparatively much more regular and smooth.
6.3. Surface remeshing results. With the conformal parameterization and the Teichmüller adaptive mesh on the parameter domain, the regular mesh can be projected onto the surface to obtain a remeshed surface. Figure 8(A) shows the remeshed surface by projecting the base mesh onto the surface. Figure 8(B) is obtained by projecting the subdivision mesh of the base mesh onto the surface. Note that the surface remesh is not smooth. Jumps can be seen across the edges of the base mesh. Figure 8(C) shows our remeshing result of our proposed method. The surface mesh is much more regular and smooth. Figure 9 and Figure 10 show the remeshing results of the foot surface at different viewing angles. The left column   shows the original mesh at different angles. The right column shows the remeshed surface at the corresponding viewing angles. Figure 11 shows the remeshing results of the hand surface at different viewing angles. Figure 12, 13, 14 and 15 show the remeshing results of the face, Venus, mask and lion surfaces respectively. The remeshed surfaces have much better triangle qualities than their original surface meshes. Figure 16 and 17 show the histograms of the radius-ratio for surface meshes of the foot and Venus respectively. (A) show the histogram of the radius-ratio of the original mesh. The histograms of the radius-ratio of the remeshed surfaces using    have radius-ratio close to 1. It means the triangular meshes are much improved compared to the original mesh.
Furthermore, note that our proposed Teichmüller adaptive remeshing algorithm has two major components, namely, 1. the subdivision mesh and 2. the Teichmüller map (T-Map). It turns out both components are crucial. In fact, building the subdivision mesh helps to improve the triangle qualities of the remeshed surface.  original meshes, and they are worse than those obtained by our proposed algorithm. Also, deforming the base mesh using the T-Map benefits for improving the qualities of the triangulations. (C) shows the triangle qualities of the surface meshes obtained by projecting the subdivided mesh of the base mesh. Notice that the triangle qualities are not as good as our proposed algorithm.
As for the computational times, for surface meshes with less than 20k faces, the proposed parameterization algorithm takes less than 10s to compute. And the computation for the remeshing algorithm takes less than 2 minutes.
As a matter of fact, the remeshing can also be done by adjusting the vertex positions of a dense 2D regular mesh using CVT directly, and remesh the surface through the conformal parameterization. This direct method can produce meshes  with satisfactory triangle quality, although the computational cost is more expensive. Table 2 shows the comparison between the direct method and the Teichmüller remeshing method. Note that the triangle quality of both methods are comparable, although the direct method produces worse minimum radius-ratio. On the other hand, the Teichmüller remeshing method is more efficient than the direct method.  Our proposed algorithm naturally constructs a structured mesh with subdivision connectivity. The Tutte's embedding K tutte CV T of the Delaunay triangulation K CV T associated to the CVT generators is first constructed and is further refined by subdivision to obtain a sequence of refined meshes of various levels. With the Teichmüller map, the refined meshes are deformed with least conformality distortion, so that the deformed mesh follows the density function. The deformed triangulations of D are then projected back to the surface to get a sequence of remeshed surface with subdivision connectivity. Figure 18 shows the remeshing results with subdivision connectivity of an irregular triangulation of a foot surface. The first row shows the original triangulation at two different angles. The second row shows the remeshing result at the coarsest level. The third and forth rows show the subdivided triangulations at the level 2 and 3 respectively. Note that the refined triangulations of the 2D parameter domain are projected to the original surface using the Teichmüller map, which minimize the conformality distortion. Hence, the remeshing results follow the geometry of the original surface while maintaining good regularity of the triangulations. Figure 19 shows the remeshing results with subdivision connectivity of an irregular triangulation of a human face. Again, the first row shows the original triangulation at different angles. The second, third and forth rows show the remeshing result with subdivision connectivity at various levels. Again, it is observed that the remeshing results follow the geometry of the original surface, while maintaining the regularity of the triangulations.

7.
Conclusion. This paper proposes a method for surface remeshing based on quasi-conformal theories. The main idea is to conformally parameterize the irregular surface mesh onto a simple domain. A regular mesh on the parameter domain can then be projected onto the surface using the obtained parameterization. Obtaining conformal parameterization of a highly irregular surface mesh is challenging. In this work, we propose an effective algorithm to obtain optimal conformal parameterizations of highly irregular surface meshes using quasi-conformal Teichmüller theories. After obtaining the parameterization, a regular mesh on the parameter domain, which is adaptive to the area distortion of the parameterization, has to be built. In this work, we propose an algorithm, called the Teichmüller adaptive remeshing, to obtain an adaptive regular mesh with subdivision connectivity on the conformal parameter domain using the landmark-matching Teichmüller map. Experiments have been carried out to remesh real surface meshes, which show the efficacy of the algorithm to optimize the regularity of an irregular triangulation. Figure 19. Surface remeshing of the mask surface with subdivision connectivity at various levels.