Reconstruction of freeform surfaces for metrology

The application of freeform surfaces has increased since their complex shapes closely express a product's functional specifications and their machining is obtained with higher accuracy. In particular, optical surfaces exhibit enhanced performance especially when they take aspheric forms or more complex forms with multi-undulations. This study is mainly focused on the reconstruction of complex shapes such as freeform optical surfaces, and on the characterization of their form. The computer graphics community has proposed various algorithms for constructing a mesh based on the cloud of sample points. The mesh is a piecewise linear approximation of the surface and an interpolation of the point set. The mesh can further be processed for fitting parametric surfaces (Polyworks® or Geomagic®). The metrology community investigates direct fitting approaches. If the surface mathematical model is given, fitting is a straight forward task. Nonetheless, if the surface model is unknown, fitting is only possible through the association of polynomial Spline parametric surfaces. In this paper, a comparative study carried out on methods proposed by the computer graphics community will be presented to elucidate the advantages of these approaches. We stress the importance of the pre-processing phase as well as the significance of initial conditions. We further emphasize the importance of the meshing phase by stating that a proper mesh has two major advantages. First, it organizes the initially unstructured point set and it provides an insight of orientation, neighbourhood and curvature, and infers information on both its geometry and topology. Second, it conveys a better segmentation of the space, leading to a correct patching and association of parametric surfaces.


Introduction
Surface reconstruction is an inverse problem. It is the process of creating an approximation of a surface based on a cloud of points measured by optical and tactile probes. Many researchers worked on this subject and proposed a large number of solutions and scientific publications during the last 30 years. Delaunay and Voronoi-based meshing techniques [1] as well as implicit techniques [2] have been developed within the computer graphics community. These techniques relate to the problem of surface reconstruction where only the cloud of points is known and concern applications such as video games, arts and reverse engineering. Nonetheless, when studied closely and tested, these methods are ineffective unless two important initial conditions are fulfilled. The cloud of points must verify a "good" sampling criterion and has to be free of noise [3]. The first Delaunay and Voronoi-based algorithms such as the sculpting algorithm of Boissonnat and the alpha shapes of Edelsbrunner and Mucke did not impose any of these conditions [4,5]. Hence, the resulting meshes were not strictly manifold. Amenta et al [6] presents proofs and guarantees for a "good" quality mesh reconstruction for the first time. Implicit techniques are based on a different philosophy and create an implicit representation of a surface. With the absence of model parameters, these approaches are not valid for metrology since no uncertainty budget can be established later on. Moreover, implicit surfaces techniques have no guarantees for surface reconstruction, meaning that no quality assessment can be done.
Freeform surfaces have seen enlarged applications since their complex shapes closely express a product's functional specifications and their machining is obtained with higher accuracy [7]. A freeform is a surface hardly defined by an equation or any generic form and is characterized by its asymmetric profile. According to ISO 17450-1 [8] freeform surfaces are complex geometrical features which have no invariance degree. However, these surfaces have proven their complex design functionalities [9]. This paper deals with datasets of an aspheric optical lens. The surface reconstruction algorithms tested on these datasets can have an extended usage on freeform surfaces. The datasets are measured at two metrology institutes by two different measurement strategies. LNE used a point-to-point measurement strategy on an XY grid. TNO (Netherlands Organisation for Applied Scientific Research) has measured the same optical lens on its high precision measuring machine Nanomefos with a spiral-like measurement strategy.
The 3D metrology of freeform surfaces is not yet normalized as no unified approach exists [10]. It differs from one application to another. Aspheric lenses inspection requires form characterization in some cases [11] [12] and simple profile characterization in some other cases. In form characterization the global surface is of interest. Whereas in profile characterization, only points along a cross-sectional plane at a specific location are considered. If not enough points can be found within the intersection of the cutting plane and the point set or if the level of the measurement noise is large, which is typically the case, surface reconstruction becomes mandatory. Either ways, form and/or dimensional characteristics are checked for conformity with respect to the design specifications.
This paper is focused on the study of meshes based on the Delaunay and Voronoi combinatorial structures. A mesh is a linear interpolation of the points that builds a structure on them. The initially unstructured point set becomes organized and structured. The mesh is a linear approximation of the underlying surface, which gives an insight of its topology and geometry. Discrete differential metrics can be calculated and used for further processing such as filtration, partition and association. Association is of various nature depending on the type of input provided with the cloud of points. If an analytical implicit model of the surface, such as the case of aspheric surfaces, is given, the association of that model to the structured points can be performed [11]. Whereas if no model is given, parametric models are fitted.

