Merge pull request #442 from sgiraudot/CGAL_Tutorials-Surface_reconstruction-GF

CGAL Tutorial surface reconstruction

Advancing_front_surface_reconstruction/doc/Advancing_front_surface_reconstruction/Advancing_front_surface_reconstruction.txt

@@ -47,6 +47,9 @@ not provide any guarantee on the topology of the surface.
 
 We describe next the algorithm and provide examples.
 
+\note A \ref tuto_reconstruction "detailed tutorial on surface reconstruction"
+is provided with a guide to choose the most appropriate method along
+with pre- and post-processing.
 
 \section AFSR_Definitions Definitions and the Algorithm
 

Documentation/doc/Documentation/Doxyfile.in

@@ -11,7 +11,8 @@ HTML_HEADER         = ${CGAL_DOC_HEADER}
 LAYOUT_FILE         = ${CGAL_DOC_RESOURCE_DIR}/DoxygenLayout.xml
 GENERATE_TAGFILE    = ${CGAL_DOC_TAG_GEN_DIR}/Manual.tag
 EXAMPLE_PATH        = ${CGAL_Convex_hull_2_EXAMPLE_DIR} \
-                      ${CGAL_Kernel_23_EXAMPLE_DIR}
+                      ${CGAL_Kernel_23_EXAMPLE_DIR} \
+                      ${CGAL_Poisson_surface_reconstruction_3_EXAMPLE_DIR}
 FILTER_PATTERNS     = *.txt=${CMAKE_BINARY_DIR}/pkglist_filter
 
 HTML_EXTRA_FILES += ${CGAL_DOC_RESOURCE_DIR}/hacks.js \

Documentation/doc/Documentation/Tutorials/Tutorial_reconstruction.txt

