Optimization of Mesh Quality with Preserving Surface Characteristics

In this paper, we present a new technique called Trapezium Drawing to improve surface mesh quality while maintaining essential surface characteristics. In contrast to previous methods we do not tend to preserve mesh vertices on the original discrete surface. Instead our approach allows to keep a resulting mesh close to the surface approximated by the initial mesh. All operations are performed directly on the surface. As a result our technique is robust and runs at interactive speeds. It can be applied to triangular and quadrilateral meshes iteratively. Various quantitative measures are presented to demonstrate effectiveness of proposed technique.

Key words: Surface Mesh Quality, Surface Characteristics, Triangles, Quadrilaterals

1. Introduction

Improvement of surface mesh quality is important problem for numerical simulations, solid mesh generation and computer graphics applications.

There are two main ways for mesh optimization: modifications of mesh topology (remeshing) by inserting/deleting mesh nodes or edge flipping [1][2][3][4] and node movement methods commonly called mesh smoothing. Let us note that changing mesh topology may cause some problems in simulations requiring solution transfer from the original mesh to the resulting one. Moreover many remeshing algorithms can be applied only to triangular meshes. Therefore in this paper we focuses on node movement technique for surface mesh improvement.

To improve mesh quality in the plane a number of smoothing techniques have been developed ranging from simple Laplacian smoothing (see, for instance, [5] and references therein) to more sophisticated algorithms. Among them there are physically based methods [6][7] where nodes are moved under the influence of some forces so that the shape of incident elements is improved. Instead of local mesh optimization by moving each node on the basis some geometric characteristics (as is done in Laplacian smoothing, angle-based [8] and physically based methods) the optimization-based techniques allow to improve all original mesh. In these techniques so called cost function [9] is optimized. As such function aspect ratio [10] or distortion metrics [11] [12] can be used.

It is necessary to note that the good shape of mesh elements is not the only criteria for mesh quality when surface meshes are considered. It is also essential to minimize changes in surface characteristics like normals and curvatures. As it has been pointed out in [13] preservation of such characteristics is important for preventing drastic changes in the volume enclosed by the surfaces and in forces like surface tension that depend on surface properties

1.1 Previous works

Several techniques to improve surface mesh quality have been developed over the last decade. Most of them are based on an idea to constrain node movement to the underlying discrete surface. A simple way is to reposition each node in a locally derived tangent plane and project it back to the surface [14][15]. More robust algorithms use global parameterization of the original mesh, and then improvement in the parameter domain [16][17][18]. One drawback of these methods is high computational cost since they involve the solution of a large set of equations. Moreover, global parameterization may distort the complicated 3D structure. The alternative to global parameterization has been proposed in [13]. The nodes of the mesh are moved in a series of local parametric spaces derived from individual mesh elements.

Let us note, however, that all these methods allow to keep mesh nodes on the original mesh but not on the surface approximated by this mesh. As a simple example consider a sphere and a mesh with the nodes situated on this sphere. Applying algorithms described above we will obtain new nodes situated on the original mesh but not on the original sphere. Therefore, unlike initial mesh resulting mesh will not be discrete approximation of the original sphere. Furthermore, several iterations will cause considerable shrinking the initial surface.

1.2 Contribution and overview

We looked at the problem of surface mesh improvement from another point of view. We do not pose a problem to preserve new nodes on the original mesh that is on the discrete surface. Instead we propose method called Trapezium Drawing (TD) to keep resulting mesh very close to the surface approximated by the initial mesh. This method can be called explicit because all operations are performed directly on the surface. That is the reason why our technique is robust and fast. But for all that we do not sacrifice any quality in results. Moreover the method can be applied iteratively. Even after several passes there is no degenerating mesh quality and all main features of the surface are preserved.

Fig. 1 demonstrates the main difference between our approach to the problem and previous ones.

Fig. 1. (a) Original pyramid; (b) Optimized pyramid with the apex vertex constrained to the original surface; (c) Pyramid optimized with TD algorithm

