Available online at www.sciencedirect.com

SciVerse ScienceDirect

Procedía ¡engineering 29 (2012) 589 - 596

Procedía Engineering

www.elsevier.com/Iocate/procedia

2012 International Workshop on Information and Electronics Engineering (IWIEE)

An Adaptive Progressive Mesh Reconstruction Algorithm for

Spatial Discrete Points

Shujun Zhang

aCollege of Information Science & Technology, Qingdao University of Science & Technology,No. 99 Songling Road, Laoshan

District, Qingdao 266061, China bState Key Laboratory of Virtual Reality Technology and Systems, Beihang University, Beijing 100191, China

Abstract

Triangulation is an important procedure of model generation. Stereo vision can restore the depth of points but only obtain the coordinates of discrete points, which need to be formed into mesh model using surface reconstruction method. This paper proposed an adaptive progressive mesh reconstruction algorithm for spatial discrete points. Through switch of 3D and 2D Delaunay triangulation and adaptive search of the k nearest neighbors, the algorithm can restore the topology of the points and reconstruct the model rapidly. Experimental results show that the algorithm has the ability of adaptive modeling to discrete points of different scale and shape with high efficiency and low memory occupancy.

© 2011 Published by Elsevier Ltd. Selection and/or peer-review under responsibility of Harbin University of Science and Technology

Keywords: surface reconstruction; Delaunay triangulation; discrete points; mesh model

1. Introduction

Mesh surface reconstruction of spatial points is also called triangulation. It is widely used in scientific visualization, reverse engineering, computer vision and other related fields. Its purpose is to reconstruct the topology of each discrete point and obtain a linear surface (that is, polyhedral mesh) in accordance with the original structure. Because triangle is the simplest polygon and each polygon can be divided into a quantity of triangles, triangular mesh is most frequently used in 3D reconstruction. Therefore, triangulation is the key step and research hotspot in this area. Numerous triangulation methods fall into

* Corresponding author. E-mail address: lindazsj@163.com.

1877-7058 © 2011 Published by Elsevier Ltd. doi:10.1016/j.proeng.2012.01.009

two major categories: Delaunay based [1] and non-Delaunay based [2]. Delaunay based triangulation method is traditional with features of uniqueness and local optimization but still needs to be improved.

There are three commonly used 3D reconstruction methods to arbitrary topological structure [3]: isoline tracking method, carving-based method and region increasing based method.

Isoline tracking method was proposed by Hoppe [4]. Triangular surface model is acquired by first defining an index of point density and then calculating the mesh using marching cubes. It estimates the tangent plane on each sample point through computing their k nearest neighbours and linearly approaches the local shape of the curved surface waiting to be modelled. This type of method demands that the sample points should be dense and well-distributed, but it cannot deal with the data with thin boundary or shape features.

Carving-based method first carries out Delaunay triangulation and then extracts the triangle related to the shape of the point cloud according to a certain rule. The representative methods are a-shape based method [5], Voronoi neighbourhood based method [6] and UmbrellaFilter algorithm [7]. The shortcoming of this kind of method is high complexity. UmbrellaFilter algorithm needs to solve a linear programming (LP) problem in order to guarantee the lightness of the final topology. If the sample surface differs much with the reconstructed surface, the LP problem will be hard to calculate even unanswered.

Region increasing based method usually selects an initial triangle from point cloud and adds it into the surface set as the seed together with the three edges inserting into the edge set. The two sets are both updated according to a dynamic increasing principle from local to global. The representatives are Ball Pivoting Algorithm (BPA) [8] and Delaunay-based Region-growing Algorithm (DBRB) [9]. This type of method has innate strength on reconstruction of open surface mesh but it seems not very strict using only sampling uniformity degree to determine the scale of the influence region. The computation method of the influence region probably becomes the next bottle neck of the efficiency of the algorithm. Recently a fast progressive surface reconstruction algorithm was brought forth [10]. However, the result will be greatly affected by the parameter k and the paper does not clarify how to choose k.

In summary, we propose a new adaptive progressive mesh reconstruction algorithm for spatial discrete points which integrates the superiority of isoline tracking method and region growing method so as to improve the quality, efficiency and stability of the mesh modeling process.

2. Principle of the Algorithm

