Author: Simon Giraudot

在GIS应用中常用的传感器（如激光雷达LIDAR）会生成密集的点云数据。这类应用通常需要使用更高级的数据结构，例如不规则三角网（Triangulated Irregular Network, TIN），它可以作为数字高程模型（Digital Elevation Model, DEM）的基础，特别是用于生成数字地形模型（Digital Terrain Model, DTM）。点云数据还可以通过分类信息进行扩充，将点分为地面、植被和建筑物等用户自定义的类别。

不同文献中对某些数据结构的定义可能有所不同。在本教程中，我们采用以下术语定义：

TIN：不规则三角网，是一种基于点在水平面上投影连接形成的二维三角剖分结构。
DSM：数字表面模型（Digital Surface Model），表示包含建筑物和植被在内的完整扫描表面模型。我们使用TIN来存储DSM。
DTM：数字地形模型，表示不包含建筑物和植被等地物的裸地表面模型。我们同时使用TIN和栅格来存储DTM。
DEM：数字高程模型，是一个更广泛的术语，包括DSM和DTM。

本教程将演示以下处理流程：首先从输入的点云数据计算出以TIN形式存储的DSM。然后，过滤掉对应于建筑物立面或植被噪声的过大面片。保留对应于地面的大型连通分量。填充空洞并对获得的DEM进行重新网格化。从中生成栅格DEM和等高线。最后，执行有监督的三类别分类，将点分为植被、建筑物和地面。

不规则三角网（TIN）

CGAL提供了多种三角剖分数据结构和算法。TIN可以通过结合二维Delaunay三角剖分和投影特征来生成：使用点在选定平面（通常是XY平面）上的二维位置计算三角剖分结构，同时保留点的三维位置用于可视化和测量。

因此，TIN数据结构可以简单地定义如下：

 
using Kernel = CGAL::Exact_predicates_inexact_constructions_kernel;
using Projection_traits = CGAL::Projection_traits_xy_3<Kernel>;
using Point_2 = Kernel::Point_2;
using Point_3 = Kernel::Point_3;
using Segment_3 = Kernel::Segment_3;
 
// Triangulated Irregular Network
using TIN = CGAL::Delaunay_triangulation_2<Projection_traits>;
 

数字表面模型（DSM）

可以使用流操作符轻松地将多种格式（XYZ、OFF、PLY、LAS）的点云数据加载到CGAL::Point_set_3结构中。然后可以直接生成存储为TIN的DSM：

 
  // Read points
  std::ifstream ifile (fname, std::ios_base::binary);
  CGAL::Point_set_3<Point_3> points;
  ifile >> points;
  std::cerr << points.size() << " point(s) read" << std::endl;
 
  // Create DSM
  TIN dsm (points.points().begin(), points.points().end());
 

由于CGAL的Delaunay三角剖分是FaceGraph的一个模型，因此可以直接将生成的TIN转换为CGAL::Surface_mesh等网格结构，并保存为这种结构支持的任何格式：

 
  using Mesh = CGAL::Surface_mesh<Point_3>;
 
  Mesh dsm_mesh;
  CGAL::copy_face_graph (dsm, dsm_mesh);
  std::ofstream dsm_ofile ("dsm.ply", std::ios_base::binary);
  CGAL::IO::set_binary_mode (dsm_ofile);
  CGAL::IO::write_PLY (dsm_ofile, dsm_mesh);
  dsm_ofile.close();
 

\ref TutorialGIS_Reference）上计算得到的DSM示例。

Figure 0.1 输入点云和输出DSM。

数字地形模型（DTM）

生成的DSM被用作DTM计算的基础，即在过滤掉非地面点后形成另一个TIN表示的地面。

作为示例，我们提供一个简单的DTM估计方法，包含以下步骤：

对面片高度进行阈值处理，去除高程突变
将其他面片聚类为连通分量
过滤掉小于用户定义阈值的所有分量

该算法依赖两个参数：高度阈值（对应建筑物的最小高度）和周长阈值（对应建筑物在二维投影上的最大尺寸）。

带信息的TIN

得益于CGAL Delaunay三角剖分灵活的API，我们的TIN可以在顶点和/或面上添加信息。在本例中，每个顶点记录输入点云中对应点的索引（这将用于后续过滤地面点），每个面则被赋予其连通分量的索引。

 
  auto idx_to_point_with_info
    = [&](const Point_set::Index& idx) -> std::pair<Point_3, Point_set::Index>
      {
        return std::make_pair (points.point(idx), idx);
      };
 
  TIN_with_info tin_with_info
    (boost::make_transform_iterator (points.begin(), idx_to_point_with_info),
     boost::make_transform_iterator (points.end(), idx_to_point_with_info));
 

识别连通分量

通过洪水填充算法识别连通分量：从一个种子面开始，将所有相邻面加入当前连通分量，除非它们的高度超过用户定义的阈值。

 
  double spacing = CGAL::compute_average_spacing<Concurrency_tag>(points, 6);
  spacing *= 2;
 
  auto face_height
    = [&](const TIN_with_info::Face_handle fh) -> double
      {
        double out = 0.;
        for (int i = 0; i < 3; ++ i)
          out = (std::max) (out, CGAL::abs(fh->vertex(i)->point().z() - fh->vertex((i+1)%3)->point().z()));
        return out;
      };
 
  // Initialize faces info
  for (TIN_with_info::Face_handle fh : tin_with_info.all_face_handles())
    if (tin_with_info.is_infinite(fh) || face_height(fh) > spacing) // Filtered faces are given info() = -2
      fh->info() = -2;
    else // Pending faces are given info() = -1;
      fh->info() = -1;
 
  // Flooding algorithm
  std::vector<int> component_size;
  for (TIN_with_info::Face_handle fh : tin_with_info.finite_face_handles())
  {
    if (fh->info() != -1)
      continue;
 
    std::queue<TIN_with_info::Face_handle> todo;
    todo.push(fh);
 
    int size = 0;
    while (!todo.empty())
    {
      TIN_with_info::Face_handle current = todo.front();
      todo.pop();
 
      if (current->info() != -1)
        continue;
      current->info() = int(component_size.size());
      ++ size;
 
      for (int i = 0; i < 3; ++ i)
        todo.push (current->neighbor(i));
    }
 
    component_size.push_back (size);
  }
 
  std::cerr << component_size.size() << " connected component(s) found" << std::endl;
 

