Update tutorial + example

Classification/examples/Classification/gis_tutorial_example.cpp

@@ -18,6 +18,8 @@
 #include <boost/graph/adjacency_list.hpp>
 #include <CGAL/boost/graph/split_graph_into_polylines.h>
 
+#include <CGAL/IO/WKT.h>
+
 #include <CGAL/Constrained_Delaunay_triangulation_2.h>
 #include <CGAL/Constrained_triangulation_plus_2.h>
 
@@ -69,7 +71,7 @@ using Concurrency_tag = CGAL::Parallel_tag;
 using Concurrency_tag = CGAL::Sequential_tag;
 #endif
 
-///////////////////////////////////////////////////////////////////
+///////////////////////////////////////////////////////////////////
 //! [Contouring functions]
 
 bool face_has_isovalue (TIN_with_info::Face_handle fh, double isovalue)
@@ -77,6 +79,8 @@ bool face_has_isovalue (TIN_with_info::Face_handle fh, double isovalue)
   bool above = false, below = false;
   for (int i = 0; i < 3; ++ i)
   {
+    // Face has isovalue if one of its vertices is above and another
+    // one below
     if (fh->vertex(i)->point().z() > isovalue)
       above = true;
     if (fh->vertex(i)->point().z() < isovalue)
@@ -97,6 +101,7 @@ Segment_3 isocontour_in_face (TIN_with_info::Face_handle fh, double isovalue)
     Point_3 p0 = fh->vertex((i+1) % 3)->point();
     Point_3 p1 = fh->vertex((i+2) % 3)->point();
 
+    // Check if the isovalue crosses segment (p0,p1)
     if ((p0.z() - isovalue) * (p1.z() - isovalue) > 0)
       continue;
 
@@ -107,9 +112,9 @@ Segment_3 isocontour_in_face (TIN_with_info::Face_handle fh, double isovalue)
       std::swap (zbottom, ztop);
       std::swap (p0, p1);
     }
-    
+
+    // Compute position of segment vertex
     double ratio = (isovalue - zbottom) / (ztop - zbottom);
-    
     Point_3 p = CGAL::barycenter (p0, (1 - ratio), p1,ratio);
     
     if (source_found)