The core idea of our algorithm is to project the 3D discrete points into their approximate surface for 2D Delaunay triangulation, and then feedback the topology reflected by the 2D point set to the 3D point set, complete the mesh jointing, therefore, obtain the reconstructed mesh.

Through making the best of the local topology and geometric information implicated by the neighbouring point set, this algorithm rapidly implements local topology reconstruction of each point on the basis of 2D Delaunay triangulation technique. After automatic rectification of the illegal connection relationship of the local data points, it joints the local triangle net into a standard integrated 2D manifold mesh in a way of increasing expansion. The algorithm can detect the empty hole in the reconstruction procedure and judge open or close topology contained in the discrete data. It supports non-uniform point distribution.

The input of the algorithm is the spatial coordinate of the points set: p = V } and output is a Delaunay triangle net M and normal vector of each triangle. We use Vi to describe one point. It can be divided into three steps: pre-processing, mesh generation and mesh optimization.

Step 1. Pre-processing.

(1) To each Vi , find its ki nearest neighboring points and denote them as Neigh(Vi). Fit the

approximate normal vector of Vi. If the number of the nearest neighboring points of Vi is less than 3, then Vi is set to be noise point, otherwise, it is set to be isolate point. (2) Set M as a null set.

Step 2. Mesh generation: traverse each point of the point set in sequence, to edge point or isolate point V, delete the inner point and noise point from Nrigh{Vi ), then we obtain a point set Qi. If the number

of points in Qi is over 2, then the following processing should be carried out.

(1) Project Qi onto the normal surface of Vi, conduct triangulation and a triangle net is acquired. As for those triangles that are not connected with Vi, delete it from the 2D triangle net.

(2) Feedback the relationship among the points in 2D triangle net to 3D space so as to obtain the local triangle net Mi of Vi.

(3) Insert Mi to M and at the same time adjust the property of each point in Vi and Qi.

(4) Traverse the point set again until the spatial triangle net M changes no more which means no more new triangles enter.

Step 3. Mesh optimization. Compute the normal vector of each triangle and then make normalization of all the normal vectors.

3. Elaboration of the Algorithm

3. )Prr-prccrssieg

The spatial bounding box of the discrete points should first be calculated from the input point set. Using spatial partition method, the bounding box is divided into a quantity of cubes according to the distribution of the points. Therefore, the nearest k points of Vi need to be searched only from the cube that the point lies in and its 27 neighbouring cubes. To each point in the set, find its k (4< k <20) nearest points, save them into the neighbour set of this point and fit the approaching normal vector using inverse iteration method. Meanwhile, the approximate normal tangent plane taking this point as point of tangency is formed. Then the pre-processing ends.

During this procedure, the computation of k neighbours is the key problem. Since this neighbourhood means the k points who are nearest with the reference point in the discrete point set, the most straightforward idea is to calculate the distance between this point and each point then find the k points in comparison.

It will directly influence the efficiency of the algorithm how to determinate the number ki of the

closest points of data point Vi In order to guarantee the locality of the regional triangulation, ki shall

not be too large. But it shall not be too small either; otherwise the local triangulation will not be integrate and will make the inner points become edge points. Therefore, we propose an adaptive method to determinate the ki value. First, set every ki the same initial value between 5 and 20. Then, its value is

dynamically adjusted according to the curvature at the point Vi during the reconstruction period. On the

region where the surface is relatively smooth, ki can be increased; on other regions where curvature is high, ki can be properly decreased. Due to the fact that we usually do not know the surface of the

discrete points before reconstruction, still less the curvature of a certain point, we can estimate the even degree of the local surface when computing the similar normal vector. If the sum of deviation Error is small, it means the surface at this point is comparably smooth so ki will automatically grow. Otherwise,

kt is accordingly shortened. Through this dynamic adjustment method, the reconstruction algorithm can be adaptive, automatic and fit for any topologic data set.

As to point V, its similar tangent plane is computed through its k-neighborhood. This process can be

described as the following: Given the neighborhood Q of V ■

Q = Qq I i = 1,k} c P (1)

Solve the similar normal vector N of V so as to make the expression (1) the smallest, that is

error = Y\(Qi - V)• N]2 (2)

Applying the least square method, we can obtain a matrix C of 3 rows and 3 columns:

c=£ (a - V)t •(Q, - v) (3)

