A Variational Principle for Improving 2D Triangle Meshes based on Hyperbolic Volume
In this paper, we consider the problem of improving 2D triangle meshes tessellating planar regions. We propose a new variational principle for improving 2D triangle meshes where the energy functional is a convex function over the angle structures whose maximizer is unique and consists only of equilateral triangles. This energy functional is related to hyperbolic volume of ideal 3-simplex. Even with extra constraints on the angles for embedding the mesh into the plane and preserving the boundary, the energy functional remains well-behaved. We devise an efficient algorithm for maximizing the energy functional over these extra constraints. We apply our algorithm to various datasets and compare its performance with that of CVT. The experimental results show that our algorithm produces the meshes with both the angles and the aspect ratios of triangles lying in tighter intervals.
In this paper, we consider the problem of improving 2D triangle meshes tessellating planar regions. The applications in scientific computing require quality meshes. The quality here refers to the shape and size of the elements: a triangle is of good shape if it is close to the equilateral triangle, i.e., its inner angles are close to . The most popular approach for generating quality meshes is Delaunay refinement . Delaunay refinement algorithms commonly perform one local change at a time, until the criteria of the shape and size of elements is satisfied. They are greedy approaches which may produce bad shaped triangles, especially near the boundary. Quite a few methods have been proposed to deal with this issue. Among them, centroidal Voronoi tessellation (CVT) is widely used where the vertices are iteratively moved to the barycenters of the corresponding Voronoi cells. There is a variational principle associated with this so-call Lloyd iteration where the energy functional measures the difference between the site points and the barycenters of the corresponding Voronoi cells .
This paper proposes a new variational principle for improving 2D triangle meshes where the energy functional is a convex function over the angle structures whose maximizer is unique and consists only of equilateral triangles. This energy functional is related to hyperbolic volume of ideal 3-simplex. Of course, one needs to impose the extra constraints on the angles to embed the mesh into the plane and preserve the boundary. Nevertheless, we show the energy functional is still well-behaved even with these constraints. The space of angle structures is much bigger than the space of coordinates, which provides more freedom for algorithms to search for maximizer and thus find better (local) maximizer and generate better meshes, as demonstrated in Section 5. We devise an algorithm based on interior-point method for maximizing the energy functional over extra constraints. We apply our algorithm to various datasets and compare its performance with that of CVT. The experimental results show that our algorithm produces the meshes with both the angles and the aspect ratios of triangles lying in tighter intervals.
Previous work: There are great amounts of research work on quality mesh generation and many meshing strategies have been proposed and studied. Here we summarize those most relevant to our work. Readers are referred to  and the references therein for more related work on quality mesh generation. The most popular approach is Delaunay refinement, which iteratively inserts Steiner points to improve the quality of the mesh until the initial criteria is satisfied for each triangle. Delaunay refinement approach is pioneered by Chew , and later improved and extended by many others[9, 10]. Shewchuk  shows that it terminates with a finite number of Steiner points and with bounds on the angles. Delaunay refinement can be improved either by carefully designing the order of inserting Steiner points or choosing the positions of Steiner points other than the circumcenters of triangles . However, it remains a greedy approach, which may make globally bad decisions which are not reversible. To address this issue, variational approaches are proposed where an energy functional is chosen so that the low levels of this energy correspond to the meshes with good quality. The widely used energy functional is the one used in CVT which sums the differences between the site points and the barycenters of the corresponding Voronoi cells . This is based on the observation that in 2D, evenly distributed points lead to well-shaped triangles in Delaunay triangulation . Du et al.  proposed the Lloyd iteration to transforms an initial ordinary Voronoi diagram into a centroidal Voronoi diagram. Finally, Tournois et al.  proposed to interleave Delaunay refinement and CVT for generating and improving 2D triangle meshes.
2 The Energy Functional
In this section, we describe the energy functional and discuss its properties. We start with a single triangle . Let be the inner angles of . Assign the following energy.
where is Lobachevsky function:
This energy assigned to is in fact the volume of an ideal hyperbolic 3-simplex. Consider the upper half space: with the hyperbolic metric . Place the triangle on the plane , and add the fourth vertex at infinity. Then all four vertices are at infinity, and thus form an ideal hyperbolic 3-simples, denoted . Place the vertex at a position so that the dihedral angles along three edges meeting at are the inner angles of the triangle , as shown in Figure 2. The fact that is the volume of the ideal 3-simple follows from the following lemma.
Lemma 2.1 (e.g., )
Consider an ideal hyperbolic 3-simplex, that is a simplex with all four vertices at infinity. If , , are the dihedral angles along three edges meeting at a common vertex, then , and
Remark that it does not matter which particular vertex we choose, since it follows easily that the dihedral angles along the opposite edges of are equal so that we have the same three dihedral angles , , incident to any vertex.
There are many nice properties of the energy function . Here we state two of them which are most relevant to our setting. The following lemma says that reaches maximum when is equilateral.
Lemma 2.2 (e.g., )
The volume of a hyperbolic 3-simplex reaches the maximum when .
Since , is a function of two inner angles and which parametrize the space of all Euclidean triangles up to similar transformations. The following lemma tells is a strictly concave function over the space of Euclidean triangles up to similar transformations.
Lemma 2.3 ()
The Hessian of is
and is negative definite.
Now we are ready to define the energy functional for a triangle mesh . Denote the sets of vertices, edges and triangles of by , and . We identify a vertex in with an index, i.e., , where is the number of vertices in . We denote by the edge with vertices and , by , the triangle with vertices , and . The triangles in are Euclidean. Let denote the inner angle at vertex in triangle . Let denote the cardinality of a set . Define the energy functional as the sum of the energy over all triangles in , i.e.,
Let denote all possible angle structures given the combinatorial structure of , i.e.,
If the combinatorial structure of is fixed, the energy is a function over the angle structures . The following lemma follows easily from Lemma 2.3.
is a concave function over the angle structures .
We remark that (1) From Lemma 2.2, reaches the maximum over when all triangles in become equilateral, (2) In the meshing application considered in this paper, additional constraints are necessary to impose on the angles as we will discuss in the next section. Nevertheless, the triangles in the mesh become more well-shaped (closer to equilateral) when the energy increases.
3 Angle Structures and Embeddings
In the paper, we consider the problem of improving a triangle mesh which tessellates a planar region. Therefore the triangle mesh is embedded in the plane where each vertex has a coordinate in the plane. The quality of triangles are improved by adjusting the coordinates of vertices. On the other hand, our energy functional is defined over the angles. In our method, the coordinates of the vertices are adjusted by changing the angles. Therefore, it is necessary to relate the embeddings of and its angle structures. We relate them using the metric. The metric of a triangle mesh specifies the length for each edge in the mesh, or equivalently, is a function such that the triangle inequalities hold for each triangle, i.e., for triangle .
It is obvious that an embedding of a triangle mesh, or equivalently, the coordinates , induces a metric for the triangle mesh where , for any , and a metric of a triangle mesh induces an angle structure where for any . Furthermore, the induced angle structure is an invariant of rigid transformations (translations, rotations and reflections) and uniform scaling of the embedding of the triangle mesh. In fact, one can recover the embedding from the induced angle structure as follows. Observe that if the inner angles of a triangle are given, one can calculate the coordinate of the third vertex from the coordinates of the other two. First, pick a triangle in and embed it into the plane with given inner angles, i.e., compute the coordinates for its three vertices. Up to a rigid transformation and a uniform scaling, these coordinates are uniquely determined. Then consider its neighboring triangles, each of which has two of its vertices already embedded into the plane. Based on the previous observation, the coordinate of its third vertex is uniquely determined from the angles and can be easily computed. Once embed these neighboring triangles in the plane, consider their neighboring triangles and repeat the above procedure until all triangles get embedded into the plane. This leads to an algorithm to layout a triangle mesh based on its angle structure. See Algorithm 3.
On the other hand, not every angle structure is induced from an embedding of on the plane. To see what are the constraints that the angle structure induced from an embedded triangle mesh satisfies, consider the one-ring neighbor of a vertex either in the interior or on the boundary, as shown in Figure 3. We use the following two quantities to describe the constraints. One quantity is the the angle sum at vertex given an angle structure :
The other quantity is the so-called holonomy at vertex given an angle structure 111Precisely, it is the holonomy of a sequence of the triangles incident to vertex . See :
Consider an interior vertex of a triangle mesh on the plane. See Figure 3(a). First, the angle sum at interior vertex has to be , i.e.,
Second, notice that if let be the length of the edge opposite to vertex in triangle , then . By law of sines, we have for any interior vertex
One can show that if its angle structure satisfies equation (8, 9) for each interior vertex , the triangle mesh is locally flat and can then be immersed into the plane. An immersion is locally a one-to-one map. However, note that an immersion is not necessary an embedding, which may have global self-intersections . It is relatively hard to impose or even describe the constraints on the angles to circumvent global self-intersections. Fortunately, we have not observed such global self-intersections in our experiments. This may be due to the fact that we also preserve the boundary as described below, which makes them very rare.
In addition, as our purpose is to improve the quality of a mesh tessellating a fixed planar region, we want to preserve the boundary. For simplicity, assume the boundary has one connected component. See section 4 for how we deal with multiple connected components on the boundary. Consider a vertex on the boundary. Let and are two edges on the boundary incident to vertex . See Figure 3(b). First, the angle sum at a boundary vertex need to be preserved, i.e.,
where is the angle between and containing the interior of the planar region. Second, the ratio between the length of two consecutive edges on the boundary need to be preserved. One can write this ratio in terms of the holonomy of the boundary vertex
where are the length of respectively. Note if the angle structure satisfies equation (10, 11) for each boundary vertex , then the boundary is preserved, up to a rigid transformation and a uniform scaling.
In summary, we have two types of extra constraints imposed on the angle structures so that the triangle mesh can be immersed into the plane with the shape of the boundary preserved. One type is the angle sum imposed on each vertex. Given with specifying the angle sum at vertex , define a subset of as
The other type is the holonomy condition imposed on each vertex. Given with specifying the holonomy at vertex , define another subset of as
Since the total sum of all the inner angles has to equal , there are only number of independent equality constraints in defining the subset . In addition, the sum off the holonomy of all the vertices has to be since each inner angle appears twice in this sum, once positive and once negative. So there are also number of independent equality constraints in defining the subset .
Observe that the equality constraints on the angle sum is linear in angles and thus is always a convex subset of , which leads to the following lemma.
is a concave function over the subset of the angle structures for any . Furthermore, if is an extremal point of , then the holonomy at each interior vertex is automatically .
The proof for an interior vertex having holonomy given an extremal point of can be find for example in . However, in order to preserve the shape of the boundary, one need to impose the extra holonomy conditions to all boundary vertices, which unfortunately does not hold automatically. These equality constraints on the holonomy are nonlinear in angles. Thus the maximization of the energy functional becomes non-convex.
4 The Algorithms
In this section, we describe an algorithm which takes input a triangle mesh tessellating a region on the plane, and outputs another triangle mesh tessellating the same region with the shape of the triangles improved to closer to the equilateral triangles. We sketch our remeshing algorithm in pseudo-code as follows:
When there are multiple connected components on the boundary, the constraints described in Section 3 only preserve the shape for each component but not their relative position. To deal with this issue, in the first step, we connect the different components using the shortest paths and then cut the mesh along them into a topological disk with only one boundary component. In this step, any vertex on a cutting path may be split into several copies in each of which takes as its initial coordinates. In the final step of outputting the new mesh, these copies of vertex in have the same new coordinates as the shape of the boundary is preserved. So just take one of them as the new coordinate for vertex in . Figure 4 illustrates the algorithm. As we can see, that our algorithm improves the quality of the triangles, especially those near the boundary.
|input mesh||cut paths||output mesh|
In this subsection, we describe an algorithm to find an angle structure which maximizes the energy functional over the subset of angle structures. This is the key part of the algorithm. The subset is nonlinear in angles. Thus it is an optimization problem over nonlinear constraints. We use the interior-point method and follow the implementation of the Matlab routing fmincon to solve this optimization problem.
Since there are number of nonlinear equality constraints in defining the subset , it is hard for the interior-point method to search for feasible solutions in a subset of this high codimension. To address this issue, we break the optimization procedure into two steps. In the first step. we maximize over . By Lemma 3.1, is concave over and has a unique maximum. This step can be done very efficiently using the routing . In the second step, we minimize another energy functional over where given a vector ,
for any . measures how much the angle structure violates the nonlinear equality constraints, and reaches the minimum when is in . This minimization is the most time consuming step of the algorithm. In the implementation, we supply the gradient and the Hessian matrix of both energy functional and to the routing fmincon, which significantly improves the efficiency of the algorithm. The pseudo-code of the algorithm argmax is follows:
Note that in Algorithm 2, the second step of minimizing is necessary for preserving the boundary but may decrease and deteriorate the quality of triangles. In fact, the initial angle structure is a minimizer of . To see the performance of this step, we tested the following procedure to maximize over . For any , define
Notice that . We relax the constraints of nonlinear equality to inequalities and then maximize over an enlarged subset for some , and then reduce within along its negative gradient to . Repeat the above steps with . In this way, we make sure the final is a local maximum of in . We observe that this procedure produces a mesh with the same quality as argmax but needs significantly more computation time. This shows that the second step of minimizing somehow also respects the energy functional well, which is worth further investigation.
4.2 Layout Mesh
The procedure of layout a triangle mesh from its angle structure is described in Section 3. Since the shape of the boundary is preserved, we embed the vertices on the boundary using their original coordinates. Now for any triangle with one edge on the boundary, the coordinates of the third vertex can then be computed uniquely from the angles. Then propagate and compute the coordinates for the other vertices. We sketch the algorithm layout in pseudo-code as follows:
4.3 Cut Mesh
In this subsection, we describe a method to implement the cut algorithm. The basic idea is to find some paths to connect the different connected components together and then cut the mesh along these paths. The pseudo-code of algorithm cut is as follows:
Complexity: The complexity of the algorithm is dominated by optimization step. The complexity of algorithm layout and algorithm cut are and respectively. See section 5 for the performance of the algorithm over various datasets.
In this section, we apply our meshing algorithm to various datasets and show its performance. All input meshes are obtained by the Delaunay refinement algorithm triangle by Shewchuk .
Figure 5 and 6 show two meshes and compare our algorithm with the standard centroidal Voronoi tessellation (CVT) where the boundary is fixed. Following , we measure the quality of elements based on its inner angles and aspect ratios (circumradii to shortest edge ratios). As we can see, our method performs better in sense while CVT may do better in sense. Namely, both the angles and the aspect ratios in the meshes generated by our method lie in tighter intervals than CVT, while the meshes generated by CVT may have more triangles close to the equilateral one.
|Delaunay refinement||Our method||CVT|
|Delaunay refinement||Our method||CVT|
Table 1 shows the timings of the procedure argmax over input meshes with increasing number of vertices. The other procedures of the algorithm are negligible in terms of timing. Figure 7 shows the meshes with three different resolutions for each model. As we can see, the procedure of minimizing consumes most of computation time, which also increases faster as the number of vertices increases.
|in model H||393||756||1476||3584||7072|
|in model Hole3||340||660||1468||3136||6135|
In this paper, we have proposed a new variational principle for improving 2D triangle meshes based on hyperbolic volume, devised an efficient algorithm to maximize the energy functional over nonlinear constraints and to improve the quality of meshes, and applied our algorithm to various datasets and compared its performance to CVT. Here we point out a couple of possible directions for future work. First, notice that the combinatorial structures of input meshes are fixed in the current framework. However, one can also make them Delaunay 222The sum of the opposite angles is less than for each interior edge by flipping edges to improve the quality of meshes. It is interesting to see how the energy functional changes while flipping edges. Second, notice that the energy functional can be extended to higher dimensional spaces by considering higher dimensional ideal simplex. Thus one can follow the similar framework to improve higher dimensional meshes.
The authors acknowledge The National Basic Research Program of China (973 Program 2012CB825501); Tsinghua National Laboratory for Information Science and Technologyï¼TNListï¼Cross-discipline Foundation.
-  Boris Springborn Alexander Bobenko, Ulrich Pinkall. Discrete conformal maps and ideal hyperbolic polyhedra. arXiv:1005.2698, 2010.
-  Siu-Wing Cheng, Tamal Krishna Dey, and Jonathan Richard Shewchuk. Delaunay Mesh Generation. CRC Press, 2012.
-  L. P. Chew. Constrained delaunay triangulations. In Proceedings of the third annual symposium on Computational geometry, SCG ’87, pages 215–222, New York, NY, USA, 1987. ACM.
-  Qiang Du, Vance Faber, and Max Gunzburger. Centroidal voronoi tessellations: Applications and algorithms. SIAM Rev., 41(4):637–676, December 1999.
-  David Eppstein. Global optimization of mesh quality. Tutorial at the 10th International Meshing Roundtable, 2001.
-  Morris W. Hirsch. Differential Topology. Springer, 1997.
-  John Milnor. Hyperbolic geometry, the first 150 years. Bull.AMS, 6:9–24, 1982.
-  Igor Rivin. Euclidean structures on simplicial surfaces and hyperbolic volume. Ann. Math., 139:553–580, 1994.
-  Jim Ruppert. A delaunay refinement algorithm for quality 2-dimensional mesh generation. J. Algorithms, 18(3):548–585, May 1995.
-  Jonathan Richard Shewchuk. Delaunay refinement algorithms for triangular mesh generation. Comput. Geom. Theory Appl., 22(1-3):21–74, May 2002.
-  Jane Tournois, Camille Wormser, Pierre Alliez, and Mathieu Desbrun. Interleaving delaunay refinement and optimization for practical isotropic tetrahedron mesh generation. ACM Trans. Graph., 28(3):75:1–75:9, July 2009.
-  Alper Üngör. Off-centers: A new type of steiner points for computing size-optimal quality-guaranteed delaunay triangulations. Comput. Geom. Theory Appl., 42(2):109–118, February 2009.