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数据结构可以简单地定义如下:

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>;

2、数字表面模型(DSM)

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

// Read pointsstd::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 DSMTIN 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();

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

3、数字地形模型(DTM)

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

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

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

3.1 带有信息的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));

3.2 识别连接组件

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

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 infofor (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() = -2fh->info() = -2;else // Pending faces are given info() = -1;fh->info() = -1;// Flooding algorithmstd::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 grayif (ff.first->info() < 0)color_map[ff.second] = CGAL::IO::Color(128, 128, 128);else{// Random color seeded by the component IDCGAL::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();

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

3.3 数据清洗

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

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 thresholdbool 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);

3.4 孔洞填充和网格重建

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

// Copy and keep track of overly large facesMesh 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 DTMstd::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 removalCGAL::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 DTMstd::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 holesstd::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 bigest (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::Emptyset_iterator(), CGAL::Emptyset_iterator());// Save DTM with holes filledstd::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();

为了产生不受二维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();

4 光栅化

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

CGAL::Bbox_3 bbox = CGAL::bbox_3 (points.points().begin(), points.points().end());// Generate raster image 1920-pixels largestd::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 simplicitystd::ofstream raster_ofile ("raster.ppm", std::ios_base::binary);// PPM headerraster_ofile << "P6" << std::endl // magic number<< width << " " << height << std::endl // dimensions of the image<< 255 << std::endl; // maximum color value// Use rainbow color ramp outputColor_ramp color_ramp;// Keeping track of location from one point to its neighbor allows// for fast locate in DTTIN::Face_handle location;// Query each pixel of the imagefor (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 2Dlocation = dtm_clean.locate (query, location);// Points outside the convex hull will be colored blackstd::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 1double 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();

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

5 等高线生成

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

5.1绘制等高线图

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

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 belowif (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 vertexdouble 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 isovaluesfor (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 graphusing 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 graphMap_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);}

5.2分割成折线

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

// Split segments into polylinesstd::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 filestd::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 polylinesif (polylines.back().size() < 50)polylines.pop_back();}
};

5.3 等高线简化

由于输出的噪声很大，用户可能希望简化折线。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 constraintsCTP ctp;for (const std::vector<Point_3>& poly : polylines)ctp.insert_constraint (poly.begin(), poly.end());// Simplification algorithm with limit on distancePS::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 filestd::ofstream simplified_ofile ("simplified.wkt");simplified_ofile.precision(18);CGAL::IO::write_multi_linestring_WKT (simplified_ofile, polylines);simplified_ofile.close();

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

6、点云分类

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

// Get training from inputPoint_set::Property_map<int> training_map;bool training_found;std::tie (training_map, training_found) = points.property_map<int>("training");if (training_found){std::cerr << "Classifying ground/vegetation/building" << std::endl;// Create labelsClassification::Label_set labels ({ "ground", "vegetation", "building" });// Generate featuresClassification::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 parallelfeatures.begin_parallel_additions();generator.generate_point_based_features (features);features.end_parallel_additions();
#elsegenerator.generate_point_based_features (features);
#endif// Train a random forest classifierClassification::ETHZ::Random_forest_classifier classifier (labels, features);classifier.train (points.range(training_map));// Classify with graphcut regularizationPoint_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 graph0.5f, // graphcut weight12, // Subdivide to speed-up processlabel_map);// Evaluatestd::cerr << "Mean IoU on training data = "<< Classification::Evaluation(labels,points.range(training_map),points.range(label_map)).mean_intersection_over_union() << std::endl;// Save the classified point setstd::ofstream classified_ofile ("classified.ply");CGAL::IO::set_binary_mode (classified_ofile);classified_ofile << points;classified_ofile.close();}

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

7、完整代码示例

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

