devbox
从前慢现在也慢
这个屌丝很懒，什么也没留下！
关注作者
热门标签
热门文章
当前位置:   article > 正文

CGAL 点云数据生成DSM、DTM、等高线和数据分类

作者：从前慢现在也慢 | 2024-05-02 05:21:34

赞

踩

CGAL 点云数据生成DSM、DTM、等高线和数据分类

原文链接 CGAL 点云数据生成DSM、DTM、等高线和数据分类 - 知乎

在GIS应用软件中使用的许多传感器(如激光雷达)都会产生密集的点云。这类应用软件通常利用更高级的数据结构：如：不规则三角格网 (TIN)是生成数字高程模型 (DEM) 的基础，也可以利用TIN生成数字地形模型 (DTM)。对点云数据进行分类，提取地面、植被和建筑点（或其他用户定义的标签）等分类数据，从而使得获取的信息更加丰富。

因空间数据源获取方式不同，数据结构的定义也不尽一致。CGAL中使用以下术语定义这些数据结构：

  • TIN：不规则三角格网，一种 2D 三角剖分结构，根据 3D 点在水平面上的投影关系连接它们，使用众多三角形构造地物表面。
  • DSM：数字表面模型，包括建筑和植被的整个扫描表面的模型。CGAL中使用一个TIN来存储DSM。
  • DTM：数字地形模型，一个没有建筑物或植被等物体的裸地面模型。CGAL中使用TIN和光栅来存储DTM。
  • DEM：数字高程模型，是一个更通用的术语，包括DSM和DTM。

1、不规则三角网(TIN)

提供多种三角测量数据结构和算法。可以通过结合二维Delaunay三角剖分和投影特征来生成TIN。三角剖分结构是利用点在选定平面(通常是xy平面)上的二维位置来计算的，而点的三维位置则保留来进行可视化和测量。

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

  1. using Kernel = CGAL::Exact_predicates_inexact_constructions_kernel;
  2. using Projection_traits = CGAL::Projection_traits_xy_3<Kernel>;
  3. using Point_2 = Kernel::Point_2;
  4. using Point_3 = Kernel::Point_3;
  5. using Segment_3 = Kernel::Segment_3;
  6.  
  7. // Triangulated Irregular Network
  8. using TIN = CGAL::Delaunay_triangulation_2<Projection_traits>;

2、数字表面模型(DSM)

许多格式的点云(XYZ, OFF, PLY, LAS)都可以使用流操作符轻松加载到CGAL::Point_set_3结构中，然后生成存DSM储在TIN中。DSM表示如下：

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

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

  1. using Mesh = CGAL::Surface_mesh<Point_3>;
  2.  
  3.   Mesh dsm_mesh;
  4.   CGAL::copy_face_graph (dsm, dsm_mesh);
  5.   std::ofstream dsm_ofile ("dsm.ply", std::ios_base::binary);
  6.   CGAL::IO::set_binary_mode (dsm_ofile);
  7.   CGAL::IO::write_PLY (dsm_ofile, dsm_mesh);
  8.   dsm_ofile.close();

在旧金山数据集上以这种方式计算DSM的一个例子，如下图所示：

3、数字地形模型(DTM)

DSM 是 计算DTM的基础，将DSM上的非地形点进行数据清洗后，通过计算可以生成DTM，即另一个用TIN表示的DTM。作为一个例子，可以通过以下步骤计算得到一个 DTM ：

  • 设置高度阈值消除高度的剧烈变化
  • 设置周长阈值聚类周边的地物作为连接的组件
  • 过滤所有小于定义阈值的组件

该算法依赖于 2 个参数：高度阈值对应于建筑物的最小高度，以及周长阈值对应于 2D 投影上建筑物的最大尺寸。

3.1 带有信息的TIN

CGAL Delaunay 提供了一组三角测量的 API，用于计算顶点、三角面的空间相互关系，可以获取比TIN本身更多的信息。在例子中，计算每个顶点跟踪的输入点云中对应点索引，并且每个面都被赋予其连接组件的索引。

  1. auto idx_to_point_with_info
  2.     = [&](const Point_set::Index& idx) -> std::pair<Point_3, Point_set::Index>
  3.       {
  4.         return std::make_pair (points.point(idx), idx);
  5.       };
  6.  
  7.   TIN_with_info tin_with_info
  8.     (boost::make_transform_iterator (points.begin(), idx_to_point_with_info),
  9.      boost::make_transform_iterator (points.end(), idx_to_point_with_info));

3.2 识别连接组件

使用泛洪算法计算不规则三角网构成的连接组件，以种子点进行泛洪，在未超过设定的阈值时，将附近的顶点坐标信息加入到同一个连接组件TIN中。

  1. double spacing = CGAL::compute_average_spacing<Concurrency_tag>(points, 6);
  2.   spacing *= 2;
  3.  
  4.   auto face_height
  5.     = [&](const TIN_with_info::Face_handle fh) -> double
  6.       {
  7.         double out = 0.;
  8.         for (int i = 0; i < 3; ++ i)
  9.           out = (std::max) (out, CGAL::abs(fh->vertex(i)->point().z() - fh->vertex((i+1)%3)->point().z()));
  10.         return out;
  11.       };
  12.  
  13.   // Initialize faces info
  14.   for (TIN_with_info::Face_handle fh : tin_with_info.all_face_handles())
  15.     if (tin_with_info.is_infinite(fh) || face_height(fh) > spacing) // Filtered faces are given info() = -2
  16.       fh->info() = -2;
  17.     else // Pending faces are given info() = -1;
  18.       fh->info() = -1;
  19.  
  20.   // Flooding algorithm
  21.   std::vector<int> component_size;
  22.   for (TIN_with_info::Face_handle fh : tin_with_info.finite_face_handles())
  23.   {
  24.     if (fh->info() != -1)
  25.       continue;
  26.  
  27.     std::queue<TIN_with_info::Face_handle> todo;
  28.     todo.push(fh);
  29.  
  30.     int size = 0;
  31.     while (!todo.empty())
  32.     {
  33.       TIN_with_info::Face_handle current = todo.front();
  34.       todo.pop();
  35.  
  36.       if (current->info() != -1)
  37.         continue;
  38.       current->info() = int(component_size.size());
  39.       ++ size;
  40.  
  41.       for (int i = 0; i < 3; ++ i)
  42.         todo.push (current->neighbor(i));
  43.     }
  44.  
  45.     component_size.push_back (size);
  46.   }
  47.  
  48.   std::cerr << component_size.size() << " connected component(s) found" << std::endl;

这个包含了连接部件信息的TIN可以保存为彩色网格:

  1. Mesh tin_colored_mesh;
  2.  
  3.   Mesh::Property_map<Mesh::Face_index, CGAL::IO::Color>
  4.     color_map = tin_colored_mesh.add_property_map<Mesh::Face_index, CGAL::IO::Color>("f:color").first;
  5.  
  6.   CGAL::copy_face_graph (tin_with_info, tin_colored_mesh,
  7.                          CGAL::parameters::face_to_face_output_iterator
  8.                          (boost::make_function_output_iterator
  9.                           ([&](const std::pair<TIN_with_info::Face_handle, Mesh::Face_index>& ff)
  10.                            {
  11.                              // Color unassigned faces gray
  12.                              if (ff.first->info() < 0)
  13.                                color_map[ff.second] = CGAL::IO::Color(128, 128, 128);
  14.                              else
  15.                              {
  16.                                // Random color seeded by the component ID
  17.                                CGAL::Random r (ff.first->info());
  18.                                color_map[ff.second] = CGAL::IO::Color (r.get_int(64, 192),
  19.                                                                    r.get_int(64, 192),
  20.                                                                    r.get_int(64, 192));
  21.                              }
  22.                            })));
  23.  
  24.   std::ofstream tin_colored_ofile ("colored_tin.ply", std::ios_base::binary);
  25.   CGAL::IO::set_binary_mode (tin_colored_ofile);
  26.   CGAL::IO::write_PLY (tin_colored_ofile, tin_colored_mesh);
  27.   tin_colored_ofile.close();

