mirror of https://github.com/CGAL/cgal
321 lines
12 KiB
Plaintext
321 lines
12 KiB
Plaintext
namespace CGAL {
|
|
/*!
|
|
\example Classification/gis_tutorial_example.cpp
|
|
*/
|
|
|
|
/*!
|
|
|
|
\page tuto_gis GIS (Geographic Information System)
|
|
\cgalAutoToc
|
|
|
|
\author Simon Giraudot
|
|
|
|
Many sensors used in GIS applications (for example LIDAR), generate
|
|
dense point clouds. Such applications often take advantage of more
|
|
advanced data structures: for example, Triangulated Irregular Networks
|
|
(TIN) that can be the base for Digital Evelation Models (DEM) and in
|
|
particular for the generation of Digital Terrain Models (DTM). Point
|
|
clouds can also be enriched by classification information that
|
|
segments the points into ground, vegetation and building points (or
|
|
other user-defined labels).
|
|
|
|
The definitions of some data structures may vary according to
|
|
different sources. We will use the following terms within this
|
|
tutorial:
|
|
|
|
- TIN: _Triangulated Irregular Network_, a 2D triangulation structure
|
|
that connects 3D points based on their projections on the horizontal
|
|
plane.
|
|
|
|
- DSM: _Digital Surface Model_, a model of the whole scanned surface
|
|
including buildings and vegetation. We use a TIN to store the DSM.
|
|
|
|
- DTM: _Digital Terrain Model_, a model of the bare ground surface
|
|
without objects like buildings or vegetation. We both use a TIN and a
|
|
raster to store a DTM.
|
|
|
|
- DEM: _Digital Elevation Model_, a more general term that includes
|
|
both DSM and DTM.
|
|
|
|
This tutorial illustrates the following scenario. From an input point
|
|
cloud, we first compute a DSM stored as a TIN. Then, we filter out
|
|
overly large facets that correspond either to building facades or to
|
|
vegetation noise. Large components corresponding to the ground are
|
|
kept. Holes are filled and the obtained DEM is remeshed. A raster DEM
|
|
and a set of contour polylines are generated from it. Finally,
|
|
supervised 3-label classification is performed to segment vegetation,
|
|
building and group points.
|
|
|
|
\section TutorialGIS_TIN Triangulated Irregular Network (TIN)
|
|
|
|
\cgal provides several triangulation data structures and algorithms. A
|
|
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), 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:
|
|
|
|
\snippet Classification/gis_tutorial_example.cpp TIN DS
|
|
|
|
\section TutorialGIS_DSM Digital Surface Model (DSM)
|
|
|
|
Point clouds in many formats (XYZ, OFF, PLY, LAS) can be easily loaded
|
|
into a `CGAL::Point_set_3` structure, using the stream
|
|
operator. Generating a DSM stored in TIN is then straightforward:
|
|
|
|
\snippet Classification/gis_tutorial_example.cpp Init DSM
|
|
|
|
Because the Delaunay triangulation of \cgal is a model of
|
|
`FaceGraph`, it is straightforward to convert the generated TIN 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 DSM
|
|
|
|
An example of a DSM computed this way on the San Fransisco data set
|
|
(see \ref TutorialGIS_Reference) is given in
|
|
\cgalFigureRef{TutorialGISFigTIN}.
|
|
|
|
\cgalFigureBegin{TutorialGISFigTIN, tin.jpg}
|
|
Input point cloud and output DSM.
|
|
\cgalFigureEnd
|
|
|
|
\section TutorialGIS_DTM Digital Terrain Model (DTM)
|
|
|
|
The generated DSM is used as a base for DTM computation, that is a
|
|
representation of the ground as another TIN after filtering non-ground
|
|
points.
|
|
|
|
We propose, as an example, a simple DTM estimation decomposed in the
|
|
following steps:
|
|
|
|
1. Thresholding the height of the facets to remove brutal changes of
|
|
elevation
|
|
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
|
|
corresponds to the minimum height of a building, and a perimeter
|
|
threshold that corresponds to the maximum size of a building on the 2D
|
|
projection.
|
|
|
|
\subsection TutorialGIS_DTM_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 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
|
|
|
|
\subsection TutorialGIS_DTM_components Identifying Connected Components
|
|
|
|
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 exceed the user-defined threshold.
|
|
|
|
\snippet Classification/gis_tutorial_example.cpp Components
|
|
|
|
This TIN enriched with connected component information can be saved as
|
|
a colored mesh:
|
|
|
|
\snippet Classification/gis_tutorial_example.cpp Save TIN with info
|
|
|
|
An example of a TIN colored by connected components is
|
|
given in \cgalFigureRef{TutorialGISFigComponents}.
|
|
|
|
\cgalFigureBegin{TutorialGISFigComponents, components.jpg}
|
|
TIN colored by connected components. Faces above height threshold are
|
|
not assigned to any component and are displayed in gray.
|
|
\cgalFigureEnd
|
|
|
|
\subsection TutorialGIS_DTM_filtering Filtering
|
|
|
|
Components smaller than the largest building can be removed this way:
|
|
|
|
\snippet Classification/gis_tutorial_example.cpp Filtering
|
|
|
|
\subsection Tutorial_DTM_hole_filling Hole Filling and Remeshing
|
|
|
|
Because simply removing vertices in the large areas covered by
|
|
buildings results in large Delaunay faces that offer a poor 3D
|
|
representation of the DTM, an additional step can help producing
|
|
better shaped meshes: faces larger than a threshold are removed and
|
|
filled with a hole filling algorithm that triangulates, refines and
|
|
fairs the holes.
|
|
|
|
The following snippet copies the TIN into a mesh while filtering out
|
|
overly large faces, then identifies the holes and fills them all
|
|
except for the largest one (which is the outer hull).
|
|
|
|
\snippet Classification/gis_tutorial_example.cpp Hole filling
|
|
|
|
Isotropic remeshing can also be performed as a final step in order to
|
|
produce a more regular mesh that is not constrained by the shape of 2D
|
|
Delaunay faces.
|
|
|
|
\snippet Classification/gis_tutorial_example.cpp Remeshing
|
|
|
|
\cgalFigureRef{TutorialGISFigProc} shows how these different steps
|
|
affect the output mesh and \cgalFigureRef{TutorialGISFigDTM} shows the
|
|
DTM isotropic mesh.
|
|
|
|
\cgalFigureBegin{TutorialGISFigProc, dtm_proc.jpg}
|
|
Raw DTM and postprocessing
|
|
\cgalFigureEnd
|
|
|
|
\cgalFigureBegin{TutorialGISFigDTM, dtm.jpg}
|
|
Final DTM.
|
|
\cgalFigureEnd
|
|
|
|
|
|
\section TutorialGIS_Raster Rastering
|
|
|
|
The TIN data structure can be combined with barycentric coordinates in
|
|
order to interpolate and thus rasterize a height map at any resolution
|
|
needed the information embedded in the vertices.
|
|
|
|
Because the latest two steps (hole filling and remeshing) were
|
|
performed on a 3D mesh, the hypothesis that our DTM is a 2.5D
|
|
representation may no longer be valid. We thus first rebuild a TIN
|
|
using the vertices of the isotropic DTM mesh lastly computed.
|
|
|
|
The following snippet generates a raster image of the height in the
|
|
form of rainbow ramp PPM file (simple bitmap format):
|
|
|
|
\snippet Classification/gis_tutorial_example.cpp Rastering
|
|
|
|
An example of a raster image with a rainbow ramp representing height
|
|
is given in \cgalFigureRef{TutorialGISFigRastering}.
|
|
|
|
\cgalFigureBegin{TutorialGISFigRastering, raster.jpg}
|
|
Raster visualization of height using a rainbow ramp, ranging from
|
|
light blue for low values to dark red for high values.
|
|
\cgalFigureEnd
|
|
|
|
\section TutorialGIS_Contour Contouring
|
|
|
|
Extracting isolevels of a function defined on a TIN is another
|
|
application that can be done with \cgal. We demonstrate here how to
|
|
extract isolevels of height to build a topographic map.
|
|
|
|
\subsection TutorialGIS_Contour_Extraction Building a Contour Graph
|
|
|
|
The first step is to extract, from all faces of the triangulation, the
|
|
section of each isolevel that goes through that face, in the form of a
|
|
segment. The following functions allow to test if one isovalue does
|
|
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
|
|
[boost::adjacency_list](https://www.boost.org/doc/libs/1_72_0/libs/graph/doc/adjacency_list.html)
|
|
structure 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 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 \link split_graph_into_polylines()
|
|
`CGAL::split_graph_into_polylines()` \endlink:
|
|
|
|
\snippet Classification/gis_tutorial_example.cpp Contouring split
|
|
|
|
This function requires a visitor which is called when starting a
|
|
polyline, when adding a point to it and when ending it. Defining such
|
|
a class is straightforward in our case:
|
|
|
|
\snippet Classification/gis_tutorial_example.cpp Contouring visitor
|
|
|
|
\subsection TutorialGIS_Contour_Simplifying Simplifying
|
|
|
|
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
|
|
`CGAL::Constrained_triangulation_plus_2`, which embeds
|
|
polylines as a set of constraints:
|
|
|
|
\snippet Classification/gis_tutorial_example.cpp CDT
|
|
|
|
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.
|
|
|
|
\snippet Classification/gis_tutorial_example.cpp Contouring simplify
|
|
|
|
Examples of contours and simplifications are given in
|
|
\cgalFigureRef{TutorialGISFigContours}.
|
|
|
|
\cgalFigureBegin{TutorialGISFigContours, contours.png}
|
|
Contouring using 50 isovalues evenly spaced. Top: original contouring
|
|
using 148k vertices and simplification with a tolerance equal to the
|
|
average spacing of the input point cloud (3.4\% of the original
|
|
polyline vertices remaining). Bottom: simplification with tolerance 4
|
|
times the average spacing (1.3\% of vertices remaining) and with
|
|
tolerance 10 times the average spacing (0.9\% of vertices
|
|
remaining). Polylines are intersection-free in all cases.
|
|
\cgalFigureEnd
|
|
|
|
\section TutorialGIS_Classify Classifying
|
|
|
|
\cgal provides a package %Classification 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
|
|
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
|
|
regularized by a graph cut algorithm:
|
|
|
|
\snippet Classification/gis_tutorial_example.cpp Classification
|
|
|
|
An example of training set and resulting classification is given in
|
|
\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 manually
|
|
classified slices.
|
|
\cgalFigureEnd
|
|
|
|
\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
|
|
used). We give a full code example which achieves all the steps
|
|
described in this tutorial.
|
|
|
|
\include Classification/gis_tutorial_example.cpp
|
|
|
|
\section TutorialGIS_Reference References
|
|
|
|
This tutorial is 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_.
|
|
|
|
*/
|
|
}
|