Mesh smoothing algorithm based on exterior angles split.

Since meshes of poor quality give rise to low accuracy in finite element analysis and kinds of inconveniences in many other applications, mesh smoothing is widely used as an essential technique for the improvement of mesh quality. With respect to this issue, the main contribution of this paper is that a novel mesh smoothing method based on an exterior-angle-split process is proposed. The proposed method contains three main stages: the first stage is independent element geometric transformation performed by exterior-angle-split operations, treating elements unconnected; the second stage is to offset scaling and displacement induced by element transformation; the third stage is to determine the final positions of nodes with a weighted strategy. Theoretical proof describes the regularity of this method and many numerical experiments illustrate its convergence. Not only is this method applicable for triangular mesh, but also can be naturally extended to arbitrary polygonal surface mesh. Quality improvements of demonstrations on triangular and quadrilateral meshes show the effectiveness of this method.


Introduction
Nowadays, mesh casts an irreplaceable role in extensive fields as an effective way of discretization. In advanced manufacturing field, there is finite element analysis for structures, fluids, electromagnetics and thermodynamics; in biological field, there is medical imaging, organs 3D printing and surgical simulation; in the field of computer vision and graphics, there is Simultaneous Localization and Mapping (SLAM), graphical semantic segmentation, Geographic Information System (GIS), video special effects and so on. There are many mature mesh generation approaches, such as marching cubes [1], Delaunay [2], advancing front [3,4], frame field [5] and so on. But the mesh quality cannot meet all the requirements among these applications, and the quality influence final performances to a great extent. For instance, meshes with poor quality can lead to precision decline in finite element analysis, distortions in geometry surface modeling and bad effects even failure in graphics rendering. Such situations put forwards higher requirements for mesh quality. For this purpose, mesh quality improvement needs to be introduced and the research about this technique has also attracted a lot of attention. This paper presents a mesh smoothing method based on an exterior-angle-split (EAS) process to improve mesh quality. a1111111111 a1111111111 a1111111111 a1111111111 a1111111111 This paper is organized as follows. Section 2 reviews related previous work. Section 3 illustrates the rule of element transformation for triangular mesh and the global triangular mesh smoothing strategy. And properties of this algorithm are also discussed in this section. Section 4 applies this smoothing method to quadrilateral mesh and arbitrary polygon mesh. Results and discussions of some numerical tests are presented in Section 5. Section 6 presents the conclusion and future work of this paper.

Related work
Methods to improve mesh quality can be roughly divided into geometric methods, topological methods and combined methods. Geometric methods focus on nodes' geometric positions without changing connectivity of nodes. Among geometric methods, there are also two subclasses, optimization-based methods and smoothing methods. Relying on mathematical models, optimization-based methods [6][7][8][9] regard mesh quality as a form of energy. The quality can be characterized as an energy formula relying on mathematical models, so that the issue of quality improvement is converted into an optimization problem. Through minimizing the energy function, these methods can gain better global quality but usually a time-consuming work. While smoothing methods make a balance between mesh quality and computation efficiency, which is also what this paper concentrates on. Such easy-implement methods can meet the quality requirement of applications in most cases. Elements surrounding a node can be seen as a ring. An idea of smoothing is to relocate the central point based on its ring. The wellknown Laplacian smoothing [10,11] moves the central point to the average of its surrounding vertices' positions. On the basis of this there are also modifications [12][13][14] to avoid generating negative elements, known as constrained Laplacian smoothing. The outer outline of a ring is a polygon. Based on the polygon of a node's ring, angle-based smoothing [15] takes an iterative scheme to move the central point to the average of intersection points of adjacent interior angle bisectors of the polygon. It can obtain better angle distribution and maximize the minimum angle. Meanwhile, GETMe [16,17] starts another viewpoint of mesh smoothing based on element transformation. This viewpoint focuses on each element shape instead of the central point position of its ring. SSO [18] is also an element transformation method, which optimizes mesh quality with the assistance of the height of element.
Topological methods improve mesh quality through the ways of changing the connectivity of nodes as well as the number of nodes. Basic operations including edge split, edge collapse and edge swap can be usually seen in [19,20]. Small polyhedral reconnection [21] is proposed to eliminate poorly-shaped element, which is also can be used in surface mesh.
Combined methods have both geometric operations and topological operations. Centroidal Voronoi Tessellation is the concept of generating points at mass centroids of the corresponding Voronoi tessellations derived from a triangulation, which is applied to mesh optimization in [22][23][24]. Optimal Delaunay Triangulations [25,26] improve quality through minimizing a defined interpolation error. This paper presents a new geometric method based on exterior angle split to improve mesh quality. The proposed method achieves quality improvement mainly relying on a linear element transformation and an assembly strategy.

