From a573cb93b05fce77004a95ae499061258c6bcb7f Mon Sep 17 00:00:00 2001 From: Simon Giraudot Date: Wed, 5 Feb 2020 16:05:30 +0100 Subject: [PATCH] Update tutorial + example --- .../Classification/gis_tutorial_example.cpp | 83 +++++++++--------- .../Documentation/Tutorials/Tutorial_GIS.txt | 85 ++++++++++++------- 2 files changed, 97 insertions(+), 71 deletions(-) diff --git a/Classification/examples/Classification/gis_tutorial_example.cpp b/Classification/examples/Classification/gis_tutorial_example.cpp index efbaaac1dfc..cffff288503 100644 --- a/Classification/examples/Classification/gis_tutorial_example.cpp +++ b/Classification/examples/Classification/gis_tutorial_example.cpp @@ -18,6 +18,8 @@ #include #include +#include + #include #include @@ -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& 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 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 dsm_mesh; - CGAL::copy_face_graph (tin_with_info, dsm_mesh); - std::ofstream dsm_ofile ("dsm.ply", std::ios_base::binary); - CGAL::set_binary_mode (dsm_ofile); - CGAL::write_ply (dsm_ofile, dsm_mesh); + CGAL::Surface_mesh dem_mesh; + CGAL::copy_face_graph (tin_with_info, dem_mesh); + std::ofstream dem_ofile ("dem.ply", std::ios_base::binary); + CGAL::set_binary_mode (dem_ofile); + 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 colors { 0, 0, 0 }; - if (!tin_with_info.is_infinite(location)) { std::array 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 - std::ofstream contour_ofile ("contour.polylines.txt"); + // Output to WKT file + std::ofstream contour_ofile ("contour.wkt"); contour_ofile.precision(18); - for (const std::vector& poly : polylines) - { - contour_ofile << poly.size(); - for (const Point_3& p : poly) - contour_ofile << " " << p; - contour_ofile << std::endl; - } + CGAL::write_multi_linestring_WKT (contour_ofile, polylines); //! [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 - std::ofstream simplified_ofile ("simplified.polylines.txt"); - simplified_ofile.precision(18); - std::size_t nb_vertices = 0; + polylines.clear(); for (CTP::Constraint_id cid : ctp.constraints()) { - simplified_ofile << ctp.vertices_in_constraint (cid).size(); - nb_vertices += ctp.vertices_in_constraint (cid).size(); + polylines.push_back (std::vector()); + polylines.back().reserve (ctp.vertices_in_constraint (cid).size()); for (CTP::Vertex_handle vh : ctp.vertices_in_constraint(cid)) - simplified_ofile << " " << vh->point(); - simplified_ofile << std::endl; + polylines.back().push_back (vh->point()); } - 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& 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] /////////////////////////////////////////////////////////////////// diff --git a/Documentation/doc/Documentation/Tutorials/Tutorial_GIS.txt b/Documentation/doc/Documentation/Tutorials/Tutorial_GIS.txt index b06130d2f8b..72b8bd4f225 100644 --- a/Documentation/doc/Documentation/Tutorials/Tutorial_GIS.txt +++ b/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 -clouds can also be enriched by classification information that -segments the points into a set of semantic labels. +(TIN) that can be the basis for Digital Surface Modeling (DSM) and for +the generation of Digital Elevation Modeling (DEM). Point clouds can +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 -visualisation and measurements. +XY-plane), while the 3D positions of the points are kept for +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 -given on Figure \cgalFigureRef{TutorialGISFigTIN}. +An example of a TIN computed this way on the San Fransisco data set +(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 -the following steps: +We propose, as an example, a simple DEM estimation decomposed in the +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 -vertices and/or on faces. In our case, each vertices keeps track of -the index of the corresponding point in the input point cloud (which -will allow to filter ground points afterwards), and each faces is -given the index of its connected component. +%Triangulation API, our TIN can be enriched with information on +vertices and/or on faces. In our case, each vertex keeps track of the +index of the corresponding point in the input point cloud (which will +allow to filter ground points afterwards), and each face is given the +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 -not assigned to any component and are displayed in grey. +TIN colored by connected components. Faces above height threshold are +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` -structure from boost and keep track of a mapping from the positions of -the end points to the vertices of the graph. +later into a set of polylines. To do so, we use the +`boost::adjacency_list` structure from boost and keep track of a +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 -regularize by a Graph Cut algorithm: +set to train a %Random Forest classifier and compute a classification +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_. */ }