CGAL 从DSM到DTM filtering
CGAL 从DSM到DTM filtering
上一节通过连通区域计算并将连通信息保存到三角面片中，获取了多个连通区域，本节将设置阈值将建筑物区域移除，生成一个最初的DTM。
建筑物区域去除
设置阈值为min_size，遍历三角面片，对连通区域（>= 0）且大于阈值的区域予以保留，剩下的面片予以移除。
代码
#include<iostream>
#include <queue>#include<CGAL/Surface_mesh.h>
#include<CGAL/Surface_mesh/IO/PLY.h>
#include <CGAL/draw_surface_mesh.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/Polygon_mesh_processing/triangulate_hole.h>
#include <CGAL/Polygon_mesh_processing/border.h>
#include <CGAL/Polygon_mesh_processing/remesh.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_2 = Kernel::Point_2;
using Point_3 = Kernel::Point_3;
using Mesh = CGAL::Surface_mesh<Point_3>;using Concurrency_tag = CGAL::Sequential_tag;using TIN = CGAL::Delaunay_triangulation_2<Projection_traits>;
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 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;/////! [Filtering]int min_size = int(mesh.vertices().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;CGAL::copy_face_graph (tin_with_info, dtm_mesh);// Save original DTMstd::ofstream dtm_ofile ("dtm_filtering.ply", std::ios_base::binary);CGAL::IO::set_binary_mode (dtm_ofile);CGAL::IO::write_PLY (dtm_ofile, dtm_mesh);dtm_ofile.close();return 0;}
如下图所示，移除建筑物的区域，并在孔洞处重新剖分为large faces大面:
132 connected component(s) found
10805 vertices(s) will be removed after filtering
构建编译运行
cmake -B build -S . -DCMAKE_TOOLCHAIN_FILE=D:\vcpkg\scripts\buildsystems\vcpkg.cmake
cmake --build build --config Debug
.\build\Debug\filtering.exe
参考
- https://doc.cgal.org/latest/Manual/tuto_gis.html