下图给出了一个由连接组件着色的TIN示例。

3.3 数据清洗

比最大的建筑还小的组件可以这样移除:

  1. int min_size = int(points.size() / 2);
  2.  
  3.   std::vector<TIN_with_info::Vertex_handle> to_remove;
  4.   for (TIN_with_info::Vertex_handle vh : tin_with_info.finite_vertex_handles())
  5.   {
  6.     TIN_with_info::Face_circulator circ = tin_with_info.incident_faces (vh),
  7.       start = circ;
  8.  
  9.     // Remove a vertex if it's only adjacent to components smaller than threshold
  10.     bool keep = false;
  11.     do
  12.     {
  13.       if (circ->info() >= 0 && component_size[std::size_t(circ->info())] > min_size)
  14.       {
  15.         keep = true;
  16.         break;
  17.       }
  18.     }
  19.     while (++ circ != start);
  20.  
  21.     if (!keep)
  22.       to_remove.push_back (vh);
  23.   }
  24.  
  25.   std::cerr << to_remove.size() << " vertices(s) will be removed after filtering" << std::endl;
  26.   for (TIN_with_info::Vertex_handle vh : to_remove)
  27.     tin_with_info.remove (vh);

3.4 孔洞填充和网格重建

因为简单地移除被建筑物覆盖的大面积的顶点会导致生成大的Delaunay三角面片，从而造成DTM的质量较差，对于这些孔洞，使用孔填充算法填充对孔进行三角测量、细化和平整，以生成形状更好的网格数据模型。以下代码段将 TIN 复制到网格中，同时过滤掉过大的面，然后识别孔并填充所有孔，除了最大的孔（即外壳）。

  1. // Copy and keep track of overly large faces
  2.   Mesh dtm_mesh;
  3.  
  4.   std::vector<Mesh::Face_index> face_selection;
  5.   Mesh::Property_map<Mesh::Face_index, bool> face_selection_map
  6.    = dtm_mesh.add_property_map<Mesh::Face_index, bool>("is_selected", false).first;
  7.  
  8.   double limit = CGAL::square (5 * spacing);
  9.   CGAL::copy_face_graph (tin_with_info, dtm_mesh,
  10.                          CGAL::parameters::face_to_face_output_iterator
  11.                          (boost::make_function_output_iterator
  12.                           ([&](const std::pair<TIN_with_info::Face_handle, Mesh::Face_index>& ff)
  13.                            {
  14.                              double longest_edge = 0.;
  15.                              bool border = false;
  16.                              for (int i = 0; i < 3; ++ i)
  17.                              {
  18.                                longest_edge = (std::max)(longest_edge, CGAL::squared_distance
  19.                                                          (ff.first->vertex((i+1)%3)->point(),
  20.                                                           ff.first->vertex((i+2)%3)->point()));
  21.  
  22.                                TIN_with_info::Face_circulator circ
  23.                                  = tin_with_info.incident_faces (ff.first->vertex(i)),
  24.                                  start = circ;
  25.                                do
  26.                                {
  27.                                  if (tin_with_info.is_infinite (circ))
  28.                                  {
  29.                                    border = true;
  30.                                    break;
  31.                                  }
  32.                                }
  33.                                while (++ circ != start);
  34.  
  35.                                if (border)
  36.                                  break;
  37.                              }
  38.  
  39.                              // Select if face is too big AND it's not
  40.                              // on the border (to have closed holes)
  41.                              if (!border && longest_edge > limit)
  42.                              {
  43.                                face_selection_map[ff.second] = true;
  44.                                face_selection.push_back (ff.second);
  45.                              }
  46.                            })));
  47.  
  48.   // Save original DTM
  49.   std::ofstream dtm_ofile ("dtm.ply", std::ios_base::binary);
  50.   CGAL::IO::set_binary_mode (dtm_ofile);
  51.   CGAL::IO::write_PLY (dtm_ofile, dtm_mesh);
  52.   dtm_ofile.close();
  53.  
  54.   std::cerr << face_selection.size() << " face(s) are selected for removal" << std::endl;
  55.  
  56.   // Expand face selection to keep a well formed 2-manifold mesh after removal
  57.   CGAL::expand_face_selection_for_removal (face_selection, dtm_mesh, face_selection_map);
  58.   face_selection.clear();
  59.   for (Mesh::Face_index fi : faces(dtm_mesh))
  60.     if (face_selection_map[fi])
  61.       face_selection.push_back(fi);
  62.  
  63.   std::cerr << face_selection.size() << " face(s) are selected for removal after expansion" << std::endl;
  64.  
  65.   for (Mesh::Face_index fi : face_selection)
  66.     CGAL::Euler::remove_face (halfedge(fi, dtm_mesh), dtm_mesh);
  67.   dtm_mesh.collect_garbage();
  68.  
  69.   if (!dtm_mesh.is_valid())
  70.     std::cerr << "Invalid mesh!" << std::endl;
  71.  
  72.   // Save filtered DTM
  73.   std::ofstream dtm_holes_ofile ("dtm_with_holes.ply", std::ios_base::binary);
  74.   CGAL::IO::set_binary_mode (dtm_holes_ofile);
  75.   CGAL::IO::write_PLY (dtm_holes_ofile, dtm_mesh);
  76.   dtm_holes_ofile.close();
  77.  
  78.   // Get all holes
  79.   std::vector<Mesh::Halfedge_index> holes;
  80.   CGAL::Polygon_mesh_processing::extract_boundary_cycles (dtm_mesh, std::back_inserter (holes));
  81.  
  82.   std::cerr << holes.size() << " hole(s) identified" << std::endl;
  83.  
  84.   // Identify outer hull (hole with maximum size)
  85.   double max_size = 0.;
  86.   Mesh::Halfedge_index outer_hull;
  87.   for (Mesh::Halfedge_index hi : holes)
  88.   {
  89.     CGAL::Bbox_3 hole_bbox;
  90.     for (Mesh::Halfedge_index haf : CGAL::halfedges_around_face(hi, dtm_mesh))
  91.     {
  92.       const Point_3& p = dtm_mesh.point(target(haf, dtm_mesh));
  93.       hole_bbox += p.bbox();
  94.     }
  95.     double size = CGAL::squared_distance (Point_2(hole_bbox.xmin(), hole_bbox.ymin()),
  96.                                           Point_2(hole_bbox.xmax(), hole_bbox.ymax()));
  97.     if (size > max_size)
  98.     {
  99.       max_size = size;
  100.       outer_hull = hi;
  101.     }
  102.   }
  103.  
  104.   // Fill all holes except the bigest (which is the outer hull of the mesh)
  105.   for (Mesh::Halfedge_index hi : holes)
  106.     if (hi != outer_hull)
  107.       CGAL::Polygon_mesh_processing::triangulate_refine_and_fair_hole
  108.         (dtm_mesh, hi, CGAL::Emptyset_iterator(), CGAL::Emptyset_iterator());
  109.  
  110.   // Save DTM with holes filled
  111.   std::ofstream dtm_filled_ofile ("dtm_filled.ply", std::ios_base::binary);
  112.   CGAL::IO::set_binary_mode (dtm_filled_ofile);
  113.   CGAL::IO::write_PLY (dtm_filled_ofile, dtm_mesh);
  114.   dtm_filled_ofile.close();

为了产生不受二维Delaunay面形状约束的更规则的网格，各向同性网格也可以作为最后一步进行。

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

4 光栅化