The rest of the paper is organized as follows. In Sections 2 and 3 we give a detailed description of TD algorithm for triangular and quadrilateral meshes respectively. Section 4 describes some quantitative measures to examine mesh quality and deviation of the resulting surface from the original one. Examples of applying TD algorithm to various meshes are presented in Section 5. We close by offering some concluding remarks in Section 6.

2. Algorithm for triangular meshes

Let us consider some node of the original mesh and all nodes, associated with this node. For each node the new position of the nodeis obtained according to the following procedure.

At first, we find the node such as vectors and compose a maximal angle as shown in Fig. 2. The new position of the nodewith regard to will be a vertex of a trapezium such as a triangle consisting of the points is isosceles (Fig. 3).

Fig. 2. Search for the node such as vectors and compose a maximal angle

Fig. 3. Search for the new location of the node with regard to the node

Let us denote by the midpoint of the segment. The new coordinates of may be found using following formulas (see Fig. 4):

After all have been found the new position of the nodeis obtained by averaging coordinates of,.

Fig. 4. Search for the coordinates of

2.1 Mathematical justification of the algorithm

Let us emphasize two aspects of the algorithm. We tend to make all edges adjacent to the considered node equal that provides improvement of mesh element quality. On the other hand, we preserve new node at the same height as its old position from the lines passing through the adjacent nodes. This ensures that resulting mesh will be very close to the surface approximated by the original mesh.

It is easy to prove that new position of the node always stays within a distance of from the original mesh, where is the length of the longest edge of elements adjacent to the considered node.

Let us make an important remark. We update positions of mesh nodes only after completion of iteration. Actually, many heuristic approaches update position of the node right after it has been found and use this new position to move next node. Such procedure may cause some problems for our algorithm. Consider mesh node and all its neighboring nodes , . Suppose positions of some nodes have been already changed. That means these nodes deviate a little from the original surface. It is clear that the error will accumulate as we use the nodes to find new positions of other nodes.

2.2 Quality control

Using any optimization technique we must be sure that there will be no invalid or badly shaped elements in the resulting mesh. The simplest way is to check whether the quality of mesh elements improved after applying technique or not. It is common to use for that minimal angle or signed minimal angle like it is done in smart Laplacian smoothing [11]. Let us note that using minimal angle as a quality measure reduces the risk of obtaining inverted elements but still cannot guarantee validity of the new mesh. To solve this problem it is necessary to introduce so called signed aspect ratio, which allow to detect inverted elements. In our work we propose the following measure, where are lengths of triangle sides and is the signed area of triangle (i.e., with respect to counter-clockwise orientation). For oriented mesh all valid triangles have positive areas and all inverted triangles have negative areas. It is obvious that for any triangle . A value 1 corresponds to an equilateral triangle, while –1 indicates an inverted equilateral triangle. When all three vertexes of the triangle are co-linear, the triangle is degenerate that yields a value of 0.

After the new position of the node has been found we need to calculate minimal aspect ratio for the triangles adjacent to this node and compare it with the minimal aspect ratio with regard to the old position of the node. If we move the node to the new position. Otherwise we keep the node at its initial position. Let us note that computational cost of calculating signed aspect ratio is the same as computational cost of calculating minimal angle. But unlike minimal angle signed aspect ratio guarantees validity of the resulting mesh.

3. Algorithm for quadrilateral meshes

As it has been described above for triangular meshes TD algorithm uses the neighboring nodes connected with the central node. To apply this algorithm to quadrilateral meshes we simply need to consider all surrounding nodes (Fig. 5a).

Fig. 5. (a) TD algorithm for quadrilateral meshes; (b) Quadrangle adjacent to the node; are the areas of the triangles; are the lengths of quadrangle diagonals

3.1 Quality control

Just as in the case of triangular meshes we need some measure to guarantee validity of the resulting mesh. It seems a good idea to decompose quadrangle in the pair of triangles and estimate the quality of these triangles. Frey in [19] proposed to use following quality measure:, where are the lengths of quadrangle sides, and, , are the areas of triangles (Figure 5b). In [20] it has been proved that this measure has all desired properties for a quadrangle quality measure: extremal and asymptotic. But let us remark that we need some addition checks to detect invalid, e.g. self-intersected, elements (see [19] for details). Therefore we propose slight modification of this measure that allows to effectively discover inverted elements. , where,, are the signed areas of triangles. It can easily be checked that for any quadrangle, in the case of square, in the case of degenerate element, and in the case of inverted element.

