There are simple and robust refinements (almost) as good as Delaunay

Original article There are simple and robust refinements (almost) as good as Delaunay Alberto Márquez a, Auxiliadora Moreno-González a, Ángel Plaza b,∗, José P. Suárez c a Department of Applied Mathematics I, University of Seville, Spain b Department of Mathematics, University of Las Palmas de Gran Canaria, 35017 Las Palmas de Gran Canaria, Spain c Department of Cartography and Graphic Engineering, University of Las Palmas de Gran Canaria, Spain Abstract A new edge-based partition for triangle meshes is presented, the Seven Triangle Quasi-Delaunay partition (7T-QD). The proposed partition joins together ideas of the Seven Triangle Longest-Edge partition (7T-LE), and the classical criteria for constructing Delaunay meshes. The new partition performs similarly compared to the Delaunay triangulation (7T-D) with the benefit of being more robust and with a cheaper cost in computation. It will be proved that in most of the cases the 7T-QD is equal to the 7T-D. In addition, numerical tests will show that the difference on the minimum angle obtained by the 7T-QD and by the 7T-D is negligible. Keywords: Longest-edge; Edge-refinement; Delaunay 1. Introduction Longest edge-based refinement has become popular in last decade in the context of mesh refinement [9,11,14].A well accepted acronym used to name this class of mesh subdivision is nT-LE, where nis the number of new triangles (T) produced after a single subdivision and LE stands for longest edge. So, we found in the literature well studied longest-edge partition as 2T-LE, 3T-LE, 4T-LE and 7T-LE, see [14,15,9] and the references therein. It should be noted that the iterative application of these partitions yields good-quality meshes, in the sense that they do not degenerate. Additionally, longest edge refinements have the advantage of its propagation, i.e., if we subdivide a triangle, we know how to subdivide its adjacent triangles in order to obtain a conforming triangulation. Of course, if we add some points to a triangle, and we want to obtain a subdivision with the best quality (in the sense that we want to maximize the minimum angle), the optimal solution is the Delaunay triangulation (see, for example [1,2]). Mesh generation algorithms based on Delaunay refinement have been effective tools both in theory and in practice in the last 20 years [3,19,5]. The first provably good Delaunay refinement algorithm is due to Chew [4]. Much attention has received this class of algorithms afterwards, in particular thanks to authors like Ruppert [16] and Shewchuk [17,18] among others. Longest-edge based algorithms have been used together with Delaunay triangulation for the quality triangulation problem by Rivara and co-workers [8,13]. ∗Corresponding author. Tel.: +34 928 45 88 27; fax: +34 928 45 87 11. E-mail address: [email protected] (Á. Plaza). A. Márquez et al. / Mathematics and Computers in Simulation 106 (2014) 84–94 85 Unfortunately, Delaunay refinement algorithms present some disadvantages from a practical point of view: on one hand they require a considerable amount of computation, and, on the other hand, they are not robust dealing with point sets that are not in general position (a point set is in general position if no three of them lie on the same straight line and no four lie on the same circle). It should be noted here that longest edge refinements always present point sets that are not in general position (we will discuss this topic in the next section). In order to solve these inconveniences, some refinements that perform only comparisons of distances have been presented (see [14,15,9]. Of course, to measure a distance between pairs of points is easier and far more robust than compare angles or doing circumcircle tests, and so in the mentioned works (and in many others) the authors avoid the theoretical advantages of Delaunay triangulations for the sake of the simplicity and robustness. In this way, many methods appear as those mentioned above [14,15,9]. In this work, we try to obtain the best of both worlds, we propose a longest edge refinement (the 7T-QD refinement) that can be obtained performing only lengths comparisons, and that in most of the cases, actually more than in 97% of the triangles, coincides with the Delaunay triangulation. Even in the cases that the 7T-QD is not a Delaunay triangulation, we are not far from that optimal refinement in the following ways: five out of seven of the triangles that are obtained with the 7T-QD refinement from a original triangle are Delaunay triangles, and the minimum angle of the other two triangles are, in the worst case, only a 20% worse than the minimum angle of the Delaunay refinement, but the refinement presented is better if we measure the average of the minimum angle of the two triangles. The structure of this paper is as follows: Section 2 gives a short background of the class of refinement methods treated in this paper, Section 3 introduces the Seven Quasi-Delaunay partition for triangles and gives a comprehensive comparison with the pure Delaunay triangulation and the Seven Triangles Longest-Edge partition. In Section 4,we provide a numerical study considering the min angle that stress the quality of 7T-QD. Finally some useful conclusions regarding the new introduced partition are offered. 2. Triangulation with naligned points By locating midpoints on the edges of the triangle we can compute some quality triangulations in the plane, see [9] and the references therein. This can be viewed as triangulation with naligned points. One of the interesting points of such family of triangulations is the low cost for obtaining the subdivision. In last decade, subdivision methods inserting one point per edge – based on the longest edge or not – have been sufficiently explored. Some of these partitions are: “red–green”, longest edge bisection (2T-LE), and four triangles longest edge bisection (4T-LE). Less attention has been given to subdivisions based on the insertion of two points per edge. Lastly, the Seven Triangle Longest-Edge partition (7T-LE) has been presented in [9]. The 7T-LE partition of a triangle tis obtained by putting two equally spaced points per edge. After cutting off three triangles at the corners, the remaining hexagon is subdivided further by joining each point of the longest-edge of tto the base points of the opposite sub-triangle. Finally, the interior quadrangle is subdivided into two sub-triangles by the shortest diagonal, see Fig. 1. Fig. 1. Scheme for the 7T-LE partition of triangle t0and new class of triangles generated, t1and tn1. 86 A. Márquez et al. / Mathematics and Computers in Simulation 106 (2014) 84–94 Fig. 2. Schemes for 7T-LE and 7T-D partitions. It has been proved in [9] that this partition is superior to 2T-LE and 4T-LE in terms of the quality of new triangles generated. Some other partitions also based in the trisection of the edges may be considered. We call 7T-D partition of a triangle tthe Delaunay triangulation of the cloud of points formed by the three vertices of tand the two equally spaced points per edge (Fig. 2). It is clear that Delaunay based triangulations in 2D are optimal, maximize the minimum angle. The main cost to compute a Delaunay triangulation is that the circle circumscribing any Delaunay triangle does not contain any other input points in its interior. This test can be achieved for example by computing a matrix determinant. Although some efficient algorithm to obtain Delaunay refinements are known, they take O(nlogn), in part, due to the circumcircle test. A known problem with the triangulation of naligned points is the case of collinearity or degenerated cases. For a set of points on the same line there is no Delaunay triangulation, in the sense that the notion of triangulation is degenerate for this case. Some recent cases have been reported by the authors in which some implementations of Delaunay algorithm fails in triangulating a given triangle with equally spaced points per edge. One case is the implementation of CGAL in Matlab R2009b. We try the triangulation of a set of points p, with coordinates X=p(1, :) and Y=p(2, :), see Fig. 3. The output of Delaunay command in Matlab for that set of points reveals a triangle with collinear coordinates, triangle 423. It should be pointed out that when introducing a point in the middle of an edge inevitable floating point errors would make the three points not perfectly collinear paving the way to the creation of degenerate triangles. A numerical technique to cope with this problem is the Simulation of Simplicity (SoS) [6,10]. Intuitively, it simulates a conceptual perturbation of the input data that eliminates all degeneracies. The basic idea of SoS is to perturb the given objects slightly, which amounts to changing the numbers that represent the objects. In our case this perturbation means changing the coordinates of the aligned points. A class of triangulation algorithms that help to overcome such problems are edge-refinement methods, based in the insertion of naligned points. In practice one can expect to introduce one or two points per edge, but a greater number of points may produce degeneracy of the triangulation. Fig. 3. A failed output triangulation using CGAL in Matlab R2009b. A. Márquez et al. / Mathematics and Computers in Simulation 106 (2014) 84–94 87 Fig. 4. A mapping diagram for normalized triangles. 3. The Seven Triangle Quasi-Delaunay partition Taking into account the superiority of the Delaunay triangulation on the maximizing the minimum angle, our goal in this paper is to introduce a simple and robust partition as similar as possible to the Delaunay triangulation of the cloud of points formed by the three vertices and two equally spaced points per edge in each triangle. At the same time, we should be able to overcome the collinearity problem above mentioned. Following our idea of subdividing triangles inserting two aligned points, we present here the Seven Triangle Quasi-Delaunay (7T-QD) partition. This new partition will be equivalent to the Delaunay partition except for an small region of triangles which amounts less that 3% of the area of the geometric diagram for the triangles. On addition, the computation of the different patterns used by the new paradigm is simpler than the circumcircle test. Finally, in most of the cases the new introduced partition will be also superior to the previous Seven Triangle Longest-Edge partition. Only for less than 0.4% of the triangles the minimum angle of the 7T-LE will be greater than the minimum angle of the new 7T-QD partition. Before giving a definition of such partition we use here a mapping diagram to normalize triangles [12], that facilitates the understanding of the method and its comparison with other partitions. The mapping diagram is as follows: for a given triangle we scale the longest edge to have unit length. This forms the base of the diagram. It follows that the set of all triangles is bounded by this horizontal segment (longest edge) defined by the points (0, 0), (1, 0), and by two bounding exterior circular arcs of unit radius, centered, respectively, at (1, 0) and at (0, 0), as shown in Fig. 4. In the figure it is showed the boundary curves for the diagram and a sample triangle represented with apex labeled by “O”. It is not difficult to see that the diagram and the triangles represented therein have symmetry with respect to x=1/2. For this reason, in the following we will consider only the left half part of the diagram. Definition 1. Let MON be a triangle of longest edge MN and apex O, with coordinates M= (0, 0), N= (1, 0) and O=(x,y), with x≤1/2. The Seven Triangle Quasi-Delaunay (7T-QD) partition is defined as follows, see Fig. 5: (1) Position two equally spaced points per edge (let us name such points as indicated in Fig. 5). (2) Construct triangles MPD,ADP,BCQ and CNQ. (3) Subdivide the interior pentagon, depending on the location of the apex, as follows: (a) Case I: |AB|<|OD|. Construct triangles ABO,ABD and BCD. (b) Case II: |AB|>|OD| and |BD|<|OC|. Construct triangles ADO,BDO and BCD. (c) In other case, construct triangles ADO,CDO and BCO. Notice that the three possibilities at dividing the interior pentagon ABCDO correspond respectively to three different regions into the mapping diagram: region I, II and III in Fig. 6. The boundary curves between these regions correspond to circular arcs of equations (x−1/3)2+y2= 1/9, and (x−6/5)2+y2= 16/25. As we did note above, by reflection around the vertical line x=1/2 in the mapping diagram, it is enough to consider the left half of the diagram, which is equivalent 88 A. Márquez et al. / Mathematics and Computers in Simulation 106 (2014) 84–94 Fig. 5. Scheme for the 7T-Quasi-Delaunay partition. to positioning the smallest angle at point (1, 0) and the second largest angle at point (0, 0). We can formulate new conditions that let us extend the partition to the right part of the diagram (not treated here for the sake of brevity). Fig. 4 shows the interior curves delimiting regions of interest. It is worth to note that the 7T-QD partition introduces seven triangles, from which four of them are formed in a fixed way independently of the targeted triangle. These four triangles are Delaunay as will be shown below. The other three triangles are constructed depending on the length of the diagonals of the interior pentagon ABCDO, see Fig. 5. Fig. 6. Regions of different patterns for the 7T-QD partition. A. Márquez et al. / Mathematics and Computers in Simulation 106 (2014) 84–94 89 Lemma 1. The corner triangles on the longest edge,MPD and CNQ in Fig. 5, are Delaunay triangles. Proof. To prove this result, the circumcircle test is going to be used. That is, given triangles MPD and CNQ, whose circumcircles do not contain any other point of the initial set. First of all, consider triangle t0=MPD and denote Ct0the circumcenter of t0. By the initial conditions in the construction of the 7T-QD partition, M= (0, 0), D= ((1/3), 0) and P=((x/3), (y/3)), so Ct0=((1/6),((x2−x+y2)/(6y))). It is obvious that points C,N,Aand Oare outside the circumcircle of t0, so it only must be proved that points Band Qare also outside it. Consider then B=1+2x 3,2y 3and Q=2+x 3,y 3. In order to show that Bis outside the circumcircle of t0, it is going to be proved that d(M, Ct0)≤d(B, Ct0). Now, d(M, Ct0)=1 62 +x2−x+y2 6y2 and d(B, Ct0)=1+2x 3−1 62 +2y 3−x2−x+y2 6y2 so d(M, Ct0)≤d(B, Ct0)if 1 62 +x2−x+y2 6y2 ≤1+2x 3−1 62 +2y 3−x2−x+y2 6y2 , 0≤1+2x 32 −1+2x 9+4y2 9−2 9(x2−x+y2), Therefore, d(M, Ct0)≤d(B, Ct0)if0≤x2+y2+2x. But, by the initial conditions x,y≥0, so Bis outside the circumcircle of t0. In the same way, to prove that Qis outside the circumcircle of t0, it must be proved that d(M, Ct0)≤d(Q, Ct0). Consider then d(Q, Ct0)=x+2 3−1 62 +y 3−x2−x+y2 6y2 . In this case, d(M, Ct0)≤d(Q, Ct0)if 1 62 +x2−x+y2 6y2 ≤x+2 3−1 62 +y 3−x2−x+y2 6y2 , 0≤x+2 32 −x+2 9+y2 9−1 9(x2−x+y2), Hence, d(M, Ct0)≤d(Q, Ct0)ifx≥−1/2, that it is true by the initial conditions, so Qis outside the circumcircle of t0and triangle t=MPD is a Delaunay triangle. Finally, by symmetry, it can be proved that triangle CNQ is a Delaunay triangle too.  90 A. Márquez et al. / Mathematics and Computers in Simulation 106 (2014) 84–94 Lemma 2. Triangles ADP and BCQ are Delaunay triangles. Proof. To prove this result, the circumcircle test is going to be used as in the previous one. First of all, consider triangle t1=ADP and denote Ct1the circumcenter of t1. By the initial conditions in the construction of the 7T-QD partition, A=(2x/3, 2y/3), D= (1/3, 0) and P=(x/3, y/3), so Ct1=(xt1,y t1)=2x2+2y2+1 6,x(x2+y2)(2x2+2y2+1) 36y. It is obvious that points Mand Oare outside the circumcircle of t1, so it must be proved that points B,C,Nand Q are also outside it. Consider then B=1+2x 3,2y 3,C =2 3,0,N=(1,0) and Q=2+x 3,y 3. It must be noted that any of these four points has the same y-coordinate than one of the three vertices of triangle t1, so in order to prove that these four points are outside the circumcircle of t1, the point that determines its radius must be chosen properly. First of all, to show that Bis outside the circumcircle of t1, it is going to be proved that d(A, Ct1)≤d(B, Ct1). Now, d(A, Ct1)=2x 3−xt12 +2y 3−yt12 and d(B, Ct1)=(1 +2x) 3−xt12 +2y 3−yt12 , so d(A, Ct1)≤d(B, Ct1)if 2x 3−2x2+2y2+1 62 ≤(1 +2x) 3−2x2+2y2+1 62 , 0≤(4x+1) 9−2x2+2y2+1 9. Therefore, d(A, Ct1)≤d(B, Ct1)if(x−1)2+y2≤1, which it evident because point O=(x,y) is inside the mapping diagram. In analogous way it is proved that Cis outside the circumcircle of t1if x2+y2≤1, that is, if point O=(x,y) is inside the mapping diagram. Similar arguments apply to demonstrate that Nis outside the circumcircle of t1, which is equivalent to d(D, Ct1)≤ d(N, Ct1); and also to show that Qis outside the circumcircle of t1, which is equivalent to d(P, Ct1)≤d(Q, Ct1). By symmetry, it can be proved that triangle BCQ is a Delaunay triangle too.  To see that 7T-QD enjoys good properties in mesh generation it is necessary to compare it with other similar partitions. The mapping diagram for the 7T-QD partition reveals that, in terms of the minimum angle, and for those triangles within the regions I and II, the partition is the same as the 7T-Delaunay partition. However, in a small subset of region III they are different. This subset appears shaded in grey color in Fig. 7. The equation of the curve separating the region in which these two partitions are different can be obtained using the circumcircle test: (x−1)2+y2=1/3. Notice that the area of the shaded region is approximately 0.0091 which means that the ratio over the total area of the mapping diagram results 0.0296796. This implies that less than 3% of the triangles yield different ways of dividing by the 7T-QD and the 7T-D. A. Márquez et al. / Mathematics and Computers in Simulation 106 (2014) 84–94 91 Fig. 7. In color the subset of region III in which 7T-QD and 7T-D differ. 4. Some numerical tests on 7T-Quasi-Delaunay partition We present here some numerical studies that show the practice behavior of the 7T-QD partition in comparison with 7T-LE and 7T-D partitions. The study consists in getting representative triangle cases from regions I, II, and III and then applying the partitions 7T-LE, 7T-QD and 7T-D to those triangles. Notice that for triangles in region I the three partitions are the same. As we have mentioned before, for triangles in region II 7T-QD coincides with 7T-D. Also in the central part of region III 7T-QD is equal than 7T-D. Therefore we will focus here the comparison for triangles in the colored subset of region III. See Fig. 7. Fig. 8(a) and (b) shows the minimum angle and its mean, for the triangles generated into the interior pentagon by the three partitions 7T-QD, 7T-D and 7T-LE and when the initial triangle is chosen into the colored region of interest. The triangles are taken inside this region from bottom to top and from left to right. Then, the resulting min angle and mean min angle are sorted from lower to higher and graphed. Note that the minimum angle for the 7T-QD is close to the minimum angle for the 7T-D partition. Fig. 9 shows in detail the distribution of the values for the minimum angles corresponding to the 7T-QD and 7T-D partitions in gray levels. Triangles near to the right side of the region present similar values for the minimum angles generated, while triangles near to the left side show the higher differences. This Fig. 8. Minimum angle comparison for 7T-QD, 7T-D and 7T-LE partitions. 92 A. Márquez et al. / Mathematics and Computers in Simulation 106 (2014) 84–94 Fig. 9. Distribution of minimum angle for 7T-QD and 7T-D partitions. is explained taking into account that 7T-QD and 7T-D partitions are the same out of the colored region and the pattern for the 7T-QD changes precisely at the left side boundary of the colored region. On the other hand, in comparison with the 7T-LE partition, the minimum angle obtained by the 7T-QD is higher. Only for a small number of triangles the 7T-LE presents a minimum angle greater than the one of the 7T-QD. The region of these triangles is given in Fig. 10. The area of this region is approximately 0.0014270 and the ratio with the total area of the mapping diagram is 0.00464707. In order to a deeper comparison we consider now the mean of the minimum angles. The results are shown in Fig. 8(b). It is surprising that in most of the cases, the new introduced partition 7T-QD is better even than the 7T-D partition. The region for these triangles is colored in black in Fig. 11. Note that in the only region in which the minimum angle generated by the 7T-LE is greater than the generated by the 7T-QD, the mean of the minimum angles is lower. It is clear that the partitions 7T-QD, 7T-D and 7T-LE perform differently depending on the triangle shapes and so, on the region of study in the diagram. But it is also clear that this difference concentrates only on a small located area that supposes the 3% of the total area in the diagram. To better clarify the behavior of min and mean min-angle we collect in Table 1 the possible cases of regions in the diagram, Regions I, II and III with the 97% and regions C1, C2 and C3 with the other 3%. To indicate that a partition is better than other according to min or mean min-angle we use symbol < and > otherwise. If partition is similar according to these values we use symbol =. Fig. 10. Region C1 in which the minimum angle obtained by the 7T-LE partition is better than the obtained by the 7T-QD.