CGAL笔记之形状重建—多边形曲面重建

最新推荐文章于 2025-04-17 19:16:05 发布

3333yyt

最新推荐文章于 2025-04-17 19:16:05 发布

阅读量1.2k

点赞数 3

分类专栏： CGAL 文章标签： c++ 图形渲染算法数据结构 3d

本文链接：https://blog.youkuaiyun.com/younggung/article/details/129338312

版权

CGAL 专栏收录该内容

13 篇文章

订阅专栏

本文介绍了使用CGAL库进行多边形曲面重建的方法，包括从点云数据中提取平面，通过RANSAC或区域生长算法，以及通过优化选择的平面组合生成流形和水密的多边形模型。该过程涉及处理噪声、异常值，并通过整数规划解决约束优化问题，以控制模型的复杂度和细节水平。

摘要生成于 C知道，由 DeepSeek-R1 满血版支持，前往体验 >

1 引言

该软件包实现了基于假设和选择的方法，用于从点云分段重建平面对象。该方法将从分段平面对象采样的无序点集作为输入。输出是一个紧凑且密闭的表面网格，用于插补输入点集。该方法假定提供了所有必要的主平面（或者可以使用形状检测中描述的形状检测方法或任何其他替代方法从输入点集中提取）。

CGAL 中现有的表面重建方法（即泊松曲面重建、推进前曲面重建和尺度-空间曲面重建）适用于表示由光滑曲面描述的对象的点集。对于建筑物等人造物体，由于重建模型的缺陷和复杂性（即巨大的网格、缺失的区域、噪声和不需要的结构），结果并不令人满意。这主要是因为这些方法倾向于密切关注表面细节。此软件包中引入的算法旨在生成简化且密闭的曲面网格。它可以恢复物体的尖锐特征，并且可以处理大量的噪声和异常值，补充其他表面重建方法。

需要混合整数规划（MIP）求解器来求解该方法制定的约束优化问题

2 算法

与传统的分段平面对象重建方法不同，该方法侧重于提取良好的几何基元或获得适当的基元排列，该方法的重点在于与基元（即平面）相交并寻求它们的适当组合以获得流形和水密多边形曲面模型。

该方法将曲面重建转换为基于假设和选择策略的二元标记问题。重建包括以下步骤：

从输入点集中提取平面（如果平面已知或通过其他方式提供，则可以跳过）;
通过与提取的平面基元相交来生成一组候选面;
通过在硬约束下进行优化来选择候选面的最佳子集，从而强制最终多边形表面在拓扑上正确。

（a）设置输入点;（b）提取的平面段;（c）成对交集生成的候选面;（d）通过优化选择的人脸;（e）重建模型。

在这里插入图片描述

3 示例

3.1 点云重建

该方法假定可以从输入点集中提取所有必需的平面。以下示例首先从输入点云中提取平面，然后重建表面模型。在第一个示例中，高效 RANSAC 方法用于提取平面。它非常快，但不是确定性的，与第二个示例中的区域增长方法相反，后者更慢，但更精确，并且始终为相同的给定参数返回相同的结果。