这个带有连通分量信息的TIN可以保存为彩色网格：

 
  Mesh tin_colored_mesh;
 
  Mesh::Property_map<Mesh::Face_index, CGAL::IO::Color>
    color_map = tin_colored_mesh.add_property_map<Mesh::Face_index, CGAL::IO::Color>("f:color").first;
 
  CGAL::copy_face_graph (tin_with_info, tin_colored_mesh,
                         CGAL::parameters::face_to_face_output_iterator
                         (boost::make_function_output_iterator
                          ([&](const std::pair<TIN_with_info::Face_handle, Mesh::Face_index>& ff)
                           {
                             // Color unassigned faces gray
                             if (ff.first->info() < 0)
                               color_map[ff.second] = CGAL::IO::Color(128, 128, 128);
                             else
                             {
                               // Random color seeded by the component ID
                               CGAL::Random r (ff.first->info());
                               color_map[ff.second] = CGAL::IO::Color (r.get_int(64, 192),
                                                                   r.get_int(64, 192),
                                                                   r.get_int(64, 192));
                             }
                           })));
 
  std::ofstream tin_colored_ofile ("colored_tin.ply", std::ios_base::binary);
  CGAL::IO::set_binary_mode (tin_colored_ofile);
  CGAL::IO::write_PLY (tin_colored_ofile, tin_colored_mesh);
  tin_colored_ofile.close();
 

fig__TutorialGISFigComponents展示了一个按连通分量着色的TIN示例。

Figure 0.2 按连通分量着色的TIN。超过高度阈值的面片不属于任何分量，显示为灰色。

过滤

可以通过以下方式移除小于最大建筑物的分量：

 
  int min_size = int(points.size() / 2);
 
  std::vector<TIN_with_info::Vertex_handle> to_remove;
  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),
      start = circ;
 
    // Remove a vertex if it's only adjacent to components smaller than threshold
    bool keep = false;
    do
    {
      if (circ->info() >= 0 && component_size[std::size_t(circ->info())] > min_size)
      {
        keep = true;
        break;
      }
    }
    while (++ circ != start);
 
    if (!keep)
      to_remove.push_back (vh);
  }
 
  std::cerr << to_remove.size() << " vertices(s) will be removed after filtering" << std::endl;
  for (TIN_with_info::Vertex_handle vh : to_remove)
    tin_with_info.remove (vh);
 

空洞填充和重新网格化

由于简单地移除建筑物覆盖的大面积区域内的顶点会导致Delaunay面片过大，无法很好地表示DTM的3D特征，因此可以通过额外的步骤来生成更好的网格形状：移除大于阈值的面片，然后使用空洞填充算法对空洞进行三角剖分、细化和平滑处理。

以下代码片段将TIN复制到网格中，同时过滤掉过大的面片，然后识别空洞并填充除最大空洞（外部轮廓）之外的所有空洞。

 
  // Copy and keep track of overly large faces
  Mesh dtm_mesh;
 
  std::vector<Mesh::Face_index> face_selection;
  Mesh::Property_map<Mesh::Face_index, bool> face_selection_map
   = dtm_mesh.add_property_map<Mesh::Face_index, bool>("is_selected", false).first;
 
  double limit = CGAL::square (5 * spacing);
  CGAL::copy_face_graph (tin_with_info, dtm_mesh,
                         CGAL::parameters::face_to_face_output_iterator
                         (boost::make_function_output_iterator
                          ([&](const std::pair<TIN_with_info::Face_handle, Mesh::Face_index>& ff)
                           {
                             double longest_edge = 0.;
                             bool border = false;
                             for (int i = 0; i < 3; ++ i)
                             {
                               longest_edge = (std::max)(longest_edge, CGAL::squared_distance
                                                         (ff.first->vertex((i+1)%3)->point(),
                                                          ff.first->vertex((i+2)%3)->point()));
 
                               TIN_with_info::Face_circulator circ
                                 = tin_with_info.incident_faces (ff.first->vertex(i)),
                                 start = circ;
                               do
                               {
                                 if (tin_with_info.is_infinite (circ))
                                 {
                                   border = true;
                                   break;
                                 }
                               }
                               while (++ circ != start);
 
                               if (border)
                                 break;
                             }
 
                             // Select if face is too big AND it's not
                             // on the border (to have closed holes)
                             if (!border && longest_edge > limit)
                             {
                               face_selection_map[ff.second] = true;
                               face_selection.push_back (ff.second);
                             }
                           })));
 
  // Save original DTM
  std::ofstream dtm_ofile ("dtm.ply", std::ios_base::binary);
  CGAL::IO::set_binary_mode (dtm_ofile);
  CGAL::IO::write_PLY (dtm_ofile, dtm_mesh);
  dtm_ofile.close();
 
  std::cerr << face_selection.size() << " face(s) are selected for removal" << std::endl;
 
  // Expand face selection to keep a well formed 2-manifold mesh after removal
  CGAL::expand_face_selection_for_removal (face_selection, dtm_mesh, face_selection_map);
  face_selection.clear();
  for (Mesh::Face_index fi : faces(dtm_mesh))
    if (face_selection_map[fi])
      face_selection.push_back(fi);
 
  std::cerr << face_selection.size() << " face(s) are selected for removal after expansion" << std::endl;
 
  for (Mesh::Face_index fi : face_selection)
    CGAL::Euler::remove_face (halfedge(fi, dtm_mesh), dtm_mesh);
  dtm_mesh.collect_garbage();
 
  if (!dtm_mesh.is_valid())
    std::cerr << "Invalid mesh!" << std::endl;
 
  // Save filtered DTM
  std::ofstream dtm_holes_ofile ("dtm_with_holes.ply", std::ios_base::binary);
  CGAL::IO::set_binary_mode (dtm_holes_ofile);
  CGAL::IO::write_PLY (dtm_holes_ofile, dtm_mesh);
  dtm_holes_ofile.close();
 
  // Get all holes
  std::vector<Mesh::Halfedge_index> holes;
  CGAL::Polygon_mesh_processing::extract_boundary_cycles (dtm_mesh, std::back_inserter (holes));
 
  std::cerr << holes.size() << " hole(s) identified" << std::endl;
 
  // Identify outer hull (hole with maximum size)
  double max_size = 0.;
  Mesh::Halfedge_index outer_hull;
  for (Mesh::Halfedge_index hi : holes)
  {
    CGAL::Bbox_3 hole_bbox;
    for (Mesh::Halfedge_index haf : CGAL::halfedges_around_face(hi, dtm_mesh))
    {
      const Point_3& p = dtm_mesh.point(target(haf, dtm_mesh));
      hole_bbox += p.bbox();
    }
    double size = CGAL::squared_distance (Point_2(hole_bbox.xmin(), hole_bbox.ymin()),
                                          Point_2(hole_bbox.xmax(), hole_bbox.ymax()));
    if (size > max_size)
    {
      max_size = size;
      outer_hull = hi;
    }
  }
 
  // Fill all holes except the biggest (which is the outer hull of the mesh)
  for (Mesh::Halfedge_index hi : holes)
    if (hi != outer_hull)
      CGAL::Polygon_mesh_processing::triangulate_refine_and_fair_hole
         (dtm_mesh, hi, CGAL::parameters::fairing_continuity(0));
 
  // Save DTM with holes filled
  std::ofstream dtm_filled_ofile ("dtm_filled.ply", std::ios_base::binary);
  CGAL::IO::set_binary_mode (dtm_filled_ofile);
  CGAL::IO::write_PLY (dtm_filled_ofile, dtm_mesh);
  dtm_filled_ofile.close();
 

