2.cgal教程 Surface Reconstruction from Point Clouds

Surface Reconstruction from Point Clouds

原文地址: https://doc.cgal.org/latest/Manual/tuto_reconstruction.html

点云表面重建是几何处理的核心主题[3]。这是一个不适定问题：有无数个表面近似于单个点云，而点云本身并不定义表面。因此，用户必须定义额外的假设和约束，并且可以通过许多不同的方式实现重建。本教程提供有关如何使用 CGAL 的不同算法有效执行表面重建的指导。

1 我应该使用哪种算法？

CGAL 提供三种不同的表面重建算法：

• 泊松曲面重建
• 推进前表面重建
• 尺度空间表面重建

由于重构是一个不适定问题，因此必须通过先验知识对其进行正则化。先验的差异导致不同的算法，并且选择这些方法中的一种或另一种取决于这些先验。例如，泊松总是生成闭合形状（限制体积）并需要法线，但不会对输入点进行插值（输出表面不会完全穿过输入点）。下表列出了输入和输出的不同属性，以帮助用户选择最适合每个问题的方法：

	泊松	前进前沿	尺度空间
是否需要法线？	是的	不	不
噪音处理了吗？	是的	通过预处理	是的
是否处理变量抽样？	是的	是的	通过预处理
输入点是否正好位于表面上？	不	是的	是的
输出始终关闭吗？	是的	不	不
输出总是平滑的吗？	是的	不	不
输出总是多种多样的吗？	是的	是的	选修的
输出总是可定向的吗？	是的	是的	选修的

图 0.1应用于相同输入（完整形状和特写）的重建方法的比较。从左到右：原始点云；泊松；前进前沿；尺度空间。

关于这些不同方法的更多信息可以在它们各自的手册页和重建部分找到。

2 管道概览

本教程旨在更全面地介绍 CGAL 为处理点云和表面重建而提供的可能性。下图显示了使用 CGAL 工具的常见重建步骤的概述（并非详尽）。

我们现在更详细地回顾其中一些步骤。

3 读取输入

CGAL 中的重建算法将容器上的一系列迭代器作为输入，并使用属性映射来访问点（以及需要时的法线）。点通常以纯文本格式存储（表示为“XYZ”格式），其中每个点由换行符分隔，每个坐标由空格分隔。其他可用的格式有“OFF”、“PLY”和“LAS”。CGAL 提供了读取此类格式的函数：

• read_XYZ()
• read_OFF()
• read_PLY()
• read_PLY_with_properties()读取其他 PLY 属性
• read_LAS()
• read_LAS_with_properties()读取其他 LAS 属性

CGAL 还提供了一个专用容器CGAL::Point_set_3来处理具有附加属性（例如法向量）的点集。在这种情况下，属性映射很容易处理，如以下部分所示。该结构还处理流运算符以读取前面描述的任何格式的点集。使用此方法会产生更短的代码，如以下示例所示：

  Point_set points;
  std::string fname = argc==1?CGAL::data_file_path("points_3/kitten.xyz") : argv[1];
  if (argc < 2)
  {
    std::cerr << "Usage: " << argv[0] << " [input.xyz/off/ply/las]" << std::endl;
    std::cerr <<"Running " << argv[0] << " data/kitten.xyz -1\n";
  }
  std::ifstream stream (fname, std::ios_base::binary);
  if (!stream)
  {
    std::cerr << "Error: cannot read file " << fname << std::endl;
    return EXIT_FAILURE;
  }
  stream >> points;
  std::cout << "Read " << points.size () << " point(s)" << std::endl;
  if (points.empty())
    return EXIT_FAILURE;

4 预处理

由于重建算法有一些点云并不总是满足的特定要求，因此可能需要进行一些预处理才能产生最佳结果。

请注意，此预处理步骤是可选的：当输入点云没有缺陷时，可以对其应用重建而无需任何预处理。

4.1 异常值去除

一些采集技术生成远离表面的点。这些点通常称为“异常值”，与重建无关。在充满异常值的点云上使用 CGAL 重建算法会产生过度扭曲的输出，因此强烈建议在执行重建之前过滤这些异常值。

  typename Point_set::iterator rout_it = CGAL::remove_outliers<CGAL::Sequential_tag>
    (points,
     24, // Number of neighbors considered for evaluation
     points.parameters().threshold_percent (5.0)); // Percentage of points to remove
  points.remove(rout_it, points.end());
  std::cout << points.number_of_removed_points()
            << " point(s) are outliers." << std::endl;
  // Applying point set processing algorithm to a CGAL::Point_set_3
  // object does not erase the points from memory but place them in
  // the garbage of the object: memory can be freed by the user.
  points.collect_garbage();

4.2 简化

