EFFICIENT RECONSTRUCTION OF NURBS SURFACES FOR SHAPE ANALYSIS AND SURFACE INSPECTION Christian Teutsch Dirk Berndt Erik Trostmann Michael Weber Fraunhofer Institute for Factory Operation and Automation IFF, Magdeburg, Germany Email:
[email protected]
Abstract: We present a new approach for the fast visualization and evaluation of large point clouds derived from different structured light 3D-scanners. By using the underlying measuring principle, we generate regular quadrilateral meshes in nearly real-time. These meshes are additionally used for a partial NURBS approximation. Based on the derived parameterized surfaces we are able to calculate and measure higher degree surface features, e.g. curvature and its direction. Finally, we evaluate the efficiency and feasibility of our approaches at the example of a variety of different point clouds.
a) b) c) Figure 1: Illustration of the generated mesh (a), together with the finally shaded model (b) and the computed Gaussian curvature on the object surface (c). 1. Introduction The industrial application of 3D measurement techniques is a common tool in quality assurance and the automated measurement and comparison of multidimensional geometries. Most 3D scanners use the structured light and triangulation principle to compute the 3D data. Therefore, the evaluation of surface features needs known interrelationships between single points in the generated point clouds. In general, polygon meshes are used to compute ordered connections in an unstructured set of points. Mesh generation is typically based on Delaunay triangulations or marching cubes methods. Another usual mechanism is the approximation of parametric surfaces to given meshes. In this case, the underlying meshes typically consist of quadrilateral polygons than of triangles. In practise, there is a wide variety of applications which use regular, quadrilateral surface meshes. The most popular is the computer aided design (CAD). Thus, the conversion from 3D points back to CAD models (reverse engineering) is still an active research area. In this work, we also consider the generation of
regular quadrilateral meshes from raw measured 3D points with respect to the fast visualization and precise geometry evaluation of large point clouds. Duplicating unique works of art, archaeological finds or work pieces for the industrial mass production are usual applications for 3D scanners. The scanning procedure can be timeconsuming, especially for large or complex objects with lots of small and important surface features. In most cases, the shape that should be scanned is unknown. To assure that all areas of the object are captured within the actual partial scan, we need a fast preview method, which assists the user to control and evaluate the scanning process. A fast mesh-based preview allows a user to check for gaps and interfering reflections. As a result the scanning procedure becomes faster and more comfortable and thus, the results are even better. Furthermore, the surface features have to be detected and evaluated, to check the geometry and optimize the scanning procedure. Based on the known scanning system and the underlying measuring principle, we present a new approach for a fast mesh generation procedure. The resulting meshes give rather precise information on point-to-point connections. In section two we have a look at related work in this research area to show advantages, disadvantages and differences to existing methods. The main approach for a fast grid and mesh construction is considered in section 3. We describe new methods for pre-processing the scanned point cloud, the computation of the quadrilateral mesh and a post-processing step to optimize the meshes. The approximation of NURBS surfaces and the calculation of higher degree surface features are discussed in section 4. 2. Related Work The analysis of scanned and unknown free-form surfaces is an active area of research. Building up an adequate scanning and measuring system is just the first step. Additionally, a surface reconstruction and in particular a meshing and triangulation of the generated point clouds is needed. Thus, the reconstruction problem is an interpolation and approximation problem for a given set of points in space. Typically, this is done by constructing a triangular mesh which interpolates in all the given points and in this way defines a piecewise linear interpolating surface. Many methods have been developed over the past years. Levoy et al. applied the entire process chain to the famous statues of Michelangelo in the Digital Michelangelo Project [15]. They present hard- and software based procedures for capturing large objects including the surface texture. In addition Bernardini et al. use a combined multi-view and photometric system to subsequently scan the surface of Michelangelos’ Florentine Pietá [5]. This approach is based on the general model acquisition pipeline, which is presented in [3]. Typically, there are several classes of triangulation algorithms. Many approaches use Delaunay triangulation based techniques or the Delaunay tetrahedrization and a filtering strategy for constructing a triangular mesh as a subset of the Delaunay triangulation. These methods are often robust and have proven reconstruction properties concerning sampling density, triangular relationships and quality guarantees [2, 6]. However they tend to be costly, because of the time-consuming computation of the Delaunay triangulation or its dual Voronoi diagram [9, 11]. Additional there are implicit methods, which usually involve the computation of a signed distance function, whose zero isosurface is considered to be the reconstructed implicit surface. The reconstructed implicit surface can be used to generate a piecewise triangle mesh, by means of an iso-surfacing algorithm such as Marching Cubes [7,14]. The reconstructed surfaces are usually very smooth, which makes this method suitable for noisy data or an esthetic visualization [8]. This method is rarely used for surface measurements.
There are also local triangulations or advancing front algorithms [16]. Bernardini et al. presented the Ball-Pivoting Algorithm (BPA) [4]. This method uses a ball, which is rolling over the entire surface. Subsequently, new points and edges, touched by the ball, are added to the triangle mesh. However this approach uses an approximation of the normal at each data point, and thus, requires multiple passes for unevenly sampled surfaces. Besides classical 3D triangulation algorithms, there are several methods based on a meshless parameterization. Floater and Reimers presented a mechanism to obtain a 3D triangulation by triangulating the parameterized points in the plane and lifting the resulting mesh into 3D [12]. Large systems of linear equations have to be solved to compute the parameterization. However, this approach is restricted by the requirement that the point cloud must be homeomorphic to a disc, and a valid boundary has to be computed for each point cloud. The most interesting surface description four our approach are b-spline surfaces. This kind of surface guarantees a certain geometric continuity depending on the order of their basis functions. Park et al. presented a method to describe a measured point cloud as a mesh of lots of small rational b-spline surfaces [17, 10]. Furthermore, b-spline representations allow a weighting of single control points, for example, based on their quality. 3. Partial Surface Reconstruction For the generation of point clouds, we developed a flexible 3D scanning system. This system consists of a multi-axis locomotor system and two structured light sensors with digital cameras and line lasers. The object is rotated and translated in two different directions in front of the sensors. For each position, both of the sensors project and measure a line on the objects surface. Capturing large and complex objects often requires multiple scans. For each scan, the laser or the object is translatory moved and the resulting set of scanlines is stored as one operation. Subsequently, the object is rotated and new translatory scans are applied. This procedure is repeated until the entire surface is sampled. Finally, the generated point cloud consists of a number of scans, which themselves consist of several scanned lines. The translatory movement is not necessary for area-based scanners, e.g. stripe projection, phase-shift and photogrammetric methods. In these cases, a single scan implicitly consists of an ordered set of scanlines. 3.1. Mesh Construction The main assumption is that nearly no point cloud is really unstructured. The underlying measuring principle and, for example, the movements of the different sensors yields a lot of additional information, that can be used for new and faster algorithms. In this work, we mainly consider the meshing of neighbouring scanlines, which seems to be a trivial task. But only using the Euclidian coordinates, allows no unique assignment. Points on neighbouring scanlines must not be neighbours, even if some of the spatial coordinates are nearly equal. This problem is additionally illustrated in Figure 2. Therefore, we have to apply a new and straightforward approach which uses the known sorting and the projection parameters. The first step is the pre-processing and optimization of the raw input data. This is a necessary procedure to minimize the influence of erroneous points.
3.2. Pre-Processing In a fast pre-processing step, single lines are evaluated and a weight is calculated for each point. For each sensor, we use the camera viewing direction and the laser direction of projection to evaluate the quality of a surface point with respect to its normal vector. Surface points that have been seen under good viewing/projection conditions get a higher weight than the others. The weight is additionally influenced by the curvature of the point on the scanline to keep edges. Furthermore, we apply a controlled noise reduction to minimize the influence of erroneous points in further measurements. Depending on the surface and the distances between sensor and surface, the point distribution is irregular. This problem is solved by resampling the scanlines based on a NURBS curve approximation with controlled precision. As a result of the automated correction methods, we get an optimized scanline, which can be post-processed with respect to other neighbouring lines. The line evaluation and the correction procedures are explicitly described in [18]. 3.3. Grid Projection Our approach aims to construct a regular quadrilateral mesh without pre-calculating triangulated polygon meshes. We use a method, which is sometimes called mapped meshing or grid projection. Therefore, a grid and its dimensions are manually defined. Actually the vertical dimension (x-direction) of the mesh is equal to the number of scanlines. The horizontal sampling rate is defined by the user. The relationship between two scanlines is a known system parameter. Depending on the scanning direction and the number of sampling steps, we compute the sorting and the distances between each two scanlines per scan. Additionally the sorting of the points on a scanline is checked, depending on the actual laser or projector position. The main procedure is the projection of the pre-defined grid onto the partial point cloud. To achieve a regular distribution and an unambiguous mapping, we use the position of the projector. This is the only unique position for each scanline, from which all points have been seen without intersections and shadowing effects. Based on the perspective projection, we assign a unique direction of projection di and the corresponding projection angle to each point per scanline (see Fig. 2).
laser 1
object
position n
dnmax d0max
position 0 plane n plane 0 dnmin laser 0
d0min
a) b) Figure 2: Measuring principle with two sensors and highlighted overlay regions (a) and calculation of the minimum/maximum directions of projection with respect to one scan (b).
The next step is to find a corresponding neighbour for each point with respect to all scanlines of the scans. We assume, that points that have been seen from the same direction of projection, are at least near or neighbours. This is valid, if the surface is slightly changing from scanline to scanline. For precise surface measurements, the sampling rate in all directions is surely high enough. In practise, neighbouring scanlines of complex objects have different point densities and distributions. Additionally, the length and shape are changing. We have to assure, that all points from the same direction of projection for all scanlines of one scan are mapped to one horizontal line in the mesh. Therefore, we use the unique projection angle. Considering all scanlines, there is a minimum and a maximum projection angle. These two angles limit the mesh in the projection space. All mesh positions are calculated based on the manually defined dimension by linear interpolation between the two angles. Finally, each point of a scanline is back projected to the mesh and for each knot a real 3D coordinate is assigned. Changing point densities and gaps in the point cloud cause knot points without valid 3D positions (see Fig. 3). To keep the mesh consistent, we assign a dummy point to these special positions. In further considerations, such knots are not evaluated. The discretization of the measured points may result in single gaps within the computed mesh. Thus, the initial basic meshes possibly have to be repaired and smoothed. Gaps are located as unused dummy knots in the grid structure. Small gaps are closed by linear interpolation. Because of the very regular structure of the meshes fast algorithms can be applied. For example, smoothing is realized with the Laplacian operator, also known from image processing (1). For each mesh point, all n surrounding points are considered. The influence of these points on the new position of the actual point is controlled with the weighting parameter λ. As a result, the mesh structure is more regular, and noise is reduced.
⎛ x old − x0old x0new = x0old + λ ∑ ⎜⎜ i n ⎝
⎞ ⎟⎟, 0 < λ < 1 ⎠
(1)
Initially unstructured point clouds are now structured and interpolated by a regular mesh. Therefore, we also reduced the spatial dimension from 3D to 2.5D by using the projector positions for a unique grid mapping. 3.4. Mesh Registration Because of uncertainties from the sensor calibration procedures, there are small deviations and distortions between overlaying scanlines from different sensors. The resulting meshes do not match exactly. Several approaches can be used to perform surface-based registration. A popular one, that we also used, is the Iterative Closest Point (ICP) algorithm. In this approach, the objective is to compute a transformation that minimizes the disparity function X in which {xj} for j= 1,..,Nx is a set of point on the surface X and yj is a point on the surface Y that corresponds to the point xj (e.g., the closest point on the surface Y). To find this closest point, one typically triangulates the surface Y and projects the point xj onto the surface. d (T ( X ), Y ) =
Nx
∑ T (x j =1
j
)− yj
2
,
y j = C (T ( x j ), Y )
(2)
This is done by computing the intersection of a line passing through the point that is perpendicular to the surface patches close to the point. For the generated regular meshes, the correspondence problem can be solved easily. The information on neighbouring points needs no additional triangulation. Furthermore, the distances between overlaying meshes are only small and thus, only a few iterations are needed. Once the point correspondence is established, the transformation that minimizes (2) can be computed using the least-squares method. The process is repeated until convergence. This algorithm is guaranteed to converge to, at least, a local minimum of (2). 4. Surface Measurements Precise measurements on the interpolated meshes require a mathematical surface description. The described meshing approach as well as triangulation methods (Delaunay, Marching Cubes, Ball-Pivoting) generate piecewise linear surfaces. Determining surface features on such models with differential geometry methods is impossible. Therefore, we need surfaces of higher degree and geometric continuity.
a) b) c) Figure 3: Approximating the basic grid (with marked unused knots) (a) by a NURBS surface (b), and an example of a real grid after back projection to the 3D space (c). The optimal solution is a parametric surface with a flexible shape description, which gives the possibility to measure at arbitrary positions. Thus, we use weighted b-spline surfaces with a non-uniform knot distribution (NURBS). We use the generated quadrilateral meshes as control net for a b-spline surface approximation. Because of some unused areas and existing gaps in the underlying mesh, we have to define weights. That’s why we use rational functions. Knot positions in the mesh with unused points get a zero weight. The influence of these points on their environment becomes zero and the consistency of the control net is kept. The weights for valid knot points are either one or defined by the pre-processing step in section 3.2. 4.1. NURBS-Approximation A non-uniform rational b-spline surface of degree (p,q) is defined by the basis functions N, the weights w and the points of the control net P (3).
⎞ ⎞ ⎛ n m ⎛ n m S (u, v) = ⎜⎜ ∑∑ N i , p (u ) N j ,q (v) wi , j Pi , j ⎟⎟ ⋅ ⎜⎜ ∑∑ N i , p (u ) N j ,q (v) wi , j ⎟⎟ ⎠ ⎠ ⎝ i =0 j =0 ⎝ i =0 j =0 The basis functions are calculated using the following equation:
-1
(3)
N i , p (u ) =
u i + p +1 − u u − ui N i , p −1 (u ) + N i +1, p −1 (u ) , ui+ p − ui u i + p +1 − u i +1
⎧1 u i ≤ u ≤ u i +1 N i , 0 (u ) = ⎨ (4) ⎩0 otherwise
Typically, 3rd degree surface are used in CAD applications. The analysis with differential geometry methods requires the calculation of at least the 1st and the 2nd derivative of the surface. Furthermore, the influence of one point of the control net is defined by the order of the surface. To increase the influence region, higher degree surfaces have to be calculated. Therefore, our implementation of the NURBS surface algorithms is variable in degree, even if this becomes time-consuming. Additionally, we want to compensate the quadrilateral structure of the polygons in the initial control mesh to get a square structure. This is realized on the basis of the non-uniform knot vectors. The result is illustrated in Figure 3. 4.2. Calculation of Surface Features The analysis of parameterized geometrical surfaces is usually based on Gauss’ fundamental theorem surface theory. The parametric form of surface is defined with the three Cartesian coordinates x,y and z. A surface is defined by the relation r=r(u,v) with the real parameters u,v and the relations x=x(u,v), y=y(u,v) and z=z(u,v). For the determination of surface features and the other metric properties of a surface, there are three types of so-called fundamental forms. The most important are the first and second (since the third can be expressed in terms of these). The first Gaussian fundamental theorem for curved surface is explicitly given by the Riemannian metric (5). The second fundamental form is the symmetric bilinear form with respect to the tangent space of the first, and is given by relation (6). ds 2 = Edu 2 + Fdudv + Gdv 2
(5)
− dNdr = Ldu 2 + 2 Mdudv + Ndv 2
(6)
The Gaussian curvature K and the mean curvature H are one of the most important surface features. Depending on their values, one can derive information about the shape. For example, a positive curvature indicates bumps and negative values pits. The curvatures are defined based on the parameters of the first (E,F,G) and second (L,N,M) fundamental forms: K = ( LN − M 2 ) /( EG − F 2 ), H = (LG - 2FM + EN) /(2(EG - F 2 )) Two exemplary illustrations for the industrial application of high degree surface measurement are given in Figure 4. The position and head height of a rivet and the contours of a casting have been detected, segmented and measured on the basis of surface curvatures.
a) b) Figure 4: Mean curvature of a rivet on a plate (a) and the Gaussian curvature of a casting (b).
(7)
4.3. Summary and Future Work We presented a new and straightforward approach for the fast interpolation of scanned 3D surfaces by regular meshes. Based on the scan device properties and the measuring principle, we described a method to project the generated 3D points to a regular grid. The resulting meshes are used for a fast visualization to control the scanning process. Additionally, the meshes are registered to minimize and balance the uncertainties and noise caused by different sensor calibrations. Since, our approach bases on the structure of neighbouring scanlines, the mesh construction procedure can also be applied to other scan systems, which base on structured light. The performance and feasibility of the presented methods are evaluated in table 1 with respect to Figure 5. For the analyses we used a standard PC (Pentium4 with 512MB RAM). model total points scanlines meshes polygons total time time/mesh
woman 2.057.610 2936 8 871.196 2573ms 322ms
boot 1.246.994 1796 8 559.248 1683ms 210ms
pepper 2.135.745 2740 16 1.283.491 5804ms 363ms
casting 1.656.360 3160 8 1.216.437 2444ms 306ms
duck 1.628.515 2316 8 849.448 3341ms 418ms
can 1.681.967 2264 10 868.860 2713ms 271ms
Table 1: Performance evaluation for the construction of the basic and initial meshes with respect to the point clouds illustrated in Figure 5. Furthermore, we applied a piecewise NURBS surface approximation to the point cloud by using the initial meshes as control polygons. This allows precise and robust computation of surface features (e.g. different kinds of curvature) to evaluate the object. This work will continue, by applying fast methods to close small gaps in the initial mesh by nonlinear surface interpolation. Additional algorithms, that connect the quadrilateral meshes of several scans to one mesh with preserving the structure, are needed. Furthermore, the data exchange to other systems often requires one triangulation for all meshes. This will be realized by fast and optimized triangulation procedures using the given mesh neighbourhoods. References: [1] [2] [3] [4]
[5]
[6]
Alexa, M., J. Behr, D. Cohen-Or, S. Fleishman, D. Levin and C. T. Silva, “Point set surfaces,” IEEE Visualization 2001, S. 21–28, Okt. 2001. Amenta, N., Sunghee Choi, and Ravi Kolluri. The power crust. In Proc. of 6th ACM Symposium on Solid Modeling, pp. 249–260. ACM Press, 2001. Bernardini, F. and H. Rushmeier, “The 3D model acquisition pipeline,” Computer Graphics Forum 21, pp. 149–172, June 2002. Bernardini, F., J. Mittleman, H. Rushmeier, C. Silva, G. Taubin: “The Ball-Pivoting Algorithm for Surface Reconstruction”. IEEE Transactions on Visualization and Computer Graphics, pp. 349-359, IEEE Press, 1999 Bernardini, F., I. Martin, J. Mittleman, H. Rushmeier and G. Taubin, “Building a Digital Model of Michelangelo’s Florentine Pietà,” IEEE Computer Graphics & Applications 22, pp. 59–67, Jan. 2002. Boissonnat, J.-D. “Geometric structures for three-dimensional shape representation”. ACM Trans. Graph., 3(4):266–286, 1984.
[7]
[8] [9]
[10] [11] [12] [13]
[14]
[15]
[16] [17]
[18]
Carr, J. C., R. K. Beatson, J.B. Cherrie, T. J. Mitchell, W. R. Fright, B. C. McCallum, and T. R. Evans. “Reconstruction and representation of 3d objects with radial basis functions”. In ACM SIGGRAPH, pp. 67–76, 2001 Curless, B. and M. Levoy. “A volumetric method for building complex models from range images”. Computer Graphics, vol. 30:303–312, 1996. Dey, T. K., J. Giesen, and J. Hudson. “Delaunay based shape reconstruction from large data”. In Proc. IEEE Symposium in Parallel and Large Data Visualization and Graphics (PVG2001), pp. 19–27, 2001. Eck, M. and H. Hoppe, “Automatic Reconstruction of B-spline Surfaces of Arbitrary Topological Type,” Computer Graphics 30, pp. 325–334, Aug. 1996. Edelsbrunner, H. and Ernst P. Muecke. „Three-dimensional alpha shapes”. ACM Transactions on Graphics, 13(1):43–72, 1994. Floater, M. S. and M. Reimers. “Meshless parameterization and surface reconstruction”. Comp. Aided Geom. Design, 18:77–92, 2001. Gopi, M., S. Krishnan, and C. T. Silva. “Surface Reconstruction based on lower Dimensional Localized Delaunay Triangulation”. In Computer Graphics Forum (Eurographics 2000), volume 19(3), 2000. Hoppe, H., T. DeRose, T. Duchamp, J. McDonald, and W. Stuetzle. “Surface Reconstruction from Unorganized Points”. In ACM SIGGRAPH 1992, pp. 71–78. ACM SIGGRAPH, 1992. Levoy, M., K. Pulli, B. Curless, S. Rusinkiewicz, D. Koller, L. Pereira, M. Ginzton, S. Anderson, J. Davis, J. Ginsberg, J. Shade und D. Fulk, “The Digital Michelangelo Project: 3D Scanning of Large Statues,” in Proc. SIGGRAPH 2000, K. Akeley, ed., pp. 131–144, ACM SIGGRAPH, Addison Wesley, 2000. Lorensen, W. E. and H. E. Cline, “Marching Cubes: A High Resolution 3D Surface Construction Algorithm,” Computer Graphics 21, pp. 163–169, July 1987. Park, I. K., I. D. Yun and S. U. Lee, “Constructing nurbs surface model from scattered and unorganized range data,”. In Proc. of the 2nd International Conference on 3-D Imaging and Modeling (Ottawa, Canada, October 04–08, 1999), Annual Conference Proceedings, S. 312–340, National Research Council of Canada, 1999. Teutsch, C, T. Isenberg, E. Trostmann, M. Weber, D. Berndt and T. Strothotte. Evaluation and Correction of Laser-Scanned Point Clouds. In J.-A. Beraldin, S. F. ElHakim, A. Grün, and J. S. Walton, ed., Proc.of Videometrics VIII (Electronic Imaging 2005, January 16-20, 2005, San Jose, California, USA), vol. 5665, pp. 172-183, Bellingham, Washington, 2005. SPIE/IS&T.
a)
b)
c)
d)
e) f) Figure 5: Illustration of registered meshes in different colors at the example of the point clouds from the models of woman (a), boot (a), pepper(c), casting(d), can(e) and a duck (f).