@@ -0,0 +1,295 @@
+namespace CGAL {
+/*!
+\example Poisson_surface_reconstruction_3/tutorial_example.cpp
+*/
+
+/*!
+
+\page tuto_reconstruction Surface Reconstruction from Point Clouds
+\cgalAutoToc
+
+\author Simon Giraudot
+
+Surface reconstruction from point clouds is a core topic in geometry
+processing \cgalCite{cgal:btsag-asosr-16}. It is an ill-posed problem:
+there is an infinite number of surfaces that approximate a single
+point cloud and a point cloud does not define a surface in
+itself. Thus additional assumptions and constraints must be defined by
+the user and reconstruction can be achieved in many different
+ways. This tutorial provides guidance on how to use the different
+algorithms of \cgal to effectively perform surface reconstruction.
+
+\section TutorialsReconstruction_algorithms Which algorithm should I use?
+
+\cgal offers three different algorithms for surface reconstruction:
+
+- \ref Chapter_Poisson_Surface_Reconstruction "Poisson Surface Reconstruction"
+- \ref Chapter_Advancing_Front_Surface_Reconstruction "Advancing Front Surface Reconstruction"
+- \ref Chapter_Scale_space_reconstruction "Scale Space Surface Reconstruction"
+
+
+Because reconstruction is an ill-posed problem, it must be regularized
+via prior knowledge. Differences in prior lead to different
+algorithms, and choosing one or the other of these methods is
+dependent on these priors. For example, Poisson always generates
+closed shapes (bounding a volume) and requires normals but does not
+interpolate input points (the output surface does not pass exactly
+through the input points). The following table lists different
+properties of the input and output to help the user choose the method
+best suited to each problem:
+
+<center>
+|                                          | Poisson | Advancing front  | Scale space      |
+|------------------------------------------|:-------:|:----------------:|:----------------:|
+| Are normals required?                    | Yes     | No               | No               |
+| Is noise handled?                        | Yes     | By preprocessing | Yes              |
+| Is variable sampling handled?            | Yes     | Yes              | By preprocessing |
+| Are input points exactly on the surface? | No      | Yes              | Yes              |
+| Is the output always closed?             | Yes     | No               | No               |
+| Is the output always smooth?             | Yes     | No               | No               |
+| Is the output always manifold?           | Yes     | Yes              | Optional         |
+| Is the output always orientable?         | Yes     | Yes              | Optional         |
+</center>
+
+\cgalFigureBegin{TutorialsReconstructionFigComparisons, compare_reconstructions.png}
+Comparison of reconstruction methods applied to the same input (full shape and close-up). From left to right: original point cloud; Poisson; advancing front; scale space.
+\cgalFigureEnd
+
+More information on these different methods can be found on their
+respective manual pages and in Section \ref TutorialsReconstruction_reconstruction.
+
+\section TutorialsReconstruction_overview Pipeline Overview
+
+This tutorial aims at providing a more comprehensive view of the
+possibilities offered by \cgal for dealing with point clouds, for
+surface reconstruction purposes. The following diagram shows an
+overview (not exhaustive) of common reconstruction steps using \cgal
+tools.
+
+\cgalFigureBegin{TutorialsReconstructionFigPipeline, reconstruction.svg}
+Pipeline Overview
+\cgalFigureEnd
+
+We now review some of these steps in more detail.
+
+
+\section TutorialsReconstruction_input Reading Input
+
+The reconstruction algorithms in \cgal take a range of iterators on a
+container as input and use property maps to access the points (and the
+normals when they are required). Points are typically stored in plain
+text format (denoted as 'XYZ' format) where each point is separated by
+a newline character and each coordinate separated by a white
+space. Other formats available are 'OFF', 'PLY' and 'LAS'. \cgal
+provides functions to read such formats:
+
+- `read_xyz_points()`
+- `read_off_points()`
+- `read_ply_points()`
+- `read_ply_points_with_properties()` to read additional PLY properties
+- `read_las_points()`
+- `read_las_points_with_properties()` to read additional LAS properties
+
+\cgal also provides a dedicated container `CGAL::Point_set_3` to
+handle point sets with additional properties such as normal
+vectors. In this case, property maps are easily handled as shown in
+the following sections. This structure also handles the stream
+operator to read point sets in any of the formats previously
+described. Using this method yields substantially shorter code, as can
+be seen on the following example:
+
+\snippet Poisson_surface_reconstruction_3/tutorial_example.cpp Reading input
+
+\section TutorialsReconstruction_preprocessing Preprocessing
+
+Because reconstruction algorithms have some specific requirements that
+point clouds do not always meet, some preprocessing might be necessary
+to yield the best results.
+
+Note that this _preprocessing_ step is optional: when the input point
+cloud has no imperfections, reconstruction can be applied to it
+without any preprocessing.
+
+
+\cgalFigureBegin{TutorialsReconstructionFigPreprocessing, reconstruction_preproc.png}
+Comparison of advancing front reconstruction output using different
+preprocessing on the same input. The smooth point cloud was generated
+using _jet smoothing_; the simplified point cloud was generated using
+_grid simplification_.
+\cgalFigureEnd
+
+
+\subsection TutorialsReconstruction_preprocessing_outliers Outlier removal
+
+Some acquisition techniques generate points which are far away from
+the surface. These points, commonly referred to as "outliers", have no
+relevance for reconstruction. Using the \cgal reconstruction
+algorithms on outlier-ridden point clouds produce overly distorted
+output, it is therefore strongly advised to filter these outliers
+_before_ performing reconstruction.
+
+\snippet Poisson_surface_reconstruction_3/tutorial_example.cpp Outlier removal
+
+\subsection TutorialsReconstruction_preprocessing_simplification Simplification
+
+Some laser scanners generate points with widely variable
+sampling. Typically, lines of scan are very densely sampled but the
+gap between two lines of scan is much larger, leading to an overly
+massive point cloud with large variations of sampling density. This
+type of input point cloud might generate imperfect output using
+algorithms which, in general, only handle small variations of sampling
+density.
+
+\cgal provides several simplification algorithms. In addition to
+reducing the size of the input point cloud and therefore decreasing
+computation time, some of them can help making the input more
+uniform. This is the case of the function `grid_simplify_point_set()`
+which defines a grid of a user-specified size and keeps one point per
+occupied cell.
+
+\snippet Poisson_surface_reconstruction_3/tutorial_example.cpp Simplification
+
+
+\subsection TutorialsReconstruction_preprocessing_smoothing Smoothing
+
+Although reconstructions via 'Poisson' or 'Scale space' handle noise
+internally, one may want to get tighter control over the smoothing
+step. For example, a slightly noisy point cloud can benefit from some
+reliable smoothing algorithms and be reconstructed via 'Advancing
+front' which provides relevant properties (oriented mesh with
+boundaries).
+
+Two functions are provided to smooth a noisy point cloud with a good
+approximation (i.e. without degrading curvature, for example):
+
+- `jet_smooth_point_set()`
+- `bilateral_smooth_point_set()`
+
+These functions directly modify the container:
+
+\snippet Poisson_surface_reconstruction_3/tutorial_example.cpp Smoothing
+
+
+\subsection TutorialsReconstruction_preprocessing_normal Normal Estimation and Orientation
+
+\ref Chapter_Poisson_Surface_Reconstruction "Poisson Surface Reconstruction"
+requires points with oriented normal vectors. To apply the algorithm
+to a raw point cloud, normals must be estimated first, for example
+with one of these two functions:
+
+- `pca_estimate_normals()`
+- `jet_estimate_normals()`
+
+PCA is faster but jet is more accurate in the presence of high
+curvatures. These function only estimates the _direction_ of the
+normals, not their orientation (the orientation of the vectors might
+not be locally consistent). To properly orient the normals, the
+function `mst_orient_normals()` can be used. Notice that it can also
+be used directly on input normals if their orientation is not
+consistent.
+
+\snippet Poisson_surface_reconstruction_3/tutorial_example.cpp Normal estimation
+
+\section TutorialsReconstruction_reconstruction Reconstruction
+
+\subsection TutorialsReconstruction_reconstruction_poisson Poisson
+
+Poisson reconstruction consists in computing an implicit function
+whose gradient matches the input normal vector field: this indicator
+function has opposite signs inside and outside of the inferred shape
+(hence the need for closed shapes). This method thus requires normals
+and produces smooth closed surfaces. It is not appropriate if the
+surface is expected to interpolate the input points. On the contrary,
+it performs well if the aim is to approximate a noisy point cloud with
+a smooth surface.
+
+\snippet Poisson_surface_reconstruction_3/tutorial_example.cpp Poisson reconstruction
+
+\subsection TutorialsReconstruction_reconstruction_advancing Advancing Front
+
+Advancing front is a Delaunay-based approach which interpolates a
+subset of the input points. It generates triples of point indices
+which describe the triangular facets of the reconstruction: it uses a
+priority queue to sequentially pick the Delaunay facet the most likely
+to be part of the surface, based on a size criterion (to favor the
+small facets) and an angle criterion (to favor smoothness). Its main
+virtue is to generate oriented manifold surfaces with boundaries:
+contrary to Poisson, it does not require normals and is not bound to
+reconstruct closed shapes. However, it requires preprocessing if the
+point cloud is noisy.
+
+The \ref Chapter_Advancing_Front_Surface_Reconstruction "Advancing Front"
+package provides several ways of constructing the function. Here is
+a simple example:
+
+\snippet Poisson_surface_reconstruction_3/tutorial_example.cpp Advancing front reconstruction
+
+\subsection TutorialsReconstruction_reconstruction_scale_space Scale Space
+
+Scale space reconstruction aims at producing a surface which
+interpolates the input points (interpolant) while offering some
+robustness to noise. More specifically, it first applies several times
+a smoothing filter (such as Jet Smoothing) to the input point set to
+produce a scale space; then, the smoothest scale is meshed (using for
+example the Advancing Front mesher); finally, the resulting
+connectivity between smoothed points is propagated to the original raw
+input point set. This method is the right choice if the input point
+cloud is noisy but the user still wants the surface to pass exactly
+through the points.
+
+\snippet Poisson_surface_reconstruction_3/tutorial_example.cpp Scale space reconstruction
+
+\section TutorialsReconstruction_postprocessing Output and Postprocessing
+
+Each of these methods produce a triangle mesh stored in different
+ways. If this output mesh is hampered by defects such as holes or
+self-intersections, \cgal provide several algorithms to post-process
+it (hole filling, remeshing, etc.) in the package \ref
+Chapter_PolygonMeshProcessing "Polygon Mesh Processing".
+
+We do not discuss these functions here as there are many
+postprocessing possibilities whose relevance strongly depends on the
+user's expectations on the output mesh.
+
+The mesh (postprocessed or not) can easily be saved in the PLY format
+(here, using the binary variant):
+
+\snippet Poisson_surface_reconstruction_3/tutorial_example.cpp Output poisson
+
+A polygon soup can also be saved in the OFF format by iterating on the
+points and faces:
+
+\snippet Poisson_surface_reconstruction_3/tutorial_example.cpp Output scale space
+
+Finally, if the polygon soup can be converted into a polygon mesh, it
+can also be saved directly in the OFF format using the stream
+operator:
+
+\snippet Poisson_surface_reconstruction_3/tutorial_example.cpp Output advancing front
+
+\section TutorialsReconstruction_recap Full Code Example
+
+All the code snippets used in this tutorial can be assembled to create
+a full algorithm pipeline (provided the correct _includes_ are
+used). We give a full code example which achieves all the steps
+described in this tutorial. The reconstruction method can be selected
+by the user at runtime with the second argument.
+
+\include Poisson_surface_reconstruction_3/tutorial_example.cpp
+
+\section TutorialsReconstruction_pipeline Full Pipeline Images
+
+The following figure an example of a full reconstruction pipeline
+applied to a bear statue (courtesy _EPFL Computer Graphics and
+Geometry Laboratory_ \cgalCite{cgal:e-esmr}). Two mesh processing
+algorithms (hole filling and isotropic remeshing) are also applied
+(refer to the chapter \ref Chapter_PolygonMeshProcessing "Polygon Mesh Processing"
+for more information).
+
+\cgalFigureBegin{TutorialsReconstructionFigFull, reconstruction_pipeline.png}
+Full reconstruction pipeline (with close-ups).
+\cgalFigureEnd
+
+
+*/
+}