Smoothing algorithm for triangular mesh
This section presents an exterior-angle-split algorithm for triangular mesh smoothing. Element geometric transformation is to make each element equilateral. And the final result is obtained with a weighted assembly strategy. The regularity for element transformations can be proved on the theory level and numerical experiments illustrate its convergence after finite times of iterations. As for triangular mesh, this paper selects the ratio of two times incircle radius (2r) and circumcircle radius (R) of the triangle as quality criteria, denoted as γ.

Transformation of single element
EAS method regards elements unconnected in the first two stages, therefore, each element transformation is independent to each other. A transformation for one element does not influence its adjacent elements at all, this will be specified in Section 3.3.
A local Cartesian right-handed coordinate system is built to describe the transformation where z-axis coincides with outer normal vector of the element. Nodes of an element are recorded in a counterclockwise direction, denoted as p m , m�{0,1,2}. A counterclockwise directed edge is denoted as p mod(m−1,3) p m , and a clockwise directed edge can be denoted as p mod(m+1,3) p m . An interior angle can be denoted as α m , m�{0,1,2}, whose subscript is the same as the subscript of the corresponding node p m . And the superscript stands for the current stage. For example, p ð0Þ m is the original position and p ð1Þ m is the position after the first stage. The first stage EAS method consists of two substages, a counterclockwise transformation substage and a clockwise transformation substage. Take the counterclockwise substage as the first substage for instance as depicted in Fig 1A. First, exterior angles can be constructed by the extension lines of counterclockwise directed edges p ð0Þ modðmÀ 1;3Þ p ð0Þ m . For a clearer description, we define a virtual point q ð0Þ m at the infinite position of extension line, hence the extension line can be denoted as p ð0Þ m q ð0Þ m . Next, a parameter t, 0<t<1, is introduced to define a split line which splits each exterior angle into two angles. One of the angles denoted as b ð0Þ m is expressed as t � ffq ð0Þ m p ð0Þ m p ð0Þ modðmþ1;3Þ . According to geometric relations, b ð0Þ m can also be expressed as t � ðffa ð0Þ modðmÀ 1;3Þ þ ffa ð0Þ modðmþ1;3Þ Þ. Then, three intersection points of split lines can be calculated with arithmetic. Finally, endpoints p ð0Þ m of the directed edges are moved to their corresponding intersection points p ð1=2Þ m . The second substage is a reperform of the first substage but in clockwise direction as depicted in Fig 1B, where p ð1Þ m can be obtained. The successive two substages are in the opposite direction to each other to rotation offset. It doesn't matter that which direction is the first, and it is just customary to take the counterclockwise as the first one. Although an element can obtain a better shape after EAS transformation as shown in Fig 2, necessary shrinking operation which is the second stage of EAS method have to be taken to offset the scaling and displacement without shape change. The second stage selects either the ratio between the perimeter or the ratio between the square root of area before and after transformation as shrinking coefficient μ. Eq 1 specifies the shrinking operation with shape preserving and barycenter preserving, where bc (0) , bc (1) are the barycenter before and after the transformation of the first stage, p ð1Þ m ; p ð2Þ m are nodes' positions before and after shrinking. The second stage is to preserve the positive influence generated by the first stage and get rid of the negative influence. It drags the element back to its original barycenter and maintains its original size, meanwhile, a better shape is preserved. Fig 3 shows the shape-changing process of a triangular element of poor quality.

