Volume 43, Issue 6, June 2011, Pages 639-650
doi:10.1016/j.cad.2011.02.012 | How to Cite or Link Using DOI
Permissions & Reprints
A new mesh-growing algorithm for fast surface reconstruction
a, a, aLuca Di Angelo, Paolo Di Stefanoand Luigi Giaccari
a Università degli Studi de L‘Aquila, Italy
Received 7 September 2010;
accepted 22 February 2011.
Available online 1 March 2011.
This paper presents a new high-performance method for triangular mesh generation based on a mesh-growing approach. Starting from a seed triangle, the algorithm grows the triangular mesh by selecting a new point based on the Gabriel 2—Simplex criterion. This criterion can be considered to be a good approximation of the 2D Delaunay if the point cloud is well-sampled and not too rough. The performance of the proposed method is compared with that of the Cocone family and that of Ball Pivoting
as regards the tessellation rate and the quality of the surface being generated from some benchmark point clouds and artificially noised test cases. The results are analysed and critically discussed. Highlights
? A new surface reconstruction method based on a mesh growing approach is presented. ? Surface reconstruction is carried out by using the Gabriel two—Simplex Criterion. ? For well-sampled and not
rough point clouds it is similar to the 2D Delaunay approach. ? The proposed method is competitive for
tessellation rate, quality and defectiveness.
Keywords: Surface reconstruction; Triangular meshes; Reverse engineering
2.1. The implicit methods
2.2. The explicit methods
2.2.1. Voronoi/Delaunay-based methods
2.2.2. Mesh-growing methods
Gabriel 2—Simplex criterion
4.1. Selection of the seed triangle
5.1. Typical benchmarks
5.2. Noisy point clouds
The process of constructing a triangular mesh from a point cloud which has been acquired from a 3D real object is known as a surface tessellation or reconstruction. The applications of the triangular meshing of 3D point clouds are wide-ranging and include Reverse Engineering, Collaborative Design, Inspection, Computer Vision, Dissemination of Museum artefacts, Medicine, Special Effects, Games and Virtual Worlds. For many of the already-mentioned applications, robustness, reliability and triangulation speed are important factors that are required in any tessellation method. Equally critical is the capability to tessellate according to the nature of the surface which is to be reconstructed.
Over the last few years 3D scanning systems have been offering high resolutions with a measuring accuracy as high as 10 μm, which, on the one hand, makes it possible to capture the smallest surface features but, on the other hand, generates very large data sets. The typical triangulation speed of the methods presented in literature may not be enough for the tessellation of a cloud with over one million
points. Furthermore, in analysing the peak memory usage of these algorithms, it becomes clear that most of them cannot be run on personal computers.
In order to overcome these limitations, this paper proposes a new algorithm. It is based on the concept of the Gabriel 2—Simplex (G2S). This concept allows us to extend 2D tessellation approaches to 3D triangulations. As long as point clouds are locally flat, the method being proposed guarantees a correct triangulation. Unlike some other methods presented in literature, the one being put forth does not require that the surface defined by the point cloud should enclose a volume or have a strictly uniform sampling. This method is tested for the tessellation of some benchmark point clouds and artificially noised test cases. The results derived from those experiments are critically discussed further down. 2. Related works
Numerous algorithms to tessellate point clouds are presented in literature. The more recent exhaustive overviews of 3D surface reconstruction methods are presented in  and . Surface reconstruction
algorithms are generally divided into two categories: implicit and explicit.
2.1. The implicit methods
Implicit methods seek to reconstruct an implicit function , where is either the whole point cloud or only a part of it. The final triangulation is obtained by extracting the iso-surface for . The most common methods define the implicit function as:
– the sum of radial basis functions (RBF) centred at the points
– a set of constraints that force the function to assume given values at the sample points and also force its upward gradient to match the assigned normal at the sample points (Moving Least Squares)  and
– a Poisson problem
These implicit methods carry out a watertight surface reconstruction also in the case of sparse and noisy data. However, these methods may generally require many computations and, sometimes, they may even need the surface normal at each data point. Moreover, the triangles of the final surface may not pass through all the points, and this in turn may imply the loss of some details of the shape of the original model. The more points we use to compute the implicit function, the higher will be the fitting accuracy and the longer will be the computation time.
2.2. The explicit methods
Explicit methods attempt to triangulate the points directly. In contrast with the implicit methods, they are less robust against noise although they are generally faster. The most common explicit methods can be classified into two main groups: Voronoi/Delaunay-based and mesh-growing approaches. 2.2.1. Voronoi/Delaunay-based methods
This first group includes algorithms that compute a volume tetrahedralisation by means of a 3D Delaunay triangulation of the sample points. The most important methods presented in the related literature (, , , , ,  and ) essentially differ in the way they remove the tetrahedra and build the external triangular mesh. In the Crust method, proposed by Amenta et al. in , the set of
candidate triangles (called Crust) are those having three vertices which are not poles. Since the Crust
is still not a manifold, in order to extract a topologically correct surface, a walking strategy is employed
2in candidate triangles. The worst time complexity of the Crust algorithm is Θ(m), where m is the
number of points and poles. In order to reduce running time and memory consumption, while still providing the same theoretical guarantees, Amenta et al. in  proposed an improvement of the Crust
algorithm, which they called Cocone. The candidate triangles are those triangles inside the Cocone
region that is the complement of a double cone which has its apex at the point under analysis () and an assigned opening angle and whose axis is the normal at . The Cocone’s worst time complexity is
2Θ(n), where n is the number of points. As pointed out by Chang et al. in , the theoretical guarantee
of obtaining a correct reconstruction with the Crust and Cocone methods is only possible as long as the
point cloud is well sampled. Dey and Goswami proposed other Cocone versions, one suited to provide
a watertight closed surface reconstruction, which they called Tight Cocone  and one suited for
noisy data called Robust Cocone . Dey et al. in  proposed the Super Cocone algorithm so as to
manage large amounts of data. In that method, the entire set of sample points is partitioned into smaller clusters using an octree subdivision, and the Cocone algorithm is applied to each cluster
separately. A further evolution of the Crust is the Power Crust proposed by Amenta et al. in . This
algorithm can generate a watertight mesh for any point cloud which is sampled enough but it could be non-manifold. Furthermore, sharp edges or noisy data can result in defectiveness. Gopi et al. in 
proposed a different approach which, for each sample point, provides the projection of the neighbouring points onto the approximating tangent plane and the tessellation of the projections by
means of a 2-D Delaunay triangulation. The 2D edges obtained are then applied to 3D space. With a view to enhancing the performance of the Delaunay-based methods, Cohen-Steiner and Da in 
proposed the Greedy algorithm, which selects the triangles sequentially. The selection of candidate
triangles is carried out by the plausibility grade, which is an empirical function of the circumradius and
the dihedral angle between adjacent triangles. As stated by Cazals and Giesen in , this method
―may fail to interpolate all points or to provide a closed surface essentially due to the presence of slivers‖.
2.2.2. Mesh-growing methods
Mesh-growing approaches generate the tessellated surface from a seed triangle and grow the meshed area pushing the fronts ahead by using some criteria. Over the last few years a number of algorithms have been presented. Bernardini et al.  introduced the Ball Pivoting algorithm, by which the front
grows as a ball of user-defined radius pivots around the front edge. When the ball touches three points, a new triangle is formed. Generally speaking, this method affords a correct triangulation for any uniform data point, which is also noisy, but which does not have any concave sharp features. Huang and Menq in  proposed an algorithm which, for each front edge, projects the k-nearest points of two endpoints onto the plane that is defined by the triangle adjacent to the front edge. Any point generating triangles whose edges intersect edges of already-existing triangles is discarded. Of the points which are retained, the point showing the minimum sum of distances from the front-edge endpoints is chosen. Nonetheless, and as pointed out by Lin et al. in , this algorithm presents some shortcomings. In
order to overcome them, Lin et al. in  introduced the Intrinsic Property Driven (IPD) algorithm,
which improves the way of searching for the point to be triangulated. As stated by Chang et al. in , all
the methods based on mesh-growing approaches are fast, efficient and simple to implement but they, however, may fall short whenever two surfaces are either close together or near sharp features. More recently, Li et al. in  proposed a method based on a Priority Driven approach that evaluates shape
changes from an estimation of the original surface that is made at the front of the mesh-growing area. The experimental results in  evidenced that the triangulation speed of the method was higher than that of the Ball Pivoting and the Cocone. However, no reckoning was made of the defectiveness it
3. Gabriel 2—Simplex criterion
The algorithm being proposed is based on the criterion, henceforth referred to as the Gabriel
2—Simplex criterion (G2S): A triangle is a G2S if its smallest circumscribing ball is empty . We
should bear in mind the fact that if all the points are coplanar, the obtained triangulation is a classical 2D Delaunay. Therefore, the flatter the surface appears to be locally, the more the criterion tends to be a 2D Delaunay and the better is the guarantee of a good reconstruction. More generally, as pointed out by Dyer et al. in , any mesh in which every triangle verifies the G2S criterion is a Delaunay mesh.
This criterion had already been used by Ruiz et al. in  for the tessellation of parametric surfaces.
The method cannot be used directly for the triangulation of point clouds since it is based on information which needs to be extracted from the parametric surface (pre-image in a parametric space, boundary, and local curvatures). Along the same lines, Cohen-Steiner and Da in  used a similar criterion to
G2S in order to choose the external triangular facets of tetrahedra coming from a 3D Delaunay triangulation. However, in the algorithm here proposed, the G2S is applied directly to a point cloud so
as to increasingly identify triangles pertaining to the external surface. The G2S takes into consideration
spheres not in order to identify tetrahedra, as the Delaunay 3D criterion does, but rather to generate triangles, as the Delaunay 2D criterion does. Fig. 1 shows examples of two G2S triangles. The criterion
being applied does not ensure that the G2S triangles belong to the external surface. In fact, the G2S
triangles can either lie on the external surface or be transversal to it. That is due to the nature of the method which searches for the nearest points in the 3D area inside a sphere. This problem may arise in thin parts.
Full-size image (35K)
Examples of G2S triangles; the two spheres respectively circumscribing the vertices of the red and the yellow triangles are empty.
View Within Article
For a given data set, the G2S triangulation is not unique and depends on the seed triangle which is chosen, as opposed to the 2D and 3D classical Delaunay criterion. Such is the case of the four points () shown in Fig. 2(a). It is possible to verify that all three triangles coloured blue in Fig. 2
( and ) fulfil the G2S criterion. If
the method builds the triangle T first, then the G2S triangle T is generated (Fig. 2(b)). However, 1,2,41,3,4
if the method builds the triangle T first, then the triangle being generated is 1,2,3
because its smallest circumscribing sphere contains only the that has just been triangulated (Fig. 2(c)).
Full-size image (41K)
An ambiguous data set for the G2S criterion.
View Within Article
4. Algorithm description
The algorithm herein proposed performs the tessellation of a 3D point cloud with a mesh-growing method. The algorithm takes a 3D point cloud as input and returns an oriented triangulated surface. The algorithm can be summarised in the following steps:
Step 1 Import of the point cloud;
Step 2 Building of a specific data structure to speed up the search for the nearest point; Step 3 Search for the seed triangle;
Step 4 Triangulation based on the G2S criterion.
The first step involves importing the point cloud which pertains to a continuous surface in the form of the coordinates x,y,z. Each point of the point cloud is kept in a hash table data structure for both point indexing and nearest-neighbour search. In the algorithm being presented the data structure which has been used is an improvement of those proposed by Hoppe et al. in  and Turk et al. in . Next, the
seed triangle is selected by using a new approach which guarantees that it pertains to the external surface and not to sharp features of the point cloud. The edges of the seed triangle constitute the initial advancing front of the mesh-growing method. In order to push the fronts ahead in a topologically consistent way, a new strategy is introduced, which is based on the previously-defined G2S criterion.
The procedure ends when all the free edges have been processed.
4.1. Selection of the seed triangle
The selection of the seed triangle is an important aspect of any tessellation method based on a mesh-growing approach. Some methods are presented in the related literature. Huang and Menq in  proposed an algorithm that constructs the seed triangle from the point having the maximum value for the z coordinate and its two closest neighbours. Very often the points having the maximum value for the z coordinate are either outliers or points which result from noise in the measuring process. Cohen-Steiner and Da in  selected the seed triangle having the smallest circumradius among those coming from a previous Delaunay triangulation. Li et al. in  proposed a complicated method that
searches for the seed triangle in the flattest region selected from a set of randomly chosen areas. The method brought forth by Bernardini et al.  builds the seed triangle by randomly choosing a point and
its two closest neighbours. This triangle is assumed to be a seed triangle if we are able to verify both the consistency of its normal with the estimation of the normal at three vertices and the Ball Pivoting
The method herein introduced chooses the seed triangle from among G2S triangles which lie on the
external surface of the object.
The seed triangle is selected according to the following sub-steps:
Step 3.1 A point of the cloud is randomly chosen.
Step 3.2 Its nearest neighbour point is searched for, and an edge is formed between these two points.
Step 3.3 A range search is performed inside a sphere which is centred at the midpoint of the edge and whose radius is k-times the length of the edge (k can be either a user-defined parameter or a
Step 3.4 For every point in the range, the method manages to build a triangle that is formed by the edge and the concerning point. The triangle is verified to be G2S. Then, in order to ensure that the
triangle lies on the surface, we need to verify that an infinite cylinder passing through the triangle vertices and having its axis parallel to the normal of the triangle is empty on the outer side of the triangle. In other words, the points contained in this cylinder must be either all above or all under the triangle (
Full-size image (45K)
The cylinder test to select a good seed triangle.
View Within Article
These steps are repeated until a triangle is found which satisfies the two above-mentioned conditions. The cylinder test also excludes any G2S triangles which are transversal to the surface.
Generally speaking, in a mesh-growing approach, for each front edge, a point in the neighbourhood (henceforth called reference point) must be chosen so that a new triangle pertaining to the surface may be found. The front of the growing mesh is made up of free edges, which are edges pertaining to only
one triangle. The first edges getting into the front queue are the ones from the seed triangle. For every free edge () a triangle is generated according to the following sub-steps:
Step 4.1: Identification of all the candidate points near .
The candidate points are those inside a sphere having its centre on the axis of the free edge which lies
on the plane of the front triangle, in the growing direction, and having as radius the search radius (Fig. 4). If no points are found inside the search region, the free edge is classified as boundary edge.
The search radius value affects tessellation quality and time. The value of the search radius is specified as a fraction of the front edge length: a search radius equal to 1 indicates that it has the same length as the front edge. A low search radius keeps the search region near the front edge, so the method falls short in the meshing of under-sampled areas. A high search radius, on the other hand, reduces the tessellation rate and sometimes generates defectiveness. An automatic approach has been implemented to set the search radius value. It changes the search radius value from an initial
value (in the following test cases it is 1) to an assigned maximum value according to a given step value. The search radius value is increased if no points are found inside the search region.
Full-size image (40K)
The search region‘s definition terms.
View Within Article
Since the search region centre lies on the same plane as the front triangle, the search method usually favours finding candidate points in the area wherein the surface is expected to grow. In typical practical situations, this avoids further controls that may cause the algorithm to slacken. Furthermore, this search region excludes as candidate points any outliers or points which could give generate thin and slivery triangles. The method considers outliers (points that should not be triangulated) those points which, as regards the mesh growing direction, are farther away or those points which, for every front edge, fall outside the search region. Finally, this criterion generally stops the front in the presence of sharp edges.
Step 4.2: Selection of the reference point in the search region.
When more than one point is found inside the search region ( and in Fig. 5(a)), a point is
randomly chosen ( in Fig. 5(b)) and the smallest sphere circumscribing that point and the front edge‘s points is traced (the first search sphere in Fig. 5(b)). Then, if this sphere is not empty ( and
in Fig. 5(c)), another point inside the region is picked () and the procedure is repeated (the second search sphere in Fig. 5(c)). The process ends when the sphere passing through the chosen
reference point and the extremes of is empty. The triangle defined by the front edge and the reference point is a candidate triangle which needs to be verified in the next step (Fig. 5(d)). In order to
speed up the triangulation process, if only one point is found inside the search region, that point is directly assumed to identify a candidate triangle with the free edge () without verifying whether or