一些激光扫描仪生成具有广泛可变采样的点。通常，扫描线的采样非常密集，但两条扫描线之间的间隙要大得多，导致采样密度变化较大的点云过大。这种类型的输入点云可能会使用通常仅处理采样密度的微小变化的算法生成不完美的输出。

CGAL 提供了多种简化算法。除了减小输入点云的大小并因此减少计算时间之外，其中一些还可以帮助使输入更加均匀。这是该函数的情况grid_simplify_point_set()，该函数定义用户指定大小的网格并为每个占用的单元保留一个点。

  // Compute average spacing using neighborhood of 6 points
  double spacing = CGAL::compute_average_spacing<CGAL::Sequential_tag> (points, 6);
  // Simplify using a grid of size 2 * average spacing
  typename Point_set::iterator gsim_it = CGAL::grid_simplify_point_set (points, 2. * spacing);
  points.remove(gsim_it, points.end());
  std::cout << points.number_of_removed_points()
            << " point(s) removed after simplification." << std::endl;
  points.collect_garbage();

4.3 平滑

尽管通过“泊松”或“尺度空间”进行重建可以在内部处理噪声，但人们可能希望对平滑步骤进行更严格的控制。例如，稍微有噪声的点云可以受益于一些可靠的平滑算法，并通过提供相关属性（具有边界的定向网格）的“前进前沿”进行重建。

提供了两个函数来平滑具有良好近似值的噪声点云（例如，不会降低曲率）：

• jet_smooth_point_set()
• bilateral_smooth_point_set()

这些函数直接修改容器：

CGAL::jet_smooth_point_set<CGAL::Sequential_tag> (points, 24);

4.4 正常估计和方向

泊松曲面重建需要具有定向法向量的点。要将算法应用于原始点云，必须首先估计法线，例如使用以下两个函数之一：

• pca_estimate_normals()
• jet_estimate_normals()

PCA 速度更快，但在存在高曲率的情况下，jet 更准确。这些函数仅估计法线的*方向，而不是它们的方向（矢量的方向可能局部不一致）。*为了正确定向法线，可以使用以下函数：

• mst_orient_normals()
• scanline_orient_normals()

第一个使用最小生成树在越来越大的邻域中一致地传播法线的方向。在数据具有许多尖锐特征和遮挡的情况下（例如，这在机载激光雷达数据中很常见），第二种算法可能会产生更好的结果：它利用排列成扫描线的点云来估计物体的视线。每个点，从而相应地定向法线。

请注意，如果它们的方向不一致，也可以直接在输入法线上使用它们。

    CGAL::jet_estimate_normals<CGAL::Sequential_tag>
      (points, 24); // Use 24 neighbors
    // Orientation of normals, returns iterator to first unoriented point
    typename Point_set::iterator unoriented_points_begin =
      CGAL::mst_orient_normals(points, 24); // Use 24 neighbors
    points.remove (unoriented_points_begin, points.end());

5 重建

5.1 泊松

泊松重建包括计算一个隐式函数，其梯度与输入法线向量场匹配：该指示函数在推断形状的内部和外部具有相反的符号（因此需要闭合形状）。因此，该方法需要法线并产生平滑的封闭表面。如果希望表面对输入点进行插值，则这是不合适的。相反，如果目标是用光滑表面近似噪声点云，则它表现良好。

    CGAL::Surface_mesh<Point_3> output_mesh;
    CGAL::poisson_surface_reconstruction_delaunay
      (points.begin(), points.end(),
       points.point_map(), points.normal_map(),
       output_mesh, spacing);

5.2 前进阵线

Advancing front 是一种基于 Delaunay 的方法，它对输入点的子集进行插值。它生成描述重建的三角面的点索引的三元组：它使用优先级队列根据尺寸标准（有利于小面）和角度标准（有利于平滑度）。它的主要优点是生成带有边界的定向流形表面：与泊松相反，它不需要法线，也不必重建封闭形状。然而，如果点云有噪声，则需要进行预处理。

Advancing Front包提供了多种构造该函数的方法。这是一个简单的例子：

    typedef std::array<std::size_t, 3> Facet; // Triple of indices
    std::vector<Facet> facets;
    // The function is called using directly the points raw iterators
    CGAL::advancing_front_surface_reconstruction(points.points().begin(),
                                                 points.points().end(),
                                                 std::back_inserter(facets));
    std::cout << facets.size ()
              << " facet(s) generated by reconstruction." << std::endl;

5.3 尺度空间

