CGAL 5.3.1 - Scale-Space Surface Reconstruction
User Manual

Authors
Thijs van Lankveld

This CGAL package implements a surface reconstruction method which takes as input an unordered point set and computes a triangulated surface mesh interpolating the point set. We assume that the input points were sampled from the surface of an object. The method can also process point sets sampled from the interior of the object, although we cannot provide guarantees on the output. This method can handle a decent amount of noise and outliers. The point set may greatly undersample the object in occluded regions, although no surface will be reconstructed to fill these regions.

knot2_pts_surf.png
Figure 62.1 Left: 5760 points on a synthetic knot data set. Right: reconstructed surface mesh.

A triangulated surface mesh is generated by first computing the point set at a coarse scale, then constructing a mesh of the point set at this scale, and finally reverting the points of the mesh back to their original scale.

Note
A detailed tutorial on surface reconstruction is provided with a guide to choose the most appropriate method along with pre- and post-processing.

Scale-Space

The problem of surface reconstruction from a point sample is ill-defined. Because the points only provide a sample of the surface, we cannot determine with exact certainty how the surface behaves between the points. We may get more information about the surface by sampling the surface further, but even then the problem remains on a smaller scale.

For real-world data sets we also have to contend with noise and outliers, small errors in the point location and points that do not sample the surface respectively. These are generally caused by faults in the measurement procedure. Unless we actually know the surface sampled, we cannot determine whether a distinctive point is an outlier, is caused by noise, or samples a small feature of the surface.

To overcome these obstacles, we model the point set at different scales using a scale-space [1]. For our purpose, a scale is an abstraction level of the point set. At a lower, finer scale the point set describes smaller features of the object. Conversely at a higher, coarser scale the point set only describes larger features. A scale-space describes the point set at a dynamic scale. It is a four-dimensional space with the scale as fourth dimension. Increasing the scale is similar to smoothing the underlying surface described by the point set. Figure 62.2 shows an example point set at different scales. The overall effect of modeling this scale-space, is that the surface reconstruction problem is more well-defined at higher scales.

scale_space_pts.png
Figure 62.2 The black point set contains a smooth region, a region with large features, and a region with either small features or noisy points. The blue points show the point cloud at a higher scale.

Small changes in scale result in small changes in the points set that depends on the local shape. As the scale increases, the geometry approaches a number of points based on the connected components. Given a set of points at a fixed scale, the point set at a higher scale can be determined unambiguously, but this is not the case for the point set at a lower scale, as demonstrated by Figure 62.3

line.png
Figure 62.3 The black points and crosses show two different point sets. The blue boxes show both point sets at a higher scale. Note that given just the blue boxes, we cannot determine its point set at a lower scale.

The scale-space surface reconstruction method works in two phases. Firstly, a scale-space is constructed from the point set and an intersection of this space at a higher scale \( s \) is computed. Then, a triangulated surface mesh of the point set at scale \( s \) is computed. Neighboring triangles in the mesh have consistent orientation, as expressed using the right-hand rule on the ordered vertices of the triangle.

Generally, a reconstruction interpolating the original point set is desired. This can be achieved by reverting the points of the mesh to their original scale, i.e. their original locations. Figure 62.4 shows example reconstructions from the point set at its original scale and at a higher scale with reverted mesh.

scale_space_surf_1.png
scale_space_surf_2.png

Figure 62.4 The surfaces constructed from the point set at the original scale (left) and at a higher scale (right). The surface points in the right figure are reverted back to the original scale. The dashed paths shows the inward-facing shell where it does not overlap the outward-facing shell. Note how the scale-space reconstruction assumes the points on the right side of the object sample small features.

Both the smoothing operator and the mesh reconstruction assume that points near each other belong to the same part of the object. This is usually expressed in the notion of balls with a fixed size, the neighborhood radius. If such a ball contains multiple points, these points are near each other and will influence each other while increasing the scale. If such a ball is empty, it lies outside the object. Note that outside is based on regions empty of points, not on whether a volume is enclosed by the surface.

The neighborhood size can be estimated through statistical analysis. We use a kD-tree to estimate the mean distance to the n-th nearest neighbor and we use this distance as an approximator for the resolution.