最后还可以进行各向同性重网格化，以生成更规则的网格，不受2D Delaunay面片形状的约束。

 
  CGAL::Polygon_mesh_processing::isotropic_remeshing (faces(dtm_mesh), spacing, dtm_mesh);
 
  std::ofstream dtm_remeshed_ofile ("dtm_remeshed.ply", std::ios_base::binary);
  CGAL::IO::set_binary_mode (dtm_remeshed_ofile);
  CGAL::IO::write_PLY (dtm_remeshed_ofile, dtm_mesh);
  dtm_remeshed_ofile.close();
 

\ref fig__TutorialGISFigDTM展示了DTM各向同性网格。

Figure 0.3 原始DTM和后处理

Figure 0.4 最终DTM。

栅格化

TIN数据结构可以与重心坐标结合使用，以便在任何需要的分辨率下对顶点中嵌入的信息进行插值和栅格化高程图。

由于最后两个步骤（空洞填充和重网格化）是在3D网格上执行的，我们的DTM作为2.5D表示的假设可能不再有效。因此，我们首先使用最后计算的各向同性DTM网格的顶点重建TIN。

以下代码片段生成一个使用彩虹渐变的PPM格式（简单位图格式）高程栅格图像：

 
  CGAL::Bbox_3 bbox = CGAL::bbox_3 (points.points().begin(), points.points().end());
 
  // Generate raster image 1920-pixels large
  std::size_t width = 1920;
  std::size_t height = std::size_t((bbox.ymax() - bbox.ymin()) * 1920 / (bbox.xmax() - bbox.xmin()));
 
  std::cerr << "Rastering with resolution " << width << "x" << height << std::endl;
 
  // Use PPM format (Portable PixMap) for simplicity
  std::ofstream raster_ofile ("raster.ppm", std::ios_base::binary);
 
  // 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::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)
    {
      Point_3 query (bbox.xmin() + x * (bbox.xmax() - bbox.xmin()) / double(width),
                     bbox.ymin() + (height-y) * (bbox.ymax() - bbox.ymin()) / double(height),
                     0); // not relevant for location in 2D
 
      location = dtm_clean.locate (query, location);
 
      // Points outside the convex hull will be colored black
      std::array<unsigned char, 3> colors { 0, 0, 0 };
      if (!dtm_clean.is_infinite(location))
      {
        std::array<double, 3> barycentric_coordinates
          = CGAL::Polygon_mesh_processing::barycentric_coordinates
          (Point_2 (location->vertex(0)->point().x(), location->vertex(0)->point().y()),
           Point_2 (location->vertex(1)->point().x(), location->vertex(1)->point().y()),
           Point_2 (location->vertex(2)->point().x(), location->vertex(2)->point().y()),
           Point_2 (query.x(), query.y()),
           Kernel());
 
        double height_at_query
          = (barycentric_coordinates[0] * location->vertex(0)->point().z()
             + 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 (reinterpret_cast<char*>(&colors), 3);
    }
 
  raster_ofile.close();
 

fig__TutorialGISFigRastering展示了一个使用彩虹渐变表示高程的栅格图像示例。

Figure 0.5 使用彩虹渐变可视化高程，从浅蓝色（低值）到深红色（高值）。

等高线生成

CGAL还可以用于提取TIN上定义函数的等值线。这里我们演示如何提取高程的等值线来构建地形图。

构建等高线图

第一步是从三角剖分的所有面中提取每个等值线穿过该面的部分，形式为线段。以下函数用于测试一个等值是否穿过面片，并提取等值线：

 
bool face_has_isovalue (TIN::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)
      below = true;
  }
 
  return (above && below);
}
 
Segment_3 isocontour_in_face (TIN::Face_handle fh, double isovalue)
{
  Point_3 source;
  Point_3 target;
  bool source_found = false;
 
  for (int i = 0; i < 3; ++ i)
  {
    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;
 
    double zbottom = p0.z();
    double ztop = p1.z();
    if (zbottom > ztop)
    {
      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)
      target = p;
    else
    {
      source = p;
      source_found = true;
    }
  }
 
  return Segment_3 (source, target);
}
 
 