尺度空间重建旨在生成一个对输入点进行插值（插值）的表面，同时提供一定的噪声鲁棒性。更具体地说，它首先对输入点集应用多次平滑滤波器（例如Jet Smoothing）以产生尺度空间；然后，对最平滑的比例进行网格划分（例如使用 Advancing Front 网格划分器）；最后，平滑点之间的连接性被传播到原始输入点集。如果输入点云有噪声但用户仍然希望表面精确地通过点，则此方法是正确的选择。

    CGAL::Scale_space_surface_reconstruction_3<Kernel> reconstruct
      (points.points().begin(), points.points().end());
    // Smooth using 4 iterations of Jet Smoothing
    reconstruct.increase_scale (4, CGAL::Scale_space_reconstruction_3::Jet_smoother<Kernel>());
    // Mesh with the Advancing Front mesher with a maximum facet length of 0.5
    reconstruct.reconstruct_surface (CGAL::Scale_space_reconstruction_3::Advancing_front_mesher<Kernel>(0.5));

6 输出和后处理

这些方法中的每一种都会生成以不同方式存储的三角形网格。如果此输出网格受到诸如孔或自相交之类的缺陷的阻碍，CGAL 会在包“多边形网格处理”中提供多种算法对其进行后处理（孔填充、重新划分网格等）。

我们在这里不讨论这些函数，因为存在许多后处理可能性，其相关性很大程度上取决于用户对输出网格的期望。

网格（无论是否经过后处理）可以轻松地以 PLY 格式保存（此处使用二进制变量）：

    std::ofstream f ("out_poisson.ply", std::ios_base::binary);
    CGAL::IO::set_binary_mode (f);
    CGAL::IO::write_PLY(f, output_mesh);
    f.close ();

多边形汤还可以通过迭代点和面以 OFF 格式保存：

    std::ofstream f ("out_sp.off");
    f << "OFF" << std::endl << points.size () << " "
      << reconstruct.number_of_facets() << " 0" << std::endl;
    for (Point_set::Index idx : points)
      f << points.point (idx) << std::endl;
    for (const auto& facet : CGAL::make_range (reconstruct.facets_begin(), reconstruct.facets_end()))
      f << "3 "<< facet[0] << " " << facet[1] << " " << facet[2] << std::endl;
    f.close ();

最后，如果多边形汤可以转换为多边形网格，也可以使用流运算符直接以OFF格式保存：

    // copy points for random access
    std::vector<Point_3> vertices;
    vertices.reserve (points.size());
    std::copy (points.points().begin(), points.points().end(), std::back_inserter (vertices));
    CGAL::Surface_mesh<Point_3> output_mesh;
    CGAL::Polygon_mesh_processing::polygon_soup_to_polygon_mesh (vertices, facets, output_mesh);
    std::ofstream f ("out_af.off");
    f << output_mesh;
    f.close ();

7 完整代码示例

本教程中使用的所有代码片段都可以组合起来创建完整的算法管道（前提是使用正确的*包含）。*我们提供了一个完整的代码示例，它实现了本教程中描述的所有步骤。用户可以在运行时使用第二个参数选择重建方法。