The approximate value of N can be calculated as the unitization of eigenvector which is corresponding to the smallest feature value of C. We make use of inverse iteration method to solve this feature value and eigenvector. The equation of tangent plane is obtained using point norm form equation after the feature normal vector is computed. Then the tangent plane is stored into vertices. at(i). tangent plane. The accuracy of tangent plane estimation plays a very important role in the projected location of the point cloud and further affects the topology of triangulation and quality of the reconstructed mesh.

3.2 Mesh Generation

Mesh generation is the core of the 3D reconstruction. Its principle is to project the spatial discrete points into the approximate tangent plane of the corresponding point to carry out local 2D Delaunay triangulation, then feed the topology of the two dimension point set back to the 3D point set and finish the mesh jointing. The procedure is shown in Fig.1.

Fig. 1 Projection of 2D local triangulation to 3D space

A. Delaunay triangulation of planar point set

After calculating the approximation vector, we can change the spatial triangulation of V and its

closest point set Q to 2D triangulation of projection points on its tangent plane. The 2D Delaunay triangulation algorithm we apply is as follows. To a 2D point set pSet, if its scale if less than 3, then the processing ends. Otherwise, set the formed triangle set as empty and add three appendent points in order to construct a huge enough triangle HjgrTriaegle. Assign the three points three different number less than 0 as their indexes so as to distinguish them from the original point set. Add the HjgrTriaegle into the formed triangle set and construct a 2D Delaunay triangle net using a progressive method. As to each unprocessed pset.at(i) in the point set, traverse the triangle set currently formed and structure all the triangles whose circumcircle includes this point into a convex hull (the vertex is clockwise arrayed). Connect this point and every two neighboring vertexes into a triangle tempTri, insert it into triaegleSet until there are no point unprocessed in the set. Now, delete the triangles whose vertexes include the appendent points. Finally, delete the triangle that is not directly connected with the vertex vertices.at(i).

B. Judgment of Spatial Triangle Overlay

Overlay of two triangles is relative to a certain vertex. We project all the vertexes of two triangles into the approximate cutting plane of their public vertex and obtain another two triangles on the same plane. Through judgment whether these two triangles on the same plane is overlapped, we can determinate whether the two spatial triangles are overlapped. Suppose the public vertex of two projected triangles is A, we can denote two vertexes of one triangle B and C, two vertexes of the other triangle D and E, as shown in Fig. 2.

We use the following method to check whether two coplanar triangles are overlapped. If they really lap over, then a certain edge of one triangle whose vertexes include this public vertex is sure to divide the other triangle into two parts, as shown in Fig.2. The sum of inclination of edge AB and edge AD (¿BAD) and that of edge of AD and edge AC (¿DAC) is equal to the inclination of edge AB and edge AC (ZBAC). Therefore, through checking whether there are edges with this relationship in the two triangles, we can determinate the overlay problem. It is worth noting that, if the two triangles have only one public edge or completely overlapped, then they need to further judgment.

C. Jointing of 3D Mesh

Every new generated local spatial triangle mesh in our algorithm should be integrated to currently formed global mesh. This is progressive extension. This method takes each triangle from the local mesh in sequence and check whether it lap over the triangles in the global mesh. If not, then this triangle is inserted into the global mesh. The property of corresponding spatial discrete point will be adaptively changed when the triangle adds into the global mesh.

3.3 Mesh Optimization

Fig.2 Overlay of two coplanar triangles

Normal vector of each triangle is calculated after global spatial triangular net is generated. Since some of the normal vectors point inside, some point outside, mesh optimization step focuses on the unification

of normal vector, which means adjusting outward all the normal vectors of the triangles in global mesh. This is necessary for mesh model rendering and display. To a triangle AABC, its normal vector is described as:

N = (B - A)x (C - A)/||(B - A)x (C - A)| (4)

If N points outward, we need to change the storage sequence of two vertexes in the triangle and set N to the opposite direction.

4. Experimental Results and Analysis

We made experiments to test the proposed algorithm. Spatial discrete points are loaded from files for 3D reconstruction. The generated mesh can be evaluated from the difference between the original model. Here we select eight point files with different scale as the testing data. Tab.1 gives the number of points, the k parameters and number of triangles after the algorithm runs.

Tab.1 Experimental data and results

Point set file Number of discrete points k value Number of triangles

sphere 146 8 288