利用这些函数，我们可以创建一个线段图，后续将其处理成一组折线。为此，我们使用boost::adjacency_list结构，并维护端点位置到图顶点的映射。

以下代码计算点云最小和最大高度之间均匀分布的50个等值，并创建包含所有等值线的图：

 
  std::array<double, 50> isovalues; // Contour 50 isovalues
  for (std::size_t i = 0; i < isovalues.size(); ++ i)
    isovalues[i] = bbox.zmin() + ((i+1) * (bbox.zmax() - bbox.zmin()) / (isovalues.size() - 2));
 
  // First find on each face if they are crossed by some isovalues and
  // extract segments in a graph
  using Segment_graph = boost::adjacency_list<boost::listS, boost::vecS, boost::undirectedS, Point_3>;
  Segment_graph graph;
  using Map_p2v = std::map<Point_3, Segment_graph::vertex_descriptor>;
  Map_p2v map_p2v;
  for (TIN::Face_handle vh : dtm_clean.finite_face_handles())
    for (double iv : isovalues)
      if (face_has_isovalue (vh, iv))
      {
        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()));
          if (inserted)
          {
            iter->second = boost::add_vertex (graph);
            graph[iter->second] = p;
          }
        }
        boost::add_edge (map_p2v[segment.source()], map_p2v[segment.target()], graph);
      }
 

分割为折线

创建图后，使用函数CGAL::split_graph_into_polylines() 可以轻松地将其分割为折线：

 
  // Split segments into polylines
  std::vector<std::vector<Point_3> > polylines;
  Polylines_visitor<Segment_graph> visitor (graph, polylines);
  CGAL::split_graph_into_polylines (graph, visitor);
 
  std::cerr << polylines.size() << " polylines computed, with "
            << map_p2v.size() << " vertices in total" << std::endl;
 
  // Output to WKT file
  std::ofstream contour_ofile ("contour.wkt");
  contour_ofile.precision(18);
  CGAL::IO::write_multi_linestring_WKT (contour_ofile, polylines);
  contour_ofile.close();
 

该函数需要一个访问器，在开始折线、向其添加点和结束时被调用。在我们的例子中，定义这样的类很简单：

 
template <typename Graph>
class Polylines_visitor
{
private:
  std::vector<std::vector<Point_3> >& polylines;
  Graph& graph;
 
public:
 
  Polylines_visitor (Graph& graph, std::vector<std::vector<Point_3> >& polylines)
    : polylines (polylines), graph(graph) { }
 
  void start_new_polyline()
  {
    polylines.push_back (std::vector<Point_3>());
  }
 
  void add_node (typename Graph::vertex_descriptor vd)
  {
    polylines.back().push_back (graph[vd]);
  }
 
  void end_polyline()
  {
    // filter small polylines
    if (polylines.back().size() < 50)
      polylines.pop_back();
  }
};
 

简化

由于输出可能比较嘈杂，用户可能希望简化折线。CGAL提供了一个折线简化算法，可以保证简化后的折线不会相交。该算法利用CGAL::Constrained_triangulation_plus_2，将折线作为约束集合嵌入：

 
namespace PS = CGAL::Polyline_simplification_2;
using CDT_vertex_base = PS::Vertex_base_2<Projection_traits>;
using CDT_face_base = CGAL::Constrained_triangulation_face_base_2<Projection_traits>;
using CDT_TDS = CGAL::Triangulation_data_structure_2<CDT_vertex_base, CDT_face_base>;
using CDT = CGAL::Constrained_Delaunay_triangulation_2<Projection_traits, CDT_TDS>;
using CTP = CGAL::Constrained_triangulation_plus_2<CDT>;
 