Delaunay and Voronoi structures
The Delaunay triangulation and Voronoi diagram are dual representations that infer complete knowledge about spatial neighbourhood. They build a structure that allows determining a point's exclusive neighbours. The Voronoi diagram is built based on this notion of proximity and cuts space into neighbourhood cells (3). Each point x in the space is assigned to one single point p in the dataset P as being in its exclusive neighbourhood. The Delaunay triangulation is a simplicial complex that connects the points that are exclusive neighbours among themselves. It is also known as a triangulation of the convex hull of the dataset. Triangulation is a broad word referring to tetrahedralization in 3D and triangulation in 2D, but we use it for simplicity. The simplices forming the 3D Delaunay complex are tetrahedra as opposed to triangles in two dimensional Delaunay complexes. Figure 1 shows how a Delaunay simplex is circumscribed in a circle empty of other points in the dataset. The general course of the surface reconstruction algorithms hereafter starts with the 3D Delaunay triangulation of the point set and then extracts a triangular surface mesh from it. Only triangles that approximate the underlying surface are selected in the output mesh. These triangles constitute the set of restricted Delaunay simplices. The Delaunay triangulation grants the most regular triangulation, as it builds the most equiangular triangles possible [13]. The duality between Voronoi and Delaunay is done according to their simplices dimensions. By taking n=3 as the dimension of the ambient space, a 0-simplex (vertex) of the Delaunay triangulation corresponds to a n-simplex (Voronoi cell) of the Voronoi diagram. Similarly, a 1-simplex (edge), 2-simplex (triangle) and nsimplex (tetrahedron) in the Delaunay triangulation, correspond to a (n-1)-simplex (Voronoi facet), a (n-2)-simplex (Voronoi edge) and a 0-simplex (Voronoi vertex) in the Voronoi diagram, respectively.

Expected mesh quality and algorithm requirements
The mesh of the dataset is a linear interpolation of the points, and is based on the Delaunay structure. The common approach for 3D datasets is to build a 3D Delaunay triangulation and extract triangular facets that are a linear approximation of the underlying surface. A good mesh quality is a mesh that approximates well the underlying surface and that connects the points correctly in accordance with the topology of the surface. The quality of the mesh sought is based on well-defined criteria. First, the reconstructed surface should be topologically equivalent to the underlying surface of the points set. This means that the reconstructed surface should be homeomorphic to a 2-manifold and to ℜ². In topology we usually refer to basic topological elements such as disks. A mesh is a discrete 2-manifold only if each vertex umbrella is homeomorphic to a 2-disk [13]. The umbrella of a point on a mesh is the set of triangles (and their vertices) that are incident to that point. Second, the reconstructed surface should also be geometrically equivalent to the underlying surface (figure 2). Both facet normals and surface normals should be similar in direction. Both representations should be geometrically close in terms of Hausdorff distance. The chord error between the facets and the surface should not exceed a given tolerance and be proportional to the sampling density. Surface reconstruction is the process of extracting the restricted Delaunay facets that approximate the first order of the surface out of the 3D Delaunay complex. Algorithmic complexity is also a major issue to our application. The reconstruction algorithm is to be optimized with regard to space and time complexities. Space, because very large data need to be processed, and time because the characterization of the parts is an on-line process that should not exceed -at least-the measurement time. In fact, the Delaunay triangulation algorithm is incremental so it can be applied incrementally with a sequential insertion scheme as the points are being measured. Furthermore, in the case of an on-line metrology application, an automatic reconstruction process is required. No human intervention shall be necessary.  Figure 2. Geometrical equivalence between mesh and theoretical surface

Alpha shapes
The Alpha shapes is a kick-off algorithm for surface reconstruction and is used as a data filter in applications such as in form fitting of freeform surfaces [5,11]. From the 3D Delaunay triangulation of the dataset, the α-shape is an extracted sub-complex formed by all simplices that have the property of being α-exposed. Figure 3 shows the difference between an α-exposed simplex and a non α-exposed simplex. The simplices considered are triangular facets. The algorithm passes a sphere of constant diameter α across the facets and checks for its content. Either the sphere is empty of points and the facet is α-exposed or the sphere contains points from the dataset and the facet is non α-exposed. All the facets that are α-exposed are selected and constitute the simplices of the output mesh. Alpha shapes are not designed to build manifold meshes but only an approximation of the shape inherent to the points. It can be shown that a non-uniform point distribution causes the result to be non manifold. Figure 3. 2D example, a constant diameter circle. Non α-exposed simplex (left), α-exposed simplex (right)

