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() = -2
      fh->info() = -2;
    else // Pending faces are given info() = -1;
      fh->info() = -1;
  // Flooding algorithm
  std::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 threshold
    bool 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 faces
  Mesh dtm_mesh;
  CGAL::copy_face_graph (tin_with_info, dtm_mesh);

  // Save original DTM
  std::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

参考

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