g2 1 3200 8 6400

g2 3 4795 8 9590

head 12097 12 24002

bunny 35947 15 71828

locust 64226 18 124750

horse 83516 16 165936

jug 193592 20 384552

As Tab.1 shows, to the discrete point sets with different scales, the proposed algorithm can effectively reconstruct the mesh surface and from the number of triangles we can see, generated models have comparably fine accuracy. Fig.3 gives the model result of sphere and bunny.

From Fig.3 we can see, the reconstructed surfaces are very uniform and smooth. Our algorithm can effectively reconstruct the discrete points and restore the original model.

We also apply the algorithm to handle the discrete point cloud obtained using image-based reconstruction method. Fig.4 shows one frame of the original video and reconstructed wireframe model. The reconstructed model will completely restore the information of the object after texture mapping.

Fig.4 Original image frame and point cloud reconstructed wireframe model

To analyze of the experimental results, we give the following graphs. According to the data from Tab. 1, we conclude the relationship between the number of reconstructed triangles Nt and the number of discrete points Np in Fig. 5(a). We also test the relationship between the time for 3D reconstruction (denoted as T) and Np in Fig. 5(b).

Np/thousani

-C ■H

MP/thr

Fig.5 (a) Line chart of Nt and Np ; (b) Line chart of T and Np

From Fig.5(a) we can see that, the two parameters have linear relationship. The number of reconstructed triangles is about twice of the number of discrete points. This is in accordance with the feature that the generated mesh is 2D manifold.

From Fig.5(b) it is found that when the scale of the point set is comparably small (e.g, less than 50 thousand) the time and the number of points have linear relationship, while when Np increases, the time

of reconstruction becomes unsteady, and also has great influence with the parameter k. Therefore, through adaptive determination of k value, the efficiency of reconstruction algorithm is effectively controlled.

5. Conclusions

This paper proposed a new adaptive progressive rapid surface reconstruction algorithm based on local 2D Delaunay triangulation. This method can automatically determinate the number of neighboring point set and fits for almost all point set with different scale and shape. It has high efficiency. The algorithm can be applied into point cloud generated from stereo vision to obtain smooth and accurate 3D mesh model. This technique has wide application prospect in many fields including 3D information recover, surface measuring and visualization.

Acknowledgments

This paper is supported by the National Natural Science Foundation (No. 60903064 and No. 61040047) and Open Foundation of the State Key Laboratory of Virtual Reality Technology and Systems (No. BUAA-VR-11KF).

References

[1] Okabe, Atsuyuki, Boots Barry, Sugihara Kokichi. Spatial Tessellations: Concepts and Applications of Voronoi Diagrams. John Wiley & Sons Ltd. 2009

[2] M.de Berg, M.van Kreveld. Computational Geometry Algorithms and Applications, Second Edition. Beijing: Tsinghua University Press, 2005.

[3] Zhenming Zhang. Research on Delaunay-based Surface Reconstruction Algorithms from 3D Unorganized Point Clouds. Nangjing: Nangjing University of Aeronautics & Astronautics, dissertation of Master Degree, 2006.

[4] H Hoppe, Spherical parametrization and remeshing, ACM Transactions on Graphics (TOG), 2003

[5] Dd H. Edelsbrunner, E.P. Mucke. Three- dimensional alpha shapes. ACM Transactions on Graphics, 13:43 - 72, 1994.

[6] Amenta N, Bern M, Kamvysselis M. A New Voronoi-based Surface Reconstruction Algorithm. Orlando: Proceedings of SIGGRAPH, ACM Press, 1998, pp. 415- 420

[7] Adamy U, Giesen J, John M. Surface reconstruction using umbrella filters. Comput Geometry 2002, (21):63-86.

[8] Fausto Bernardini, Joshua Mittleman, Holly Rushmeier, et al. The Ball-Pivoting Algorithm for Surface Reconstruction, IEEE Transactions on Visualization and Computer Graphics, v.5 n.4, p.349-359, October 1999

[9] CC Kuo. A Delaunay-based region-growing approach to surface reconstruction from unorganized points, Computer-Aided Design, 2005

[10] Wang Qing, Wang Rongqing, Bao Hujun, et al. A Fast Progressive Surface Reconstruction Algorithm for Unorganized Points. Journal of Software, 2000, 11(9): 1221-1227.