#include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
#include <CGAL/IO/read_points.h>
#include <CGAL/property_map.h>
#include <CGAL/Surface_mesh.h>
#include <CGAL/Shape_detection/Efficient_RANSAC.h>
#include <CGAL/Polygonal_surface_reconstruction.h>
#ifdef CGAL_USE_SCIP  // defined (or not) by CMake scripts, do not define by hand
#include <CGAL/SCIP_mixed_integer_program_traits.h>
typedef CGAL::SCIP_mixed_integer_program_traits<double>                        MIP_Solver;
#elif defined(CGAL_USE_GLPK)  // defined (or not) by CMake scripts, do not define by hand
#include <CGAL/GLPK_mixed_integer_program_traits.h>
typedef CGAL::GLPK_mixed_integer_program_traits<double>                        MIP_Solver;
#endif
#if defined(CGAL_USE_GLPK) || defined(CGAL_USE_SCIP)
#include <CGAL/Timer.h>
#include <fstream>
typedef CGAL::Exact_predicates_inexact_constructions_kernel                        Kernel;
typedef Kernel::Point_3                                                                                                Point;
typedef Kernel::Vector_3                                                                                        Vector;
// Point with normal, and plane index
typedef boost::tuple<Point, Vector, int>                                                        PNI;
typedef std::vector<PNI>                                                                                        Point_vector;
typedef CGAL::Nth_of_tuple_property_map<0, PNI>                                                Point_map;
typedef CGAL::Nth_of_tuple_property_map<1, PNI>                                                Normal_map;
typedef CGAL::Nth_of_tuple_property_map<2, PNI>                                                Plane_index_map;
typedef CGAL::Shape_detection::Efficient_RANSAC_traits<Kernel, Point_vector, Point_map, Normal_map>     Traits;
typedef CGAL::Shape_detection::Efficient_RANSAC<Traits>             Efficient_ransac;
typedef CGAL::Shape_detection::Plane<Traits>                                                Plane;
typedef CGAL::Shape_detection::Point_to_shape_index_map<Traits>     Point_to_shape_index_map;
typedef        CGAL::Polygonal_surface_reconstruction<Kernel>                                Polygonal_surface_reconstruction;
typedef CGAL::Surface_mesh<Point>                                                                        Surface_mesh;
/*
* This example first extracts planes from the input point cloud
* (using RANSAC with default parameters) and then reconstructs
* the surface model from the planes.
*/
int main()
{
  Point_vector points;
  // Loads point set from a file.
  const std::string input_file(CGAL::data_file_path("points_3/cube.pwn"));
    std::ifstream input_stream(input_file.c_str());
  if (input_stream.fail()) {
    std::cerr << "failed open file \'" <<input_file << "\'" << std::endl;
    return EXIT_FAILURE;
  }
  input_stream.close();
  std::cout << "Loading point cloud: " << input_file << "...";
  CGAL::Timer t;
  t.start();
  if (!CGAL::IO::read_points(input_file.c_str(), std::back_inserter(points),
                             CGAL::parameters::point_map(Point_map()).normal_map(Normal_map())))
  {
    std::cerr << "Error: cannot read file " << input_file << std::endl;
    return EXIT_FAILURE;
  }
  else
    std::cout << " Done. " << points.size() << " points. Time: " << t.time() << " sec." << std::endl;
  // Shape detection
  Efficient_ransac ransac;
  ransac.set_input(points);
  ransac.add_shape_factory<Plane>();
  std::cout << "Extracting planes...";
  t.reset();
  ransac.detect();
  Efficient_ransac::Plane_range planes = ransac.planes();
  std::size_t num_planes = planes.size();
  std::cout << " Done. " << num_planes << " planes extracted. Time: " << t.time() << " sec." << std::endl;
  // Stores the plane index of each point as the third element of the tuple.
  Point_to_shape_index_map shape_index_map(points, planes);
  for (std::size_t i = 0; i < points.size(); ++i) {
    // Uses the get function from the property map that accesses the 3rd element of the tuple.
    int plane_index = get(shape_index_map, i);
    points[i].get<2>() = plane_index;
  }
  std::cout << "Generating candidate faces...";
  t.reset();
  Polygonal_surface_reconstruction algo(
    points,
    Point_map(),
    Normal_map(),
    Plane_index_map()
  );
  std::cout << " Done. Time: " << t.time() << " sec." << std::endl;
  Surface_mesh model;
  std::cout << "Reconstructing...";
  t.reset();
  if (!algo.reconstruct<MIP_Solver>(model)) {
    std::cerr << " Failed: " << algo.error_message() << std::endl;
    return EXIT_FAILURE;
  }
  const std::string& output_file("data/cube_result.off");
  if (CGAL::IO::write_OFF(output_file, model))
    std::cout << " Done. Saved to " << output_file << ". Time: " << t.time() << " sec." << std::endl;
  else {
    std::cerr << " Failed saving file." << std::endl;
    return EXIT_FAILURE;
  }
  // Also stores the candidate faces as a surface mesh to a file
  Surface_mesh candidate_faces;
  algo.output_candidate_faces(candidate_faces);
  const std::string& candidate_faces_file("data/cube_candidate_faces.off");
  std::ofstream candidate_stream(candidate_faces_file.c_str());
  if (CGAL::IO::write_OFF(candidate_stream, candidate_faces))
    std::cout << "Candidate faces saved to " << candidate_faces_file << "." << std::endl;
  return EXIT_SUCCESS;
}
#else
int main(int, char**)
{
    std::cerr << "This test requires either GLPK or SCIP.\n";
    return EXIT_SUCCESS;
}
#endif  // defined(CGAL_USE_GLPK) || defined(CGAL_USE_SCIP)