#include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
#include <CGAL/Point_set_3.h>
#include <CGAL/Point_set_3/IO.h>
#include <CGAL/remove_outliers.h>
#include <CGAL/grid_simplify_point_set.h>
#include <CGAL/jet_smooth_point_set.h>
#include <CGAL/jet_estimate_normals.h>
#include <CGAL/mst_orient_normals.h>
#include <CGAL/poisson_surface_reconstruction.h>
#include <CGAL/Advancing_front_surface_reconstruction.h>
#include <CGAL/Scale_space_surface_reconstruction_3.h>
#include <CGAL/Scale_space_reconstruction_3/Jet_smoother.h>
#include <CGAL/Scale_space_reconstruction_3/Advancing_front_mesher.h>
#include <CGAL/Surface_mesh.h>
#include <CGAL/Polygon_mesh_processing/polygon_soup_to_polygon_mesh.h>
#include <cstdlib>
#include <vector>
#include <fstream>
// types
typedef CGAL::Exact_predicates_inexact_constructions_kernel Kernel;
typedef Kernel::FT FT;
typedef Kernel::Point_3 Point_3;
typedef Kernel::Vector_3 Vector_3;
typedef Kernel::Sphere_3 Sphere_3;
typedef CGAL::Point_set_3<Point_3, Vector_3> Point_set;
int main(int argc, char*argv[])
{
  Point_set points;
  std::string fname = argc==1?CGAL::data_file_path("points_3/kitten.xyz") : argv[1];
  if (argc < 2)
  {
    std::cerr << "Usage: " << argv[0] << " [input.xyz/off/ply/las]" << std::endl;
    std::cerr <<"Running " << argv[0] << " data/kitten.xyz -1\n";
  }
  std::ifstream stream (fname, std::ios_base::binary);
  if (!stream)
  {
    std::cerr << "Error: cannot read file " << fname << std::endl;
    return EXIT_FAILURE;
  }
  stream >> points;
  std::cout << "Read " << points.size () << " point(s)" << std::endl;
  if (points.empty())
    return EXIT_FAILURE;
  typename Point_set::iterator rout_it = CGAL::remove_outliers<CGAL::Sequential_tag>
    (points,
     24, // Number of neighbors considered for evaluation
     points.parameters().threshold_percent (5.0)); // Percentage of points to remove
  points.remove(rout_it, points.end());
  std::cout << points.number_of_removed_points()
            << " point(s) are outliers." << std::endl;
  // Applying point set processing algorithm to a CGAL::Point_set_3
  // object does not erase the points from memory but place them in
  // the garbage of the object: memory can be freed by the user.
  points.collect_garbage();
  // Compute average spacing using neighborhood of 6 points
  double spacing = CGAL::compute_average_spacing<CGAL::Sequential_tag> (points, 6);
  // Simplify using a grid of size 2 * average spacing
  typename Point_set::iterator gsim_it = CGAL::grid_simplify_point_set (points, 2. * spacing);
  points.remove(gsim_it, points.end());
  std::cout << points.number_of_removed_points()
            << " point(s) removed after simplification." << std::endl;
  points.collect_garbage();
  CGAL::jet_smooth_point_set<CGAL::Sequential_tag> (points, 24);
  int reconstruction_choice
    = argc==1? -1 : (argc < 3 ? 0 : atoi(argv[2]));
  if (reconstruction_choice == 0 || reconstruction_choice==-1) // Poisson
  {
    CGAL::jet_estimate_normals<CGAL::Sequential_tag>
      (points, 24); // Use 24 neighbors
    // Orientation of normals, returns iterator to first unoriented point
    typename Point_set::iterator unoriented_points_begin =
      CGAL::mst_orient_normals(points, 24); // Use 24 neighbors
    points.remove (unoriented_points_begin, points.end());
    CGAL::Surface_mesh<Point_3> output_mesh;
    CGAL::poisson_surface_reconstruction_delaunay
      (points.begin(), points.end(),
       points.point_map(), points.normal_map(),
       output_mesh, spacing);
    std::ofstream f ("out_poisson.ply", std::ios_base::binary);
    CGAL::IO::set_binary_mode (f);
    CGAL::IO::write_PLY(f, output_mesh);
    f.close ();
  }
  if (reconstruction_choice == 1 || reconstruction_choice==-1) // Advancing front
  {
    typedef std::array<std::size_t, 3> Facet; // Triple of indices
    std::vector<Facet> facets;
    // The function is called using directly the points raw iterators
    CGAL::advancing_front_surface_reconstruction(points.points().begin(),
                                                 points.points().end(),
                                                 std::back_inserter(facets));
    std::cout << facets.size ()
              << " facet(s) generated by reconstruction." << std::endl;
    // copy points for random access
    std::vector<Point_3> vertices;
    vertices.reserve (points.size());
    std::copy (points.points().begin(), points.points().end(), std::back_inserter (vertices));
    CGAL::Surface_mesh<Point_3> output_mesh;
    CGAL::Polygon_mesh_processing::polygon_soup_to_polygon_mesh (vertices, facets, output_mesh);
    std::ofstream f ("out_af.off");
    f << output_mesh;
    f.close ();
  }
  if (reconstruction_choice == 2 || reconstruction_choice==-1) // Scale space
  {
    CGAL::Scale_space_surface_reconstruction_3<Kernel> reconstruct
      (points.points().begin(), points.points().end());
    // Smooth using 4 iterations of Jet Smoothing
    reconstruct.increase_scale (4, CGAL::Scale_space_reconstruction_3::Jet_smoother<Kernel>());
    // Mesh with the Advancing Front mesher with a maximum facet length of 0.5
    reconstruct.reconstruct_surface (CGAL::Scale_space_reconstruction_3::Advancing_front_mesher<Kernel>(0.5));
    std::ofstream f ("out_sp.off");
    f << "OFF" << std::endl << points.size () << " "
      << reconstruct.number_of_facets() << " 0" << std::endl;
    for (Point_set::Index idx : points)
      f << points.point (idx) << std::endl;
    for (const auto& facet : CGAL::make_range (reconstruct.facets_begin(), reconstruct.facets_end()))
      f << "3 "<< facet[0] << " " << facet[1] << " " << facet[2] << std::endl;
    f.close ();
  }
  else // Handle error
  {
    std::cerr << "Error: invalid reconstruction id: " << reconstruction_choice << std::endl;
    return EXIT_FAILURE;
  }
  return EXIT_SUCCESS;
}

8 完整的管道图像

下图显示了应用于熊雕像的完整重建流程（由EPFL 计算机图形和几何实验室提供 [5]）。还应用了两种网格处理算法（孔填充和各向同性重新网格化）（有关更多信息，请参阅多边形网格处理一章）。