以下代码基于到原始折线的平方距离进行简化，当没有更多顶点可以在不超过平均间距4倍的情况下移除时停止。

 
  // Construct constrained Delaunay triangulation with polylines as constraints
  CTP ctp;
  for (const std::vector<Point_3>& poly : polylines)
    ctp.insert_constraint (poly.begin(), poly.end());
 
  // Simplification algorithm with limit on distance
  PS::simplify (ctp, PS::Squared_distance_cost(), PS::Stop_above_cost_threshold (16 * spacing * spacing));
 
  polylines.clear();
  for (CTP::Constraint_id cid : ctp.constraints())
  {
    polylines.push_back (std::vector<Point_3>());
    polylines.back().reserve (ctp.vertices_in_constraint (cid).size());
    for (CTP::Vertex_handle vh : ctp.vertices_in_constraint(cid))
      polylines.back().push_back (vh->point());
  }
 
  std::size_t nb_vertices
    = std::accumulate (polylines.begin(), polylines.end(), std::size_t(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::IO::write_multi_linestring_WKT (simplified_ofile, polylines);
  simplified_ofile.close();
 

fig__TutorialGISFigContours展示了等高线及其简化示例。

Figure 0.6 使用50个均匀间隔的等值生成等高线。上图：原始等高线（148k个顶点）和使用等于输入点云平均间距的容差进行简化（保留3.4的原始折线顶点）。下图：使用4倍平均间距的容差（保留1.3的顶点）和10倍平均间距的容差（保留0.9的顶点）进行简化。所有情况下的折线都不相交。

分类

CGAL提供了分类包，可用于将点云分割为用户定义的标签集。目前CGAL中可用的最先进分类器是来自ETHZ的随机森林。作为一个有监督分类器，它需要训练集。

以下代码片段展示了如何使用手动选择的训练集来训练随机森林分类器，并计算通过图割算法正则化的分类结果：

 
  // Get training from input
  std::optional<Point_set::Property_map<int>> training_map = points.property_map<int>("training");
 
  if (training_map.has_value())
  {
    std::cerr << "Classifying ground/vegetation/building" << std::endl;
 
    // Create labels
    Classification::Label_set labels ({ "ground", "vegetation", "building" });
 
    // Generate features
    Classification::Feature_set features;
    Classification::Point_set_feature_generator<Kernel, Point_set, Point_set::Point_map>
      generator (points, points.point_map(), 5); // 5 scales
 
#ifdef CGAL_LINKED_WITH_TBB
    // If TBB is used, features can be computed in parallel
    features.begin_parallel_additions();
    generator.generate_point_based_features (features);
    features.end_parallel_additions();
#else
    generator.generate_point_based_features (features);
#endif
 
    // Train a random forest classifier
    Classification::ETHZ::Random_forest_classifier classifier (labels, features);
    classifier.train (points.range(training_map.value()));
 
    // Classify with graphcut regularization
    Point_set::Property_map<int> label_map = points.add_property_map<int>("labels").first;
    Classification::classify_with_graphcut<Concurrency_tag>
      (points, points.point_map(), labels, classifier,
       generator.neighborhood().k_neighbor_query(12), // regularize on 12-neighbors graph
       0.5f, // graphcut weight
       12, // Subdivide to speed-up process
       label_map);
 
    // Evaluate
    std::cerr << "Mean IoU on training data = "
              << Classification::Evaluation(labels,
                                            points.range(training_map.value()),
                                            points.range(label_map)).mean_intersection_over_union() << std::endl;
 
    // Save the classified point set
    std::ofstream classified_ofile ("classification_gis_tutorial.ply");
    CGAL::IO::set_binary_mode (classified_ofile);
    classified_ofile << points;
    classified_ofile.close();
  }
 

fig__TutorialGISFigClassif展示了训练集和分类结果示例。

Figure 0.7 上图：手动分类的点云切片。下图：在3个手动分类切片上训练后，通过图割正则化的分类结果。

完整代码示例

本教程中使用的所有代码片段可以组合成一个完整的GIS处理流程（需要包含正确的头文件）。我们提供了一个实现本教程所有步骤的完整代码示例。

#include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
#include <CGAL/Projection_traits_xy_3.h>
 
#include <CGAL/Delaunay_triangulation_2.h>
#include <CGAL/Triangulation_vertex_base_with_info_2.h>
#include <CGAL/Triangulation_face_base_with_info_2.h>
 
#include <CGAL/boost/graph/graph_traits_Delaunay_triangulation_2.h>
#include <CGAL/boost/graph/copy_face_graph.h>
 
#include <CGAL/Point_set_3.h>
#include <CGAL/Point_set_3/IO.h>
#include <CGAL/compute_average_spacing.h>
 
#include <CGAL/Surface_mesh.h>
#include <CGAL/Polygon_mesh_processing/locate.h>
 
#include <CGAL/Polygon_mesh_processing/triangulate_hole.h>
#include <CGAL/Polygon_mesh_processing/border.h>
#include <CGAL/Polygon_mesh_processing/remesh.h>
 
#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>
 
#include <CGAL/Polyline_simplification_2/simplify.h>
#include <CGAL/Polyline_simplification_2/Squared_distance_cost.h>
 
#include <CGAL/Classification.h>
 
#include <CGAL/Random.h>
 
#include <fstream>
#include <queue>
 
#include "include/Color_ramp.h"
 
 
using Kernel = CGAL::Exact_predicates_inexact_constructions_kernel;
using Projection_traits = CGAL::Projection_traits_xy_3<Kernel>;
using Point_2 = Kernel::Point_2;
using Point_3 = Kernel::Point_3;
using Segment_3 = Kernel::Segment_3;
 
// Triangulated Irregular Network
using TIN = CGAL::Delaunay_triangulation_2<Projection_traits>;
 
 
 
// Triangulated Irregular Network (with info)
using Point_set = CGAL::Point_set_3<Point_3>;
using Vbi = CGAL::Triangulation_vertex_base_with_info_2 <Point_set::Index, Projection_traits>;
using Fbi = CGAL::Triangulation_face_base_with_info_2<int, Projection_traits>;
using TDS = CGAL::Triangulation_data_structure_2<Vbi, Fbi>;
using TIN_with_info = CGAL::Delaunay_triangulation_2<Projection_traits, TDS>;
 
 
namespace Classification = CGAL::Classification;
 
#ifdef CGAL_LINKED_WITH_TBB
using Concurrency_tag = CGAL::Parallel_tag;
#else
using Concurrency_tag = CGAL::Sequential_tag;
#endif
 
 
bool face_has_isovalue (TIN::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)
      below = true;
  }
 
  return (above && below);
}
 
Segment_3 isocontour_in_face (TIN::Face_handle fh, double isovalue)
{
  Point_3 source;
  Point_3 target;
  bool source_found = false;
 
  for (int i = 0; i < 3; ++ i)
  {
    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;
 
    double zbottom = p0.z();
    double ztop = p1.z();
    if (zbottom > ztop)
    {
      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)
      target = p;
    else
    {
      source = p;
      source_found = true;
    }
  }
 
  return Segment_3 (source, target);
}
 
 
 
 
template <typename Graph>
class Polylines_visitor
{
private:
  std::vector<std::vector<Point_3> >& polylines;
  Graph& graph;
 
public:
 
  Polylines_visitor (Graph& graph, std::vector<std::vector<Point_3> >& polylines)
    : polylines (polylines), graph(graph) { }
 
  void start_new_polyline()
  {
    polylines.push_back (std::vector<Point_3>());
  }
 
  void add_node (typename Graph::vertex_descriptor vd)
  {
    polylines.back().push_back (graph[vd]);
  }
 
  void end_polyline()
  {
    // filter small polylines
    if (polylines.back().size() < 50)
      polylines.pop_back();
  }
};
 
 
 