Regularity
Through observing changes of interior angles, the underlying laws can be discovered. We use the symbols the same as Fig 1 and Fig 2 to clarify the property of EAS transformation. The interior angle after the first stage of EAS method is denoted as a ð1Þ m . Since the second stage of EAS method does not change the shape as well as the interior angle, the value of a ð1Þ m is preserved after the second stage. As illustrated in Fig 1A,  is a According to the theorem that an exterior angle is equal to the sum of two nonadjacent interior angles in a triangle, ffp ð1=2Þ Clearly, the element in row i column j is equal to element in row j column i, making K a symmetric matrix. According to the theorem of spectral decomposition, K has three eigenvalues and three orthogonal eigenvectors. And its eigenvalues λ m , m�{0,1,2}, are 3t 2 −3t+1, 3t 2 −3t +1,1, respectively, which means K is a positive definite matrix. After performing the normalization of eigenvectors, Z m ¼ r m = ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi r m � r m p , the orthogonal matrix F = [η 0 η 1 η 2 ] consisting of normalized eigenvectors is obtained. The spectral decomposition of symmetric Matrix K can be expressed as Eq 6, where I n denotes a 3×3 matrix with the only non-zero coefficient in row n column n, and (I n ) n,n = 1.
As is known to all, the sum of interior angles of a triangle is 180˚. Eq 7 implies that a series of EAS transformation can average interior angles of an element to 60˚i.e. each element tends to be a regular triangle.

Assembly
In the first two stages, each element is treated unconnected, and transformations and shrinking operations are conducted independently for each isolated element, therefore, one element's operations will not deteriorate others' quality. By the end of the second stage, an integral mesh may be transformed into a set of isolated elements. Section 3.3 elaborates the assembly strategy for isolated elements to gain the final result which is also the third stage of EAS method.
A node with the valence of v on the original mesh means that it is included by its v adjacent elements. Naturally, after two stages of EAS method, v replicas of this original node exist its v unconnected included elements as illustrated in Fig 4. In order to determine the final position of nodes, a weighted strategy is selected to relocate the position of each node based on the positions of its replicas as well as the sizes and quality changes of its included elements, as represented in Eq 8.
In Eq 8, p ð3Þ i is the final position of node with a global index i, p ij is the position of its corresponding replica with local index j, s ij is the area of original element, and Δγ ij is the quality change compared to the original element with local index j. Through the third stage, isolated elements are assembled into an integral one again. Clearly, the time cost of the whole EAS method is linear, whose time complexity is O(n), where the n is the number of elements.
EAS method is a mutually independent work for each element, which provides access to parallel programming for simultaneous transformations and assembly on the global mesh.

Boundary treatments
EAS method takes boundary nodes into account. Fig 5 illustrates treatments for boundary nodes, which is either fixing boundary nodes at their original positions or merely allowing them to move along the direction of boundary edges.

Convergence
In terms of quality versus computational cost, the number of smoothing times have to be finite. Such is an ideal situation that mesh quality is improved and stabilizes at a higher level with the less increasing smoothing time. Although it is hard to obtain a rigorous mathematical proof for the convergence of this linear system, plenty of numerical experiments are carried out to testify the convergence of EAS algorithm. Through plenty of numerical tests, EAS can reach a high performance when t has the value of 0.5. Typical demonstrations of average quality and small angle percentage (0~20˚) changes are illustrated in Fig 6 where EAS method shows its performance and stability.