3.2 使用用户提供的平面进行重建

以下示例显示了使用用户提供的平面线段进行重建。

#include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
#include <CGAL/IO/read_points.h>
#include <CGAL/property_map.h>
#include <CGAL/Surface_mesh.h>
#include <CGAL/Shape_detection/Region_growing/Region_growing.h>
#include <CGAL/Shape_detection/Region_growing/Region_growing_on_point_set.h>
#include <CGAL/Polygonal_surface_reconstruction.h>
#ifdef CGAL_USE_SCIP  // defined (or not) by CMake scripts, do not define by hand
#include <CGAL/SCIP_mixed_integer_program_traits.h>
typedef CGAL::SCIP_mixed_integer_program_traits<double> MIP_Solver;
#elif defined(CGAL_USE_GLPK)  // defined (or not) by CMake scripts, do not define by hand
#include <CGAL/GLPK_mixed_integer_program_traits.h>
typedef CGAL::GLPK_mixed_integer_program_traits<double>        MIP_Solver;
#endif
#if defined(CGAL_USE_GLPK) || defined(CGAL_USE_SCIP)
#include <fstream>
#include <CGAL/Timer.h>
typedef CGAL::Exact_predicates_inexact_constructions_kernel        Kernel;
typedef Kernel::FT       FT;
typedef Kernel::Point_3         Point;
typedef Kernel::Vector_3 Vector;
// Point with normal, and plane index.
typedef boost::tuple<Point, Vector, int> PNI;
typedef std::vector<PNI> Point_vector;
typedef CGAL::Nth_of_tuple_property_map<0, PNI>        Point_map;
typedef CGAL::Nth_of_tuple_property_map<1, PNI>        Normal_map;
typedef CGAL::Nth_of_tuple_property_map<2, PNI>        Plane_index_map;
typedef CGAL::Shape_detection::Point_set::
Sphere_neighbor_query<Kernel, Point_vector, Point_map> Neighbor_query;
typedef CGAL::Shape_detection::Point_set::
Least_squares_plane_fit_region<Kernel, Point_vector, Point_map, Normal_map> Region_type;
typedef CGAL::Shape_detection::
Region_growing<Point_vector, Neighbor_query, Region_type> Region_growing;
typedef CGAL::Surface_mesh<Point>        Surface_mesh;
typedef        CGAL::Polygonal_surface_reconstruction<Kernel> Polygonal_surface_reconstruction;
class Index_map {
public:
  using key_type = std::size_t;
  using value_type = int;
  using reference = value_type;
  using category = boost::readable_property_map_tag;
  Index_map() { }
  template<typename PointRange>
  Index_map(const PointRange& points,
            const std::vector< std::vector<std::size_t> >& regions)
    : m_indices(new std::vector<int>(points.size(), -1))
  {
    for (std::size_t i = 0; i < regions.size(); ++i)
      for (const std::size_t idx : regions[i])
        (*m_indices)[idx] = static_cast<int>(i);
  }
  inline friend value_type get(const Index_map& index_map,
                               const key_type key)
  {
    const auto& indices = *(index_map.m_indices);
    return indices[key];
  }
private:
  std::shared_ptr< std::vector<int> > m_indices;
};
/*
* This example first extracts planes from the input point cloud
* (using region growing) and then reconstructs
* the surface model from the planes.
*/
int main()
{
  Point_vector points;
  // Load point set from a file.
  const std::string input_file(CGAL::data_file_path("points_3/cube.pwn"));
    std::ifstream input_stream(input_file.c_str());
  if (input_stream.fail()) {
    std::cerr << "Failed open file \'" << input_file << "\'" << std::endl;
    return EXIT_FAILURE;
  }
  input_stream.close();
  std::cout << "Loading point cloud: " << input_file << "...";
  CGAL::Timer t;
  t.start();
  if (!CGAL::IO::read_points(input_file.c_str(), std::back_inserter(points),
                             CGAL::parameters::point_map(Point_map()).normal_map(Normal_map()))) {
    std::cerr << "Error: cannot read file " << input_file << std::endl;
    return EXIT_FAILURE;
  }
  else
    std::cout << " Done. " << points.size() << " points. Time: "
    << t.time() << " sec." << std::endl;
  // Shape detection.
  // Default parameter values for the data file cube.pwn.
  const FT          search_sphere_radius  = FT(2) / FT(100);
  const FT          max_distance_to_plane = FT(2) / FT(1000);
  const FT          max_accepted_angle    = FT(25);
  const std::size_t min_region_size       = 200;
  // Create instances of the classes Neighbor_query and Region_type.
  Neighbor_query neighbor_query(
    points,
    search_sphere_radius);
  Region_type region_type(
    points,
    max_distance_to_plane, max_accepted_angle, min_region_size);
  // Create an instance of the region growing class.
  Region_growing region_growing(
    points, neighbor_query, region_type);
  std::cout << "Extracting planes...";
  std::vector< std::vector<std::size_t> > regions;
  t.reset();
  region_growing.detect(std::back_inserter(regions));
  std::cout << " Done. " << regions.size() << " planes extracted. Time: "
  << t.time() << " sec." << std::endl;
  // Stores the plane index of each point as the third element of the tuple.
  Index_map index_map(points, regions);
  for (std::size_t i = 0; i < points.size(); ++i) {
    // Uses the get function from the property map that accesses the 3rd element of the tuple.
    const int plane_index = get(index_map, i);
    points[i].get<2>() = plane_index;
  }
  // Reconstruction.
  std::cout << "Generating candidate faces...";
  t.reset();
  Polygonal_surface_reconstruction algo(
    points,
    Point_map(),
    Normal_map(),
    Plane_index_map()
  );
  std::cout << " Done. Time: " << t.time() << " sec." << std::endl;
  Surface_mesh model;
  std::cout << "Reconstructing...";
  t.reset();
  if (!algo.reconstruct<MIP_Solver>(model)) {
    std::cerr << "Failed: " << algo.error_message() << std::endl;
    return EXIT_FAILURE;
  }
  std::cout << " Done. Time: " << t.time() << " sec." << std::endl;
  std::cout << "Saving...";
  t.reset();
  const std::string& output_file("data/cube_result.off");
  if (CGAL::IO::write_OFF(output_file, model))
    std::cout << " Done. Saved to " << output_file << ". Time: " << t.time() << " sec." << std::endl;
  else {
    std::cerr << " Failed saving file." << std::endl;
    return EXIT_FAILURE;
  }
  return EXIT_SUCCESS;
}
#else
int main(int, char**)
{
    std::cerr << "This test requires either GLPK or SCIP.\n";
    return EXIT_SUCCESS;
}
#endif  // defined(CGAL_USE_GLPK) || defined(CGAL_USE_SCIP)