The point set at the initial scale is equivalent to the input point set. In theory the scale can be varied continuously, but in practice the scale is increased in discrete iterations for efficiency reasons. The scale is increased one iteration by transforming each point using a smoothing operator. CGAL provides two smoothing operators for scale space reconstruction:

  • CGAL::Scale_space_reconstruction_3::Weighted_PCA_smoother (default) uses the density-weighted principal component analysis (PCA) of the local neighborhood. If the point set was sampled from a surface for which any unwanted deformation and sampling noise is smaller than the neighborhood size, the scale is coarse enough for mesh reconstruction after a few iterations of increasing the scale.
  • CGAL::Scale_space_reconstruction_3::Jet_smoother uses the function CGAL::jet_smooth_point_set() that projects the point to a local smooth (jet) surface. This smoothing is less aggressive and should be used if the point set is not very noisy and if a higher precision is sough.

Users can define their own smoothing operators by following the concept CGAL::Scale_space_reconstruction_3::Smoother.

Meshing is achieved by interpolating the smoothed point cloud and propagating back the connectivity to the original point cloud. CGAL provides two meshing operators for scale space reconstruction:

The method provides access to intermediate results and users can adjust these to better suit their needs. The (intermediate) results are the estimate of the resolution, the scale, and the final collection of surface triangles.

Note that Jet_smoother and Advancing_front_mesher are methods that rely upon other CGAL packages and that do not correspond to the original scale space algorithm. They can be used as alternative operators for the cases of point clouds with low noise that sample surfaces with boundaries. For more information on their parameters and effect, please refer to their respective manual pages (Smoothing and Chapter_Advancing_Front_Surface_Reconstruction). The rest of this documentation focuses on the default behavior using Weighted_PCA_smoother and Alpha_shape_mesher.

Parameter Settings

The scale-space surface reconstruction method has two main global parameters: the radius of the neighborhood ball and the number iterations of increasing the scale. If no appropriate neighborhood radius is known, this can be estimated using another two parameters: the mean number of neighbors and the number of samples. We have empirically determined values for each of these parameters that work well on a broad spectrum of data sets. However, we advise to carefully fine-tune these parameters for each type of data set.

The scale-space reconstruction operates locally. Points within a local neighborhood influence each other when increasing the scale. Similarly, points near each other in space are more likely to be near each other on the reconstructed surface. Therefore, it is important to have a good indication of which points are near each other. This proximity is expressed in the neighborhood ball radius parameter.

The neighborhood ball radius is related to the point density and the object thickness. In the ideal case, the radius is chosen such that any ball \( B \) with this radius centered on any point of surface of the object will contain several points that should be connected by a local patch in the surface to reconstruct. At the same time, \( B \) should not contain any points that should not be part of this local patch, such as points on the other side of the object.

The neighborhood radius can be estimated automatically from the point set. In this case, the radius is estimated such that on average it contains a given number of points. The number of points should be chosen such that the conditions on the neighborhood radius described above are preserved. Generally, an average of 30 neighbors gives good results. This number may be increased if parts of the object have a very low point density, it may be decreased if thin features of the object should be reconstructed better.

The radius estimation is based on a number of sample points. This number of samples is related to how regularly the point set covers the object. A large sample size may be required if the object is covered very irregularly. However, a larger number of samples may require more computation time. We generally use 200 samples to estimate the neighborhood radius.

The scale-space method tries to make the surface reconstruction problem less ill-posed by increasing the scale to an appropriately coarse scale. At this coarse scale the point set samples a smoothed surface, which is easier to reconstruct. The number of scale increase iterations determines how much this surface is smoothed.

The number of iterations is related to the noise in the point set, the acuteness of surface features, and the thickness of the object. Point sets with a lot of noise and objects with sharp or small features will require more iterations. At the same time, processing too many iterations can degenerate a volume into a plane. These degeneracies may cause the reconstructed surface to connect points on opposite sides of the object. Generally, 4 iterations of increasing the scale are appropriate.

Output Surface

The surface mesh constructed at scale \( s \) is non-self-intersecting. The interiors of any two triangles cannot pairwise intersect in a line segment. However, the surface does not need to be 2-manifold. An edge may be incident to more than two triangles and triangles may overlap exactly if large regions on both sides of the triangle are empty of points. Note that we count overlapping triangles with opposite orientations as separate triangles. In many cases where the points sample the surface of an object, the computed surface will contain both an outward-facing and a similar inward-facing surface, with a thin volume between them.

