CGAL 从DSM到DTM-建筑物区域提取

CGAL 从DSM到DTM-建筑物区域提取

生成的DSM被用作DTM计算的基础，即地面表示为过滤掉非地面点后的另一个TIN。主要是去除一些建筑物和植被非地形点。

建筑物立面及连通区域提取

建筑物立面的特征是三角形面片的高度变化剧烈。
通过遍历每一个三角面片，计算顶点之间最大的高差face_height，与事先设置的阈值spacing比较，大于阈值或包含无限远点infinite_vertex标记为建筑物立面。

除去建筑物立面，对剩下的三角面片使用洪水算法确定连通区域，并标记对应的区域id，取值范围：[0，component_size.size()]。

代码

#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/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/compute_average_spacing.h>using Kernel = CGAL::Exact_predicates_inexact_constructions_kernel;using Projection_traits = CGAL::Projection_traits_xy_3<Kernel>;using Point_3 = Kernel::Point_3;
using Mesh = CGAL::Surface_mesh<Point_3>;using Concurrency_tag = CGAL::Sequential_tag;using Vbi = CGAL::Triangulation_vertex_base_with_info_2 <Mesh::Vertex_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>;int main() {Mesh mesh;CGAL::IO::read_PLY("./data/dsm.ply", mesh);auto idx_to_point_with_info= [&](const Mesh::Vertex_index& idx) -> std::pair<Point_3, Mesh::Vertex_index>{return std::make_pair ( mesh.point(idx), idx);};TIN_with_info tin_with_info(boost::make_transform_iterator (mesh.vertices().begin(), idx_to_point_with_info),boost::make_transform_iterator (mesh.vertices().end(), idx_to_point_with_info));double spacing = CGAL::compute_average_spacing<Concurrency_tag>(mesh.points(), 6);spacing *= 2;auto face_height= [&](const TIN_with_info::Face_handle fh) -> double{double out = 0.;for (int i = 0; i < 3; ++ i)out = (std::max) (out, CGAL::abs(fh->vertex(i)->point().z() - fh->vertex((i+1)%3)->point().z()));return out;};// Initialize faces info 初始化三角面附件信息（int类型）for (TIN_with_info::Face_handle fh : tin_with_info.all_face_handles())if (tin_with_info.is_infinite(fh) || face_height(fh) > spacing) // Filtered faces are given info() = -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){// 灰色标记建筑物立面及待处理的其他非地面点if (ff.first->info() < 0){color_map[ff.second] = CGAL::IO::Color(128, 128, 128);// put(color_map, ff.second, CGAL::IO::Color(128, 128, 128));}else{// 其他标记随机颜色CGAL::Random r (ff.first->info());color_map[ff.second] = CGAL::IO::Color (r.get_int(64, 192),r.get_int(64, 192),r.get_int(64, 192));/*put(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();return 0;
}

如下图所示，灰色标记为建筑物立面，其他联通区域标记随机颜色：

构建编译运行

cmake -B build -S . -DCMAKE_TOOLCHAIN_FILE=D:\vcpkg\scripts\buildsystems\vcpkg.cmake
cmake --build build --config Debug
.\build\Debug\coloredTIN.exe

参考

1. https://doc.cgal.org/latest/Manual/tuto_gis.html