TIN数据结构可以与重心坐标相结合，以便插值并栅格化任何分辨率的高度图，需要嵌入顶点的信息。由于最新的两个步骤(孔填充和网格划分)是在3D网格上进行的，因此DTM是2.5D表示的假设可能不再有效。因此，我们首先使用最后计算的各向同性DTM网格的顶点重建TIN。以下使用代码段使用简单位图格式（PPM）生成栅格DEM。

  1. CGAL::Bbox_3 bbox = CGAL::bbox_3 (points.points().begin(), points.points().end());
  2.  
  3.   // Generate raster image 1920-pixels large
  4.   std::size_t width = 1920;
  5.   std::size_t height = std::size_t((bbox.ymax() - bbox.ymin()) * 1920 / (bbox.xmax() - bbox.xmin()));
  6.  
  7.   std::cerr << "Rastering with resolution " << width << "x" << height << std::endl;
  8.  
  9.   // Use PPM format (Portable PixMap) for simplicity
  10.   std::ofstream raster_ofile ("raster.ppm", std::ios_base::binary);
  11.  
  12.   // PPM header
  13.   raster_ofile << "P6" << std::endl // magic number
  14.                << width << " " << height << std::endl // dimensions of the image
  15.                << 255 << std::endl; // maximum color value
  16.  
  17.   // Use rainbow color ramp output
  18.   Color_ramp color_ramp;
  19.  
  20.   // Keeping track of location from one point to its neighbor allows
  21.   // for fast locate in DT
  22.   TIN::Face_handle location;
  23.  
  24.   // Query each pixel of the image
  25.   for (std::size_t y = 0; y < height; ++ y)
  26.     for (std::size_t x = 0; x < width; ++ x)
  27.     {
  28.       Point_3 query (bbox.xmin() + x * (bbox.xmax() - bbox.xmin()) / double(width),
  29.                      bbox.ymin() + (height-y) * (bbox.ymax() - bbox.ymin()) / double(height),
  30.                      0); // not relevant for location in 2D
  31.  
  32.       location = dtm_clean.locate (query, location);
  33.  
  34.       // Points outside the convex hull will be colored black
  35.       std::array<unsigned char, 3> colors { 0, 0, 0 };
  36.       if (!dtm_clean.is_infinite(location))
  37.       {
  38.         std::array<double, 3> barycentric_coordinates
  39.           = CGAL::Polygon_mesh_processing::barycentric_coordinates
  40.           (Point_2 (location->vertex(0)->point().x(), location->vertex(0)->point().y()),
  41.            Point_2 (location->vertex(1)->point().x(), location->vertex(1)->point().y()),
  42.            Point_2 (location->vertex(2)->point().x(), location->vertex(2)->point().y()),
  43.            Point_2 (query.x(), query.y()),
  44.            Kernel());
  45.  
  46.         double height_at_query
  47.           = (barycentric_coordinates[0] * location->vertex(0)->point().z()
  48.              + barycentric_coordinates[1] * location->vertex(1)->point().z()
  49.              + barycentric_coordinates[2] * location->vertex(2)->point().z());
  50.  
  51.         // Color ramp generates a color depending on a value from 0 to 1
  52.         double height_ratio = (height_at_query - bbox.zmin()) / (bbox.zmax() - bbox.zmin());
  53.         colors = color_ramp.get(height_ratio);
  54.       }
  55.       raster_ofile.write (reinterpret_cast<char*>(&colors), 3);
  56.     }
  57.  
  58.   raster_ofile.close();

下图给出了一个用彩虹斜坡表示高度的光栅图像示例。

5 等高线生成

提取TIN上定义的函数的等值是可以用CGAL完成的另一个应用程序。我们在这里演示如何提取等高线来构建地形图。

5.1绘制等高线图

第一步是，从三角剖分的所有面中，以分段的形式提取每个等值面经过该面的截面。下面的函数允许测试一个等值值是否穿过一个面，然后提取它:

  1. bool face_has_isovalue (TIN::Face_handle fh, double isovalue)
  2. {
  3.   bool above = false, below = false;
  4.   for (int i = 0; i < 3; ++ i)
  5.   {
  6.     // Face has isovalue if one of its vertices is above and another
  7.     // one below
  8.     if (fh->vertex(i)->point().z() > isovalue)
  9.       above = true;
  10.     if (fh->vertex(i)->point().z() < isovalue)
  11.       below = true;
  12.   }
  13.  
  14.   return (above && below);
  15. }
  16.  
  17. Segment_3 isocontour_in_face (TIN::Face_handle fh, double isovalue)
  18. {
  19.   Point_3 source;
  20.   Point_3 target;
  21.   bool source_found = false;
  22.  
  23.   for (int i = 0; i < 3; ++ i)
  24.   {
  25.     Point_3 p0 = fh->vertex((i+1) % 3)->point();
  26.     Point_3 p1 = fh->vertex((i+2) % 3)->point();
  27.  
  28.     // Check if the isovalue crosses segment (p0,p1)
  29.     if ((p0.z() - isovalue) * (p1.z() - isovalue) > 0)
  30.       continue;
  31.  
  32.     double zbottom = p0.z();
  33.     double ztop = p1.z();
  34.     if (zbottom > ztop)
  35.     {
  36.       std::swap (zbottom, ztop);
  37.       std::swap (p0, p1);
  38.     }
  39.  
  40.     // Compute position of segment vertex
  41.     double ratio = (isovalue - zbottom) / (ztop - zbottom);
  42.     Point_3 p = CGAL::barycenter (p0, (1 - ratio), p1,ratio);
  43.  
  44.     if (source_found)
  45.       target = p;
  46.     else
  47.     {
  48.       source = p;
  49.       source_found = true;
  50.     }
  51.   }
  52.  
  53.   return Segment_3 (source, target);
  54. }

 通过这些函数，我们可以创建一个线段图，然后将其处理成一组折线。为此，我们使用`boost::adjacency_list`结构，并跟踪从端点到顶点的映射。下以下代码计算最小和最大高度差的50 个等值，并创建一个包含所有等值的图形：

  1. std::array<double, 50> isovalues; // Contour 50 isovalues
  2.   for (std::size_t i = 0; i < isovalues.size(); ++ i)
  3.     isovalues[i] = bbox.zmin() + ((i+1) * (bbox.zmax() - bbox.zmin()) / (isovalues.size() - 2));
  4.  
  5.   // First find on each face if they are crossed by some isovalues and
  6.   // extract segments in a graph
  7.   using Segment_graph = boost::adjacency_list<boost::listS, boost::vecS, boost::undirectedS, Point_3>;
  8.   Segment_graph graph;
  9.   using Map_p2v = std::map<Point_3, Segment_graph::vertex_descriptor>;
  10.   Map_p2v map_p2v;
  11.   for (TIN::Face_handle vh : dtm_clean.finite_face_handles())
  12.     for (double iv : isovalues)
  13.       if (face_has_isovalue (vh, iv))
  14.       {
  15.         Segment_3 segment = isocontour_in_face (vh, iv);
  16.         for (const Point_3& p : { segment.source(), segment.target() })
  17.         {
  18.           // Only insert end points of segments once to get a well connected graph
  19.           Map_p2v::iterator iter;
  20.           bool inserted;
  21.           std::tie (iter, inserted) = map_p2v.insert (std::make_pair (p, Segment_graph::vertex_descriptor()));
  22.           if (inserted)
  23.           {
  24.             iter->second = boost::add_vertex (graph);
  25.             graph[iter->second] = p;
  26.           }
  27.         }
  28.         boost::add_edge (map_p2v[segment.source()], map_p2v[segment.target()], graph);
  29.       }

5.2分割成折线