Similarly to the case of triangular meshes after the new location of the nodehas been found we calculate minimal aspect ratio for the quadrangles adjacent to this node and compare it with the minimal aspect ratio for the old node position. If we move the node to the new position. Otherwise we keep the node at its initial position.

4. Quantitative measures

We use various quantitative measures to demonstrate that after applying TD algorithm mesh quality is improved and the deviation of the resulting mesh from the original surface is small.

4.1 Signed aspect ratio

To measure the geometric properties of the obtained mesh we use signed aspect ratios described in Sections 2.2 and 3.1. Let us remind that the closer the value of to 1, the closer element to equilateral one. corresponds inverted elements while indicates degenerate elements. The histograms of aspect ratio for original and resulting meshes are presented to demonstrate the improvement of mesh quality.

4.2 Change of the volume

Having a triangular or quadrilateral mesh it is easy to calculate the interior volume. This can be done by summing the volumes of all oriented pyramids centered at a point in space (the origin, for example) and with a triangle or quadrangle of the mesh as a base. We computed the volumes enclosed by the original and resulting surfaces and calculated the difference between them.

4.3 Changes in surface normals

For numerical simulation it is very important to minimize changes in the surface normals while improving mesh quality. We calculated differences between normals at the nodes of the resulting mesh and corresponding nodes of the original mesh to show that these changes are small.

Fig. 6. Dihedral angle used for approximation of mean curvature

4.4 Changes in surface curvatures

It is also important to keep changes in surface curvatures small. Roughly speaking, the degree of surface “bending” in space is defined by mean curvatures. Therefore we use notion of mean curvature to prove that there are no large deviations between curvatures of the original and resulting surfaces. It is easy to see that polyhedral surface “bends” along its edges.

Denote by the dihedral angle between two faces adjacent to the edge of the mesh (Fig. 6). We calculate this angle in such a way that if is convex edge, if is concave edge, and if is plane edge.

The difference in dihedral angles of the original and resulting meshes is calculated to demonstrate the change of the surface shape.

5. Experimental results

In order to show the effectiveness of the proposed algorithm consider two aspects of the problem of surface mesh improvement. Firstly let us demonstrate the visual effect of applying TD algorithm. From the results shown in Fig. 7, we can see that even after fifth iteration the model of Mannequin preserves all its characteristic features. The statistics for the mesh of Mannequin after first and fifth passes are presented in Tables 1 and 2. We calculated maximum and average changes in dihedral angles (,(degree)), maximum and average changes in normals (, (degree)), and difference of the volume (). We can see that average changes in normals and dihedral angles are small. Let us note also that without any special constraints [21] TD algorithm gives very good results in terms of volume preserving. Fig. 8 demonstrates that TD algorithm does not destroy special features of the model such as anisotropy and local refinement. We can see that local refinement near wolf’s mouse and eyes is preserved even after fifth iteration. All statistics for the mesh of Wolf can be found in Tables 3 and 4.

Fig. 7. (a) Mesh of Mannequin; (b) The mesh optimized with TD algorithm after 5^th iteration

1^st pass

5^th pass

Table 1. Statistics on the model of Mannequin

Initial

1^st pass

5^th pass

Table 2. Histograms of signed aspect ratio for the mesh of Mannequin (percentage)

Fig. 8. (a) Model of Wolf; (b) The model processed with TD algorithm after 5^th iteration; (c) Original mesh of Wolf; (d) The mesh optimized with TD algorithm after 5^th iteration

1^st pass

5^th pass

Table 3. Statistics on the model of Wolf

Initial

1^st pass

5^th pass

Table 4. Histograms of signed aspect ratio for the mesh of Wolf

Now let us demonstrate more precisely how TD algorithm influences on the shape and size of mesh elements.