namespace PS = CGAL::Polyline_simplification_2;
using CDT_vertex_base = PS::Vertex_base_2<Projection_traits>;
using CDT_face_base = CGAL::Constrained_triangulation_face_base_2<Projection_traits>;
using CDT_TDS = CGAL::Triangulation_data_structure_2<CDT_vertex_base, CDT_face_base>;
using CDT = CGAL::Constrained_Delaunay_triangulation_2<Projection_traits, CDT_TDS>;
using CTP = CGAL::Constrained_triangulation_plus_2<CDT>;
 
 
int main (int argc, char** argv)
{
  const std::string fname = argc != 2 ? CGAL::data_file_path("points_3/b9_training.ply") : argv[1];
  if (argc != 2)
  {
    std::cerr << "Usage: " << argv[0] << " points.ply" << std::endl;
    std::cerr << "Running with default value " << fname << "\n";
  }
 
 
  // Read points
  std::ifstream ifile (fname, std::ios_base::binary);
  CGAL::Point_set_3<Point_3> points;
  ifile >> points;
  std::cerr << points.size() << " point(s) read" << std::endl;
 
  // Create DSM
  TIN dsm (points.points().begin(), points.points().end());
 
 
 
  using Mesh = CGAL::Surface_mesh<Point_3>;
 
  Mesh dsm_mesh;
  CGAL::copy_face_graph (dsm, dsm_mesh);
  std::ofstream dsm_ofile ("dsm.ply", std::ios_base::binary);
  CGAL::IO::set_binary_mode (dsm_ofile);
  CGAL::IO::write_PLY (dsm_ofile, dsm_mesh);
  dsm_ofile.close();
 
 
 
  auto idx_to_point_with_info
    = [&](const Point_set::Index& idx) -> std::pair<Point_3, Point_set::Index>
      {
        return std::make_pair (points.point(idx), idx);
      };
 
  TIN_with_info tin_with_info
    (boost::make_transform_iterator (points.begin(), idx_to_point_with_info),
     boost::make_transform_iterator (points.end(), idx_to_point_with_info));
 
 
 
  double spacing = CGAL::compute_average_spacing<Concurrency_tag>(points, 6);
  spacing *= 2;
 
  auto face_height
    = [&](const TIN_with_info::Face_handle fh) -> double
      {
        double out = 0.;
        for (int i = 0; i < 3; ++ i)
          out = (std::max) (out, CGAL::abs(fh->vertex(i)->point().z() - fh->vertex((i+1)%3)->point().z()));
        return out;
      };
 
  // Initialize faces info
  for (TIN_with_info::Face_handle fh : tin_with_info.all_face_handles())
    if (tin_with_info.is_infinite(fh) || face_height(fh) > spacing) // Filtered faces are given info() = -2
      fh->info() = -2;
    else // Pending faces are given info() = -1;
      fh->info() = -1;
 
  // Flooding algorithm
  std::vector<int> component_size;
  for (TIN_with_info::Face_handle fh : tin_with_info.finite_face_handles())
  {
    if (fh->info() != -1)
      continue;
 
    std::queue<TIN_with_info::Face_handle> todo;
    todo.push(fh);
 
    int size = 0;
    while (!todo.empty())
    {
      TIN_with_info::Face_handle current = todo.front();
      todo.pop();
 
      if (current->info() != -1)
        continue;
      current->info() = int(component_size.size());
      ++ size;
 
      for (int i = 0; i < 3; ++ i)
        todo.push (current->neighbor(i));
    }
 
    component_size.push_back (size);
  }
 
  std::cerr << component_size.size() << " connected component(s) found" << std::endl;
 
 
 
  Mesh tin_colored_mesh;
 
  Mesh::Property_map<Mesh::Face_index, CGAL::IO::Color>
    color_map = tin_colored_mesh.add_property_map<Mesh::Face_index, CGAL::IO::Color>("f:color").first;
 
  CGAL::copy_face_graph (tin_with_info, tin_colored_mesh,
                         CGAL::parameters::face_to_face_output_iterator
                         (boost::make_function_output_iterator
                          ([&](const std::pair<TIN_with_info::Face_handle, Mesh::Face_index>& ff)
                           {
                             // Color unassigned faces gray
                             if (ff.first->info() < 0)
                               color_map[ff.second] = CGAL::IO::Color(128, 128, 128);
                             else
                             {
                               // Random color seeded by the component ID
                               CGAL::Random r (ff.first->info());
                               color_map[ff.second] = CGAL::IO::Color (r.get_int(64, 192),
                                                                   r.get_int(64, 192),
                                                                   r.get_int(64, 192));
                             }
                           })));
 
  std::ofstream tin_colored_ofile ("colored_tin.ply", std::ios_base::binary);
  CGAL::IO::set_binary_mode (tin_colored_ofile);
  CGAL::IO::write_PLY (tin_colored_ofile, tin_colored_mesh);
  tin_colored_ofile.close();
 
 
 
  int min_size = int(points.size() / 2);
 
  std::vector<TIN_with_info::Vertex_handle> to_remove;
  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),
      start = circ;
 
    // Remove a vertex if it's only adjacent to components smaller than threshold
    bool keep = false;
    do
    {
      if (circ->info() >= 0 && component_size[std::size_t(circ->info())] > min_size)
      {
        keep = true;
        break;
      }
    }
    while (++ circ != start);
 
    if (!keep)
      to_remove.push_back (vh);
  }
 
  std::cerr << to_remove.size() << " vertices(s) will be removed after filtering" << std::endl;
  for (TIN_with_info::Vertex_handle vh : to_remove)
    tin_with_info.remove (vh);
 
 
 
  // Copy and keep track of overly large faces
  Mesh dtm_mesh;
 
  std::vector<Mesh::Face_index> face_selection;
  Mesh::Property_map<Mesh::Face_index, bool> face_selection_map
   = dtm_mesh.add_property_map<Mesh::Face_index, bool>("is_selected", false).first;
 
  double limit = CGAL::square (5 * spacing);
  CGAL::copy_face_graph (tin_with_info, dtm_mesh,
                         CGAL::parameters::face_to_face_output_iterator
                         (boost::make_function_output_iterator
                          ([&](const std::pair<TIN_with_info::Face_handle, Mesh::Face_index>& ff)
                           {
                             double longest_edge = 0.;
                             bool border = false;
                             for (int i = 0; i < 3; ++ i)
                             {
                               longest_edge = (std::max)(longest_edge, CGAL::squared_distance
                                                         (ff.first->vertex((i+1)%3)->point(),
                                                          ff.first->vertex((i+2)%3)->point()));
 
                               TIN_with_info::Face_circulator circ
                                 = tin_with_info.incident_faces (ff.first->vertex(i)),
                                 start = circ;
                               do
                               {
                                 if (tin_with_info.is_infinite (circ))
                                 {
                                   border = true;
                                   break;
                                 }
                               }
                               while (++ circ != start);
 
                               if (border)
                                 break;
                             }
 
                             // Select if face is too big AND it's not
                             // on the border (to have closed holes)
                             if (!border && longest_edge > limit)
                             {
                               face_selection_map[ff.second] = true;
                               face_selection.push_back (ff.second);
                             }
                           })));
 
  // Save original DTM
  std::ofstream dtm_ofile ("dtm.ply", std::ios_base::binary);
  CGAL::IO::set_binary_mode (dtm_ofile);
  CGAL::IO::write_PLY (dtm_ofile, dtm_mesh);
  dtm_ofile.close();
 
  std::cerr << face_selection.size() << " face(s) are selected for removal" << std::endl;
 
  // Expand face selection to keep a well formed 2-manifold mesh after removal
  CGAL::expand_face_selection_for_removal (face_selection, dtm_mesh, face_selection_map);
  face_selection.clear();
  for (Mesh::Face_index fi : faces(dtm_mesh))
    if (face_selection_map[fi])
      face_selection.push_back(fi);
 
  std::cerr << face_selection.size() << " face(s) are selected for removal after expansion" << std::endl;
 
  for (Mesh::Face_index fi : face_selection)
    CGAL::Euler::remove_face (halfedge(fi, dtm_mesh), dtm_mesh);
  dtm_mesh.collect_garbage();
 
  if (!dtm_mesh.is_valid())
    std::cerr << "Invalid mesh!" << std::endl;
 
  // Save filtered DTM
  std::ofstream dtm_holes_ofile ("dtm_with_holes.ply", std::ios_base::binary);
  CGAL::IO::set_binary_mode (dtm_holes_ofile);
  CGAL::IO::write_PLY (dtm_holes_ofile, dtm_mesh);
  dtm_holes_ofile.close();
 
  // Get all holes
  std::vector<Mesh::Halfedge_index> holes;
  CGAL::Polygon_mesh_processing::extract_boundary_cycles (dtm_mesh, std::back_inserter (holes));
 
  std::cerr << holes.size() << " hole(s) identified" << std::endl;
 
  // Identify outer hull (hole with maximum size)
  double max_size = 0.;
  Mesh::Halfedge_index outer_hull;
  for (Mesh::Halfedge_index hi : holes)
  {
    CGAL::Bbox_3 hole_bbox;
    for (Mesh::Halfedge_index haf : CGAL::halfedges_around_face(hi, dtm_mesh))
    {
      const Point_3& p = dtm_mesh.point(target(haf, dtm_mesh));
      hole_bbox += p.bbox();
    }
    double size = CGAL::squared_distance (Point_2(hole_bbox.xmin(), hole_bbox.ymin()),
                                          Point_2(hole_bbox.xmax(), hole_bbox.ymax()));
    if (size > max_size)
    {
      max_size = size;
      outer_hull = hi;
    }
  }
 
  // Fill all holes except the biggest (which is the outer hull of the mesh)
  for (Mesh::Halfedge_index hi : holes)
    if (hi != outer_hull)
      CGAL::Polygon_mesh_processing::triangulate_refine_and_fair_hole
         (dtm_mesh, hi, CGAL::parameters::fairing_continuity(0));
 
  // Save DTM with holes filled
  std::ofstream dtm_filled_ofile ("dtm_filled.ply", std::ios_base::binary);
  CGAL::IO::set_binary_mode (dtm_filled_ofile);
  CGAL::IO::write_PLY (dtm_filled_ofile, dtm_mesh);
  dtm_filled_ofile.close();
 
 
 
  CGAL::Polygon_mesh_processing::isotropic_remeshing (faces(dtm_mesh), spacing, dtm_mesh);
 
  std::ofstream dtm_remeshed_ofile ("dtm_remeshed.ply", std::ios_base::binary);
  CGAL::IO::set_binary_mode (dtm_remeshed_ofile);
  CGAL::IO::write_PLY (dtm_remeshed_ofile, dtm_mesh);
  dtm_remeshed_ofile.close();
 
 
  TIN dtm_clean (dtm_mesh.points().begin(), dtm_mesh.points().end());
 
 
  CGAL::Bbox_3 bbox = CGAL::bbox_3 (points.points().begin(), points.points().end());
 
  // Generate raster image 1920-pixels large
  std::size_t width = 1920;
  std::size_t height = std::size_t((bbox.ymax() - bbox.ymin()) * 1920 / (bbox.xmax() - bbox.xmin()));
 
  std::cerr << "Rastering with resolution " << width << "x" << height << std::endl;
 
  // Use PPM format (Portable PixMap) for simplicity
  std::ofstream raster_ofile ("raster.ppm", std::ios_base::binary);
 
  // 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::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)
    {
      Point_3 query (bbox.xmin() + x * (bbox.xmax() - bbox.xmin()) / double(width),
                     bbox.ymin() + (height-y) * (bbox.ymax() - bbox.ymin()) / double(height),
                     0); // not relevant for location in 2D
 
      location = dtm_clean.locate (query, location);
 
      // Points outside the convex hull will be colored black
      std::array<unsigned char, 3> colors { 0, 0, 0 };
      if (!dtm_clean.is_infinite(location))
      {
        std::array<double, 3> barycentric_coordinates
          = CGAL::Polygon_mesh_processing::barycentric_coordinates
          (Point_2 (location->vertex(0)->point().x(), location->vertex(0)->point().y()),
           Point_2 (location->vertex(1)->point().x(), location->vertex(1)->point().y()),
           Point_2 (location->vertex(2)->point().x(), location->vertex(2)->point().y()),
           Point_2 (query.x(), query.y()),
           Kernel());
 
        double height_at_query
          = (barycentric_coordinates[0] * location->vertex(0)->point().z()
             + 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 (reinterpret_cast<char*>(&colors), 3);
    }
 
  raster_ofile.close();
 
 
  // Smooth heights with 5 successive Gaussian filters
  double gaussian_variance = 4 * spacing * spacing;
  for (TIN::Vertex_handle vh : dtm_clean.finite_vertex_handles())
  {
    double z = vh->point().z();
    double total_weight = 1;
 
    TIN::Vertex_circulator circ = dtm_clean.incident_vertices (vh),
      start = circ;
 
    do
    {
      if (!dtm_clean.is_infinite(circ))
      {
        double sq_dist = CGAL::squared_distance (vh->point(), circ->point());
 
        double weight = std::exp(- sq_dist / gaussian_variance);
        z += weight * circ->point().z();
        total_weight += weight;
      }
    }
    while (++ circ != start);
 
    z /= total_weight;
 
    vh->point() = Point_3 (vh->point().x(), vh->point().y(), z);
  }
 
 
  std::array<double, 50> isovalues; // Contour 50 isovalues
  for (std::size_t i = 0; i < isovalues.size(); ++ i)
    isovalues[i] = bbox.zmin() + ((i+1) * (bbox.zmax() - bbox.zmin()) / (isovalues.size() - 2));
 
  // First find on each face if they are crossed by some isovalues and
  // extract segments in a graph
  using Segment_graph = boost::adjacency_list<boost::listS, boost::vecS, boost::undirectedS, Point_3>;
  Segment_graph graph;
  using Map_p2v = std::map<Point_3, Segment_graph::vertex_descriptor>;
  Map_p2v map_p2v;
  for (TIN::Face_handle vh : dtm_clean.finite_face_handles())
    for (double iv : isovalues)
      if (face_has_isovalue (vh, iv))
      {
        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()));
          if (inserted)
          {
            iter->second = boost::add_vertex (graph);
            graph[iter->second] = p;
          }
        }
        boost::add_edge (map_p2v[segment.source()], map_p2v[segment.target()], graph);
      }
 
 
 
  // Split segments into polylines
  std::vector<std::vector<Point_3> > polylines;
  Polylines_visitor<Segment_graph> visitor (graph, polylines);
  CGAL::split_graph_into_polylines (graph, visitor);
 
  std::cerr << polylines.size() << " polylines computed, with "
            << map_p2v.size() << " vertices in total" << std::endl;
 
  // Output to WKT file
  std::ofstream contour_ofile ("contour.wkt");
  contour_ofile.precision(18);
  CGAL::IO::write_multi_linestring_WKT (contour_ofile, polylines);
  contour_ofile.close();
 
 
 
  // Construct constrained Delaunay triangulation with polylines as constraints
  CTP ctp;
  for (const std::vector<Point_3>& poly : polylines)
    ctp.insert_constraint (poly.begin(), poly.end());
 
  // Simplification algorithm with limit on distance
  PS::simplify (ctp, PS::Squared_distance_cost(), PS::Stop_above_cost_threshold (16 * spacing * spacing));
 
  polylines.clear();
  for (CTP::Constraint_id cid : ctp.constraints())
  {
    polylines.push_back (std::vector<Point_3>());
    polylines.back().reserve (ctp.vertices_in_constraint (cid).size());
    for (CTP::Vertex_handle vh : ctp.vertices_in_constraint(cid))
      polylines.back().push_back (vh->point());
  }
 
  std::size_t nb_vertices
    = std::accumulate (polylines.begin(), polylines.end(), std::size_t(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::IO::write_multi_linestring_WKT (simplified_ofile, polylines);
  simplified_ofile.close();
 
 
 
  // Get training from input
  std::optional<Point_set::Property_map<int>> training_map = points.property_map<int>("training");
 
  if (training_map.has_value())
  {
    std::cerr << "Classifying ground/vegetation/building" << std::endl;
 
    // Create labels
    Classification::Label_set labels ({ "ground", "vegetation", "building" });
 
    // Generate features
    Classification::Feature_set features;
    Classification::Point_set_feature_generator<Kernel, Point_set, Point_set::Point_map>
      generator (points, points.point_map(), 5); // 5 scales
 
#ifdef CGAL_LINKED_WITH_TBB
    // If TBB is used, features can be computed in parallel
    features.begin_parallel_additions();
    generator.generate_point_based_features (features);
    features.end_parallel_additions();
#else
    generator.generate_point_based_features (features);
#endif
 
    // Train a random forest classifier
    Classification::ETHZ::Random_forest_classifier classifier (labels, features);
    classifier.train (points.range(training_map.value()));
 
    // Classify with graphcut regularization
    Point_set::Property_map<int> label_map = points.add_property_map<int>("labels").first;
    Classification::classify_with_graphcut<Concurrency_tag>
      (points, points.point_map(), labels, classifier,
       generator.neighborhood().k_neighbor_query(12), // regularize on 12-neighbors graph
       0.5f, // graphcut weight
       12, // Subdivide to speed-up process
       label_map);
 
    // Evaluate
    std::cerr << "Mean IoU on training data = "
              << Classification::Evaluation(labels,
                                            points.range(training_map.value()),
                                            points.range(label_map)).mean_intersection_over_union() << std::endl;
 
    // Save the classified point set
    std::ofstream classified_ofile ("classification_gis_tutorial.ply");
    CGAL::IO::set_binary_mode (classified_ofile);
    classified_ofile << points;
    classified_ofile.close();
  }
 
 
 
  return EXIT_SUCCESS;
}

参考文献

本教程基于以下CGAL包：

本教程使用的数据集来自https://www.usgs.gov/数据库，采用美国政府公共领域许可。