创建图形后，使用函数`CGAL::split_graph_into_polylines()`就可以很容易地将其分解为折线:

  1. // Split segments into polylines
  2.   std::vector<std::vector<Point_3> > polylines;
  3.   Polylines_visitor<Segment_graph> visitor (graph, polylines);
  4.   CGAL::split_graph_into_polylines (graph, visitor);
  5.  
  6.   std::cerr << polylines.size() << " polylines computed, with "
  7.             << map_p2v.size() << " vertices in total" << std::endl;
  8.  
  9.   // Output to WKT file
  10.   std::ofstream contour_ofile ("contour.wkt");
  11.   contour_ofile.precision(18);
  12.   CGAL::IO::write_multi_linestring_WKT (contour_ofile, polylines);
  13.   contour_ofile.close();

 该函数需要提供一个可供访问的对象，以便在开始折线、向其添加点以及结束它时调用。可以定义一个简单的类使用它。

  1. template <typename Graph>
  2. class Polylines_visitor
  3. {
  4. private:
  5.   std::vector<std::vector<Point_3> >& polylines;
  6.   Graph& graph;
  7.  
  8. public:
  9.  
  10.   Polylines_visitor (Graph& graph, std::vector<std::vector<Point_3> >& polylines)
  11.     : polylines (polylines), graph(graph) { }
  12.  
  13.   void start_new_polyline()
  14.   {
  15.     polylines.push_back (std::vector<Point_3>());
  16.   }
  17.  
  18.   void add_node (typename Graph::vertex_descriptor vd)
  19.   {
  20.     polylines.back().push_back (graph[vd]);
  21.   }
  22.  
  23.   void end_polyline()
  24.   {
  25.     // filter small polylines
  26.     if (polylines.back().size() < 50)
  27.       polylines.pop_back();
  28.   }
  29. };

 5.3 等高线简化

由于输出的噪声很大，用户可能希望简化折线。CGAL提供了一种折线简化算法，该算法保证简化后两条折线不相交。该算法利用了`CGAL::Constrained_triangulation_plus_2`，它将折线嵌入为一组约束:

  1. namespace PS = CGAL::Polyline_simplification_2;
  2. using CDT_vertex_base = PS::Vertex_base_2<Projection_traits>;
  3. using CDT_face_base = CGAL::Constrained_triangulation_face_base_2<Projection_traits>;
  4. using CDT_TDS = CGAL::Triangulation_data_structure_2<CDT_vertex_base, CDT_face_base>;
  5. using CDT = CGAL::Constrained_Delaunay_triangulation_2<Projection_traits, CDT_TDS>;
  6. using CTP = CGAL::Constrained_triangulation_plus_2<CDT>;

下面的代码根据折线到原始折线的平方距离简化了折线集，当没有更多的顶点可以移除时停止，而不超过平均间距的4倍。

  1. // Construct constrained Delaunay triangulation with polylines as constraints
  2.   CTP ctp;
  3.   for (const std::vector<Point_3>& poly : polylines)
  4.     ctp.insert_constraint (poly.begin(), poly.end());
  5.  
  6.   // Simplification algorithm with limit on distance
  7.   PS::simplify (ctp, PS::Squared_distance_cost(), PS::Stop_above_cost_threshold (16 * spacing * spacing));
  8.  
  9.   polylines.clear();
  10.   for (CTP::Constraint_id cid : ctp.constraints())
  11.   {
  12.     polylines.push_back (std::vector<Point_3>());
  13.     polylines.back().reserve (ctp.vertices_in_constraint (cid).size());
  14.     for (CTP::Vertex_handle vh : ctp.vertices_in_constraint(cid))
  15.       polylines.back().push_back (vh->point());
  16.   }
  17.  
  18.   std::size_t nb_vertices
  19.     = std::accumulate (polylines.begin(), polylines.end(), std::size_t(0),
  20.                        [](std::size_t size, const std::vector<Point_3>& poly) -> std::size_t
  21.                        { return size + poly.size(); });
  22.  
  23.   std::cerr << nb_vertices
  24.             << " vertices remaining after simplification ("
  25.             << 100. * (nb_vertices / double(map_p2v.size())) << "%)" << std::endl;
  26.  
  27.   // Output to WKT file
  28.   std::ofstream simplified_ofile ("simplified.wkt");
  29.   simplified_ofile.precision(18);
  30.   CGAL::IO::write_multi_linestring_WKT (simplified_ofile, polylines);
  31.   simplified_ofile.close();

 图中给出了轮廓和简化的例子。等高线使用50等值均匀间隔。顶部:使用148k个顶点的原始轮廓和简化容差等于输入点云的平均间距(原始折线顶点的3.4%剩余)。下图:简化，公差为平均间距的4倍(剩余顶点的1.3%)，公差为平均间距的10倍(剩余顶点的0.9%)。折线在所有情况下都是无交叉的。

6、点云分类

提供了一个包分类，可用于将点云分割为用户定义的标签集。目前CGAL中最先进的分类器是随机森林。由于它是一个监督分类器，因此需要一个训练集。下面的代码片段展示了如何使用一些手动选择的训练集来训练随机森林分类器，并计算由图切割算法正则化的分类:

  1. // Get training from input
  2.   Point_set::Property_map<int> training_map;
  3.   bool training_found;
  4.   std::tie (training_map, training_found) = points.property_map<int>("training");
  5.  
  6.   if (training_found)
  7.   {
  8.     std::cerr << "Classifying ground/vegetation/building" << std::endl;
  9.  
  10.     // Create labels
  11.     Classification::Label_set labels ({ "ground", "vegetation", "building" });
  12.  
  13.     // Generate features
  14.     Classification::Feature_set features;
  15.     Classification::Point_set_feature_generator<Kernel, Point_set, Point_set::Point_map>
  16.       generator (points, points.point_map(), 5); // 5 scales
  17.  
  18. #ifdef CGAL_LINKED_WITH_TBB
  19.     // If TBB is used, features can be computed in parallel
  20.     features.begin_parallel_additions();
  21.     generator.generate_point_based_features (features);
  22.     features.end_parallel_additions();
  23. #else
  24.     generator.generate_point_based_features (features);
  25. #endif
  26.  
  27.     // Train a random forest classifier
  28.     Classification::ETHZ::Random_forest_classifier classifier (labels, features);
  29.     classifier.train (points.range(training_map));
  30.  
  31.     // Classify with graphcut regularization
  32.     Point_set::Property_map<int> label_map = points.add_property_map<int>("labels").first;
  33.     Classification::classify_with_graphcut<Concurrency_tag>
  34.       (points, points.point_map(), labels, classifier,
  35.        generator.neighborhood().k_neighbor_query(12), // regularize on 12-neighbors graph
  36.        0.5f, // graphcut weight
  37.        12, // Subdivide to speed-up process
  38.        label_map);
  39.  
  40.     // Evaluate
  41.     std::cerr << "Mean IoU on training data = "
  42.               << Classification::Evaluation(labels,
  43.                                             points.range(training_map),
  44.                                             points.range(label_map)).mean_intersection_over_union() << std::endl;
  45.  
  46.     // Save the classified point set
  47.     std::ofstream classified_ofile ("classified.ply");
  48.     CGAL::IO::set_binary_mode (classified_ofile);
  49.     classified_ofile << points;
  50.     classified_ofile.close();
  51.   }

 下图给出了训练集和分类结果的示例。顶部:用手分类的点云切片。下图:在3个人工分类切片上训练后，通过图切割规范化的分类。

7、完整代码示例

