Quality improvement of surface triangular mesh using a modified Laplacian smoothing approach avoiding intersection

We present a systematic procedure to improve the qualities of triangular molecular surface meshes and at the same time preserve the manifoldness. The procedure utilizes an algorithm to remove redundant points having three or four valences and another algorithm to smooth the mesh using a modified version of Laplacian method without causing intersecting triangles. This approach can be effectively applied to any manifold surface meshes with arbitrary complex geometry. In this paper, the tested meshes are biomolecular surface meshes exhibiting typically highly irregular geometry. The results show that the qualities of the surface meshes are greatly improved and the manifoldness of the surface meshes are preserved. Compared with the original meshes, these improved molecular surface meshes can be directly applied to boundary element simulations and generation of body-fitted volume meshes using Tetgen. The procedure has been incorporated into our triangular molecular surface mesh generator, TMSmesh 2.0. It can be also used as a standalone program and works together with any other surface triangular mesh generator to obtain qualified manifold mesh. The package is downloadable at https://doi.org/10.6084/m9.figshare.5346169.v1 and can be run online at http://www.xyzgate.com.


Introduction
Surface mesh generation arises in many applications, such as numerical simulation, computer visualization and geometry processing. Most computational applications involve triangulation of a complex surface geometry, especially in computational biology. Molecular surface plays an important role in computational biology, such as protein folding, structure prediction, docking and implicit solvent modeling. Recent developments in realistic mathematical modeling and numerical simulation of biomolecular systems raise new demands for qualified, stable, and efficient surface meshing, especially in implicit-solvent modeling [1]. Many triangulated meshes are generated by scanning devices or by isosurfacing implicit representations. However, it is not easy to generate high-quality mesh for complex surface geometry in such a1111111111 a1111111111 a1111111111 a1111111111 a1111111111 processes, especially if automated. Low-quality triangular meshes can undermine the order of accuracy or even cause non-convergence in numerical computations. In this case, it is desirable to build high-quality meshes from those low-quality meshes before performing any numerical simulation. [2] Mesh quality refers to the faithfulness, manifoldness and uniformness. Faithfulness is measured by how accurately the surface mesh preserves the original geometry and topology, such as surface area, volume and curvature of the referenced surface. Faithfulness related to the accuracy of the numerical simulation and geometry processing. Manifoldness of a surface mesh means that each point on the surface has a neighborhood which is homeomorphic to a disk in a real plane. Non-manifold surface mesh brings difficulties in generating corresponding surface conforming volume mesh that is an essential prerequisite for successful finite element method computations. One of the typical non-manifold errors in surface meshes is intersecting triangles. It is an important issue that preserving the manifoldness of the original surface meshes in the process of improvements of mesh qualities. Uniformness includes the triangle shape, regularity, complexity and so on. The quality of triangles in surface mesh is crucial for robustness of boundary element method (BEM) and finite element method (FEM) computations. In some numerical analysis, such as finite element analysis, regular or qualified mesh is always required. The mesh size is related to the complexity of the mesh. Generally, more sampling points capture more surface details. But more sampling points increase the complexity of the mesh, which increases the difficulty of numerical simulation. Therefore, meshes can be improved with respect to any number of quality metrics including shape, size, manifoldness, solution error or combinations of these. The improvement process is called remeshing. [3] There are two fundamental remeshing approaches, parameterization techniques [4][5][6][7][8][9][10] and mesh adaptation strategies [11][12][13][14][15]. In the context of the parameterizationbased remeshing, the initial surface is parameterized onto a surface in 2D, the 2D surface is meshed by the standard mesh generation method and the new mesh is projected back to the original surface. The existing parameterization method is classified to linear method [16][17][18][19][20][21], non-linear method [22][23][24] and hybrid method [25,26]. The parameterization-based methods can work even for coarse resolution remeshing, but the computation is expensive and the result is sensitive to the specific parameterization. The mesh adaptation method treats the surface as a set of points in space and operate directly on the surface. The strategies use local mesh modifications to improve the quality of the input surface mesh or adapt the mesh to a given mesh size criterion. Typical techniques for local mesh modifications combine vertex smoothing, vertex insertion/deletion, vertex shift, edge split, edge flip, edge collapse and so on. [10,[27][28][29][30] Recently, we have developed a method, TMSmesh [31,32] and its improved version, TMSmesh 2.0 [33], to generate molecular surface mesh for large biomolecules. TMSmesh 2.0 succeeds to generate manifold surface meshes for biomolecules comprised of more than one million atoms efficiently. The generated surface mesh preserves the original molecular surface features and properties (topology, surface area, enclosed volume, local curvature and so on) [34]. However, the mesh quality in the aspects of uniformness and smoothness are not guaranteed. There are some redundant points and some singular triangles in the meshes generated by TMSmesh 2.0. In boundary element and finite element simulation for some large biomolecules, the uniformness of the surface mesh is required to be improved to make the numerical simulation more robust and accurate. The faithfulness of the surface mesh to the original geometric features are critical in getting reliable simulation results from BEM and FEM computation in continuum molecular simulations. The parameterization-based methods contain the process of resampling, so these automated operation can not guarantee the fidelity of improved meshes to the original surface. And the complex geometry and topology of the large biomolecules will increase the difficulty of parameterization dramatically. There are some packages using mesh adaption strategies to improve mesh qualities, but few of them can preserve manifold meshes for complex surfaces. Therefore, in this work, we develop a package named SMOPT to improve the mesh quality directly on the surface and perform a series of local modifications on the mesh. The qualities of the new resulted meshes are improved greatly over the original surface meshes and avoid causing non-manifold errors, such as intersecting triangles. The program SMOPT has been incorporated into TMSmesh 2.0. It can be also used as a standalone program and works together with any other surface triangular mesh generator to obtain qualified manifold mesh.
In the next Section, the method used to improve the surface mesh quality is introduced. Some examples, analysis and applications are presented in the results section. The final section, Conclusion, gives some concluding remarks.

Materials and methods
In this section, we describe the mesh generator TMSmesh 2.0 and the improvement procedure SMOPT to improve the quality of the triangular mesh generated by TMSmesh 2.0.

TMSmesh 2.0
TMSmesh 2.0 is an algorithm for manifold triangular meshing of Gaussian molecular surfaces. [31][32][33] The Gaussian surface is defined as a level set of the summation of Gaussian kernel functions: x i and r i are the location and radius of the ith atom, N is the number of atoms in the molecule. d is the decay rate of the Gaussian kernel. When d decreases, the molecules will show more geometrical details. c is the isovalue and it controls the volume enclosed by the Gaussian surface [34]. The algorithm of TMSmesh 2.0 contains two stages, the first stage is an adaptive estimation and division process. The Gaussian surface is approximated by piecewise trilinear surface within controllable error. The second stage is to partition each piece of trilinear surface into single-valued patches along x, y, z directions by tracing along the fold curves. Then each single-valued patch is triangulated by the ear clipping algorithm. TMSmesh 2.0 succeeds to generate surface meshes for biomolecules comprised of more than one million atoms, and all the generated meshes are manifold mesh preserving the original detailed geometry of molecular surface. The sampled points in TMSmesh 2.0 are distributed according to the curvature and geometric features of molecular surface adaptively. In the regions with dramatic changing shapes, the density of the sampled points is high. It may result in some very close points in the mesh, which will bring difficulties to generate body-fitted volume mesh. And the very close points may cause singularity in numerical simulations for continuum modeling of biomolecules. Therefore, some methods are required to delete the redundant points and improve the mesh quality without changing the geometric features of molecules.

Algorithm to delete redundant points
The redundant points are chosen from the vertices whose valence is three or four. Here, valence is the number of the adjacent vertices connected by edges. Fig 1 shows the redundant condition that the valence of the point is three. When the black point is or almost on the plane formed by its three red connected neighboured points, the black point is considered as a redundant one. In this situation, the black point should be deleted. Meanwhile, the elements and the edges whose vertices contain the black point should also be deleted. Fig 2 shows the case of redundant vertex with a valence of four. When the black point is almost on one of two diagonals of the polygon formed by the four red connected neighboured points, the black point is considered as a redundant one. In this case, the black vertex, the elements and the edges whose vertices contain the black vertex should be deleted, and the diagonal on which the black vertex located should be added as a new edge. Then we can get two new triangles to replace the original four triangles containing the redundant vertex. Algorithm 1 shows the whole process to delete the redundant vertices.
The valence three and four cases are easy to handle. When the valence is three, we can delete the point directly. When the valance is four, we delete the point and re-triangulate the resulted polygon with four vertices. For this two simple cases, the deletion and re-triangulate processes will not lead to intersections, because these operations do not add any new planes. But for the cases with higher number of valence, after deleting the redundant vertex, we should re-triangulate the resulted polygon composed of the neighborhood triangles of the deleted vertex. It is hard to guarantee no intersection, especially when this polygon is not flat. Therefore, in these cases, we prefer to change the position of the vertices with more than four valence rather than delete.
Algorithm 1 Delete the redundant points.

Mesh smoothing
An inherent problem of the automated triangulated surface mesh generation is that the resulted polygonal surface appears faceted and may contains singular triangles with very small angles. The original Gaussian surface of the biomolecule is smooth. But the generated mesh can make the approximated surface rough. And the singular triangles may cause difficulties in BEM computation and corresponding body-fitted volume mesh generation. Therefore, smoothing method is necessary for improving the qualities of surface meshes.
Laplacian smoothing is an algorithm to smooth a triangular surface mesh. [35,36] The Laplacian smoothing technique changes the position of nodes without modifying the topology of the mesh. In the process of Laplacian smoothing, each vertex is move to a new position defined by taking the average position of its neighbors: where N i is the number of the adjacent vertices of the ith node, q j is the position of the jth adjacent vertex. To control the rate of smoothing, the current position q i is also included in the calculation of the new position p i : where β is the parameter to control the rate of smoothing. Usually, β takes a value between [0, 1]. However, if the geometrical characters are very complex in a region, the Laplacian smoothing may cause self-intersection. To ensure manifoldness, we proposed the modified version of Laplacian smoothing as follows. For each vertex, initially, β is chosen to be 1 and the vertex is moved to a new position through Eq (4). Then we check whether this movement causes intersection near this vertex. And if it does result in intersecting triangles, redo this movement by Eq (4) using half β. This process is repeated until the new position of current vertex does not bring intersecting triangles. And in each time of checking intersecting triangles, we only need to check whether there is an edge intersecting a triangle in the vicinity of the current moved vertex.
The key to detect self-intersection is to decide whether an edge is intersecting with a triangle. Fig 3 shows the case that an edge intersects with a triangle. In Fig 3, P 1 P 2 P 3 is a triangle and Q 1 Q 2 is an edge. If Q 2 is in the shadow of P 1 P 2 P 3 blocking light from the source Q 1 , the edge Q 1 Q 2 intersect with the triangle P 1 P 2 P 3 . In other words, if Q 1 Q 2 intersect with P 1 P 2 P 3 , following four conditions should be satisfied: • the point Q 2 and P 3 are on the same side of the plane Q 1 P 1 P 2 • the point Q 2 and P 1 are on the same side of the plane Q 1 P 2 P 3 • the point Q 2 and P 2 are on the same side of the plane Q 1 P 3 P 1 • the point Q 1 and Q 2 are on different sides of the plane P 1 P 2 P 3 . Algorithm 2 shows the whole process to detect whether an edge intersects with a triangle.

Algorithm 2
The algorithm to decide whether an edge intersects with a triangle. Input: three vertices of a triangle P 1 , P 2 , P 3 and two endpoints of an edge Q 1 , Q 2ñ ;ñ 1 )< 0 then return False (False means no intersection and True means that the edge intersects with the triangle) elsẽ  For large molecules, the amount of points and triangles are huge. If we do smoothing for all points every round, the computation is expensive. To reduce the computation complexity, we use a screening method. Every round doing the modified Laplacian smoothing, we just choose the vertices and the neighbored vertices of the low-quality triangles (the triangles with low edge ratio or containing tiny angles). The detailed algorithm is shown in algorithm 3. Firstly, we do modified Laplacian smoothing for all the points once. Then we find out all the triangles with low-quality. The vertices and the neighbored vertices of these low-quality triangles are considered as the target points which need to be smoothed next round. And we do smoothing to these points several rounds, usually three or four rounds. After this, we will find out the new irregular triangles from the targets we got last round. From the irregular triangles, we can get the new target points which should be smoothed next round. Repeating in this manner, the points which need to be smoothed are getting fewer and fewer. In this way, we don't need to do modified Laplacian smoothing for all the points in the mesh every round, which reduces the computational cost greatly. When the quality of the mesh satisfies preset goal, the smoothing process is ended. Usually, we control the number of rounds of doing modified Laplacian smoothing for irregular triangles. In our program, the default rounds is set to 100. However, in most cases, the smoothing is convergent before reaching 100 rounds.
Algorithm 3 The algorithm of mesh smoothing using modified Laplacian method.