One important factor affecting the accuracy and efficiency of the finite element analysis is size quality. There are two cases that must be considered: uniform and graded meshes. For uniform meshes it is necessary to keep difference among sizes of the elements as small as possible. Graded meshes contain several layers that should not be diffused after applying some technique and size change at the transitional zone should be smooth. In Fig. 9 we can see that even after several iterations there is no “diffusion” of layers, size change is very smooth and uniform zones remain uniform.

Fig. 9: (a) Quadrilateral mesh of Rocker Arm; (b) The mesh optimized with TD algorithm after 5^th iteration

As is well known Laplacian smoothing may degrade mesh quality if the algorithm iterates more than a few times. On the contrary, as it can be seen from Table 5 TD algorithm tends to improve the aspect ratio as more passes are performed.

Initial

1^st pass

5^th pass

Table 5. Histograms of signed aspect ratio for quadrilateral mesh of Rocker Arm

1^st pass

5^th pass

Table 6. Statistics on the quadrilateral mesh of Rocker Arm

The rest of the statistics for quadrilateral mesh of Rocker Arm is presented in Tables 6.

5.1 About applying TD to plane meshes

Although TD algorithm is meant for working with surface meshes it also gives good results for plane meshes. Actually in many cases, the results are even superior to the mostly used techniques.

As is well known Laplacian and smart Laplacian smoothing are the most commonly used techniques to improve the quality of plane meshes because of their low computational cost. But Laplacian smoothing may produce badly shaped or even inverted elements while smart Laplacian tends to “diffuse” layers of graded meshes. In [8] there has been proposed effective method called angle-based approach with low computational cost. Quality of a mesh optimized with this method is much better than after Laplacian and smart Laplacian smoothing and chance to obtain inverted elements is reduced. But in [22] it has been pointed out that this is true mostly for the meshes with regular connectivity. When mesh contains very distorted elements angle-based approach also may fail. To solve this problem Surazsky and Gostman in [22] proposed to use weights. This scheme reduces the risk of generating inverted elements but still cannot guarantee that the obtained mesh will be valid. In fact, both angle-based and weighted angle-based approaches may produce the meshes with even worse quality than Laplacian smoothing.

Now let us give some interesting example of applying TD algorithm. As is well known, data can be obtained from different sources of information, in particular from existing contour maps; an example is given in Figure 10a. Despite a flurry of activity in the generation and visualization of terrain data from scattered data points this matter remains a difficult and computationally expensive problem. Traditionally, visualization of elevation data consists of three major steps: tessellation, modeling and rendering. The most commonly used tessellation is Delaunay triangulation. Unfortunately, this method has some serious drawbacks resulting in a confusing image because clearly visible “traces” of triangulation can be observed as it can be seen in Figure 10a. To attenuate this effect, we process the mesh of Bandai Mountain, Japan, with all foregoing techniques and TD algorithm. The results are following: Laplacian smoothing produced 211 inverted triangles; smart Laplacian – no inverted elements; angle-based approach – 3673; weighted angle-based approach – 1449; TD algorithm – no inverted triangles during all performed iterations.

As it has been mentioned in Section 2.1 we update positions of mesh nodes only after completion of iteration. That prevents accumulation of the error. However, for plane meshes we can change node position right after it has been found that allows to guarantee validity of the resulting mesh.

Let us note that since signed aspect ratio is used in smart Laplacian smoothing and TD algorithm both techniques ensure that there will be no inverted elements not only after first iteration but also after several iterations. However, these techniques give different results. Figure 10bc demonstrates that smart Laplacian diffuses layers of the original graded mesh while TD algorithm tends to preserve these layers.

5.2 Comparison of Laplacian smoothing, TD and BDS algorithms

To compare our approach with the previous ones we chose simple method illustrated the main idea of the preceding works. At each node of the original mesh we define local tangent plane, project neighboring nodes to it and find new location of the concerned node using smart Laplacian smoothing. After that we project obtained node back to the discrete surface. To be precise, let us call this approach “Back to the Discrete Surface” (BDS).