本教程中使用的所有代码片段都可以组装起来创建一个完整的GIS管道(如果使用了正确的include)。我们给出了一个完整的代码示例，它实现了本教程中描述的所有步骤。

  1. #include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
  2. #include <CGAL/Projection_traits_xy_3.h>
  3.  
  4. #include <CGAL/Delaunay_triangulation_2.h>
  5. #include <CGAL/Triangulation_vertex_base_with_info_2.h>
  6. #include <CGAL/Triangulation_face_base_with_info_2.h>
  7.  
  8. #include <CGAL/boost/graph/graph_traits_Delaunay_triangulation_2.h>
  9. #include <CGAL/boost/graph/copy_face_graph.h>
  10.  
  11. #include <CGAL/Point_set_3.h>
  12. #include <CGAL/Point_set_3/IO.h>
  13. #include <CGAL/compute_average_spacing.h>
  14.  
  15. #include <CGAL/Surface_mesh.h>
  16. #include <CGAL/Polygon_mesh_processing/locate.h>
  17.  
  18. #include <CGAL/Polygon_mesh_processing/triangulate_hole.h>
  19. #include <CGAL/Polygon_mesh_processing/border.h>
  20. #include <CGAL/Polygon_mesh_processing/remesh.h>
  21.  
  22. #include <boost/graph/adjacency_list.hpp>
  23. #include <CGAL/boost/graph/split_graph_into_polylines.h>
  24.  
  25. #include <CGAL/IO/WKT.h>
  26.  
  27. #include <CGAL/Constrained_Delaunay_triangulation_2.h>
  28. #include <CGAL/Constrained_triangulation_plus_2.h>
  29.  
  30. #include <CGAL/Polyline_simplification_2/simplify.h>
  31. #include <CGAL/Polyline_simplification_2/Squared_distance_cost.h>
  32.  
  33. #include <CGAL/Classification.h>
  34.  
  35. #include <CGAL/Random.h>
  36.  
  37. #include <fstream>
  38. #include <queue>
  39.  
  40. #include "include/Color_ramp.h"
  41.  
  42.  
  43. using Kernel = CGAL::Exact_predicates_inexact_constructions_kernel;
  44. using Projection_traits = CGAL::Projection_traits_xy_3<Kernel>;
  45. using Point_2 = Kernel::Point_2;
  46. using Point_3 = Kernel::Point_3;
  47. using Segment_3 = Kernel::Segment_3;
  48.  
  49. // Triangulated Irregular Network
  50. using TIN = CGAL::Delaunay_triangulation_2<Projection_traits>;
  51.  
  52.  
  53.  
  54. // Triangulated Irregular Network (with info)
  55. using Point_set = CGAL::Point_set_3<Point_3>;
  56. using Vbi = CGAL::Triangulation_vertex_base_with_info_2 <Point_set::Index, Projection_traits>;
  57. using Fbi = CGAL::Triangulation_face_base_with_info_2<int, Projection_traits>;
  58. using TDS = CGAL::Triangulation_data_structure_2<Vbi, Fbi>;
  59. using TIN_with_info = CGAL::Delaunay_triangulation_2<Projection_traits, TDS>;
  60.  
  61.  
  62. namespace Classification = CGAL::Classification;
  63.  
  64. #ifdef CGAL_LINKED_WITH_TBB
  65. using Concurrency_tag = CGAL::Parallel_tag;
  66. #else
  67. using Concurrency_tag = CGAL::Sequential_tag;
  68. #endif
  69.  
  70.  
  71. bool face_has_isovalue (TIN::Face_handle fh, double isovalue)
  72. {
  73.   bool above = false, below = false;
  74.   for (int i = 0; i < 3; ++ i)
  75.   {
  76.     // Face has isovalue if one of its vertices is above and another
  77.     // one below
  78.     if (fh->vertex(i)->point().z() > isovalue)
  79.       above = true;
  80.     if (fh->vertex(i)->point().z() < isovalue)
  81.       below = true;
  82.   }
  83.  
  84.   return (above && below);
  85. }
  86.  
  87. Segment_3 isocontour_in_face (TIN::Face_handle fh, double isovalue)
  88. {
  89.   Point_3 source;
  90.   Point_3 target;
  91.   bool source_found = false;
  92.  
  93.   for (int i = 0; i < 3; ++ i)
  94.   {
  95.     Point_3 p0 = fh->vertex((i+1) % 3)->point();
  96.     Point_3 p1 = fh->vertex((i+2) % 3)->point();
  97.  
  98.     // Check if the isovalue crosses segment (p0,p1)
  99.     if ((p0.z() - isovalue) * (p1.z() - isovalue) > 0)
  100.       continue;
  101.  
  102.     double zbottom = p0.z();
  103.     double ztop = p1.z();
  104.     if (zbottom > ztop)
  105.     {
  106.       std::swap (zbottom, ztop);
  107.       std::swap (p0, p1);
  108.     }
  109.  
  110.     // Compute position of segment vertex
  111.     double ratio = (isovalue - zbottom) / (ztop - zbottom);
  112.     Point_3 p = CGAL::barycenter (p0, (1 - ratio), p1,ratio);
  113.  
  114.     if (source_found)
  115.       target = p;
  116.     else
  117.     {
  118.       source = p;
  119.       source_found = true;
  120.     }
  121.   }
  122.  
  123.   return Segment_3 (source, target);
  124. }
  125.  
  126.  
  127.  
  128.  
  129. template <typename Graph>
  130. class Polylines_visitor
  131. {
  132. private:
  133.   std::vector<std::vector<Point_3> >& polylines;
  134.   Graph& graph;
  135.  
  136. public:
  137.  
  138.   Polylines_visitor (Graph& graph, std::vector<std::vector<Point_3> >& polylines)
  139.     : polylines (polylines), graph(graph) { }
  140.  
  141.   void start_new_polyline()
  142.   {
  143.     polylines.push_back (std::vector<Point_3>());
  144.   }
  145.  
  146.   void add_node (typename Graph::vertex_descriptor vd)
  147.   {
  148.     polylines.back().push_back (graph[vd]);
  149.   }
  150.  
  151.   void end_polyline()
  152.   {
  153.     // filter small polylines
  154.     if (polylines.back().size() < 50)
  155.       polylines.pop_back();
  156.   }
  157. };
  158.  
  159.  
  160.  
  161. namespace PS = CGAL::Polyline_simplification_2;
  162. using CDT_vertex_base = PS::Vertex_base_2<Projection_traits>;
  163. using CDT_face_base = CGAL::Constrained_triangulation_face_base_2<Projection_traits>;
  164. using CDT_TDS = CGAL::Triangulation_data_structure_2<CDT_vertex_base, CDT_face_base>;
  165. using CDT = CGAL::Constrained_Delaunay_triangulation_2<Projection_traits, CDT_TDS>;
  166. using CTP = CGAL::Constrained_triangulation_plus_2<CDT>;
  167.  
  168.  
  169. int main (int argc, char** argv)
  170. {
  171.   const std::string fname = argc != 2 ? CGAL::data_file_path("points_3/b9_training.ply") : argv[1];
  172.   if (argc != 2)
  173.   {
  174.     std::cerr << "Usage: " << argv[0] << " points.ply" << std::endl;
  175.     std::cerr << "Running with default value " << fname << "\n";
  176.   }
  177.  
  178.  
  179.   // Read points
  180.   std::ifstream ifile (fname, std::ios_base::binary);
  181.   CGAL::Point_set_3<Point_3> points;
  182.   ifile >> points;
  183.   std::cerr << points.size() << " point(s) read" << std::endl;
  184.  
  185.   // Create DSM
  186.   TIN dsm (points.points().begin(), points.points().end());
  187.  
  188.  
  189.  
  190.   using Mesh = CGAL::Surface_mesh<Point_3>;
  191.  
  192.   Mesh dsm_mesh;
  193.   CGAL::copy_face_graph (dsm, dsm_mesh);
  194.   std::ofstream dsm_ofile ("dsm.ply", std::ios_base::binary);
  195.   CGAL::IO::set_binary_mode (dsm_ofile);
  196.   CGAL::IO::write_PLY (dsm_ofile, dsm_mesh);
  197.   dsm_ofile.close();
  198.  
  199.  
  200.  
  201.   auto idx_to_point_with_info
  202.     = [&](const Point_set::Index& idx) -> std::pair<Point_3, Point_set::Index>
  203.       {
  204.         return std::make_pair (points.point(idx), idx);
  205.       };
  206.  
  207.   TIN_with_info tin_with_info
  208.     (boost::make_transform_iterator (points.begin(), idx_to_point_with_info),
  209.      boost::make_transform_iterator (points.end(), idx_to_point_with_info));
  210.  
  211.  
  212.  
  213.   double spacing = CGAL::compute_average_spacing<Concurrency_tag>(points, 6);
  214.   spacing *= 2;
  215.  
  216.   auto face_height
  217.     = [&](const TIN_with_info::Face_handle fh) -> double
  218.       {
  219.         double out = 0.;
  220.         for (int i = 0; i < 3; ++ i)
  221.           out = (std::max) (out, CGAL::abs(fh->vertex(i)->point().z() - fh->vertex((i+1)%3)->point().z()));
  222.         return out;
  223.       };
  224.  
  225.   // Initialize faces info
  226.   for (TIN_with_info::Face_handle fh : tin_with_info.all_face_handles())
  227.     if (tin_with_info.is_infinite(fh) || face_height(fh) > spacing) // Filtered faces are given info() = -2
  228.       fh->info() = -2;
  229.     else // Pending faces are given info() = -1;
  230.       fh->info() = -1;
  231.  
  232.   // Flooding algorithm
  233.   std::vector<int> component_size;
  234.   for (TIN_with_info::Face_handle fh : tin_with_info.finite_face_handles())
  235.   {
  236.     if (fh->info() != -1)
  237.       continue;
  238.  
  239.     std::queue<TIN_with_info::Face_handle> todo;
  240.     todo.push(fh);
  241.  
  242.     int size = 0;
  243.     while (!todo.empty())
  244.     {
  245.       TIN_with_info::Face_handle current = todo.front();
  246.       todo.pop();
  247.  
  248.       if (current->info() != -1)
  249.         continue;
  250.       current->info() = int(component_size.size());
  251.       ++ size;
  252.  
  253.       for (int i = 0; i < 3; ++ i)
  254.         todo.push (current->neighbor(i));
  255.     }
  256.  
  257.     component_size.push_back (size);
  258.   }
  259.  
  260.   std::cerr << component_size.size() << " connected component(s) found" << std::endl;
  261.  
  262.  
  263.  
  264.   Mesh tin_colored_mesh;
  265.  
  266.   Mesh::Property_map<Mesh::Face_index, CGAL::IO::Color>
  267.     color_map = tin_colored_mesh.add_property_map<Mesh::Face_index, CGAL::IO::Color>("f:color").first;
  268.  
  269.   CGAL::copy_face_graph (tin_with_info, tin_colored_mesh,
  270.                          CGAL::parameters::face_to_face_output_iterator
  271.                          (boost::make_function_output_iterator
  272.                           ([&](const std::pair<TIN_with_info::Face_handle, Mesh::Face_index>& ff)
  273.                            {
  274.                              // Color unassigned faces gray
  275.                              if (ff.first->info() < 0)
  276.                                color_map[ff.second] = CGAL::IO::Color(128, 128, 128);
  277.                              else
  278.                              {
  279.                                // Random color seeded by the component ID
  280.                                CGAL::Random r (ff.first->info());
  281.                                color_map[ff.second] = CGAL::IO::Color (r.get_int(64, 192),
  282.                                                                    r.get_int(64, 192),
  283.                                                                    r.get_int(64, 192));
  284.                              }
  285.                            })));
  286.  
  287.   std::ofstream tin_colored_ofile ("colored_tin.ply", std::ios_base::binary);
  288.   CGAL::IO::set_binary_mode (tin_colored_ofile);
  289.   CGAL::IO::write_PLY (tin_colored_ofile, tin_colored_mesh);
  290.   tin_colored_ofile.close();
  291.  
  292.  
  293.  
  294.   int min_size = int(points.size() / 2);
  295.  
  296.   std::vector<TIN_with_info::Vertex_handle> to_remove;
  297.   for (TIN_with_info::Vertex_handle vh : tin_with_info.finite_vertex_handles())
  298.   {
  299.     TIN_with_info::Face_circulator circ = tin_with_info.incident_faces (vh),
  300.       start = circ;
  301.  
  302.     // Remove a vertex if it's only adjacent to components smaller than threshold
  303.     bool keep = false;
  304.     do
  305.     {
  306.       if (circ->info() >= 0 && component_size[std::size_t(circ->info())] > min_size)
  307.       {
  308.         keep = true;
  309.         break;
  310.       }
  311.     }
  312.     while (++ circ != start);
  313.  
  314.     if (!keep)
  315.       to_remove.push_back (vh);
  316.   }
  317.  
  318.   std::cerr << to_remove.size() << " vertices(s) will be removed after filtering" << std::endl;
  319.   for (TIN_with_info::Vertex_handle vh : to_remove)
  320.     tin_with_info.remove (vh);
  321.  
  322.  
  323.  
  324.   // Copy and keep track of overly large faces
  325.   Mesh dtm_mesh;
  326.  
  327.   std::vector<Mesh::Face_index> face_selection;
  328.   Mesh::Property_map<Mesh::Face_index, bool> face_selection_map
  329.    = dtm_mesh.add_property_map<Mesh::Face_index, bool>("is_selected", false).first;
  330.  
  331.   double limit = CGAL::square (5 * spacing);
  332.   CGAL::copy_face_graph (tin_with_info, dtm_mesh,
  333.                          CGAL::parameters::face_to_face_output_iterator
  334.                          (boost::make_function_output_iterator
  335.                           ([&](const std::pair<TIN_with_info::Face_handle, Mesh::Face_index>& ff)
  336.                            {
  337.                              double longest_edge = 0.;
  338.                              bool border = false;
  339.                              for (int i = 0; i < 3; ++ i)
  340.                              {
  341.                                longest_edge = (std::max)(longest_edge, CGAL::squared_distance
  342.                                                          (ff.first->vertex((i+1)%3)->point(),
  343.                                                           ff.first->vertex((i+2)%3)->point()));
  344.  
  345.                                TIN_with_info::Face_circulator circ
  346.                                  = tin_with_info.incident_faces (ff.first->vertex(i)),
  347.                                  start = circ;
  348.                                do
  349.                                {
  350.                                  if (tin_with_info.is_infinite (circ))
  351.                                  {
  352.                                    border = true;
  353.                                    break;
  354.                                  }
  355.                                }
  356.                                while (++ circ != start);
  357.  
  358.                                if (border)
  359.                                  break;
  360.                              }
  361.  
  362.                              // Select if face is too big AND it's not
  363.                              // on the border (to have closed holes)
  364.                              if (!border && longest_edge > limit)
  365.                              {
  366.                                face_selection_map[ff.second] = true;
  367.                                face_selection.push_back (ff.second);
  368.                              }
  369.                            })));
  370.  
  371.   // Save original DTM
  372.   std::ofstream dtm_ofile ("dtm.ply", std::ios_base::binary);
  373.   CGAL::IO::set_binary_mode (dtm_ofile);
  374.   CGAL::IO::write_PLY (dtm_ofile, dtm_mesh);
  375.   dtm_ofile.close();
  376.  
  377.   std::cerr << face_selection.size() << " face(s) are selected for removal" << std::endl;
  378.  
  379.   // Expand face selection to keep a well formed 2-manifold mesh after removal
  380.   CGAL::expand_face_selection_for_removal (face_selection, dtm_mesh, face_selection_map);
  381.   face_selection.clear();
  382.   for (Mesh::Face_index fi : faces(dtm_mesh))
  383.     if (face_selection_map[fi])
  384.       face_selection.push_back(fi);
  385.  
  386.   std::cerr << face_selection.size() << " face(s) are selected for removal after expansion" << std::endl;
  387.  
  388.   for (Mesh::Face_index fi : face_selection)
  389.     CGAL::Euler::remove_face (halfedge(fi, dtm_mesh), dtm_mesh);
  390.   dtm_mesh.collect_garbage();
  391.  
  392.   if (!dtm_mesh.is_valid())
  393.     std::cerr << "Invalid mesh!" << std::endl;
  394.  
  395.   // Save filtered DTM
  396.   std::ofstream dtm_holes_ofile ("dtm_with_holes.ply", std::ios_base::binary);
  397.   CGAL::IO::set_binary_mode (dtm_holes_ofile);
  398.   CGAL::IO::write_PLY (dtm_holes_ofile, dtm_mesh);
  399.   dtm_holes_ofile.close();
  400.  
  401.   // Get all holes
  402.   std::vector<Mesh::Halfedge_index> holes;
  403.   CGAL::Polygon_mesh_processing::extract_boundary_cycles (dtm_mesh, std::back_inserter (holes));
  404.  
  405.   std::cerr << holes.size() << " hole(s) identified" << std::endl;
  406.  
  407.   // Identify outer hull (hole with maximum size)
  408.   double max_size = 0.;
  409.   Mesh::Halfedge_index outer_hull;
  410.   for (Mesh::Halfedge_index hi : holes)
  411.   {
  412.     CGAL::Bbox_3 hole_bbox;
  413.     for (Mesh::Halfedge_index haf : CGAL::halfedges_around_face(hi, dtm_mesh))
  414.     {
  415.       const Point_3& p = dtm_mesh.point(target(haf, dtm_mesh));
  416.       hole_bbox += p.bbox();
  417.     }
  418.     double size = CGAL::squared_distance (Point_2(hole_bbox.xmin(), hole_bbox.ymin()),
  419.                                           Point_2(hole_bbox.xmax(), hole_bbox.ymax()));
  420.     if (size > max_size)
  421.     {
  422.       max_size = size;
  423.       outer_hull = hi;
  424.     }
  425.   }
  426.  
  427.   // Fill all holes except the bigest (which is the outer hull of the mesh)
  428.   for (Mesh::Halfedge_index hi : holes)
  429.     if (hi != outer_hull)
  430.       CGAL::Polygon_mesh_processing::triangulate_refine_and_fair_hole
  431.         (dtm_mesh, hi, CGAL::Emptyset_iterator(), CGAL::Emptyset_iterator());
  432.  
  433.   // Save DTM with holes filled
  434.   std::ofstream dtm_filled_ofile ("dtm_filled.ply", std::ios_base::binary);
  435.   CGAL::IO::set_binary_mode (dtm_filled_ofile);
  436.   CGAL::IO::write_PLY (dtm_filled_ofile, dtm_mesh);
  437.   dtm_filled_ofile.close();
  438.  
  439.  
  440.  
  441.   CGAL::Polygon_mesh_processing::isotropic_remeshing (faces(dtm_mesh), spacing, dtm_mesh);
  442.  
  443.   std::ofstream dtm_remeshed_ofile ("dtm_remeshed.ply", std::ios_base::binary);
  444.   CGAL::IO::set_binary_mode (dtm_remeshed_ofile);
  445.   CGAL::IO::write_PLY (dtm_remeshed_ofile, dtm_mesh);
  446.   dtm_remeshed_ofile.close();
  447.  
  448.  
  449.   TIN dtm_clean (dtm_mesh.points().begin(), dtm_mesh.points().end());
  450.  
  451.  
  452.   CGAL::Bbox_3 bbox = CGAL::bbox_3 (points.points().begin(), points.points().end());
  453.  
  454.   // Generate raster image 1920-pixels large
  455.   std::size_t width = 1920;
  456.   std::size_t height = std::size_t((bbox.ymax() - bbox.ymin()) * 1920 / (bbox.xmax() - bbox.xmin()));
  457.  
  458.   std::cerr << "Rastering with resolution " << width << "x" << height << std::endl;
  459.  
  460.   // Use PPM format (Portable PixMap) for simplicity
  461.   std::ofstream raster_ofile ("raster.ppm", std::ios_base::binary);
  462.  
  463.   // PPM header
  464.   raster_ofile << "P6" << std::endl // magic number
  465.                << width << " " << height << std::endl // dimensions of the image
  466.                << 255 << std::endl; // maximum color value
  467.  
  468.   // Use rainbow color ramp output
  469.   Color_ramp color_ramp;
  470.  
  471.   // Keeping track of location from one point to its neighbor allows
  472.   // for fast locate in DT
  473.   TIN::Face_handle location;
  474.  
  475.   // Query each pixel of the image
  476.   for (std::size_t y = 0; y < height; ++ y)
  477.     for (std::size_t x = 0; x < width; ++ x)
  478.     {
  479.       Point_3 query (bbox.xmin() + x * (bbox.xmax() - bbox.xmin()) / double(width),
  480.                      bbox.ymin() + (height-y) * (bbox.ymax() - bbox.ymin()) / double(height),
  481.                      0); // not relevant for location in 2D
  482.  
  483.       location = dtm_clean.locate (query, location);
  484.  
  485.       // Points outside the convex hull will be colored black
  486.       std::array<unsigned char, 3> colors { 0, 0, 0 };
  487.       if (!dtm_clean.is_infinite(location))
  488.       {
  489.         std::array<double, 3> barycentric_coordinates
  490.           = CGAL::Polygon_mesh_processing::barycentric_coordinates
  491.           (Point_2 (location->vertex(0)->point().x(), location->vertex(0)->point().y()),
  492.            Point_2 (location->vertex(1)->point().x(), location->vertex(1)->point().y()),
  493.            Point_2 (location->vertex(2)->point().x(), location->vertex(2)->point().y()),
  494.            Point_2 (query.x(), query.y()),
  495.            Kernel());
  496.  
  497.         double height_at_query
  498.           = (barycentric_coordinates[0] * location->vertex(0)->point().z()
  499.              + barycentric_coordinates[1] * location->vertex(1)->point().z()
  500.              + barycentric_coordinates[2] * location->vertex(2)->point().z());
  501.  
  502.         // Color ramp generates a color depending on a value from 0 to 1
  503.         double height_ratio = (height_at_query - bbox.zmin()) / (bbox.zmax() - bbox.zmin());
  504.         colors = color_ramp.get(height_ratio);
  505.       }
  506.       raster_ofile.write (reinterpret_cast<char*>(&colors), 3);
  507.     }
  508.  
  509.   raster_ofile.close();
  510.  
  511.  
  512.   // Smooth heights with 5 successive Gaussian filters
  513.   double gaussian_variance = 4 * spacing * spacing;
  514.   for (TIN::Vertex_handle vh : dtm_clean.finite_vertex_handles())
  515.   {
  516.     double z = vh->point().z();
  517.     double total_weight = 1;
  518.  
  519.     TIN::Vertex_circulator circ = dtm_clean.incident_vertices (vh),
  520.       start = circ;
  521.  
  522.     do
  523.     {
  524.       if (!dtm_clean.is_infinite(circ))
  525.       {
  526.         double sq_dist = CGAL::squared_distance (vh->point(), circ->point());
  527.  
  528.         double weight = std::exp(- sq_dist / gaussian_variance);
  529.         z += weight * circ->point().z();
  530.         total_weight += weight;
  531.       }
  532.     }
  533.     while (++ circ != start);
  534.  
  535.     z /= total_weight;
  536.  
  537.     vh->point() = Point_3 (vh->point().x(), vh->point().y(), z);
  538.   }
  539.  
  540.  
  541.   std::array<double, 50> isovalues; // Contour 50 isovalues
  542.   for (std::size_t i = 0; i < isovalues.size(); ++ i)
  543.     isovalues[i] = bbox.zmin() + ((i+1) * (bbox.zmax() - bbox.zmin()) / (isovalues.size() - 2));
  544.  
  545.   // First find on each face if they are crossed by some isovalues and
  546.   // extract segments in a graph
  547.   using Segment_graph = boost::adjacency_list<boost::listS, boost::vecS, boost::undirectedS, Point_3>;
  548.   Segment_graph graph;
  549.   using Map_p2v = std::map<Point_3, Segment_graph::vertex_descriptor>;
  550.   Map_p2v map_p2v;
  551.   for (TIN::Face_handle vh : dtm_clean.finite_face_handles())
  552.     for (double iv : isovalues)
  553.       if (face_has_isovalue (vh, iv))
  554.       {
  555.         Segment_3 segment = isocontour_in_face (vh, iv);
  556.         for (const Point_3& p : { segment.source(), segment.target() })
  557.         {
  558.           // Only insert end points of segments once to get a well connected graph
  559.           Map_p2v::iterator iter;
  560.           bool inserted;
  561.           std::tie (iter, inserted) = map_p2v.insert (std::make_pair (p, Segment_graph::vertex_descriptor()));
  562.           if (inserted)
  563.           {
  564.             iter->second = boost::add_vertex (graph);
  565.             graph[iter->second] = p;
  566.           }
  567.         }
  568.         boost::add_edge (map_p2v[segment.source()], map_p2v[segment.target()], graph);
  569.       }
  570.  
  571.  
  572.  
  573.   // Split segments into polylines
  574.   std::vector<std::vector<Point_3> > polylines;
  575.   Polylines_visitor<Segment_graph> visitor (graph, polylines);
  576.   CGAL::split_graph_into_polylines (graph, visitor);
  577.  
  578.   std::cerr << polylines.size() << " polylines computed, with "
  579.             << map_p2v.size() << " vertices in total" << std::endl;
  580.  
  581.   // Output to WKT file
  582.   std::ofstream contour_ofile ("contour.wkt");
  583.   contour_ofile.precision(18);
  584.   CGAL::IO::write_multi_linestring_WKT (contour_ofile, polylines);
  585.   contour_ofile.close();
  586.  
  587.  
  588.  
  589.   // Construct constrained Delaunay triangulation with polylines as constraints
  590.   CTP ctp;
  591.   for (const std::vector<Point_3>& poly : polylines)
  592.     ctp.insert_constraint (poly.begin(), poly.end());
  593.  
  594.   // Simplification algorithm with limit on distance
  595.   PS::simplify (ctp, PS::Squared_distance_cost(), PS::Stop_above_cost_threshold (16 * spacing * spacing));
  596.  
  597.   polylines.clear();
  598.   for (CTP::Constraint_id cid : ctp.constraints())
  599.   {
  600.     polylines.push_back (std::vector<Point_3>());
  601.     polylines.back().reserve (ctp.vertices_in_constraint (cid).size());
  602.     for (CTP::Vertex_handle vh : ctp.vertices_in_constraint(cid))
  603.       polylines.back().push_back (vh->point());
  604.   }
  605.  
  606.   std::size_t nb_vertices
  607.     = std::accumulate (polylines.begin(), polylines.end(), std::size_t(0),
  608.                        [](std::size_t size, const std::vector<Point_3>& poly) -> std::size_t
  609.                        { return size + poly.size(); });
  610.  
  611.   std::cerr << nb_vertices
  612.             << " vertices remaining after simplification ("
  613.             << 100. * (nb_vertices / double(map_p2v.size())) << "%)" << std::endl;
  614.  
  615.   // Output to WKT file
  616.   std::ofstream simplified_ofile ("simplified.wkt");
  617.   simplified_ofile.precision(18);
  618.   CGAL::IO::write_multi_linestring_WKT (simplified_ofile, polylines);
  619.   simplified_ofile.close();
  620.  
  621.  
  622.  
  623.   // Get training from input
  624.   Point_set::Property_map<int> training_map;
  625.   bool training_found;
  626.   std::tie (training_map, training_found) = points.property_map<int>("training");
  627.  
  628.   if (training_found)
  629.   {
  630.     std::cerr << "Classifying ground/vegetation/building" << std::endl;
  631.  
  632.     // Create labels
  633.     Classification::Label_set labels ({ "ground", "vegetation", "building" });
  634.  
  635.     // Generate features
  636.     Classification::Feature_set features;
  637.     Classification::Point_set_feature_generator<Kernel, Point_set, Point_set::Point_map>
  638.       generator (points, points.point_map(), 5); // 5 scales
  639.  
  640. #ifdef CGAL_LINKED_WITH_TBB
  641.     // If TBB is used, features can be computed in parallel
  642.     features.begin_parallel_additions();
  643.     generator.generate_point_based_features (features);
  644.     features.end_parallel_additions();
  645. #else
  646.     generator.generate_point_based_features (features);
  647. #endif
  648.  
  649.     // Train a random forest classifier
  650.     Classification::ETHZ::Random_forest_classifier classifier (labels, features);
  651.     classifier.train (points.range(training_map));
  652.  
  653.     // Classify with graphcut regularization
  654.     Point_set::Property_map<int> label_map = points.add_property_map<int>("labels").first;
  655.     Classification::classify_with_graphcut<Concurrency_tag>
  656.       (points, points.point_map(), labels, classifier,
  657.        generator.neighborhood().k_neighbor_query(12), // regularize on 12-neighbors graph
  658.        0.5f, // graphcut weight
  659.        12, // Subdivide to speed-up process
  660.        label_map);
  661.  
  662.     // Evaluate
  663.     std::cerr << "Mean IoU on training data = "
  664.               << Classification::Evaluation(labels,
  665.                                             points.range(training_map),
  666.                                             points.range(label_map)).mean_intersection_over_union() << std::endl;
  667.  
  668.     // Save the classified point set
  669.     std::ofstream classified_ofile ("classified.ply");
  670.     CGAL::IO::set_binary_mode (classified_ofile);
  671.     classified_ofile << points;
  672.     classified_ofile.close();
  673.   }
  674.  
  675.  
  676.  
  677.   return EXIT_SUCCESS;
  678. }
  679.  
  680. ```
  681. # 8、附：Color_ramp.h
  682.  
  683. ```cpp
  684. #ifndef COLOR_RAMP_H
  685. #define COLOR_RAMP_H
  686.  
  687. class Color_ramp
  688. {
  689.   typedef std::array<unsigned char, 3> Color;
  690.   typedef std::pair<double, Color> Step;
  691.  
  692.   std::vector<Step> m_steps;
  693.  
  694. public:
  695.  
  696.   Color_ramp()
  697.   {
  698.     m_steps.push_back (std::make_pair (0, Color{ 192, 192, 255}));
  699.     m_steps.push_back (std::make_pair (0.2, Color{ 0, 0, 255}));
  700.     m_steps.push_back (std::make_pair (0.4, Color{ 0, 255, 0}));
  701.     m_steps.push_back (std::make_pair (0.6, Color{ 255, 255, 0}));
  702.     m_steps.push_back (std::make_pair (0.8, Color{ 255, 0, 0}));
  703.     m_steps.push_back (std::make_pair (1.0, Color{ 128, 0, 0}));
  704.   }
  705.  
  706.   Color get (double value) const
  707.   {
  708.     std::size_t idx = 0;
  709.     while (m_steps[idx+1].first < value)
  710.       ++ idx;
  711.  
  712.     double v0 = m_steps[idx].first;
  713.     double v1 = m_steps[idx+1].first;
  714.     const Color& c0 = m_steps[idx].second;
  715.     const Color& c1 = m_steps[idx+1].second;
  716.  
  717.     double ratio = (value - v0) / (v1 - v0);
  718.  
  719.     Color out;
  720.     for (std::size_t i = 0; i < 3; ++ i)
  721.       out[i] = static_cast<unsigned char>((1 - ratio) * c0[i] + ratio * c1[i]);
  722.  
  723.     return out;
  724.   }
  725.  
  726. };
  727.  
  728. #endif // COLOR_RAMP_H

声明：本文内容由网友自发贡献，不代表【wpsshop博客】立场，版权归原作者所有，本站不承担相应法律责任。如您发现有侵权的内容，请联系我们。转载请注明出处：https://www.wpsshop.cn/w/从前慢现在也慢/article/detail/522390
推荐阅读
相关标签

Copyright © 2003-2013 www.wpsshop.cn 版权所有，并保留所有权利。

  

闽ICP备14008679号