烟台网站建设公司,asp网站模板源码免费无限下载,湖北网站设计,仿站在线工具GIS (Geographic Information System) 地理信息系统 原文地址: https://doc.cgal.org/latest/Manual/tuto_gis.html GIS 应用中使用的许多传感器#xff08;例如激光雷达#xff09;都会生成密集的点云。此类应用程序通常利用更先进的数据结构#xff1a;例如#xff0c;不…GIS (Geographic Information System) 地理信息系统 原文地址: https://doc.cgal.org/latest/Manual/tuto_gis.html GIS 应用中使用的许多传感器例如激光雷达都会生成密集的点云。此类应用程序通常利用更先进的数据结构例如不规则三角网络 (TIN) 可以作为数字高程模型 (DEM) 的基础特别是用于生成数字地形模型 (DTM) 的基础。点云还可以通过将点分割为地面、植被和建筑点或其他用户定义的标签的分类信息来丰富。
一些数据结构的定义可能根据不同的来源而有所不同。我们将在本教程中使用以下术语
TIN不规则三角网络一种 2D 三角测量结构根据 3D 点在水平面上的投影来连接它们。DSM数字表面模型包括建筑物和植被的整个扫描表面的模型。我们使用 TIN 来存储 DSM。DTM数字地形模型没有建筑物或植被等物体的裸露地面模型。我们都使用 TIN 和栅格来存储 DTM。DEM数字高程模型一个更通用的术语包括 DSM 和 DTM。
本教程说明了以下场景。根据输入点云我们首先计算存储为 TIN 的 DSM。然后我们过滤掉与建筑物外墙或植被噪声相对应的过大的面。保留与地面相对应的大组件。孔被填充获得的 DEM 被重新网格化。由此生成栅格 DEM 和一组等高线折线。最后执行监督三标签分类来分割植被、建筑物和组点。
1 不规则三角网 (TIN)
CGAL 提供了多种三角测量数据结构和算法。TIN 可以通过将 2D Delaunay 三角剖分与投影特征相结合来生成三角剖分结构是使用沿选定平面通常为 XY 平面的点的 2D 位置计算的而点的 3D 位置则保留可视化和测量。
因此TIN 数据结构可以简单地定义如下
using Kernel CGAL::Exact_predicates_inexact_constructions_kernel;
using Projection_traits CGAL::Projection_traits_xy_3Kernel;
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_2Projection_traits;2 数字表面模型 (DSM)
CGAL::Point_set_3使用流运算符可以轻松地将多种格式XYZ、OFF、PLY、LAS的点云加载到结构中。生成存储在 TIN 中的 DSM 就很简单 // Read pointsstd::ifstream ifile (fname, std::ios_base::binary);CGAL::Point_set_3Point_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_meshPoint_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();图 0.1给出了在 San Fransisco 数据集上以这种方式计算的 DSM 示例请参阅参考资料 。
3 数字地形模型 (DTM)
生成的 DSM 用作 DTM 计算的基础即过滤非地面点后将地面表示为另一个 TIN。
作为一个例子我们提出了一个简单的 DTM 估计分解为以下步骤
设定面的高度阈值以消除高度的剧烈变化将其他方面聚类成连接的组件过滤所有小于用户定义阈值的分量
该算法依赖于 2 个参数与建筑物的最小高度相对应的高度阈值以及与 2D 投影上建筑物的最大尺寸相对应的周长阈值。
3.1 含信息的 TIN
由于它利用了灵活的 CGAL Delaunay 三角剖分 API因此我们的 TIN 可以通过有关顶点和/或面的信息来丰富。在我们的例子中每个顶点都会跟踪输入点云中对应点的索引这将允许随后过滤地面点并且每个面都被赋予其连接组件的索引。 auto idx_to_point_with_info [](const Point_set::Index idx) - std::pairPoint_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 识别连接的组件
连接组件通过泛洪算法进行识别从种子面开始所有入射面都插入到当前连接组件中除非它们的高度超过用户定义的阈值。 double spacing CGAL::compute_average_spacingConcurrency_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((i1)%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::vectorint component_size;for (TIN_with_info::Face_handle fh : tin_with_info.finite_face_handles()){if (fh-info() ! -1)continue;std::queueTIN_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_mapMesh::Face_index, CGAL::IO::Colorcolor_map tin_colored_mesh.add_property_mapMesh::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::pairTIN_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();图 0.2给出了由连接分量着色的 TIN 示例。 3.3 过滤
小于最大建筑物的组件可以通过以下方式移除 int min_size int(points.size() / 2);std::vectorTIN_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 its 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 3D 表示因此额外的步骤可以帮助生成形状更好的网格删除大于阈值的面并使用孔填充算法进行填充对孔进行三角测量、细化和平整。
以下代码片段将 TIN 复制到网格中同时过滤掉过大的面然后识别孔洞并填充除最大的孔洞即外壳之外的所有孔洞。 // Copy and keep track of overly large facesMesh dtm_mesh;std::vectorMesh::Face_index face_selection;Mesh::Property_mapMesh::Face_index, bool face_selection_map dtm_mesh.add_property_mapMesh::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::pairTIN_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((i1)%3)-point(),ff.first-vertex((i2)%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 its 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::vectorMesh::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::parameters::fairing_continuity(0));// 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();
各向同性重新网格化也可以作为最后一步执行以生成不受 2D Delaunay 面形状约束的更规则的网格。 CGAL::Polygon_mesh_processing::isotropic_remeshing (faces(dtm_mesh), spacing, dtm_mesh);std::ofstream dtm_remeshed_ofile (dtm_remeshed.ply, std::ios_base::binary);CGAL::IO::set_binary_mode (dtm_remeshed_ofile);CGAL::IO::write_PLY (dtm_remeshed_ofile, dtm_mesh);dtm_remeshed_ofile.close();
图 0.3显示了这些不同步骤如何影响输出网格图 0.4显示了 DTM 各向同性网格。 4 光栅化
TIN 数据结构可以与重心坐标相结合以便在需要嵌入顶点信息的任何分辨率下对高度图进行插值和栅格化。
由于最新的两个步骤孔填充和重新划分网格是在 3D 网格上执行的因此我们的 DTM 是 2.5D 表示的假设可能不再有效。因此我们首先使用最后计算的各向同性 DTM 网格的顶点重建 TIN。
以下代码片段以彩虹渐变 PPM 文件简单位图格式的形式生成高度的光栅图像 CGAL::Bbox_3 bbox CGAL::bbox_3 (points.points().begin(), points.points().end());// Generate raster image 1920-pixels 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::arrayunsigned char, 3 colors { 0, 0, 0 };if (!dtm_clean.is_infinite(location)){std::arraydouble, 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_castchar*(colors), 3);}raster_ofile.close();图 0.5给出了带有表示高度的彩虹斜坡的光栅图像的示例。 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((i1) % 3)-point();Point_3 p1 fh-vertex((i2) % 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::arraydouble, 50 isovalues; // Contour 50 isovaluesfor (std::size_t i 0; i isovalues.size(); i)isovalues[i] bbox.zmin() ((i1) * (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_listboost::listS, boost::vecS, boost::undirectedS, Point_3;Segment_graph graph;using Map_p2v std::mapPoint_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::vectorstd::vectorPoint_3 polylines;Polylines_visitorSegment_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::vectorstd::vectorPoint_3 polylines;Graph graph;
public:Polylines_visitor (Graph graph, std::vectorstd::vectorPoint_3 polylines): polylines (polylines), graph(graph) { }void start_new_polyline(){polylines.push_back (std::vectorPoint_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_2Projection_traits;
using CDT_face_base CGAL::Constrained_triangulation_face_base_2Projection_traits;
using CDT_TDS CGAL::Triangulation_data_structure_2CDT_vertex_base, CDT_face_base;
using CDT CGAL::Constrained_Delaunay_triangulation_2Projection_traits, CDT_TDS;
using CTP CGAL::Constrained_triangulation_plus_2CDT;
以下代码根据到原始多段线的平方距离简化多段线集当在不超过平均间距 4 倍的情况下无法删除更多顶点时停止。 // Construct constrained Delaunay triangulation with polylines as constraintsCTP ctp;for (const std::vectorPoint_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::vectorPoint_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::vectorPoint_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();图 0.6给出了等高线和简化的示例。 图 0.6使用 50 个均匀分布的等值线绘制轮廓。顶部使用 148k 顶点和简化的原始轮廓其公差等于输入点云的平均间距剩余原始多段线顶点的 3.4%。底部公差为平均间距的 4 倍剩余顶点的 1.3%和公差为平均间距的 10 倍剩余顶点的 0.9%的简化。折线在所有情况下都是不相交的。
6 分类
CGAL 提供了一个包 Classification可用于将点云分割为用户定义的标签集。目前 CGAL 中最先进的分类器是 ETHZ 的随机森林。由于它是一个监督分类器因此需要一个训练集。
以下代码片段展示了如何使用一些手动选择的训练集来训练随机森林分类器并计算通过图割算法正则化的分类 // Get training from inputPoint_set::Property_mapint training_map;bool training_found;std::tie (training_map, training_found) points.property_mapint(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_generatorKernel, Point_set, Point_set::Point_mapgenerator (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_mapint label_map points.add_property_mapint(labels).first;Classification::classify_with_graphcutConcurrency_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 (classification_gis_tutorial.ply);CGAL::IO::set_binary_mode (classified_ofile);classified_ofile points;classified_ofile.close();}图 0.7给出了训练集和分类结果的示例。 7 完整代码示例
本教程中使用的所有代码片段都可以组合起来创建完整的 GIS 管道前提是使用正确的*包含内容。*我们提供了一个完整的代码示例它实现了本教程中描述的所有步骤。
#include CGAL/Exact_predicates_inexact_constructions_kernel.h
#include CGAL/Projection_traits_xy_3.h
#include CGAL/Delaunay_triangulation_2.h
#include CGAL/Triangulation_vertex_base_with_info_2.h
#include CGAL/Triangulation_face_base_with_info_2.h
#include CGAL/boost/graph/graph_traits_Delaunay_triangulation_2.h
#include CGAL/boost/graph/copy_face_graph.h
#include CGAL/Point_set_3.h
#include CGAL/Point_set_3/IO.h
#include CGAL/compute_average_spacing.h
#include CGAL/Surface_mesh.h
#include CGAL/Polygon_mesh_processing/locate.h
#include CGAL/Polygon_mesh_processing/triangulate_hole.h
#include CGAL/Polygon_mesh_processing/border.h
#include CGAL/Polygon_mesh_processing/remesh.h
#include boost/graph/adjacency_list.hpp
#include CGAL/boost/graph/split_graph_into_polylines.h
#include CGAL/IO/WKT.h
#include CGAL/Constrained_Delaunay_triangulation_2.h
#include CGAL/Constrained_triangulation_plus_2.h
#include CGAL/Polyline_simplification_2/simplify.h
#include CGAL/Polyline_simplification_2/Squared_distance_cost.h
#include CGAL/Classification.h
#include CGAL/Random.h
#include fstream
#include queue
#include include/Color_ramp.h
using Kernel CGAL::Exact_predicates_inexact_constructions_kernel;
using Projection_traits CGAL::Projection_traits_xy_3Kernel;
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_2Projection_traits;
// Triangulated Irregular Network (with info)
using Point_set CGAL::Point_set_3Point_3;
using Vbi CGAL::Triangulation_vertex_base_with_info_2 Point_set::Index, Projection_traits;
using Fbi CGAL::Triangulation_face_base_with_info_2int, Projection_traits;
using TDS CGAL::Triangulation_data_structure_2Vbi, Fbi;
using TIN_with_info CGAL::Delaunay_triangulation_2Projection_traits, TDS;
namespace Classification CGAL::Classification;
#ifdef CGAL_LINKED_WITH_TBB
using Concurrency_tag CGAL::Parallel_tag;
#else
using Concurrency_tag CGAL::Sequential_tag;
#endif
bool face_has_isovalue (TIN::Face_handle fh, double isovalue)
{bool above false, below false;for (int i 0; i 3; i){// Face has isovalue if one of its vertices is above and another// one 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((i1) % 3)-point();Point_3 p1 fh-vertex((i2) % 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::vectorstd::vectorPoint_3 polylines;Graph graph;
public:Polylines_visitor (Graph graph, std::vectorstd::vectorPoint_3 polylines): polylines (polylines), graph(graph) { }void start_new_polyline(){polylines.push_back (std::vectorPoint_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_2Projection_traits;
using CDT_face_base CGAL::Constrained_triangulation_face_base_2Projection_traits;
using CDT_TDS CGAL::Triangulation_data_structure_2CDT_vertex_base, CDT_face_base;
using CDT CGAL::Constrained_Delaunay_triangulation_2Projection_traits, CDT_TDS;
using CTP CGAL::Constrained_triangulation_plus_2CDT;
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_3Point_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_meshPoint_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::pairPoint_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_spacingConcurrency_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((i1)%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::vectorint component_size;for (TIN_with_info::Face_handle fh : tin_with_info.finite_face_handles()){if (fh-info() ! -1)continue;std::queueTIN_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_mapMesh::Face_index, CGAL::IO::Colorcolor_map tin_colored_mesh.add_property_mapMesh::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::pairTIN_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::vectorTIN_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 its 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::vectorMesh::Face_index face_selection;Mesh::Property_mapMesh::Face_index, bool face_selection_map dtm_mesh.add_property_mapMesh::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::pairTIN_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((i1)%3)-point(),ff.first-vertex((i2)%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 its 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::vectorMesh::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::parameters::fairing_continuity(0));// 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::arrayunsigned char, 3 colors { 0, 0, 0 };if (!dtm_clean.is_infinite(location)){std::arraydouble, 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_castchar*(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::arraydouble, 50 isovalues; // Contour 50 isovaluesfor (std::size_t i 0; i isovalues.size(); i)isovalues[i] bbox.zmin() ((i1) * (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_listboost::listS, boost::vecS, boost::undirectedS, Point_3;Segment_graph graph;using Map_p2v std::mapPoint_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::vectorstd::vectorPoint_3 polylines;Polylines_visitorSegment_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::vectorPoint_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::vectorPoint_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::vectorPoint_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_mapint training_map;bool training_found;std::tie (training_map, training_found) points.property_mapint(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_generatorKernel, Point_set, Point_set::Point_mapgenerator (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_mapint label_map points.add_property_mapint(labels).first;Classification::classify_with_graphcutConcurrency_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 (classification_gis_tutorial.ply);CGAL::IO::set_binary_mode (classified_ofile);classified_ofile points;classified_ofile.close();}return EXIT_SUCCESS;
}8 参考
本教程基于以下 CGAL 包
2D 三角测量参考3D 点集参考点集处理参考表面网格参考CGAL 和 Boost 图库参考多边形网格处理参考2D 折线简化参考分类参考
本教程中使用的数据集来自https://www.usgs.gov/数据库已获得美国政府公共领域许可。