Fig. 10. (a) Fragment of the initial mesh of Bandai Mountain; (b) The mesh optimized with TD algorithm after 1^st iteration; (c) The mesh optimized with smart Laplacian smoothing after 1^st iteration

We applied usual Laplacian smoothing, BDS and TD algorithms to the triangular mesh of the cubic surface. Since we have analytic surface curvatures and normals can be easily calculated at each point of the surface. Curvatures for resulting meshes can be computed numerically [23]. We calculated maximum and average changes in Gaussian and Mean curvatures (,,,), maximum and average changes in normals (, (degree)), maximum and average deviations from the original surface (,). To measure the geometric quality of the resulting meshes we again used signed aspect ratio described in Section 2.2.

Laplacian smoothing and BDS approach seem to improve shapes of mesh elements quicker than TD algorithm (Table 7). However, the data from the Table 8 verify that TD algorithm gives the best results in terms of preserving surface characteristics.

Initial

Laplacian smoothing

BDS

TD

Table 7. Histograms of signed aspect ratio for the mesh of cubic surface processed with Laplacian smoothing, BDS approach and TD algorithm

Laplacian smoothing

BDS

TD

8.48812

4.2132

3.00218

0.566505

0.504979

0.395472

2.007724

1.434133

1.075598

0.144465

0.124508

0.095078

4.89858

4.73702

3.46709

1.051849

0.809531

0.656461

0.007234

0.00361

0.00181

0.001524

0.000491

0.000326

Table 8. Statistics for the mesh of cubic surface processed with Laplacian smoothing, BDS approach and TD algorithm

6. Conclusion

In this paper, we introduced novel approach to improve surface mesh quality while maintaining essential surface characteristics. In contrast to existing methods we did not tend to keep mesh vertices on the original discrete surface. All operations are performed directly on the surface. It allows us to take an advantage in speed without sacrificing any quality in results. Moreover as it has been shown the results are even superior to the previous ones in terms of preserving surface curvatures and normals.

The procedure has been successfully tested on a number of complex triangular and quadrilateral surfaces and plane meshes. Various quantitative measures proved that proposed technique do not cause considerable changes in surface characteristics while improving mesh quality. The algorithm can be applied iteratively that allows the user to attain resulting mesh more suitable for his application.

TD algorithm can be considered as universal approach since it can be applied for any surface or plane triangular or quadrilateral mesh in a consistent manner.

An interesting and important problem is to define the number of iterations. While geometric mesh quality is improved the more passes are performed each iteration causes more damage to the original surface. Therefore we need to find some compromise between geometric quality improvement and degree of deviation from the original surface. It seems an interesting idea to construct some error-metrics (for instance, normal-based and curvature-based) for measuring the deviation between original and resulting meshes. We can define some threshold values for these metric and stop iterations when these values are reached. We leave careful study of this problem to future work.

7. Acknowledgments

Authors would like to thank Professor A. G. Kushnirenko (Scientific Institute of System Analysis of the Russian Academy of Sciences, Russia) for his support and encouragement, which were very important during our work at this paper.

References

[1] H. Hoppe, T. Duchamp, J. McDonald, W. Stuetzle, Mesh Optimization, Computer Graphics (SIGGRAPH 93), pp. 19-26, 1993.

[2] H. L. de Cougny, Refinement and Coarsening of Surface Meshes, Engineering with Computers, 14(3):214, 1998.

[3] L. Alboul, R. van Damme, Polyhedral metrics in surface reconstruction: tight triangulations, The Mathematics of Surfaces VII, T. Goodman and R.Martin (eds.), Clarendon Press, Oxford, pp. 309-336, 1997.

[4] N. Dyn, K. Hormann, S.-J. Kim, D. Levin, Optimizing 3D Triangulations Using Discrete Curvature Analysis, Mathematical Methods for Curves and Surfaces, Vanderbilt Univ. Press Innovations in Applied Mathematics Series, pp. 135-146, 2001.

[5] D. A. Field, Laplacian Smoothing and Delauney Triangulations, J. Communications in Applied Numerical Methods, Vol. 4 pp. 709-712, 1998.