Cocone
The Cocone algorithm was developed by Amenta et al for 2D curves and expanded to surfaces in 3D [3,6]. Cocone is the first algorithm to provide geometrical and topological guarantees which ensure the output triangular mesh is a good approximation of the surface to be reconstructed. These guarantees are conditioned by an ε-sampling condition on the dataset which is defined with respect to the local feature size at a given location on the surface. The condition states that any point x taken on the surface should have at least one sample point within the range of ε times the perpendicular distance from x to the medial axis of the surface. This perpendicular distance is called the local feature size of x and denoted LFS(x). Therefore the ε condition requires that the medial axis of the surface exists and the surface is closed (2). Ultimately, it also implies that the sampling is denser in regions where the curvature is high, and inversely. .
The robustness of the Cocone algorithm highly depends on the ε-sampling condition. In the absence of it, the guarantees are dropped. The algorithm is based on the Delaunay triangulation and its dual representation, the Voronoi diagram of the point set. Once the Delaunay triangulation has been constructed, the algorithm works on the triangular facets constituting the (n-1) simplices of the triangulation. The idea derives from the fact that the triangular facets of the output mesh are a first order approximation of the underlying surface. Each facet is thus an approximation of the tangent plane to the surface at the facet's location up to a chord error. It follows that the normal to the facet is an approximation of the surface normal at the same location. With this property, Cocone searches for the restricted Delaunay facets by stating that these are within a given tolerance zone delimited by a double cone. So if, within the conical zone, a facet should exist, it means that this facet belongs to the restricted Delaunay complex and that it is part of the output mesh. The Cocone is defined by two geometrical elements and an angle: the axis, the apex and the opening angle. The cone axis is defined as being orthogonal to the approximate tangent plane. Let's assume that a facet f is a good approximation of the tangent plane to the surface at a certain location. The dual Voronoi element to a Delaunay facet is an edge that is orthogonal to the facet. Therefore, the dual Voronoi edge to f is in theory both, the approximation of the normal to the surface, and the axis of the Cocone of f. Since only some Voronoi edges are in the correct direction, the algorithm seeks those ones which are directionally close to the vector that connects a query point p to its pole. The poles of a Voronoi cell are the farthest two Voronoi vertices the cell. Amenta et al [6] have proven that these poles approximate the normal directions inwards and outwards at p. Consequently, for a given facet and its dual Voronoi edge, the Cocone test is applied. It consists of simulating the plot of a double cone at each vertex of the facet (figure 4). The normal approximation of a facet is approximated by its three vertices normals. If the Cocone condition is verified for all three vertices of a facet, the facet is selected in the output mesh.