Extension to 3D surface
Given that most of the real-world models are closed surface meshes, EAS method can be extended to surface mesh smoothing. Similar to boundary edges and boundary nodes on planar meshes, there are some constrained edges and nodes called feature edges and feature nodes on closed surface meshes, respectively, which reflect some shape features such as sharpedges, corners, grooves and so on. Netgen [27] recognizes feature edges by the means of

PLOS ONE
Mesh smoothing algorithm based on exterior angles split

PLOS ONE
Mesh smoothing algorithm based on exterior angles split calculating the dihedral angle of two adjacent elements sharing one common edge, which is also used in [28].
EAS method adopts the same treatments for feature nodes as adopted for boundary nodes on planar meshes which is either fixing feature nodes at their original locations or merely allowing them to move along the direction of feature edges. For higher precision, a feature node is allowed to move along the local interpolation curve constructed by its adjacent nodes.
In order to preserve original geometric information, the original surface mesh is always recorded as M 0 before smoothing. Except for feature nodes, each node can obtain a mean normal based on its ring using the approach referred by [29], denoted as m(p i ) in Fig 7. After one time of smoothing, the final position of node p 0 i is obtained through the way of projecting p i derived from assembly back to the original mesh M 0 along m(p i ), which is also used in Mesquite [30]. On the other hand, if high precision is required, Coon patches [28] or least-squares method can help to construct a local parametric surface piecewisely based on adjacent elements of the original mesh M 0 . For the nodes updated after one time of smoothing, final positions p 0 i can be obtained by projecting p i back to the local parametric surface along m(p i ) to proceed subsequent iterations.

Smoothing algorithm for quadrilateral mesh
EAS smoothing algorithm can be also applied to quadrilateral mesh smoothing. As is similar to the case of triangular mesh, EAS transformation is to make each element equilateral and the final result is obtained with a weighted assembly strategy. The regularity for element transformations can be proved on the theory level and numerical experiments illustrate its

PLOS ONE
Similarly, interior angles at different stage are denoted as a ð0Þ m ; a ð1=2Þ m ; a ð1Þ m ; m�f0; 1; 2; 3g, whose subscript is the same as the subscript of the corresponding node and the superscript stands for the current stage. Interior angles after the first substage of the first stage can be obtained as Eq 9 with the corresponding transition matrix K + in Eq 10. fa ð1Þ g ¼ Kfa ð0Þ g ¼ Eq 11 shows the transition matrix for the first stage. Through spectral decomposition eigenvalues are computed as 4t 2 −4t+1, 2t 2 −2t+1, 2t 2 −2t+1, 1, respectively. Like the case of triangular mesh, λ 3 is the only dominant item for EAS transformation with lim l!1 λ 0 l = λ 1 l = λ 2 l = 0.
After eigenvector normalization, a series of EAS transformation can average interior angles for quadrilateral mesh from a theoretical point of view as shown in Eq 12.
Meanwhile, the second and third stages of EAS method for quadrilateral mesh can be

Smoothing algorithm for arbitrary polygonal mesh
Considering arbitrary polygonal mesh like pentagonal mesh, EAS smoothing algorithm is also able to achieved as shown in Fig 8B. In addition, EAS method has natural applications for mixed polygonal mesh with few adjustments.

Results and discussions
This section illustrates the performance of EAS smoothing algorithm described in the previous sections. Several examples from different areas are selected for demonstration and the results

