从点云重建表面 Surface Reconstruction from Point Clouds

https://doc.cgal.org/latest/Manual/tuto_reconstruction.html

从点云重建表面是几何处理的核心主题 [3]。这是一个不适定(ill-posed)的问题（不同时满足解存在、唯一、稳定这3个条件）：有无数个表面近似于单个点云，而点云本身并不定义表面。因此，用户必须定义额外的假设和约束，并且可以通过许多不同的方式实现重建。本教程提供了有关如何使用 CGAL 的不同算法来有效执行表面重建的指导。

1 我应该使用哪种算法？

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

• 泊松曲面重建(PSR,Poisson Surface Reconstruction)
• 推进前表面重建(AF,Advancing Front Surface Reconstruction)
• 尺度空间表面重建(SS,Scale Space Surface Reconstruction)

因为重建是一个不适定(ill-posed)的问题，所以必须通过先验知识来规范它。先验的差异会导致不同的算法，例如，泊松始终生成闭合形状（约束体积）并需要法线，但不插值输入点（输出表面不完全通过输入点）。下表列出了输入和输出的不同属性，以帮助用户选择最适合每个问题的方法：

	Psr	Af	ss
法线Are normals required?	Yes	No	No
噪声Is noise handled?	Yes	By preprocessing	Yes
采样Is variable sampling handled?	Yes	Yes	By preprocessing
点在表面Are input points exactly on the surface?	No	Yes	Yes
闭合Is the output always closed?	Yes	No	No
平滑Is the output always smooth?	Yes	No	No
流形？Is the output always manifold?	Yes	Yes	Optional
有向？Is the output always orientable?	Yes	Yes	Optional

有关这些不同方法的更多信息，请参阅其各自的手册页和重建部分。

2 Pipeline Overview

https://doc.cgal.org/latest/Manual/tuto_reconstruction.html#fig__TutorialsReconstructionFigPipeline

输入(Input)是 xyz,off,ply,las/laz类型文件，经过Outlier removal，简化，平滑，法线估计等预处理(Preprocessing)，进入重建(Reconstruction)，重建得到的结果进行后处理(Postprocessing)，最后将结果写入(Output)文件(off, ply)

3 读取输入

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

• read_XYZ()
• read_OFF()
• read_PLY()
• read_PLY_with_properties() to read additional PLY properties
• read_LAS()
• read_LAS_with_properties() to read additional LAS properties

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 Outlier removal 异常值删除

某些采集技术生成的点远离表面。这些点通常被称为"异常值(outliers)"，与重建无关。在异常值缠绕的点云上使用 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 freeed by the user.
  points.collect_garbage();

4.2 简化

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

CGAL 提供了几种简化算法。除了减小输入点云的大小,减少计算时间外，其中一些可以帮助使输入更加均匀。函数grid_simplify_point_set()定义了用户指定大小的网格，并为每个占用的像元(occupied cell)保留一个点。

// 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 平滑

虽然通过"泊松（Poisson）"或"刻度空间(Scale space)"进行的重建会在内部处理噪声，但人们可能希望对平滑步骤进行更严格的控制。例如，稍微嘈杂的点云可以从一些可靠的平滑算法中受益，并通过提供相关属性（oriented mesh with boundaries,具有边界的定向网格）的"前进前沿(Advancing front)"进行重建。

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

• jet_smooth_point_set()
• bilateral_smooth_point_set()

这些函数直接修改容器：

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

4.4 法线估计和方向

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

• pca_estimate_normals()
• jet_estimate_normals()

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

• 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 前进前线

前进前沿是一种基于 Delaunay 的方法，它插值了输入点的子集。它生成三个点索引，用于描述重建的三角形面：它使用优先级队列，根据大小准则（有利于小刻面）和角度准则（有利于平滑度），按顺序选择最有可能是曲面一部分的 Delaunay 小平面。它的主要优点是生成具有边界的定向流形曲面：与泊松相反，它不需要法线，也不绑定到重建闭合形状。但是，如果点云有噪声，则需要进行预处理。

"前进前沿"包提供了几种构造函数的方法。下面是一个简单的示例：

 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 缩放空间

尺度空间重建旨在产生一个表面，该表面可以插值输入点（插值），同时为噪声提供一定的鲁棒性。更具体地说，它首先将平滑滤波器（如喷射平滑）应用于输入点集以产生比例空间;然后，对最平滑的比例进行网格划分（例如使用前进前网格划分器）;最后，平滑点之间产生的连通性将传播到原始原始输入点集。如果输入点云很嘈杂，但用户仍然希望表面完全通过这些点，则此方法是正确的选择。

    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 << 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 ();