Other Methods
Other methods that were developed are also based on the ε-sampling condition. The Natural Neighbours Interpolation of Boissonnat and Cazals [14] follows, by some means, the same principle. It also seeks the restricted Delaunay facets and the dual Voronoi edges that approximate the normal direction. However, the detection of Voronoi edges is different. The algorithm calculates an implicit distance function on the points. This function would be the association of a scalar to a three dimensional point. The function returns a zero at the points of the dataset, a negative value inside the surface and a positive value outside. The Voronoi edges that have mixed vertex signs are edges that cross the implicit surface and are considered as a good approximation of surface normals. Figure 5 illustrates the logic behind this algorithm and shows Voronoi edges that are in the configuration sought and their dual restricted Delaunay elements. The Wrap algorithm of Edelsbrunner calculates a distance function on the points function and then builds a flow relation previously shown that there exists a duality between critical points and the Delaunay and Voronoi simplices and their intersections. Based on clusters. Some clusters are stable eliminated. The Wrap algorithm does not provide clusters level. The works of [17] and [ demonstrated that Wrap can supply sampling condition is satisfied.
Dey et al [19] continued working on the Cocone algorithm and developed what is known by enhanced Cocone algorithms. They datasets are processed via an Octree subdivision of the points space. algorithm with provable guarantees for noisy data. In theory, the noise is limited to a given limit proportional to that of the sampling the same non noisy dataset. The reason is that RobustCocone that do not verify the noise level condition [21] puts forward a local approach triangulation. In fact, even the Delaunay triangulation is localized. Points are clustered into cells. The Delaunay triangulation is built piecewise, cell by cell, by considering the points of of a portion of the neighbouring cells.
Boissonnat and Ghosh proposed an algorithm based on tangential Delaunay complexes Instead of natural neighbours, tangential neighbours are computed and the approx restricted simplices follows this kind of localized neighbourhood. The algorithm assumes the knowledge of normal directions and that the local feature size at any point is positively defined. these additional constraining conditions

Comparison on a simple 2D
The sculpture algorithm as well as the alpha shapes and Cocone are compared in 2D. Matlab codes are generated on the same simple dataset showing the difference in algorithmic behav shows a random dataset representing a curve in 2D the set. The algorithms are all applied on the same dataset.

2D example of the Natural Neighbours Interpolation algorithm
The Wrap algorithm of Edelsbrunner [15] follows a different approach to what preceded. It calculates a distance function on the points, identifies critical points through the gradient of that then builds a flow relation following the steepest ascent of the gradient. Gro that there exists a duality between critical points and the Delaunay and Voronoi . Based on this fact, flow relations congregate Delaunay simplices stable manifolds and some others are unstable manifold The Wrap algorithm does not provide global guarantees, only locally at the stable manifold ] and [18] have linked Wrap to the ε-sampling theory that Wrap can supply global topological and geometrical guarantees, again only if the continued working on the Cocone algorithm and developed what is known by They discussed a version of the Cocone algorithm in which large datasets are processed via an Octree subdivision of the points space. Dey et al algorithm with provable guarantees for noisy data. In theory, the noise is limited to a given limit sampling. Practically, it performs poorer than the basic Cocone algorithm on dataset. The reason is that RobustCocone filters out tetrahedra associated to points that do not verify the noise level condition which is in turn very tight. The last publication on Cocone ] puts forward a local approach to algorithmically reconstruct a surface from the 3D Delaunay triangulation. In fact, even the Delaunay triangulation is localized. Points are clustered into cells. The Delaunay triangulation is built piecewise, cell by cell, by considering the points of of a portion of the neighbouring cells.
Boissonnat and Ghosh proposed an algorithm based on tangential Delaunay complexes Instead of natural neighbours, tangential neighbours are computed and the approx this kind of localized neighbourhood. The algorithm assumes the knowledge of normal directions and that the local feature size at any point is positively defined.
conditions make of this algorithm a specific one.
point set he sculpture algorithm as well as the alpha shapes and Cocone are compared in 2D. Matlab codes are generated on the same simple dataset showing the difference in algorithmic behav shows a random dataset representing a curve in 2D and figure 6b shows the Delaunay triangulation of . The algorithms are all applied on the same dataset.

2D example of the Natural Neighbours Interpolation algorithm
] follows a different approach to what preceded. It through the gradient of that he gradient. Grove [16] has that there exists a duality between critical points and the Delaunay and Voronoi congregate Delaunay simplices in manifolds that need to be , only locally at the stable manifold theory and by that have topological and geometrical guarantees, again only if the εcontinued working on the Cocone algorithm and developed what is known by a version of the Cocone algorithm in which large et al [20] proposes an algorithm with provable guarantees for noisy data. In theory, the noise is limited to a given limit than the basic Cocone algorithm on out tetrahedra associated to points The last publication on Cocone to algorithmically reconstruct a surface from the 3D Delaunay triangulation. In fact, even the Delaunay triangulation is localized. Points are clustered into cells. The Delaunay triangulation is built piecewise, cell by cell, by considering the points of a cell and the points Boissonnat and Ghosh proposed an algorithm based on tangential Delaunay complexes [22]. Instead of natural neighbours, tangential neighbours are computed and the approximation of Delaunay this kind of localized neighbourhood. The algorithm assumes the knowledge of normal directions and that the local feature size at any point is positively defined. So he sculpture algorithm as well as the alpha shapes and Cocone are compared in 2D. Matlab codes are generated on the same simple dataset showing the difference in algorithmic behaviour. Figure 6a shows the Delaunay triangulation of  Figure 6d shows the plot of circles of diameter α=0.85 tested at each segment. Delaunay edges are traversed because it is a 2D application. The circles that are empty of other points from the dataset and that cross a Delaunay edge are called α-exposed edges. These edges constitute the α-shape of the points shown in figure 6e. It is to notice how the α-shape is not necessarily manifold. The Cocone result shown in figure 6f is manifold. Cocone reconstructs the curve only because it respects the ε-sampling condition. In the central region where the curvature is high a high point density is required, unlike low curvature regions. Based on the above comparison, the comparison in table 1 and the abundance of the Cocone citations in literature, the Cocone algorithm best fits our application. Guarantees and algorithmic complexity are the two major criteria on which the comparison is made. Another decisive criterion is the possibility to process large data. The Wrap algorithm requires that a flow relation is calculated over the Delaunay simplices and this is a relationship-based computation. Simplices cannot be treated independently. So provided that the number of Delaunay simplices is very large, the algorithm comes to a crash for excessive memory load. Unlike Wrap, the Natural Neighbours Interpolation algorithm calculates a distance function over the points and then treats the Voronoi edges separately. However, this algorithm works only on closed surfaces since it requires the knowledge of inside and outside spaces with respect to the surface. With Cocone algorithm, each facet and its dual Voronoi edge undergo the Cocone test one by one which gives it the advantage of hashing and parallelization. In addition, it does not need any information regarding inside and outside spaces. The only input, which can be generalized to an acceptable value, is the angle Θ of the Cocone.