PLOS ONE
are discussed. Taking generality and practicality into account, tests are only implemented on models with triangular elements and quadrilateral elements. Since EAS method is based on element geometric transformation, this section compares EAS with other existing state-of-theart element-transformation-based methods which have influence and is widely recognized in the field of mesh smoothing. And the optimization-based method is also chosen as a comparative item because it is usually time-consuming but recognized to reach an excellent mesh quality. Smoothing methods are programmed by C++ aided by OpenMesh [31] and the optimization-based method is implemented by Mesquite. Figures are screenshot by visualizer with a STL files and VTK files. This paper selects the ratio of two times incircle radius (2r) and circumcircle radius (R) of the triangle as quality criteria denoted as γ for triangular mesh. And for quadrilateral mesh, as proposed in Gambit, the preprocessor of Fluent [32], the angle between links of midpoints of opposite edges is regarded as quality criteria γ to measure the skewness. And there are more quality criteria in [33][34][35]. Quality criteria γ takes the value of 1 for regular element and 0 for degenerated element after normalization. Changes of node's position may produce negative elements with inversion, which is not the concern of this paper. To avoid the problem of negative elements, each movement is accepted only if no negative element emerges, or this movement is rejected. Such constraint operation is introduced to the tests for all smoothing methods mentioned in this section. To evaluate the global quality of a mesh model, this section selects average quality, minimum quality, minimum internal angle and small angles percentage (0~20˚) as measurements. All of mentioned smoothing loops execute ten times of iterations.
The first example is a planar triangular mesh with 900 nodes and 1598 elements. Several smoothing methods are taken to improve its quality and results are listed in Table 1. It can be seen that the proposed EAS method has advantages in average quality and small angles percentage. Meshes before EAS smoothing and after EAS smoothing are demonstrated in Fig 11. The second example is a practice of quadrilateral surface mesh, whose results are listed in Table 1. Meshes before EAS smoothing and after EAS smoothing are demonstrated in Fig 12. Example 3 is a model of spacecraft with 4348 nodes and 8413 triangular elements. The initial quality is poor and results are presented in Table 1. Meshes before EAS smoothing and after EAS smoothing are shown in Fig 13. Clearly, as for element transformation-based methods, EAS method still has advantages in average quality compared with other methods and GETme has better performance in minimum angle and angle distribution.
The fourth example is a triangular mesh of shoulder blade with 589 nodes and 1174 elements, as shown in Fig 14A. The smoothed mesh is presented in Fig 14B. As listed in Table 1, smoothing operations have improved mesh quality remarkably in different measurements. Among element-transformation-based methods, EAS method can obtain better average quality while SSO has better angle distribution.
The fifth example is a Stanford bunny of triangular mesh containing 43199 nodes and 86394 elements. EAS method improves the average quality and reduces small angles percentage significantly as presented in Table 1 and Fig 15. Through zooming in, it can be seen that EAS exerts positive effects for mesh quality improvement noticeably in Fig 15. And Fig 16 illustrates the distribution of element interior angles and element qualities.
Example 6 is a scanned heart model of triangular mesh with 34505 nodes and 70000 elements. The smoothed mesh and a part of the mesh before and after smoothing are illustrated in Fig 17, and the results are listed in Table 2. The last example is a spline of quadrilateral-dominant mesh as shown in Fig 18, composed of 26767 nodes and 26830 elements. Results before and after EAS smoothing are listed in Table 2.

Conclusions and future works
This paper presents an algorithm for 2D planar and surface mesh smoothing using element geometric transformation. This algorithm conducts geometric transformation independently for each element based on an exterior-angle-split process, and then assemble them into an integral mesh of high quality with a weighted strategy.
The presented method has satisfactory outcomes. From a theoretical point of view, it can be proved that EAS method can make an element equilateral. And plenty of experiments are carried out to testify the convergence of EAS algorithm from a global perspective. With simple and appropriate adjustments, this method is capable of being extended to surface meshes of arbitrary polygonal element. And numerical experiments show that EAS method has effectiveness in mesh smoothing evidently. Conclusively, the proposed EAS smoothing algorithm has very simple transformation rules and leads to significant improvement for mesh quality. EAS method is a novel geometric smoothing method, and combination with topological methods is where one of the potentials resides in the future. On the other hand, EAS method introduced in this paper remains at the level of polygonal mesh. Whereas, polyhedral mesh also has vital applications especially in industrial simulation, which puts forwards new goals and requirements for this method. Further research will naturally focus on its combination with topological methods and extension to polyhedral mesh.