The surface mesh will not have edges incident to only one triangle or holes, loops of such edges, and the triangles are all oriented away from the point set. If the point set has holes, it is likely that the surface mesh will contain overlapping triangles with opposite orientation touching this hole.

An additional treatment can be applied to force the output surface to be 2-manifold: small flat volumes are found before reconstructing and the surface is forced to only use one side of the volume while keeping consistency with the adjacent singular facets. Only one side of each singular facet is used. In addition, non-manifold edges and vertices are removed as a post-processing. The facets that have been discarded to make the surface 2-manifold are stored and accessible.

If the object is not densely sampled or has disconnected components, the reconstructed surface may have several disconnected components. The surface is either an unordered collection of triangles, or the same collection sorted per shell. A shell is a collection of connected triangles that are locally oriented towards the same side of the surface.

When reverted to the original scale, we cannot guarantee the surface is a valid embedding, because the triangles of this surface may pairwise intersect in their interior. It may also have boundary edges, although these will always be incident to two surface triangles oriented in opposite directions. However, when using appropriate parameter settings for the number of iterations and neighborhood size the surface will generally not self-intersect. The appropriate parameter settings depend on the geometry of the point set and generally need to be fine-tuned per data set, as described later.

Software Design

The main classes Scale_space_surface_reconstruction_3, Weighted_PCA_smoother and Alpha_shape_mesher contain all the functionality to estimate the neighborhood size, compute the scale-space and increase the scale, and reconstruct the surface from the point set at the current scale.

The neighborhood size is estimated using Orthogonal_k_neighbor_search. The point set is generally stored in a Orthogonal_k_neighbor_search::Tree. When the neighborhood size is estimated, this tree is searched for nearest neighbors.

The scale-space is constructed at the original scale of the points. An iteration of increasing the scale is computed using a weighted PCA procedure. As described by Digne et al. [1], unlike similar methods this procedure does not lead to an undesirable clustering effect. By default the efficient Eigen library is used for this procedure if available. Otherwise, the internal fallback Diagonalize_traits is used. It is also possible to provide your own model for the DiagonalizeTraits concept. The weighted PCA procedure is performed locally per point, so it can be performed with parallel computing (linking with Intel TBB and passing the Parallel_tag to the reconstruction class is required).

The mesh reconstruction is performed by filtering a 3D alpha shape of the point set at a fixed scale. This filtering constructs a triangle for each regular facet; each singular facet results in two triangles facing opposite directions.

The resulting surface is stored as a collection of triples. Each triple contains three indices to the point set. The scale iterations do not change the order of the points. This means these indices can be applied on the point set at its original scale to get an interpolating surface, called reverting the scale of the surface. The indices may also be applied to the point set at the current scale to get a smoothed surface. This smoothed surface will generally have shrunk somewhat compared to the reverted surface.

To reduce memory requirements, the point set is not stored at its original scale. Instead, the stored point set at the current scale is a clone of the input. This means that the original input can be used to revert the surface.

Case Sets and Performance

There is a wide variety of possible point sets, whether synthetic or measured. The scale-space surface reconstruction method attempts to make the reconstruction problem less ill-posed. However, we cannot guarantee good results for arbitrary point set. Some point sets require careful parameter setting, while for others this method is ill-suited altogether. Here, we show different cases including the ideal sampling and different problem cases.

Ideal Conditions

The scale-space surface reconstruction has lenient ideal conditions. The surface should be reasonably regularly sampled and it should not contain overly concave corners. These notions affect the neighborhood parameter in opposing ways: more regular sampling may admit a smaller neighborhood parameter, which will admit smaller concavities. The point set may contain noise that is significantly smaller than the neighborhood size. Even so, when using the method to construct a mesh interpolating the original point set, any noise will deform the mesh. Figure 62.5 shows a reconstruction under ideal conditions.

retinal_pts.png
retinal_surf.png

Figure 62.5 Left: In the ideal case, points are regularly distributed and concavities are limited. Right: reconstructed surface mesh.

Boundaries

Surfaces with boundaries are handled naturally. However, it is important to note that in these cases, the surface does not explicitly identify the boundary and loops around it. In many cases, this will lead to boundary edges being shared by neighboring triangles of the surface that are identical except for their orientation. The boundaries may be made explicit by identifying orientation discrepancies between neighboring triangles, but this functionality is not provided by this package. Figure 62.6 shows a reconstructed surface with a boundary.

mushroom_pts.png
mushroom_surf.png