Example of application of the Cocone method to measured optical lens surfaces
The idea here is to investigate surface reconstruction methods and to apply an algorithm that would build a data structure on the initially unstructured 3D point sets and that would transform the discrete point clouds into continuous piecewise linear representations. For applications such as profile metrology, it is crucial to have a continuous representation. Profile metrology is an application where dimensional characteristics are estimated on an extracted cross section of the surface. Yet, another application to surface reconstruction relates to ray-tracing simulation of optical parts. An this application requires to have a mesh representation [23].
To be representative of most common applications in 3D metrology of complex optical surfaces, an optical lens measured in two different strategies is considered. The optical lens is an asphere defined by the analytical equation (3).
The first set of measured points on an XY grid basis contains 2,782,224 points over a squared area of 5ᵡ5mm. The measurement is done using the reference profilometer developed at LNE. It ensures three degrees of freedom, x, y and z linear motions and perfectly respects the Abbe principle. Its metrology loop is optimized to be as short as possible. The part only moves in the XYplane and the probe moves only along the z axis. The second set of points contains more than 1,000,000 points covering the entire surface and is measured using the Nanomefos machine. The part is mounted on a rotating table. An optical sensor is focused at the centre of the part and while the latter undergoes rotation, the measuring head moves in the radial direction from the centre outwards.
( ) The Cocone algorithm was implemented and tested on both datasets. These datasets are known to be excessively dense. Figures 7 and 8 (left) show the failure of this algorithm for the fact that the ε-sampling condition is not fulfilled on the datasets. Figures 7 and 8 (right) show that when the sampling respects the ε condition, the reconstruction is proper. The reconstruction algorithms cited in the previous section are very sensitive to the ε-sampling condition. The implementation of the Cocone algorithm has put into evidence the limitation related to having overly dense sampled data.  Optical surfaces, whether they have aspheric or freeform shapes, are usually almost flat surfaces. Irrespective of the measurement strategy, the resulting cloud of points is a set of points that can be projected onto a plane following a bijective mapping without any superposition of points from different sides of the surface. Non bijective projections result from freeform surfaces that represent superposition of points once projected onto a plane. An example of such surfaces is turbine blades. The bijection property offers the advantage of tracing back the points to their original position without modifying either the geometry or the topology of the underlying surface. Since topology is preserved on the map, neighbourhood is conserved and points can thus be meshed in a 2D-like fashion. The application of a 2D Delaunay triangulation is the best solution as it ensures the most regular and complete triangulation possible. The data structure is built in 2D and then mapped back to 3D as simplified in figure 9. The links created among the points are preserved. Figures 10 and 11 show the resulting triangulation based on this technique on both datasets without any filtering or modification of the points. The reconstruction is proper.
With the data structure created on the points, the association of the dataset to the theoretical model of aspherical surfaces is improved. Normal and curvature estimations are more precise since they are calculated based on the true neighbourhood of a point which is the neighbourhood found on the manifold and not in the 3D neighbourhood of the point [24].

Conclusion
In this paper, a comparison of surface reconstruction algorithms is presented and lead to the choice of Cocone which fits at best our application. The chosen algorithm possesses inherent limitations related to sampling density. The tests run on real measured optical lenses highlighted the necessity of having an ε-sampled point set. If the cloud of points does not respect the ε condition, a filter is required in order to reduce the overly dense dataset. With surface reconstruction, the underlying surface can be approximated by a first order representation called the mesh. The mesh builds a data structure for initially unstructured points and most importantly creates a continuous representation of the surface. This is required whenever profile extraction and dimensional metrology are to be applied at some specified location on the surface. The mesh also offers the possibility to perform partition of the mesh vertices and thus the possibility to either apply reverse engineering or geometric feature extraction. However, the mesh is used here to structure data and to create a continuous representation of the underlying surface. A mesh is a more substantial structure to approximate normals and curvatures. Curvature estimation is an important tool for surface defects detection and for filtration which can be applied on the mesh vertices in order to reduce their number. Consequently, fitting becomes less complex. A future extension of this work will cover freeform surfaces such as freeform optics or turbine blades.