CGAL GIS 应用 - 从点云到DTM

GIS应用中使用的许多传感器(例如激光雷达)都会生成密集的点云。此类应用通常利用更高级的数据结构:例如，不规则三角网(TIN)，它可以作为数字高程模型(DEM)的基础，特别是用于生成数字地形模型(DTM)。
点云也可以通过分类信息来丰富，分类信息将点分割为地面、植被和建筑物点(或其他用户定义的标签)。

某些数据结构的定义可能因来源不同而有所不同。在本文中使用以下术语:

TIN:三角不规则网络，一种二维三角结构，根据它们在水平面上的投影连接3D点。
DSM:数字表面模型，包括建筑物和植被的整个扫描表面的模型。我们用一个TIN来存储DSM。
DTM:数字地形模型，一种没有建筑物或植被等物体的裸地表面模型。我们都使用TIN和光栅来存储DTM。
DEM:数字高程模型，一个更通用的术语，包括DSM和DTM`。

本文演示了以下场景。从输入点云，我们首先计算存储为TIN的DSM。然后，我们过滤掉与建筑立面或植被噪声相对应的过大的面片。保留与地面相对应的大型构件。对孔洞进行填充，并对得到的DEM重新网格化。由栅格DEM生成一组等高线折线。最后，对植被、建筑物和分组点进行有监督的三标签分类。

不规则三角网(TIN)

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

Digital Surface Model (DSM)

使用表示输入输出流stream操作符，可以很容易地将许多格式(XYZ、OFF、PLY、LAS)的点云加载到 CGAL::Point_set_3结构中。生成存储在TIN中的DSM很简单。

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

#include<iostream>
#include<CGAL/Surface_mesh.h>
#include<CGAL/Surface_mesh/IO/PLY.h>#include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
#include <CGAL/Projection_traits_xy_3.h>
#include <CGAL/Delaunay_triangulation_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>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;// Triangulated Irregular Network
using TIN = CGAL::Delaunay_triangulation_2<Projection_traits>;
using Mesh = CGAL::Surface_mesh<Point_3>;int main() {std::ifstream ifile("points.xyz", 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());// Write DSMMesh 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();return 0;
}

输入points.xyz：
在这里插入图片描述
输出DSM ply文件:

注：由密集点云生成三角网格Mesh，可以使用GIS等软件如QGIS TIN Mesh Creation工具，
参考：QGIS 高程点生成Mesh

官方示例：b9_training.ply点云数据生成DSM

  const std::string fname = argc != 2 ? CGAL::data_file_path("../data/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();

在这里插入图片描述

数字地形模型(DTM)

生成的DSM被用作DTM计算的基础，即地面表示为过滤掉非地面点后的另一个TIN。
本文提出一个简单的DTM估计方法，分解如下步骤:

对面片的高度进行阈值化，以消除高度的剧烈变化;
将其他面片聚类为连通分支；
过滤所有小于用户定义阈值的组件。
该算法依赖于两个参数:一个高度阈值对应于建筑物的最小高度，一个周长阈值对应于二维投影上建筑物的最大尺寸。

CGAL从DSM到DTM区域提取
coloredTIN
CGAL从DSM到DTM filtering

CGAL从DSM到DTM 孔洞填充及网格细化
DTM

栅格化

TIN数据结构可以与重心坐标相结合，以便插值并栅格化顶点中嵌入的任何信息所需的任何分辨率的高度图。

示例参见：

Mesh 网格曲面栅格化
虚幻地形高度图生成及测试

生成栅格PPM：
在这里插入图片描述
生成栅格PNG:

提取等高线

提取TIN上定义的函数的等值层是另一个可以使用CGAL完成的应用。本文演示如何提取等值高度来构建地形图。

步骤：

构建等高线图Graph；
Graph分裂为多折线Polyline;
Polyline简化。

详细示例及代码参见：

Mesh地形曲面提取等高线

在这里插入图片描述

Classifying 分类

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

// 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 ("classification_gis_tutorial.ply");CGAL::IO::set_binary_mode (classified_ofile);classified_ofile << points;classified_ofile.close();}