Documentation/doc/Documentation/Tutorials/Tutorials.txt

@@ -16,6 +16,8 @@ User Manual.
   geometric primitives, the notion of <em>traits classes</em> which
   define what primitives are used by a geometric algorithm, the
   notions of \em concept and \em model.
+- \subpage tuto_reconstruction
+
 
 \section tuto_examples Package Examples
 

Documentation/doc/biblio/cgal_manual.bib

@@ -708,6 +708,12 @@ Teillaud"
   , year = 1999
 }
 
+@Misc{ cgal:e-esmr,
+  title         = {{EPFL} statue model repository},
+  howpublished  = {{EPFL} Computer Graphics and Geometry Laboratory},
+  url           = {http://lgg.epfl.ch/statues_dataset.php}
+}
+
 @inproceedings{ cgal:eddhls-maam-95
   ,author       = {Matthias Eck and Tony DeRose and Tom Duchamp and
                   Hugues Hoppe and Michael Lounsbery and Werner Stuetzle}
@@ -1771,6 +1777,21 @@ ABSTRACT = {We present the first complete, exact and efficient C++ implementatio
   pages=""
 }
 
+@article{cgal:btsag-asosr-16,
+  TITLE = {{A Survey of Surface Reconstruction from Point Clouds}},
+  AUTHOR = {Berger, Matthew and Tagliasacchi, Andrea and Seversky, Lee and Alliez, Pierre and Guennebaud, Gael and Levine, Joshua and Sharf, Andrei and Silva, Claudio},
+  URL = {https://hal.inria.fr/hal-01348404},
+  JOURNAL = {{Computer Graphics Forum}},
+  PUBLISHER = {{Wiley}},
+  PAGES = {27},
+  YEAR = {2016},
+  DOI = {10.1111/cgf.12802},
+  PDF = {https://hal.inria.fr/hal-01348404/file/survey-author.pdf},
+  HAL_ID = {hal-01348404},
+  HAL_VERSION = {v2}
+}
+
+
 @TechReport{ cgal:pabl-cco-07,
 	author 	    = {Poudret, M. and Arnould, A. and Bertrand, Y. and Lienhardt, P.},
 	title 	    = {Cartes Combinatoires Ouvertes.},

Poisson_surface_reconstruction_3/doc/Poisson_surface_reconstruction_3/Poisson_surface_reconstruction_3.txt

@@ -28,6 +28,11 @@ function of the inferred solid (Poisson Surface Reconstruction -
 referred to as Poisson). Poisson is a two steps process: it requires
 solving for the implicit function before function evaluation.
 
+\note A \ref tuto_reconstruction "detailed tutorial on surface reconstruction"
+is provided with a guide to choose the most appropriate method along
+with pre- and post-processing.
+
+
 \section Poisson_surface_reconstruction_3Common Common Reconstruction Pipeline
 
 Surface reconstruction from point sets is often a sequential process

Poisson_surface_reconstruction_3/examples/Poisson_surface_reconstruction_3/CMakeLists.txt

@@ -31,6 +31,8 @@ if ( CGAL_FOUND )
     CGAL_target_use_Eigen(poisson_reconstruction)
     create_single_source_cgal_program( "poisson_reconstruction_function.cpp" )
     CGAL_target_use_Eigen(poisson_reconstruction_function)
+    create_single_source_cgal_program( "tutorial_example.cpp" )
+    CGAL_target_use_Eigen(tutorial_example)
   else()
     message(STATUS "NOTICE: The examples need Eigen 3.1 (or greater) will not be compiled.")
   endif()

Poisson_surface_reconstruction_3/examples/Poisson_surface_reconstruction_3/tutorial_example.cmd

@@ -0,0 +1,3 @@
+data/kitten.xyz
+data/kitten.xyz 1
+data/kitten.xyz 2

Poisson_surface_reconstruction_3/examples/Poisson_surface_reconstruction_3/tutorial_example.cpp

@@ -0,0 +1,223 @@
+#include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
+
+#include <CGAL/Point_set_3.h>
+#include <CGAL/Point_set_3/IO.h>
+
+#include <CGAL/remove_outliers.h>
+#include <CGAL/grid_simplify_point_set.h>
+#include <CGAL/jet_smooth_point_set.h>
+#include <CGAL/jet_estimate_normals.h>
+#include <CGAL/mst_orient_normals.h>
+
+#include <CGAL/poisson_surface_reconstruction.h>
+#include <CGAL/Advancing_front_surface_reconstruction.h>
+#include <CGAL/Scale_space_surface_reconstruction_3.h>
+#include <CGAL/Scale_space_reconstruction_3/Jet_smoother.h>
+#include <CGAL/Scale_space_reconstruction_3/Advancing_front_mesher.h>
+
+#include <CGAL/Surface_mesh.h>
+#include <CGAL/Polygon_mesh_processing/polygon_soup_to_polygon_mesh.h>
+
+#include <cstdlib>
+#include <vector>
+#include <fstream>
+
+// types
+typedef CGAL::Exact_predicates_inexact_constructions_kernel Kernel;
+typedef Kernel::FT FT;
+typedef Kernel::Point_3 Point_3;
+typedef Kernel::Vector_3 Vector_3;
+typedef Kernel::Sphere_3 Sphere_3;
+typedef CGAL::Point_set_3<Point_3, Vector_3> Point_set;
+
+int main(int argc, char*argv[])
+{
+  ///////////////////////////////////////////////////////////////////
+  //! [Reading input]
+
+  Point_set points;
+
+  if (argc < 2)
+  {
+    std::cerr << "Usage: " << argv[0] << " [input.xyz/off/ply/las]" << std::endl;
+    return EXIT_FAILURE;
+  }
+
+  const char* input_file = argv[1];
+  std::ifstream stream (input_file, std::ios_base::binary);
+  if (!stream)
+  {
+    std::cerr << "Error: cannot read file " << input_file << std::endl;
+    return EXIT_FAILURE;
+  }
+
+  stream >> points;
+
+  std::cout << "Read " << points.size () << " point(s)" << std::endl;
+  if (points.empty())
+    return EXIT_FAILURE;
+
+  //! [Reading input]
+  ///////////////////////////////////////////////////////////////////
+
+  ///////////////////////////////////////////////////////////////////
+  //! [Outlier removal]
+
+  CGAL::remove_outliers (points,
+                         24, // Number of neighbors considered for evaluation
+                         points.parameters().threshold_percent (5.0)); // Percentage of points to remove
+
+  std::cout << points.number_of_removed_points()
+            << " point(s) are outliers." << std::endl;
+
+  // Applying point set processing algorithm to a CGAL::Point_set_3
+  // object does not erase the points from memory but place them in
+  // the garbage of the object: memory can be freeed by the user.
+  points.collect_garbage();
+
+  //! [Outlier removal]
+  ///////////////////////////////////////////////////////////////////
+
+  ///////////////////////////////////////////////////////////////////
+  //! [Simplification]
+
+  // Compute average spacing using neighborhood of 6 points
+  double spacing = CGAL::compute_average_spacing<CGAL::Sequential_tag> (points, 6);
+
+  // Simplify using a grid of size 2 * average spacing
+  CGAL::grid_simplify_point_set (points, 2. * spacing);
+
+  std::cout << points.number_of_removed_points()
+            << " point(s) removed after simplification." << std::endl;
+
+  points.collect_garbage();
+
+  //! [Simplification]
+  ///////////////////////////////////////////////////////////////////
+
+  ///////////////////////////////////////////////////////////////////
+  //! [Smoothing]
+
+  CGAL::jet_smooth_point_set<CGAL::Sequential_tag> (points, 24);
+
+  //! [Smoothing]
+  ///////////////////////////////////////////////////////////////////
+
+  unsigned int reconstruction_choice
+    = (argc < 3 ? 0 : atoi(argv[2]));
+
+  if (reconstruction_choice == 0) // Poisson
+  {
+    ///////////////////////////////////////////////////////////////////
+    //! [Normal estimation]
+
+    CGAL::jet_estimate_normals<CGAL::Sequential_tag>
+      (points, 24); // Use 24 neighbors
+
+    // Orientation of normals, returns iterator to first unoriented point
+    typename Point_set::iterator unoriented_points_begin =
+      CGAL::mst_orient_normals(points, 24); // Use 24 neighbors
+
+    points.remove (unoriented_points_begin, points.end());
+
+    //! [Normal estimation]
+    ///////////////////////////////////////////////////////////////////
+
+    ///////////////////////////////////////////////////////////////////
+    //! [Poisson reconstruction]
+
+    CGAL::Surface_mesh<Point_3> output_mesh;
+    CGAL::poisson_surface_reconstruction_delaunay
+      (points.begin(), points.end(),
+       points.point_map(), points.normal_map(),
+       output_mesh, spacing);
+
+    //! [Poisson reconstruction]
+    ///////////////////////////////////////////////////////////////////
+
+    ///////////////////////////////////////////////////////////////////
+    //! [Output poisson]
+
+    std::ofstream f ("out.ply", std::ios_base::binary);
+    CGAL::set_binary_mode (f);
+    CGAL::write_ply (f, output_mesh);
+    f.close ();
+
+    //! [Output poisson]
+    ///////////////////////////////////////////////////////////////////
+  }
+  else if (reconstruction_choice == 1) // Advancing front
+  {
+    ///////////////////////////////////////////////////////////////////
+    //! [Advancing front reconstruction]
+
+    typedef std::array<std::size_t, 3> Facet; // Triple of indices
+
+    std::vector<Facet> facets;
+
+    // The function is called using directly the points raw iterators
+    CGAL::advancing_front_surface_reconstruction(points.points().begin(),
+                                                 points.points().end(),
+                                                 std::back_inserter(facets));
+    std::cout << facets.size ()
+              << " facet(s) generated by reconstruction." << std::endl;
+
+    //! [Advancing front reconstruction]
+    ///////////////////////////////////////////////////////////////////
+
+    ///////////////////////////////////////////////////////////////////
+    //! [Output advancing front]
+
+    // copy points for random access
+    std::vector<Point_3> vertices;
+    vertices.reserve (points.size());
+    std::copy (points.points().begin(), points.points().end(), std::back_inserter (vertices));
+
+    CGAL::Surface_mesh<Point_3> output_mesh;
+    CGAL::Polygon_mesh_processing::polygon_soup_to_polygon_mesh (vertices, facets, output_mesh);
+    std::ofstream f ("out.off");
+    f << output_mesh;
+    f.close ();
+
+    //! [Output advancing front]
+    ///////////////////////////////////////////////////////////////////
+  }
+  else if (reconstruction_choice == 2) // Scale space
+  {
+    ///////////////////////////////////////////////////////////////////
+    //! [Scale space reconstruction]
+
+    CGAL::Scale_space_surface_reconstruction_3<Kernel> reconstruct
+      (points.points().begin(), points.points().end());
+
+    // Smooth using 4 iterations of Jet Smoothing
+    reconstruct.increase_scale (4, CGAL::Scale_space_reconstruction_3::Jet_smoother<Kernel>());
+    // Mesh with the Advancing Front mesher with a maximum facet length of 0.5
+    reconstruct.reconstruct_surface (CGAL::Scale_space_reconstruction_3::Advancing_front_mesher<Kernel>(0.5));
+
+    //! [Scale space reconstruction]
+    ///////////////////////////////////////////////////////////////////
+
+    ///////////////////////////////////////////////////////////////////
+    //! [Output scale space]
+
+    std::ofstream f ("out.off");
+    f << "OFF" << std::endl << points.size () << " "
+      << reconstruct.number_of_facets() << " 0" << std::endl;
+    for (Point_set::Index idx : points)
+      f << points.point (idx) << std::endl;
+    for (const auto& facet : CGAL::make_range (reconstruct.facets_begin(), reconstruct.facets_end()))
+      f << "3 "<< facet << std::endl;
+    f.close ();
+
+    //! [Output scale space]
+    ///////////////////////////////////////////////////////////////////
+  }
+  else // Handle error
+  {
+    std::cerr << "Error: invalid reconstruction id: " << reconstruction_choice << std::endl;
+    return EXIT_FAILURE;
+  }
+
+  return EXIT_SUCCESS;
+}

Scale_space_reconstruction_3/doc/Scale_space_reconstruction_3/Scale_space_reconstruction_3.txt

@@ -17,6 +17,10 @@ Left: 5760 points on a synthetic knot data set. Right: reconstructed surface mes
 
 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 \ref tuto_reconstruction "detailed tutorial on surface reconstruction"
+is provided with a guide to choose the most appropriate method along
+with pre- and post-processing.
+
 
 \section ScaleSpaceReconstruction3secMethod Scale-Space