Figure 62.6 Left: A point sample of a surface with boundary. Right: reconstructed surface mesh.

Note that the mesher Advancing_front_mesher is a good alternative for producing well oriented surfaces with boundaries.

Sampling Regularity

As described above, ideally the point set regularly samples the surface. For example, for each point its six nearest points lie at roughly equal distance forming a rough hexagon around the point. This will generally not occur when working with realistic data. The scale-space method can reconstruct very unevenly sampled surfaces, as Figure 62.7 shows. However, as noted earlier the sampling regularity and allowed concavity form a trade-off. Figure 62.8 shows a case where we cannot simultaneously model the concave regions and sparse regions correctly.

femur_pts.png
femur_surf.png

Figure 62.7 Left: A very irregularly sampled point set. Right: reconstructed surface mesh.

lion-head_pts.png
lion-head_surf.png

Figure 62.8 Left: A very irregularly sampled point set with tight concavities. Right: reconstructed surface mesh. Note we cannot correctly model the mouth and eyes as well as the neck.

Sharp Features

The scale-space surface reconstruction method was not designed to deal specifically with surfaces with sharp features. However, it reconstructs these features adequately when using appropriate parameter settings, as demonstrated in Figure 62.9. The main prerequisite for correctly reconstructing the sharp features is that they contain sufficient sample points on the features and the features are not too thin compared to the neighborhood parameter.

fandisk_pts.png
fandisk_surf.png
fandisk_2_surf.png
fandisk_2_surf_detail.png

Figure 62.9 A point set on a surface with sharp creases and the reconstructed surface mesh seen from different sides. Note that the sharp features are mostly reconstructed correctly. However, in the thin regions circled in red the feature is not reconstructed correctly. The trade-off between regularity and concavity also means a hole appears in the third figure, bottom right.

Topology

The scale-space surface reconstruction method does not take into account the topology of the surface. This means the reconstructed surface may have a different topology than the surface sampled to obtain the point set. Figure 62.10 shows an example point set that was sampled from a surface with the topology of a sphere; the reconstructed surface has the topology of a torus.

man_pts.png
man_surf.png

Figure 62.10 Left: An irregularly sampling of a surface with the topology of a sphere. Right: reconstructed surface mesh with the topology of a torus. Note that the irregular sampling makes it impossible to reconstruct a surface with the topology of a sphere, because the method uses a global scale parameter.

Point Attributes

Many point sets contain point attributes such as normals or colors. Because the scale-space surface reconstruction constructs the connectivity between the original points, any point attributes that make sense for a surface can be propagated to the reconstruction. Figure 62.11 shows the surface reconstructed for a set of points with colors.

trophime.png
trophime_detail.png