3.3 模型复杂度控制

除了通过鼓励大平面区域来有利于干净紧凑的重建结果外，模型复杂性项还提供了对模型细节的控制，即增加该项的影响会导致重建的 3D 模型中的细节更少。通过逐渐增加模型复杂度项的影响来重建。每个子图下的值是相应优化中使用的权重。

在这里插入图片描述

以下示例演示如何通过调整模型复杂性项的权重来控制模型复杂性。此示例还演示如何重用候选生成步骤中的中间结果。

#include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
#include <CGAL/IO/read_points.h>
#include <CGAL/property_map.h>
#include <CGAL/Surface_mesh.h>
#include <CGAL/Polygonal_surface_reconstruction.h>
#ifdef CGAL_USE_SCIP // defined (or not) by CMake scripts, do not define by hand
#include <CGAL/SCIP_mixed_integer_program_traits.h>
typedef CGAL::SCIP_mixed_integer_program_traits<double>                        MIP_Solver;
#elif defined(CGAL_USE_GLPK) // defined (or not) by CMake scripts, do not define by hand
#include <CGAL/GLPK_mixed_integer_program_traits.h>
typedef CGAL::GLPK_mixed_integer_program_traits<double>                        MIP_Solver;
#endif
#if defined(CGAL_USE_GLPK) || defined(CGAL_USE_SCIP)
#include <CGAL/Timer.h>
#include <fstream>
typedef CGAL::Exact_predicates_inexact_constructions_kernel                Kernel;
typedef Kernel::Point_3                                                                                        Point;
typedef Kernel::Vector_3                                                                                Vector;
typedef        CGAL::Polygonal_surface_reconstruction<Kernel>                        Polygonal_surface_reconstruction;
typedef CGAL::Surface_mesh<Point>                                                                Surface_mesh;
// Point with normal, and plane index
typedef boost::tuple<Point, Vector, int>                                                PNI;
typedef CGAL::Nth_of_tuple_property_map<0, PNI>                                        Point_map;
typedef CGAL::Nth_of_tuple_property_map<1, PNI>                                        Normal_map;
typedef CGAL::Nth_of_tuple_property_map<2, PNI>                                        Plane_index_map;
/*
* The following example shows how to control the model complexity by
* increasing the influence of the model complexity term.
* In this example, the intermediate results from plane extraction and
* candidate generation are cached and reused.
*/
int main()
{
    const std::string& input_file(CGAL::data_file_path("points_3/building.ply"));
  std::ifstream input_stream(input_file.c_str());
  std::vector<PNI> points; // store points
  std::cout << "Loading point cloud: " << input_file << "...";
  CGAL::Timer t;
  t.start();
  if (!CGAL::IO::read_PLY_with_properties(input_stream,
                                          std::back_inserter(points),
                                          CGAL::make_ply_point_reader(Point_map()),
                                          CGAL::make_ply_normal_reader(Normal_map()),
                                          std::make_pair(Plane_index_map(), CGAL::PLY_property<int>("segment_index"))))
  {
    std::cerr << "Error: cannot read file " << input_file << std::endl;
    return EXIT_FAILURE;
  }
  else
    std::cout << " Done. " << points.size() << " points. Time: " << t.time() << " sec." << std::endl;
  std::cout << "Generating candidate faces...";
  t.reset();
  Polygonal_surface_reconstruction algo(
    points,
    Point_map(),
    Normal_map(),
    Plane_index_map()
  );
  std::cout << " Done. Time: " << t.time() << " sec." << std::endl;
  // Reconstruction with complexity control
  // Model 1: more detail
  Surface_mesh model;
  std::cout << "Reconstructing with complexity 0.05...";
  t.reset();
  if (!algo.reconstruct<MIP_Solver>(model, 0.8, 0.15, 0.05)) {
    std::cerr << " Failed: " << algo.error_message() << std::endl;
    return EXIT_FAILURE;
  }
  else {
       const std::string& output_file = "data/building_result-0.05.off";
       if (CGAL::IO::write_OFF(output_file, model)) {
           std::cout << " Done. Saved to " << output_file << ". Time: " << t.time() << " sec." << std::endl;
       }
       else {
           std::cerr << " Failed saving file." << std::endl;
           return EXIT_FAILURE;
       }
  }
  // Model 2: a little less detail
  std::cout << "Reconstructing with complexity 0.5...";
  t.reset();
  if (!algo.reconstruct<MIP_Solver>(model, 0.3, 0.2, 0.5)) {
    std::cerr << " Failed: " << algo.error_message() << std::endl;
    return EXIT_FAILURE;
  }
  else {
       const std::string& output_file = "data/building_result-0.5.off";
       if (CGAL::IO::write_OFF(output_file, model))
           std::cout << " Done. Saved to " << output_file << ". Time: " << t.time() << " sec." << std::endl;
       else {
           std::cerr << " Failed saving file." << std::endl;
           return EXIT_FAILURE;
       }
  }
  // Model 3: more less detail
  std::cout << "Reconstructing with complexity 0.7...";
  t.reset();
  if (!algo.reconstruct<MIP_Solver>(model, 0.2, 0.1, 0.7)) {
    std::cerr << " Failed: " << algo.error_message() << std::endl;
    return EXIT_FAILURE;
  }
  else {
    const std::string& output_file = "data/building_result-0.7.off";
    if (CGAL::IO::write_OFF(output_file, model)){
      std::cout << " Done. Saved to " << output_file << ". Time: " << t.time() << " sec." << std::endl;
                }
    else {
      std::cerr << " Failed saving file." << std::endl;
      return EXIT_FAILURE;
    }
  }
  return EXIT_SUCCESS;
}
#else
int main(int, char**)
{
    std::cerr << "This test requires either GLPK or SCIP.\n";
    return EXIT_SUCCESS;
}
#endif  // defined(CGAL_USE_GLPK) || defined(CGAL_USE_SCIP)