[6] I. Babushka, O. C. Zienkiewicz, J. Gago, and E. R. de A. Oliviera (eds.), Accuracy Estimates and Adaptive Refinement in Finite Element Computations, John Wiley & Sons, Chichester, pp. 281-297, 1986.

[7] K. Shimada, D. C. Gossard, Bubbled Mesh: Automated Triangular Meshing of Non-manifold Geometry by Sphere Packing, Proceedings of the ACM Third Symposium on Solid Modeling and Applications, pp. 409-419, 1995.

[8] T. Zhou, K. Shimada, An Angle-Based Approach to Two-Dimensional Mesh Smoothing, Proceedings of the 9^th International Meshing Roundtable, pp. 373-384, 2000.

[9] L. A. Freitag, On Combining Laplacian and Optimization-Based Mesh Smoothing Techniques, AMD Trends in Unstructured Mesh Generation, Vol. 220 pp. 37-43, 1997.

[10] V. Parthasarathy, S. Kodiyalam, A Constrained Optimization Approach to Finite Element Mesh Smoothing, J. Finite Element Analysis and Design, Vol. 9 pp. 309-320, 1991.

[11] S. A. Canann, J. R. Tristano, M. L. Staten, An Approach to Combined Laplacian and Optimization-Based Smoothing for Triangular, Quadrilateral, and Quad-Dominant Meshes, Proceedings of the 7^th International Meshing Roundtable, pp. 479-494, 1998.

[12] O. P. Jacquotte, G. Goussement, Structured Mesh Adaptation: Space Accuracy and Interpolation Methods, Computer Methods in Applied Mechanics and Engineering, Vol. 101 pp. 397-432, 1992.

[13] R. V. Garimella, M. J. Shashkov, P. M. Knupp, Optimization of Surface Mesh Quality Using Local Parameterization, Proceedings of the 11^th International Meshing Roundtable, pp. 41-52, 2002.

[14] P. Frey, H. Borouchaki, Geometric Surface Mesh Optimization, Computing and Visualization in Science, pp. 113-121, 1998.

[15] P. M. Knupp, Achieving Finite Element Mesh Quality via Optimization of the Jacobian Matrix Norm and Associated Quantities, International Journal for Numerical Methods in Engineering, Vol. 48 pp. 401-420, 2000.

[16] M. S. Floater, Parameterization and Smooth Approximation of Surface Triangulation, Computer Aided Geometric Design, 14:231-250,1997.

[17] K. Hormann, U. Labsik, G. Greiner, Remeshing Triangulated Surfaces with Optimal Parameterizations, Computer Aided Design, 33(11):779-788, 2001.

[18] M. Desbrun, M. Meyer, P. Alliez, Intrinsic Parameterizations of Surface Meshes, Proceedings of Eurographics 2002 Conference

[19] P. J. Frey, P. L. George, Mesh Generation, Hermes Science Publishing, Oxford & Paris, 2000.

[20] P. P. Pebay, Planar Quadrangle Quality Measures: Is There Really a Choice, Proceedings of the 11^th International Meshing Roundtable, pp. 53-62, 2002.

[21] A. Kuprat, A. Khamayseh, D. Goerge, L. Larkey, Volume Conserving Smoothing for Piecewise Linear Curves, Surfaces, and Triple Lines, Journal of Computational Physics, 172, pp. 99-118, 2001.

[22] V. Surazhsky, C. Gostman, High Quality Compatible Triangulations, Proceedings of the 11^th International Meshing Roundtable, pp. 183-192, 2002.

[23] G. Farin, Curves and Surfaces for Computer Aided Design, 4^th edition, Academic Press, Boston, 1997.

Laplacian smoothing	BDS	TD
8.48812	4.2132	3.00218
0.566505	0.504979	0.395472
2.007724	1.434133	1.075598
0.144465	0.124508	0.095078
4.89858	4.73702	3.46709
1.051849	0.809531	0.656461
0.007234	0.00361	0.00181
0.001524	0.000491	0.000326