Figure 62.11 Left: A reconstructed surface with vertex coloring using the colors of the point set. Right: detail of the surface shown as wire-frame. Images or image-based modeling: MAP (Modèles et simulations pour l'Architecture et le Patrimoine) - UMR CNRS/MCC 3495.

Performance

We have evaluated the scale-space surface reconstruction method on several data sets of various sizes. Figure 62.13 shows the parameter settings used when processing the scale-space reconstruction method.

Data set Neighbors Samples Iterations
Mushroom (s) 10 200 2
Elephant (s) 15 200 1
Knot 1 (s) 8 100 8
Retinal (s) 30 200 4
Femur 25 300 6
Kitten (s) 10 100 4
Knot 2 (s) 8 300 8
Bull (s) 100 600 1
Fandisk (s) 6 300 8
Lion-head (s) 12 400 1
Turbine (s) 60 400 1
Graill 1 10 200 4
Graill 11 10 200 4
Man (s) 100 600 1
Face 10 200 4
Carc 20 LR 10 200 4
Carc 21 LR 10 200 4
Bowl 1 10 200 4
Bowl 5 10 200 4
Fount 1 10 200 4
Fount 11 10 200 4
Trophime 10 200 4
Carc 20 10 200 4
Carc 21 10 200 4

Figure 62.12 Settings used when running the scale-space method on several data sets. Data sets indicated by "(s)" are synthetic. Neighbors indicates the number of neighbors the estimated neighborhood ball should contain on average, Samples indicates the number of sample points used to estimate the neighborhood radius, Iterations indicates the number of scale increase iterations.

Figure 62.13 shows the number of triangles of the surfaces constructed and the running time for doing so. Half of the data sets are synthetic, the others are constructed by image-based dense reconstruction. The data sets are ordered by size, the synthetic sets mainly making up the smaller ones.

The times were measured on a 64-bit Windows 7 system with 8 2.5GHz Intel processors and 16GB RAM. The scale-space constructions employ parallel computing using Intel TBB. The smaller data sets use roughly 0.1GB of memory, but the largest data sets have a peak usage of 10.5GB during the meshing step. For this benchmark we used the Exact_predicates_inexact_constructions_kernel kernel, FS enabled, and Sh disabled.

Data set Number of points Surface triangles Neighborhood (ms) Scale-space (ms) Meshing (ms) Total (ms)
Mushroom (s) 2,337 7,966 7 64 88 159
Elephant (s) 2,775 6,746 4 33 91 128
Knot 1 (s) 3,200 12,800 7 148 134 289
Retinal (s) 3,643 14,154 8 224 122 354
Femur 3,897 12,612 6 592 120 718
Kitten (s) 5,210 20,428 7 198 172 377
Knot 2 (s) 5,760 22,904 5 278 206 489
Bull (s) 6,200 11,892 20 238 149 407
Fandisk (s) 6,475 25,282 5 330 212 547
Lion-head (s) 8,356 27,128 7 125 249 382
Turbine (s) 9,210 14,084 12 293 201 507
Graill 1 12,681 45,326 14 1,166 381 1,561
Graill 11 14,617 54,550 18 1,072 454 1,544
Man (s) 17,495 38,496 25 845 421 1,291
Face 296,275 906,648 128 17,983 7,843 25,957
Carc 20 LR 1,311,850 1,311,850 154 20,652 10,822 31,632
Carc 21 LR 1,443,510 1,443,510 157 19,778 12,290 32,228
Bowl 1 2,741,328 2,841,328 351 34,688 22,559 57,608
Bowl 5 3,682,188 3,682,188 471 50,762 29,899 81,144
Fount 1 1,504,761 5,946,502 764 112,028 51,451 164,262
Fount 11 2,589,083 10,203,108 1,350 216,623 85,231 303,248
Trophime 4,572,657 16,982,194 2,413 630,502 134,427 767,389
Carc 20 5,370,681 21,190,868 2,531 304,015 179,223 485,813
Carc 21 5,800,135 23,042,576 2,916 365,731 212,840 581,535

Figure 62.13 Number of triangles of the surfaces constructed and running times in milliseconds for the scale-space method. Data sets indicated by "(s)" are synthetic.

Examples

This first example shows how to read a point set from a file, reconstruct the surface, and for each shell of the reconstruction, access the triangle vertex indices.


File Scale_space_reconstruction_3/scale_space.cpp

#include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
#include <CGAL/Scale_space_surface_reconstruction_3.h>
#include <CGAL/IO/read_points.h>
#include <CGAL/Timer.h>
#include <fstream>
#include <iostream>
typedef CGAL::Exact_predicates_inexact_constructions_kernel Kernel;
typedef Kernel::Point_3 Point;
typedef CGAL::Scale_space_surface_reconstruction_3<Kernel> Reconstruction;
typedef Reconstruction::Facet_const_iterator Facet_iterator;
int main(int argc, char** argv)
{
if (argc != 2)
{
std::cerr << "Error, no input file provided\n";
return 1;
}
std::cerr << "Reading " << std::flush;
std::vector<Point> points;
if(!CGAL::IO::read_points(argv[1], std::back_inserter(points)))
{
std::cerr << "Error: cannot read file" << std::endl;
return EXIT_FAILURE;
}
std::cerr << "done: " << points.size() << " points." << std::endl;
std::cerr << "Reconstruction ";
CGAL::Timer t;
t.start();
// Construct the mesh in a scale space.
Reconstruction reconstruct (points.begin(), points.end());
reconstruct.increase_scale(4);
reconstruct.reconstruct_surface();
std::cerr << "done in " << t.time() << " sec." << std::endl;
t.reset();
std::ofstream out ("out.off");
out << reconstruct;
std::cerr << "Writing result in " << t.time() << " sec." << std::endl;
std::cerr << "Done." << std::endl;
return EXIT_SUCCESS;
}

The different phases of the reconstruction can also be performed individually. This may be helpful in cases where the correct parameter settings are not known in advance. In Figure 62.4, the right figure shows the surface after incrementing the scale. After examining the reconstructed surface, we may decide that performing more iterations may better model the top-left part.

The next example performs two reconstructions with different smoothing of the point set by further advancing the scale space. It illustrates how the code can be used in an interactive way.


File Scale_space_reconstruction_3/scale_space_incremental.cpp

#include <CGAL/Scale_space_surface_reconstruction_3.h>
#include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
#include <CGAL/IO/read_points.h>
#include <algorithm>
#include <fstream>
#include <iostream>
typedef CGAL::Exact_predicates_inexact_constructions_kernel Kernel;
typedef CGAL::Scale_space_surface_reconstruction_3< Kernel > Reconstruction;
typedef CGAL::Scale_space_reconstruction_3::Weighted_PCA_smoother< Kernel > Smoother;
typedef CGAL::Scale_space_reconstruction_3::Alpha_shape_mesher< Kernel > Mesher;
typedef Reconstruction::Point Point;
typedef Reconstruction::Facet_const_iterator Facet_iterator;
// function for writing the reconstruction output in the off format
void dump_reconstruction(const Reconstruction& reconstruct, std::string name)
{
std::ofstream output(name.c_str());
output << "OFF " << reconstruct.number_of_points() << " "
<< reconstruct.number_of_facets() << " 0\n";
std::copy(reconstruct.points_begin(),
reconstruct.points_end(),
std::ostream_iterator<Point>(output,"\n"));
for( Facet_iterator it = reconstruct.facets_begin(); it != reconstruct.facets_end(); ++it )
output << "3 " << *it << std::endl;
}
int main(int argc, char* argv[])
{
if (argc != 2)
{
std::cerr << "Error, no input file provided\n";
return 1;
}
std::cout << "Reading " << std::flush;
std::vector<Point> points;
if(!CGAL::IO::read_points(argv[1], std::back_inserter(points)))
{
std::cerr << "Error: cannot read file" << std::endl;
return EXIT_FAILURE;
}
std::cout << "done: " << points.size() << " points." << std::endl;
// Construct the reconstruction
Reconstruction reconstruct;
// Add the points.
reconstruct.insert( points.begin(), points.end() );
// Two passes
for (std::size_t i = 0; i < 2; ++ i)
{
// Construct the smoother with parameters for
// the neighborhood squared radius estimation.
Smoother smoother( 10, 100 );
// Advance the scale-space several steps.
// This automatically estimates the scale-space.
reconstruct.increase_scale( 2, smoother );
// Reconstruct the surface from the current scale-space.
std::cout << "Neighborhood squared radius is "
<< smoother.squared_radius() << std::endl;
Mesher mesher (smoother.squared_radius());
reconstruct.reconstruct_surface(mesher);
if (i == 0)
{
std::cout << "First reconstruction done." << std::endl;
// Write the reconstruction.
dump_reconstruction(reconstruct, "reconstruction1.off");
}
else
{
std::cout << "Second reconstruction done." << std::endl;
// Write the reconstruction.
dump_reconstruction(reconstruct, "reconstruction2.off");
}
}
std::cout << "Reconstructions are ready to be examinated in your favorite viewer" << std::endl;
return EXIT_SUCCESS;
}

This example shows how to force the output surface to be 2-manifold. In this case, some facets may be discarded in order to remove non-manifold simplices. They are stored in a garbage container in an unordered fashion. They are accessed by iterating on this container.


File Scale_space_reconstruction_3/scale_space_manifold.cpp

#include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
#include <CGAL/Scale_space_surface_reconstruction_3.h>
#include <CGAL/IO/read_points.h>
#include <CGAL/Timer.h>
#include <fstream>
#include <iostream>
typedef CGAL::Exact_predicates_inexact_constructions_kernel Kernel;
typedef CGAL::Scale_space_surface_reconstruction_3< Kernel > Reconstruction;
typedef CGAL::Scale_space_reconstruction_3::Weighted_PCA_smoother< Kernel > Smoother;
typedef CGAL::Scale_space_reconstruction_3::Alpha_shape_mesher< Kernel > Mesher;
typedef Reconstruction::Point Point;
typedef Reconstruction::Facet_const_iterator Facet_iterator;
typedef Mesher::Facet_const_iterator Mesher_iterator;
typedef CGAL::Timer Timer;
int main(int argc, char* argv[])
{
if(argc != 2)
{
std::cerr << "Error, no input file provided\n";
return 1;
}
std::cerr << "Reading " << std::flush;
std::vector<Point> points;
if(!CGAL::IO::read_points(argv[1], std::back_inserter(points)))
{
std::cerr << "Error: cannot read file" << std::endl;
return EXIT_FAILURE;
}
std::cerr << "done: " << points.size() << " points." << std::endl;
Timer t;
t.start();
// Construct the mesh in a scale space.
Reconstruction reconstruct(points.begin(), points.end() );
Smoother smoother(10, 200 );
reconstruct.increase_scale(4, smoother);
Mesher mesher(smoother.squared_radius(),
false, // Do not separate shells
true // Force manifold output
);
reconstruct.reconstruct_surface(mesher);
std::cerr << "Reconstruction done in " << t.time() << " sec." << std::endl;
t.reset();
std::ofstream out("out.off");
// Write the reconstruction.
for(Facet_iterator it = reconstruct.facets_begin(); it != reconstruct.facets_end(); ++it )
out << "3 "<< *it << '\n'; // We write a '3' in front so that it can be assembled into an OFF file
std::cerr << "Writing result in " << t.time() << " sec." << std::endl;
out.close();
t.reset();
std::ofstream garbage("garbage.off");
// Write facets that were removed to force manifold output
for(Mesher_iterator it = mesher.garbage_begin(); it != mesher.garbage_end(); ++it )
garbage << "3 "<< *it << '\n'; // We write a '3' in front so that it can be assembled into an OFF file
std::cerr << "Writing garbage facets in " << t.time() << " sec." << std::endl;
std::cerr << "Done." << std::endl;
return EXIT_SUCCESS;
}

This last example shows how to use the alternative operators Jet_smoother and Advancing_front_mesher.


File Scale_space_reconstruction_3/scale_space_advancing_front.cpp

#include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
#include <CGAL/Point_set_3.h>
#include <CGAL/Scale_space_surface_reconstruction_3.h>
#include <CGAL/Scale_space_reconstruction_3/Advancing_front_mesher.h>
#include <CGAL/Scale_space_reconstruction_3/Jet_smoother.h>
#include <CGAL/Timer.h>
#include <fstream>
#include <iostream>
typedef CGAL::Exact_predicates_inexact_constructions_kernel Kernel;
typedef CGAL::Scale_space_surface_reconstruction_3<Kernel> Reconstruction;
typedef CGAL::Scale_space_reconstruction_3::Advancing_front_mesher<Kernel> Mesher;
typedef CGAL::Scale_space_reconstruction_3::Jet_smoother<Kernel> Smoother;
typedef Kernel::Point_3 Point;
typedef CGAL::Point_set_3<Point> Point_set;
typedef Reconstruction::Facet_const_iterator Facet_iterator;
int main(int argc, char** argv)
{
if (argc != 2)
{
std::cerr << "Error, no input file provided\n";
return 1;
}
// Read the data.
std::cerr << "Reading " << std::flush;
Point_set points;
if(!CGAL::IO::read_point_set(argv[1], points))
{
std::cerr << "Error: cannot read file" << std::endl;
return EXIT_FAILURE;
}
std::cerr << "done: " << points.size() << " points." << std::endl;
std::cerr << "Reconstruction ";
CGAL::Timer t;
t.start();
// Construct the mesh in a scale space.
Reconstruction reconstruct (points.points().begin(), points.points().end());
reconstruct.increase_scale<Smoother> (4);
reconstruct.reconstruct_surface (Mesher (0.5));
std::cerr << "done in " << t.time() << " sec." << std::endl;
t.reset();
std::ofstream out ("out.off");
out << "OFF" << std::endl << points.size() << " " << reconstruct.number_of_facets() << " 0" << std::endl;
for (Point_set::iterator it = points.begin(); it != points.end(); ++ it)
out << points.point(*it) << std::endl;
for (Reconstruction::Facet_iterator it = reconstruct.facets_begin();
it != reconstruct.facets_end(); ++ it)
out << "3 " << (*it)[0] << " " << (*it)[1] << " " << (*it)[2] << std::endl;
std::cerr << "Writing result in " << t.time() << " sec." << std::endl;
std::cerr << "Done." << std::endl;
return EXIT_SUCCESS;
}

Design and Implementation History

This method was developed by Julie Digne et al. in 2011 [1] and implemented by Thijs van Lankveld at Inria - Sophia Antipolis in 2014.