Results
For large molecules, some surface meshes generated by TMSmesh 2.0 can not be directly used to boundary element method (BEM) simulation or generation of body-fitted volume mesh due to singular triangles with tiny angles or very short edges. The algorithms introduced in the Method section are used to improve the quality of these meshes. We test the initial and improved meshes in BEM computations of Poisson-Boltzmann (PB) electrostatics and generating corresponding surface conforming volume meshes that are required in FEM simulations. The BEM software used is a publicly available PB solver, AFMPB [37]. In FEM simulations, the corresponding volume meshes are generated by TetGen [38]. Table 1 shows some examples which demonstrate that our improvement algorithms are effective. The number of atoms of the test molecular ranges from several thousand to several hundred thousand. The PQR Benchmark can be downloaded from https://doi.org/10.6084/m9.figshare.5349229.v1. The initial molecular surface meshes are provided by TMSmesh 2.0. They are manifold meshes but contains highly irregular and singular triangles. The algorithms introduced in the Method section can remove the redundant points and improve the qualities of triangles effectively. The meshes are smoothed 100 rounds. The improved meshes can be applicable in AFMPB [37] and TetGen [38] directly. The volume mesh generated by TetGen can used in finite element method simulations, such as Ichannel [39], SMPBS [40], mFES [41] and so on. Table 2 shows the CPU time for mesh quality improvement. All computations run on a computer with Intel 1 Xeon 1 CPU E5-4650 v2 2.4GHz and 126GB memory under 64bit Linux system. Before deleting the redundant points and doing Laplacian smoothing, we first check whether the mesh is manifold and delete the small cavities inside the molecules. If we reserve the small cavities, they tend to collapse after several times of Laplacian smoothing. The fourth column in Table 2 shows the time cost by checking manifoldness and deleting cavities. The fifth column shows the CPU time cost by deleting the redundant points and doing Laplacian smoothing.
We compared the mesh quality of the improved meshes by our SMOPT with those by ISO2mesh [42]. ISO2mesh is a free matlab/octave-based mesh generation and processing toolbox. We use the meshresample and smoothsurf modules in ISO2mesh to improve the initial meshes. To make the improved meshes from ISO2mesh comparable to those from SMOPT, in the meshresample module, we set the vertex number in output meshes from ISO2mesh close to those from SMOPT. In the smoothsurf module of ISO2mesh, the Laplacian smoothing method is chosen and the meshes are smoothed 100 rounds. In addition, we also do the same test in MeshLab [43]. MeshLab is an open source system for processing and editing 3D triangular meshes. It provides a set of tools for editing, cleaning, healing, inspecting, rendering, texturing and converting meshes. We use MeshLab to do Laplacian smoothing for the meshes generated by TMSmesh 2.0. However, after several rounds of smoothing, some points with the coordinates (NaN, NaN, NaN) appear, which prevents us to do further test. So, here we only show the results from SMOPT and ISO2mesh. Fig 4 shows the initial and improved surface   Table 3 compares mesh quality in aspect of manifoldness. There is no self-intersection in the meshes produced by SMOPT while the improved meshes by ISO2mesh can not preserve the manifoldness of the original meshes. Fig 5 shows the areas and volumes of the initial meshes and the improved meshes generated by SMOPT and ISO2mesh. The areas of the improved meshes by SMOPT has small deviations from those of the initial meshes. Many roughed regions are smoothed after improvement, which will reduce the surface area. The volume of the smoothed mesh is almost preserved and keeps a deviation within 3%. ISO2mesh preserves both the mesh area and volume. An important reason is that ISO2mesh re-samples the vertices while SMOPT doesn't. We will consider to add a re-sample module into SMOPT in the future.
We also compared the uniformness of the meshes. The distribution of the ratios of the shortest edge length to the longest edge length of each triangle and the distribution of the angles of each triangle are used to describe the uniformness of a triangulated surface mesh. A  Quality improvement of surface triangular mesh avoiding intersection ratio of 1.0 corresponds to an equilateral triangle, and a ratio close to 0 is indicative of a very poor uniformness. That is, the higher the ratio, the better the quality of the triangle. Fig 6 shows that the edge ratios of the initial meshes are uniform distributed between (0, 0.8). The edge ratios of the improved meshes by our method SMOPT and ISO2mesh are all clustered around 0.7, but the results by SMOPT are more concentrated than ISO2mesh. Additionally, there is no triangle whose edge ratio is between 0 and 0.1 in the improved meshes by our method, which means that there are few triangles with very poor quality after improvement. The angle distributions of surface triangles in the initial and improved meshes are shown in Fig 7. The angle distribution of the improved meshes are clustered around 40˚to 80˚and a small portion are close to 0˚or 180˚, which also indicates that the improved meshes have few sharp angles. The left column is the ratios of the initial mesh, the middle column is the ratios of the improved mesh by our method and the right column is the ratios of the improved mesh by ISO2mesh. https://doi.org/10.1371/journal.pone.0184206.g006

Conclusion
In this paper, we have described a systematic procedure, SMOPT, to remove redundant points in a 3D triangular surface mesh and to smooth the mesh using a modified Laplacian smoothing method. The improved meshes are shown to be of good quality and can be applied to numerical simulation successfully. The program SMOPT can be used together with or integrated into a surface mesh generator, such as TMSmesh 2.0. SMOPT can be also used as a standalone software and is downloadable at https://doi.org/10.6084/m9.figshare.5346169.v1 and can be run online at http://www.xyzgate.com.