|
CGAL 4.9 - Advancing Front Surface Reconstruction
|
|
CGAL 4.9 - Advancing Front Surface Reconstruction
|
Surface reconstruction from an unstructured point cloud amounts to generate a plausible surface that approximates well the input points. This problem is ill-posed as many surfaces can be generated. A wide range of approaches have been proposed to tackle this problem. Among them are variational methods [12][13], tensor voting [10], implicit surface [9][6], and Delaunay triangulations.
For Delaunay based algorithms the output surface is commonly generated as the union of some triangles selected in the 3D Delaunay triangulation of the input points. Such algorithms are either volume-based by generating as output the boundary of selected tetrahedra [3][4], or surface-based by selecting a set of triangles.
In most surface based Delaunay algorithms the triangles are selected independently, that is in parallel [1][2].
This chapter presents a surface-based Delaunay surface reconstruction algorithm that sequentially selects the triangles, that is it uses previously selected triangles to select a new triangle for advancing the front. At each advancing step the most plausible triangle is selected, and such that the triangles selected generates an orientable manifold triangulated surface.
Two other examples of this greedy approach are the ball pivoting algorithm and Boyer-Petitjean's algorithm [5][11]. In both algorithms a triangulated surface is incrementally grown starting from a seed triangle. Ball pivoting is fast, but the quality of the reconstruction depends on user defined parameters corresponding to the sampling density. The Boyer-Petitjean approach can handle non-uniform sampling, but fails when near co-circular points are encountered, and it does not provide any guarantee on the topology of the surface.
We describe next the algorithm and provide examples.
A detailed description of the algorithm and the underlying theory are provided in [8].
The first step of the algorithm is the construction of a 3D Delaunay triangulation of the point set. The radius of a triangle \( t \) is the radius of the smallest sphere passing through the vertices of \( t\) and enclosing no sample point. In other words, the radius \( r_t\) is the distance from any vertex of \( t\) to the Voronoi edge dual to \( t\). This triangle with three boundary edges is the initial triangulated surface, and its boundary is the advancing front. The Delaunay triangle with the smallest radius is the starting point for the greedy algorithm.
The algorithm maintains a priority queue of candidate triangles, that is of valid triangles incident to the boundary edges of the current surface. The priority is the plausibility. While the priority queue is not empty, the algorithm pops from the queue the most plausible candidate triangle and adds it to the surface. New candidate triangles are pushed to the priority queue when new boundary edges appear on the advancing front. As the algorithm creates a two-manifold surface some candidate triangles can not be selected due to topological constraints which are explained next.
Any triangle \(t\) considered as the next potential candidate shares an edge \(e\) with the front of the current reconstruction. Let \(b\) be the vertex of \(t\) opposite to \(e\). There are four configurations where \(t\) is added to the surface.
While the first three operations never induce a non-manifold edge or vertex, we only can perform gluing, if triangle \(t\) has a twin facet, that is a triangle with an edge on the front and incident to \(b\), and the third vertex on edge \(e\).
A triangle is said valid when the above operations can be applied.
Valid triangles for an edge on the front are compared through their radius. While the radius is a good criterion in the case of 2D smooth curve reconstruction [7], we need another criterion for 3D surface reconstruction, namely the dihedral angle between triangles on the surface, that is the angle between the normals of the triangles. There are three bounds namely \( \alpha_\mathrm{sliver} \), \( \beta \), and \( \delta \).
The candidate triangle of an edge \( e \) is the triangle with the smallest radius:
There may be no such triangle. In the implementation of the algorithm \( \alpha_\mathrm{sliver} \) and \( \delta\) are equal to \( 5\pi/6 \).
We denote by \( \beta_t\) the angle between the normal of a triangle \( t\) incident on a boundary edge \( e \) and the normal of the triangle on the surface incident to \( e \).
We define the plausibility grade \( p(t) \) as \( 1/r_t \), if \( \beta_t < \beta \), and \( -\beta_t \) else. The parameter \( \beta \) can be specified by the user and is set by default to \( \pi/6\).
Let's have a look at the figure below.
\( \alpha_\mathrm{sliver}\) corresponds to the red wedge. The algorithm will never select triangle t1 even if it is the only candidate triangle.
\(\beta\) corresponds to the green wedge. If there is a candidate triangle in this zone, the one with the smallest radius is the most plausible.
If there is no candidate triangle in the green wedge, the triangle with the smallest angle between its normal and the normal of t' is chosen. In the figure above this would be triangle t4.
By construction the output of the algorithm is a connected orientable manifold with or without boundary. To cope with multiple components we merely look for a new seed facet among facets disjoint from the surface. In case of noisy data or outliers, the user must filter out small surface components.
It is impossible to handle all kinds of boundaries and non uniform sampling at the same time, as a void can either be an undersampled area of the surface, or a hole.
As we do not want the algorithm to rely on a uniformity condition on the sampling it will fill holes cut off from "flat" regions of the surface. However, in many cases a boundary component cannot be closed by adding a spanning disk such that the resulting disk is well sampled. Typically, closing a boundary component due to a transversal clipping of the operation, would yield large dihedral angles at boundary edges. Moreover, if the boundary is sufficiently well sampled, the radii of the two triangles incident on a boundary edge would be very different. These heuristic facts can be used for boundary detection.
More specifically, we discard any candidate triangle \( t \), for an edge \( e \) such that \( p(t) < 0\), and \( r_t > \mathrm{radius\_ratio\_bound} \times r_{t'}\) where \( t'\) is the triangle on the surface incident on \( e \). The parameter \(\mathrm{radius\_ratio\_bound}\) is specified by the user and is set by default to 5.
For the example given in Figure 56.2, we said that if there was no triangle t3 in the green wedge, triangle t4 would be chosen as it has the smallest angle between its normal and the normal of triangle t'. However, in case its radius was \(\mathrm{radius\_ratio\_bound}\) times larger than the radius of triangle t', triangle t2 would be chosen, assuming that its radius is not \(\mathrm{radius\_ratio\_bound}\) times larger.
Note that this heuristic implies that where the sampling is too sparse with respect to curvature, it must be sufficiently uniform for our algorithm to work.
The first of the following three examples presents a free function for doing surface reconstruction. For a sequence of points the function produces a sequence of triplets of indices describing the triangles of the surface. The second example presents a class that enables to traverse the surface represented in a 2D triangulation data structure where the faces are connected with the facets of underlying 3D Delaunay triangulation. The third example shows how to get outliers and the boundaries of the surface.
The global function advancing_front_surface_reconstruction() takes an iterator range of points as input and writes for each face of the reconstructed surface a triplet of point indices into an output iterator. The following example writes the output triangulated surface to std::cout in accordance to the OFF format.
The function has an overload with an additional argument that allows to choose how to prioritize facets. It can be written in a way to avoid the generation of triangles with a perimeter larger than a given bound.
File Advancing_front_surface_reconstruction/reconstruction_fct.cpp
While the first example just writes index triples, the second example uses as output iterator a wrapper around a reference to a Surface_mesh and calls the function add_face().
File Advancing_front_surface_reconstruction/reconstruction_surface_mesh.cpp
The class Advancing_front_surface_reconstruction provides access to a 2D triangulation data structure describing the output surface. The latter can be explored by hopping from a face to its neighboring faces, and by hopping from faces of the 2D triangulation data structure to corresponding facets of the underlying 3D Delaunay triangulation.
The type of the 2D triangulation data structure describing the reconstructed surface is the nested type Advancing_front_surface_reconstruction::Triangulation_data_structure_2.
The type Advancing_front_surface_reconstruction::Triangulation_data_structure_2::Vertex is model of the concept TriangulationDataStructure_2::Vertex and has additionally the method vertex_3() that returns an Advancing_front_surface_reconstruction::Vertex_handle to the associated 3D vertex.
The type Advancing_front_surface_reconstruction::Triangulation_data_structure_2::Face is model of the concept TriangulationDataStructure_2::Face and has additionally the method facet() that returns the associated Advancing_front_surface_reconstruction::Facet, and a method is_on_surface() for testing if a face is part of the reconstructed surface.
In case the surface has boundaries, the 2D surface has one vertex which is associated to the infinite vertex of the 3D triangulation.
The underlying 3D Delaunay triangulation can be accessed as well, using the API of the class Delaunay_triangulation_3.
The following example writes the surface to std::cout in accordance to the STL (Stereo Lithography) format.
File Advancing_front_surface_reconstruction/reconstruction_class.cpp
Input points which are not on a surface are outliers. The member function outliers() returns an iterator range of those points.
Boundary edges can be traversed with the member function boundaries() It returns an iterator range type Boundary_range whose iterators have the value type Vertex_on_boundary_range. This is again an iterator range whose iterators have the value type Vertex_handle.
File Advancing_front_surface_reconstruction/boundaries.cpp
1.8.4