#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;
#endifbool 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 belowif (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 vertexdouble 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 polylinesif (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 pointsstd::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 DSMTIN 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 infofor (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() = -2fh->info() = -2;else // Pending faces are given info() = -1;fh->info() = -1;// Flooding algorithmstd::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 grayif (ff.first->info() < 0)color_map[ff.second] = CGAL::IO::Color(128, 128, 128);else{// Random color seeded by the component IDCGAL::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 thresholdbool 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 facesMesh 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 DTMstd::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 removalCGAL::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 DTMstd::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 holesstd::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 bigest (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::Emptyset_iterator(), CGAL::Emptyset_iterator());// Save DTM with holes filledstd::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 largestd::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 simplicitystd::ofstream raster_ofile ("raster.ppm", std::ios_base::binary);// PPM headerraster_ofile << "P6" << std::endl // magic number<< width << " " << height << std::endl // dimensions of the image<< 255 << std::endl; // maximum color value// Use rainbow color ramp outputColor_ramp color_ramp;// Keeping track of location from one point to its neighbor allows// for fast locate in DTTIN::Face_handle location;// Query each pixel of the imagefor (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 2Dlocation = dtm_clean.locate (query, location);// Points outside the convex hull will be colored blackstd::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 1double 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 filtersdouble 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 isovaluesfor (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 graphusing 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 graphMap_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 polylinesstd::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 filestd::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 constraintsCTP ctp;for (const std::vector<Point_3>& poly : polylines)ctp.insert_constraint (poly.begin(), poly.end());// Simplification algorithm with limit on distancePS::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 filestd::ofstream simplified_ofile ("simplified.wkt");simplified_ofile.precision(18);CGAL::IO::write_multi_linestring_WKT (simplified_ofile, polylines);simplified_ofile.close();// Get training from inputPoint_set::Property_map<int> training_map;bool training_found;std::tie (training_map, training_found) = points.property_map<int>("training");if (training_found){std::cerr << "Classifying ground/vegetation/building" << std::endl;// Create labelsClassification::Label_set labels ({ "ground", "vegetation", "building" });// Generate featuresClassification::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 parallelfeatures.begin_parallel_additions();generator.generate_point_based_features (features);features.end_parallel_additions();
#elsegenerator.generate_point_based_features (features);
#endif// Train a random forest classifierClassification::ETHZ::Random_forest_classifier classifier (labels, features);classifier.train (points.range(training_map));// Classify with graphcut regularizationPoint_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 graph0.5f, // graphcut weight12, // Subdivide to speed-up processlabel_map);// Evaluatestd::cerr << "Mean IoU on training data = "<< Classification::Evaluation(labels,points.range(training_map),points.range(label_map)).mean_intersection_over_union() << std::endl;// Save the classified point setstd::ofstream classified_ofile ("classified.ply");CGAL::IO::set_binary_mode (classified_ofile);classified_ofile << points;classified_ofile.close();}return EXIT_SUCCESS;
}```
# 8、附：Color_ramp.h```cpp
#ifndef COLOR_RAMP_H
#define COLOR_RAMP_Hclass Color_ramp
{typedef std::array<unsigned char, 3> Color;typedef std::pair<double, Color> Step;std::vector<Step> m_steps;public:Color_ramp(){m_steps.push_back (std::make_pair (0, Color{ 192, 192, 255}));m_steps.push_back (std::make_pair (0.2, Color{ 0, 0, 255}));m_steps.push_back (std::make_pair (0.4, Color{ 0, 255, 0}));m_steps.push_back (std::make_pair (0.6, Color{ 255, 255, 0}));m_steps.push_back (std::make_pair (0.8, Color{ 255, 0, 0}));m_steps.push_back (std::make_pair (1.0, Color{ 128, 0, 0}));}Color get (double value) const{std::size_t idx = 0;while (m_steps[idx+1].first < value)++ idx;double v0 = m_steps[idx].first;double v1 = m_steps[idx+1].first;const Color& c0 = m_steps[idx].second;const Color& c1 = m_steps[idx+1].second;double ratio = (value - v0) / (v1 - v0);Color out;for (std::size_t i = 0; i < 3; ++ i)out[i] = static_cast<unsigned char>((1 - ratio) * c0[i] + ratio * c1[i]);return out;}};#endif // COLOR_RAMP_H