@@ -298,7 +303,7 @@ int main (int argc, char** argv)
                          (boost::make_function_output_iterator
                           ([&](const std::pair<TIN_with_info::Face_handle, Mesh::Face_index>& ff)
                            {
-                             // Color unassigned faces grey
+                             // Color unassigned faces gray
                              if (ff.first->info() < 0)
                                color_map[ff.second] = CGAL::Color(128, 128, 128);
                              else
@@ -324,9 +329,6 @@ int main (int argc, char** argv)
   int min_size = 100000;
 
   std::vector<TIN_with_info::Vertex_handle> to_remove;
-  std::ofstream dbg ("dbg.xyz");
-  dbg.precision(18);
-
   for (TIN_with_info::Vertex_handle vh : tin_with_info.finite_vertex_handles())
   {
     TIN_with_info::Face_circulator circ = tin_with_info.incident_faces (vh),
@@ -345,10 +347,7 @@ int main (int argc, char** argv)
     while (++ circ != start);
 
     if (!keep)
-    {
       to_remove.push_back (vh);
-      dbg << vh->point() << std::endl;
-    }
   }
 
   std::cerr << to_remove.size() << " vertices(s) will be removed after filtering" << std::endl;
@@ -359,11 +358,11 @@ int main (int argc, char** argv)
   ///////////////////////////////////////////////////////////////////
 
   // Save as Mesh
-  CGAL::Surface_mesh<Point_3> dsm_mesh;
+  CGAL::Surface_mesh<Point_3> dem_mesh;
-  CGAL::copy_face_graph (tin_with_info, dsm_mesh);
+  CGAL::copy_face_graph (tin_with_info, dem_mesh);
-  std::ofstream dsm_ofile ("dsm.ply", std::ios_base::binary);
+  std::ofstream dem_ofile ("dem.ply", std::ios_base::binary);
-  CGAL::set_binary_mode (dsm_ofile);
+  CGAL::set_binary_mode (dem_ofile);
-  CGAL::write_ply (dsm_ofile, dsm_mesh);
+  CGAL::write_ply (dem_ofile, dem_mesh);
 
   ///////////////////////////////////////////////////////////////////
   //! [Rastering]
@@ -378,13 +377,20 @@ int main (int argc, char** argv)
 
   // Use PPM format (Portable PixMap) for simplicity
   std::ofstream raster_ofile ("raster.ppm", std::ios_base::binary);
-  raster_ofile << "P6" << std::endl << width << " " << height << std::endl << 255 << std::endl;
 
-  TIN_with_info::Face_handle location;
+  // PPM header
+  raster_ofile << "P6" << std::endl // magic number
+               << width << " " << height << std::endl // dimensions of the image
+               << 255 << std::endl; // maximum color value
 
   // Use rainbow color ramp output
   Color_ramp color_ramp;
-  
+
+  // Keeping track of location from one point to its neighbor allows
+  // for fast locate in DT
+  TIN_with_info::Face_handle location;
+
+  // Query each pixel of the image
   for (std::size_t y = 0; y < height; ++ y)
     for (std::size_t x = 0; x < width; ++ x)
     {
@@ -394,8 +400,8 @@ int main (int argc, char** argv)
       
       location = tin_with_info.locate (query, location);
 
+      // Points outside the convex hull will be colored black
       std::array<unsigned char, 3> colors { 0, 0, 0 };
-      
       if (!tin_with_info.is_infinite(location))
       {
         std::array<double, 3> barycentric_coordinates
@@ -411,11 +417,11 @@ int main (int argc, char** argv)
              + barycentric_coordinates[1] * location->vertex(1)->point().z()
              + barycentric_coordinates[2] * location->vertex(2)->point().z());
 
+        // Color ramp generates a color depending on a value from 0 to 1
         double height_ratio = (height_at_query - bbox.zmin()) / (bbox.zmax() - bbox.zmin());
         colors = color_ramp.get(height_ratio);
       }
       raster_ofile.write ((char*)(&colors), 3);
-
     }
   
   //! [Rastering]
@@ -470,6 +476,7 @@ int main (int argc, char** argv)
         Segment_3 segment = isocontour_in_face (vh, iv);
         for (const Point_3& p : { segment.source(), segment.target() })
         {
+          // Only insert end points of segments once to get a well connected graph
           Map_p2v::iterator iter;
           bool inserted;
           std::tie (iter, inserted) = map_p2v.insert (std::make_pair (p, Segment_graph::vertex_descriptor()));
@@ -496,16 +503,10 @@ int main (int argc, char** argv)
   std::cerr << polylines.size() << " polylines computed, with "
             << map_p2v.size() << " vertices in total" << std::endl;
 
-  // Output to polyline file
+  // Output to WKT file
-  std::ofstream contour_ofile ("contour.polylines.txt");
+  std::ofstream contour_ofile ("contour.wkt");
   contour_ofile.precision(18);
-  for (const std::vector<Point_3>& poly : polylines)
+  CGAL::write_multi_linestring_WKT (contour_ofile, polylines);
-  {
-    contour_ofile << poly.size();
-    for (const Point_3& p : poly)
-      contour_ofile << " " << p;
-    contour_ofile << std::endl;
-  }
 
   //! [Contouring split]
   ///////////////////////////////////////////////////////////////////
@@ -521,22 +522,28 @@ int main (int argc, char** argv)
   // Simplification algorithm with limit on distance
   PS::simplify (ctp, PS::Squared_distance_cost(), PS::Stop_above_cost_threshold (16 * spacing * spacing));
 
-  // Output to file
+  polylines.clear();
-  std::ofstream simplified_ofile ("simplified.polylines.txt");
-  simplified_ofile.precision(18);
-  std::size_t nb_vertices = 0;
   for (CTP::Constraint_id cid : ctp.constraints())
   {
-    simplified_ofile << ctp.vertices_in_constraint (cid).size();
+    polylines.push_back (std::vector<Point_3>());
-    nb_vertices += ctp.vertices_in_constraint (cid).size();
+    polylines.back().reserve (ctp.vertices_in_constraint (cid).size());
     for (CTP::Vertex_handle vh : ctp.vertices_in_constraint(cid))
-      simplified_ofile << " " << vh->point();
+      polylines.back().push_back (vh->point());
-    simplified_ofile << std::endl;
   }
 
-  std::cerr << nb_vertices << " vertices remaining after simplification ("
+  std::size_t nb_vertices
+    = std::accumulate (polylines.begin(), polylines.end(), 0,
+                       [](std::size_t size, const std::vector<Point_3>& poly) -> std::size_t
+                       { return size + poly.size(); });
+
+  std::cerr << nb_vertices
+            << " vertices remaining after simplification ("
             << 100. * (nb_vertices / double(map_p2v.size())) << "%)" << std::endl;
   
+  // Output to WKT file
+  std::ofstream simplified_ofile ("simplified.wkt");
+  simplified_ofile.precision(18);
+  CGAL::write_multi_linestring_WKT (simplified_ofile, polylines);
 
   //! [Contouring simplify]
   ///////////////////////////////////////////////////////////////////

Documentation/doc/Documentation/Tutorials/Tutorial_GIS.txt

@@ -13,18 +13,19 @@ namespace CGAL {
 Many sensors used in GIS applications (for example LIDAR), generate
 dense point clouds. Such applications often take advantage of more
 advanced data structure: for example, Triangulated Irregular Networks
-(TIN) that can be the basis for Digital Surface Modeling (DSM). Point
+(TIN) that can be the basis for Digital Surface Modeling (DSM) and for
-clouds can also be enriched by classification information that
+the generation of Digital Elevation Modeling (DEM). Point clouds can
-segments the points into a set of semantic labels.
+also be enriched by classification information that segments the
+points into a set of semantic labels.
 
 \section TutorialGIS_TIN TIN (Triangulated Irregular Network)
 
 \cgal provides several triangulation data structures and algorithms. A
-TIN can be generated by combining the 2D Delaunay triangulation with
+TIN can be generated by combining the 2D Delaunay %Triangulation with
 projection traits: the triangulation structure is computed using the
 2D positions of the points along a selected plane (usually, the
-XY-plane), but the 3D positions of the points are used for
+XY-plane), while the 3D positions of the points are kept for
-visualisation and measurements.
+visualization and measurements.
 
 A TIN data structure can thus simply be defined the following way:
 
@@ -36,15 +37,16 @@ operator. Generating the TIN is then straightforward:
 
 \snippet Classification/gis_tutorial_example.cpp Init TIN
 
-Because the Delaunay Triangulation of \cgal is a model of `FaceGraph`,
+Because the Delaunay %Triangulation of \cgal is a model of `FaceGraph`,
 the generated TIN can easily be converted into a mesh structure such
 as `CGAL::Surface_mesh` and be saved into whatever formats are
 supported by such structure:
 
 \snippet Classification/gis_tutorial_example.cpp Save TIN
 
-An example of a TIN computed this way on the San Fransisco data set is
+An example of a TIN computed this way on the San Fransisco data set
-given on Figure \cgalFigureRef{TutorialGISFigTIN}.
+(see \ref TutorialGIS_Reference) is given on Figure
+\cgalFigureRef{TutorialGISFigTIN}.
 
 \cgalFigureBegin{TutorialGISFigTIN, tin.jpg}
 Input point cloud and output TIN.
@@ -58,12 +60,12 @@ standard procedure is to compute a Digital Elevation Modeling, that is
 to say a representation of the ground as another TIN after filtering
 non-ground points.
 
-We propose, as an example, a simple DEM estimation that consists in
+We propose, as an example, a simple DEM estimation decomposed in the
-the following steps:
+following steps:
 
 1. Thresholding the height of the facets to remove brutal changes of
 elevation
-2. Clustering the other facets in connected components
+2. Clustering the other facets into connected components
 3. Filtering all components smaller than a user-defined threshold
 
 This algorithm relies on 2 parameters: a height threshold that
@@ -74,11 +76,11 @@ projection.
 \subsection TutorialGIS_DEM_info TIN with info
 
 Because it takes advantage of the flexible \cgal Delaunay
-Triangulation API, our TIN can be enriched with information on
+%Triangulation API, our TIN can be enriched with information on
-vertices and/or on faces. In our case, each vertices keeps track of
+vertices and/or on faces. In our case, each vertex keeps track of the
-the index of the corresponding point in the input point cloud (which
+index of the corresponding point in the input point cloud (which will
-will allow to filter ground points afterwards), and each faces is
+allow to filter ground points afterwards), and each face is given the
-given the index of its connected component.
+index of its connected component.
 
 \snippet Classification/gis_tutorial_example.cpp TIN_with_info
 
@@ -86,7 +88,7 @@ given the index of its connected component.
 
 Connected components are identified through a flooding algorithm: from
 a seed face, all incident faces are inserted in the current connected
-component unless their heights exceeds the user-defined threhold.
+component unless their heights exceed the user-defined threshold.
 
 \snippet Classification/gis_tutorial_example.cpp Components
 
@@ -99,8 +101,8 @@ An example of a TIN colored by connected components is
 given on Figure \cgalFigureRef{TutorialGISFigComponents}.
 
 \cgalFigureBegin{TutorialGISFigComponents, components.jpg}
-TIN colored by connected components. Faces above height treshold are
+TIN colored by connected components. Faces above height threshold are
-not assigned to any component and are displayed in grey.
+not assigned to any component and are displayed in gray.
 \cgalFigureEnd
 
 \subsection TutorialGIS_DEM_filtering Filtering
@@ -152,20 +154,21 @@ cross a face, and to extract it then:
 \snippet Classification/gis_tutorial_example.cpp Contouring functions
 
 From these functions, we can create a graph of segments to process
-later into a set of polylines. To do so, we use the `adjacency_list`
+later into a set of polylines. To do so, we use the
-structure from boost and keep track of a mapping from the positions of
+`boost::adjacency_list` structure from boost and keep track of a
-the end points to the vertices of the graph.
+mapping from the positions of the end points to the vertices of the
+graph.
 
 The following code computes 50 isovalues evenly distributed between
 the minimum and maximum heights of the point cloud and creates a graph
-containing all the segment sections of all isolevels:
+containing all isolevels:
 
 \snippet Classification/gis_tutorial_example.cpp Contouring extraction
 
 \subsection TutorialGIS_Contour_Splitting Splitting Into Polylines
 
 Once the graph is created, splitting it into polylines is easily
-performed using the function `CGAL::split_graph_into_polylines()`:
+performed using the function `split_graph_into_polylines()`:
 
 \snippet Classification/gis_tutorial_example.cpp Contouring split
 
@@ -188,12 +191,12 @@ Because the output is quite noisy, users may want to simplify the
 polylines. \cgal provides a polyline simplification algorithm that
 guarantees that two polylines won't intersect after
 simplification. This algorithm takes advantage of the Constrained
-Delaunay Triangulation 2 Plus, which embeds polylines as a set of
+Delaunay %Triangulation 2 Plus, which embeds polylines as a set of
 constraints:
 
 \snippet Classification/gis_tutorial_example.cpp CDT
 
-The following code simplies the polyline set based on the squared
+The following code simplifies the polyline set based on the squared
 distance to the original polylines, stopping when no more vertex can
 be removed without going farther than 4 times the average spacing.
 
@@ -209,14 +212,14 @@ vertices.
 
 \section TutorialGIS_Classify Classifying
 
-\cgal provides a Classification package which can be used to segment a
+\cgal provides a %Classification package which can be used to segment a
 point cloud into a user-defined label set. The state-of-the-art
-classifier currently available in \cgal is the Random Forest from
+classifier currently available in \cgal is the %Random Forest from
 ETHZ. As it is a supervised classifier, a training set is needed.
 
 The following snippet shows how to use some manually selected training
-set to train a Random Forest classifier and compute a classification
+set to train a %Random Forest classifier and compute a classification
-regularize by a Graph Cut algorithm:
+regularized by a Graph Cut algorithm:
 
 \snippet Classification/gis_tutorial_example.cpp Classification
 
@@ -225,10 +228,11 @@ Figure \cgalFigureRef{TutorialGISFigClassif}.
 
 \cgalFigureBegin{TutorialGISFigClassif, classif_tuto.jpg}
 Top: a slice of the point cloud classified by hand. Bottom: a
-classification regularized by Graph Cut after training on 3 slices.
+classification regularized by Graph Cut after training on 3 manually
+classified slices.
 \cgalFigureEnd
 
-\section TutorialsReconstruction_recap Full Code Example
+\section TutorialGIS_Code Full Code Example
 
 All the code snippets used in this tutorial can be assembled to create
 a full GIS pipeline (provided the correct _includes_ are
@@ -237,7 +241,22 @@ described in this tutorial.
 
 \include Classification/gis_tutorial_example.cpp
 
+\section TutorialGIS_Reference References
 
+This tutorial was based on the following \cgal packages:
+
+- \ref PkgTriangulation2Ref
+- \ref PkgPointSet3Ref
+- \ref PkgPointSetProcessing3Ref
+- \ref PkgSurface_mesh
+- \ref PkgBGLRef
+- \ref PkgPolygonMeshProcessingRef
+- \ref PkgPolylineSimplification2Ref
+- \ref PkgClassificationRef
+
+The data set used throughout this tutorial comes from the
+https://www.usgs.gov/ database, licensed under the _US Government
+Public Domain_.
 
 */
 }