CGAL笔记之形状重建——尺度-空间表面重建

---

尺度-空间表面重建

CGAL 包实现了一种曲面重建方法，该方法将无序点集作为输入，并计算插入点集的三角化曲面网格。我们假设输入点是从物体表面采样的。该方法还可以处理从对象内部采样的点集，尽管我们不能对输出提供保证。这种方法可以处理相当数量的噪声和异常值。点集可能会大大低于遮挡区域中的对象采样，尽管不会重建表面来填充这些区域。左：合成结数据集上的 5760 个点。右：重建的表面网格。三角化曲面网格是通过首先计算粗比例下的点集，然后构建该比例下的点集的网格，最后将网格的点恢复到其原始比例来生成的。

1 比例空间

从点样本进行表面重建的问题定义不明确。由于这些点仅提供表面的样本，因此我们无法确切确定表面在点之间的行为。通过进一步对表面进行采样，我们可以获得有关表面的更多信息，但即便如此，问题仍然在较小的范围内。

对于真实世界的数据集，我们还必须应对噪声和异常值，点位置中的小误差以及未分别对表面进行采样的点。这些通常是由测量过程中的故障引起的。除非我们真正知道采样的表面，否则我们无法确定一个独特的点是异常值，是由噪声引起的，还是对表面的一个小特征进行采样。

为了克服这些障碍，我们使用尺度空间对不同尺度上的点集进行建模。就我们的目的而言，尺度是点集的抽象级别。在较低、更精细的比例下，点集描述对象的较小特征。相反，在更高、更粗糙的比例下，点集仅描述较大的特征。刻度空间描述在动态刻度下设置的点。它是一个四维空间，尺度为第四维。增加比例类似于平滑点集描述的下层表面。下图 显示了在不同比例下设置的示例点。对这种尺度空间进行建模的总体效果是，曲面重建问题在更高尺度上得到了更明确的定义。黑点集包含一个平滑区域、一个具有大要素的区域以及一个具有小要素或噪点的区域。蓝点以更高的比例显示点云。

比例的微小变化会导致点集的微小变化，这取决于局部形状。随着比例的增加，几何图形会根据连接的组件接近多个点。给定一组固定比例的点，可以明确地确定在较高比例下设置的点，但对于以较低比例设置的点，情况并非如此，如图所示。黑点和十字显示两个不同的点集。蓝色框以较高的比例显示两个点集。请注意，仅给定蓝色框，我们无法确定其在较低比例下的点集。

尺度-空间表面重建方法分两个阶段进行。首先，从点集构造一个尺度空间，并计算该空间在更高尺度\（s\）下的交集。然后，计算在比例尺 s 设置的点的三角化曲面网格。网格中的相邻三角形具有一致的方向，如使用三角形有序顶点上的右手规则表示的那样。

通常，需要对原始点集进行内插的重建。这可以通过将网格的点恢复到其原始比例（即其原始位置）来实现。[图 显示了从原始比例和更高比例下使用还原网格设置的点重建的示例。

从以原始比例（左）和较高比例（右）设置的点构建的曲面。右图中的表面点将恢复为原始比例。虚线路径显示朝内的壳体，其中它不与朝外的壳体重叠。请注意比例空间重建如何假定对象样本小特征右侧的点。

平滑运算符和网格重建都假定彼此靠近的点属于对象的同一部分。这通常表现为具有固定大小的球的概念，即邻域半径。如果这样的球包含多个点，则这些点彼此靠近，并且在增加比例的同时会相互影响。如果这样的球是空的，它就位于物体之外。请注意，外部基于没有点的区域，而不是基于体积是否被表面包围。

可以通过统计分析来估计邻域大小。我们使用 k-D 树来估计到第 n 个最近邻的平均距离，并使用该距离作为分辨率的近似值。

在初始比例处设置的点等效于输入点集。理论上，尺度可以连续变化，但实际上，出于效率原因，尺度在离散迭代中增加。通过使用平滑运算符变换每个点，比例将增加一次迭代。CGAL 提供了两个用于缩放空间重建的平滑运算符：

• CGAL::Scale_space_reconstruction_3::Weighted_PCA_smoother（默认值）使用局部邻域的密度加权主成分分析 （PCA）。如果点集是从任何不需要的变形和采样噪声小于邻域大小的表面采样的，则在增加比例的几次迭代后，比例足够粗糙，可以进行网格重建。
• CGAL::Scale_space_reconstruction_3::Jet_smoother使用将点投影到局部平滑 （Jet） 表面的函数CGAL::jet_smooth_point_set()。这种平滑不太激进，如果点集不是很嘈杂，并且精度更高，则应使用这种平滑。

用户可以按照概念CGAL::Scale_space_reconstruction_3::Smoother定义自己的平滑运算符。

网格划分是通过插值平滑的点云并将连通性传播回原始点云来实现的。CGAL 提供了两种用于尺度空间重构的网格划分算子：

• CGAL::Scale_space_reconstruction_3::Alpha_shape_mesher（默认值）使用过滤算法CGAL::Alpha_shape_3生成一个或多个“外壳”。此方法适用于闭合形状（无边界）。对打开的形状进行采样的点集会导致曲面重叠。它需要一个与数据分辨率相关的固定邻域大小参数。此参数表示一个区域，如果它以表面为中心，我们可以假设它至少包含一个点。
• CGAL::Scale_space_reconstruction_3::Advancing_front_mesher使用该算法CGAL::Advancing_front_surface_reconstruction生成定向 2 流形曲面。此方法处理具有边界的形状，并让用户控制所使用的最大分面。

该方法提供对中间结果的访问，用户可以调整这些结果以更好地满足他们的需求。（中间）结果是对分辨率、比例和曲面三角形的最终集合的估计值。

Jet_smoother和 Advancing_front_mesher是依赖于其他 CGAL 包的方法，与原始缩放空间算法不对应。它们可以用作具有低噪声的点云情况的替代运算符，这些点云对具有边界的表面进行采样。

1.1 参数设置

尺度-空间表面重构方法有两个主要的全局参数：邻域球的半径和增加尺度的迭代次数。如果不知道合适的邻域半径，则可以使用另外两个参数进行估计：平均相邻要素数和样本数。我们已经根据经验确定了每个参数的值，这些值在广泛的数据集上效果很好。但是，我们建议针对每种类型的数据集仔细微调这些参数。

尺度空间重建在本地运行。增加比例时，局部邻域内的点会相互影响。同样，空间中彼此靠近的点更有可能在重建的表面上彼此靠近。因此，重要的是要很好地指示哪些点彼此靠近。此邻近性在邻域球半径参数中表示。

邻域球半径与点密度和物体厚度有关。在理想情况下，半径的选择使得任何半径以物体表面的任何点为中心的球B 将包含几个点，这些点应通过表面中的局部补丁连接以进行重建。同时，B 不应包含任何不应属于此局部修补程序的点，例如对象另一侧的点。

可以根据点集自动估计邻域半径。在这种情况下，估计半径时平均包含给定数量的点。应选择点数，以保留上述邻域半径上的条件。通常，平均 30 个邻居会产生良好的结果。如果物体的某些部分具有非常低的点密度，则此数字可能会增加，如果物体的薄特征应该更好地重建，则可能会减少。

半径估计基于多个采样点。此样本数与点集覆盖对象的频率有关。如果对象被非常不规则地覆盖，则可能需要较大的样本量。但是，更多的样本可能需要更多的计算时间。我们通常使用 200 个样本来估计邻域半径。

尺度-空间方法试图通过将尺度增加到适当的粗尺度来减少表面重建问题的不适。在此粗比例下，点集对平滑表面进行采样，该表面更易于重建。比例增加迭代次数决定了此曲面的平滑程度。

迭代次数与点集中的噪点、表面特征的锐度和对象的厚度有关。具有大量噪点的点集和具有尖锐或小特征的对象将需要更多迭代。同时，处理过多的迭代会使体积退化为平面。这些简并可能导致重建的表面连接物体相对两侧的点。通常，增加比例的 4 次迭代是合适的。

1.2 输出表面

以比例 s 构建的表面网格是非自相交的。任何两个三角形的内部不能成对地在一条线段中相交。但是，表面不需要是 2 流形。一条边可能入射到两个以上的三角形，如果三角形两侧的大区域没有点，三角形可能会完全重叠。请注意，我们将具有相反方向的重叠三角形计为单独的三角形。在许多情况下，当点对对象的表面进行采样时，计算的表面将同时包含朝外和类似的表面，它们之间有一个薄体积。

曲面网格不会有仅入射到一个或多个三角形的边，这些边的环，并且三角形的方向都远离点集。如果点集有孔，则曲面网格很可能包含与该孔接触的相反方向的重叠三角形。

可以应用额外的处理来强制输出曲面为 2 流形：在重建之前找到小的平面体积，并且表面被迫仅使用体积的一侧，同时与相邻的奇异刻面保持一致。仅使用每个奇异刻面的一侧。此外，非流形边和顶点作为后处理被移除。为使曲面 2 歧管而丢弃的刻面被存储并可访问。

如果对象未密集采样或具有断开的组件，则重建的表面可能具有多个断开连接的组件。曲面要么是三角形的无序集合，要么是按壳排序的相同集合。壳是局部朝向曲面同一侧的连接三角形的集合。

当恢复到原始比例时，我们不能保证表面是有效的嵌入，因为该表面的三角形可能成对地在其内部相交。它也可能具有边界边，尽管这些边界边将始终入射到两个方向相反的曲面三角形。但是，当对迭代次数和邻域大小使用适当的参数设置时，表面通常不会自相交。适当的参数设置取决于点集的几何形状，通常需要根据数据集进行微调

2 软件设计

主类Scale_space_surface_reconstruction_3、Weighted_PCA_smoother和Alpha_shape_mesher包含估计邻域大小、计算尺度空间和增加尺度以及从当前尺度上设置的点重建表面的所有功能。

邻域大小是使用Orthogonal_k_neighbor_search 估算的。点集通常存储在Orthogonal_k_neighbor_search::Tree 估计邻域大小时，将搜索此树以查找最近的邻域。

比例空间是在点的原始比例下构建的。使用加权 PCA 过程计算增加比例的迭代。正如Digne等人所描述的那样。 ，与类似的方法不同，此过程不会导致不希望的聚类效应。默认情况下，高效的eigen库用于此过程（如果可用）。否则，将使用内部回退。也可以为概念提供自己的模型。加权 PCA 过程按点在本地执行，因此可以通过Diagonalize_traits``DiagonalizeTraits``Parallel_tag并行计算执行（需要与英特尔 TBB 链接并传递给重建类）。

网格重建是通过过滤以固定比例设置的点的3D alpha 形状来执行的。此过滤为每个规则面构造一个三角形;每个奇异面产生两个面向相反方向的三角形。

生成的曲面存储为三元组的集合。每个三元组包含指向点集的三个索引。比例迭代不会更改点的顺序。这意味着可以将这些索引应用于以原始比例设置的点，以获得插值表面，称为还原表面的比例。索引也可以应用于在当前比例下设置的点，以获得平滑的表面。与恢复的表面相比，此平滑表面通常会有所收缩。

为了减少内存需求，点集不会以其原始比例存储。相反，在当前比例下设置的存储点是输入的克隆。这意味着原始输入可用于还原表面。

3 案例集和性能

有各种各样的可能的点集，无论是合成的还是测量的。尺度-空间曲面重构方法试图使重构问题不那么不恰当。但是，我们不能保证任意点集的良好结果。某些点集需要仔细设置参数，而对于其他点集，此方法则完全不适合。在这里，我们展示了不同的案例，包括理想样本和不同的问题案例。

3.1 理想条件

尺度-空间表面重建具有宽松的理想条件。表面应合理地定期取样，并且不应包含过于凹陷的角。这些概念以相反的方式影响邻域参数：更规则的采样可能会允许较小的邻域参数，这将允许较小的凹陷。点集可能包含明显小于邻域大小的噪声。即便如此，当使用该方法构建插入原始点集的网格时，任何噪声都会使网格变形。下图显示了理想条件下的重建。 **左：**在理想情况下，点是规则分布的，凹陷是有限的。**右：**重建的表面网格。

3.2 边界

自然处理具有边界的表面。但是，重要的是要注意，在这些情况下，曲面不会显式标识边界并围绕它进行循环。在许多情况下，这将导致边界边被曲面的相邻三角形共享，这些三角形除了方向外是相同的。可以通过识别相邻三角形之间的方向差异来显式显示边界，但此包不提供此功能。图 63.6显示了具有边界的重建表面。**左：**具有边界的表面的点样本。**右：**重建的表面网格。网格划分器Advancing_front_mesher是生成具有边界的定向良好的表面的良好替代方法。

3.3 抽样规律性

如上所述，理想情况下，点集定期对表面进行采样。例如，对于每个点，其六个最近的点位于大致相等的距离上，在该点周围形成一个粗略的六边形。处理实际数据时通常不会发生这种情况。比例空间方法可以重建非常不均匀的采样表面，如图 所示。然而，如前所述，采样规律性和允许凹形成权衡。图 中显示了我们无法同时正确建模凹区域和稀疏区域的情况。**左：**一个非常不规则的采样点集。**右：**重建的表面网格。

**左：**具有紧密凹陷的非常不规则的采样点集。**右：**重建的表面网格。请注意，我们无法正确模拟嘴巴、眼睛和颈部。

3.4 尖锐特征

尺度-空间表面重建方法并非专门用于处理具有尖锐特征的表面。但是，当使用适当的参数设置时，它会充分重建这些功能，如图所示。正确重建尖锐特征的主要前提是它们在特征上包含足够的采样点，并且与邻域参数相比，特征不会太薄。设置在具有尖锐折痕的表面上的点和从不同侧面看到的重建表面网格。请注意，尖锐的特征大多是正确重建的。但是，在以红色圈出的细区域中，要素未正确重建。规则性和凹度之间的权衡也意味着右下角的第三个图中出现了一个洞。

3.5 拓扑

尺度-空间表面重建方法不考虑表面的拓扑结构。这意味着重建的表面可能具有与为获得点集而采样的表面不同的拓扑。如图显示了一个示例点集，该点集是从具有球体拓扑的表面采样的;重建的表面具有圆环的拓扑结构。**左：**具有球体拓扑的表面的不规则采样。**右：**具有圆环拓扑结构的重建表面网格。请注意，不规则采样使得无法使用球体拓扑重建表面，因为该方法使用全局尺度参数。

3.6 点属性

许多点集包含点属性，例如法线或颜色。由于比例空间表面重建构建了原始点之间的连通性，因此任何对表面有意义的点属性都可以传播到重建中。图中显示了为一组带有颜色的点重建的表面。**左：**使用点集的颜色进行顶点着色的重建表面。**右：**显示为线框的表面细节。

4 示例

第一个示例演示如何从文件中读取点集、重建曲面以及对于重建的每个外壳，访问三角形顶点索引。

#include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
#include <CGAL/Scale_space_surface_reconstruction_3.h>
#include <CGAL/IO/read_points.h>
#include <CGAL/Timer.h>
#include <fstream>
#include <iostream>
typedef CGAL::Exact_predicates_inexact_constructions_kernel     Kernel;
typedef Kernel::Point_3                                         Point;
typedef CGAL::Scale_space_surface_reconstruction_3<Kernel>      Reconstruction;
typedef Reconstruction::Facet_const_iterator                    Facet_iterator;
int main(int argc, char** argv)
{
  // Read the data.
  std::string fname = argc==1?CGAL::data_file_path("points_3/kitten.off"):argv[1];
  std::cerr << "Reading " << std::flush;
  std::vector<Point> points;
  if(!CGAL::IO::read_points(fname, std::back_inserter(points)))
  {
    std::cerr << "Error: cannot read file" << std::endl;
    return EXIT_FAILURE;
  }
  std::cerr << "done: " << points.size() << " points." << std::endl;
  std::cerr << "Reconstruction ";
  CGAL::Timer t;
  t.start();
  // Construct the mesh in a scale space.
  Reconstruction reconstruct (points.begin(), points.end());
  reconstruct.increase_scale(4);
  reconstruct.reconstruct_surface();
  std::cerr << "done in " << t.time() << " sec." << std::endl;
  t.reset();
  std::ofstream out ("out.off");
  out << reconstruct;
  std::cerr << "Writing result in " << t.time() << " sec." << std::endl;
  std::cerr << "Done." << std::endl;
  return EXIT_SUCCESS;
}

重建的不同阶段也可以单独进行。在事先不知道正确的参数设置的情况下，这可能会有所帮助。在图中，右图显示了比例递增后的表面。在检查重建的表面后，我们可能会决定执行更多的迭代可以更好地对左上角的部分进行建模。

下一个示例通过进一步推进比例空间来执行两次重建，对设置的点进行不同的平滑处理。它说明了如何以交互方式使用代码。

#include <CGAL/Scale_space_surface_reconstruction_3.h>
#include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
#include <CGAL/IO/read_points.h>
#include <algorithm>
#include <fstream>
#include <iostream>
typedef CGAL::Exact_predicates_inexact_constructions_kernel     Kernel;
typedef CGAL::Scale_space_surface_reconstruction_3< Kernel >                Reconstruction;
typedef CGAL::Scale_space_reconstruction_3::Weighted_PCA_smoother< Kernel > Smoother;
typedef CGAL::Scale_space_reconstruction_3::Alpha_shape_mesher< Kernel >    Mesher;
typedef Reconstruction::Point                                   Point;
typedef Reconstruction::Facet_const_iterator                   Facet_iterator;
// function for writing the reconstruction output in the off format
void dump_reconstruction(const Reconstruction& reconstruct, std::string name)
{
  std::ofstream output(name.c_str());
  output << "OFF " << reconstruct.number_of_points() << " "
         << reconstruct.number_of_facets() << " 0\n";
  std::copy(reconstruct.points_begin(),
            reconstruct.points_end(),
            std::ostream_iterator<Point>(output,"\n"));
  for( Facet_iterator it = reconstruct.facets_begin(); it != reconstruct.facets_end(); ++it )
      output << "3 " << *it << std::endl;
}
int main(int argc, char* argv[])
{
    // Read the data.
    std::string fname = argc==1?CGAL::data_file_path("points_3/kitten.off"):argv[1];
    std::cout << "Reading " << std::flush;
    std::vector<Point> points;
    if(!CGAL::IO::read_points(fname, std::back_inserter(points)))
    {
      std::cerr << "Error: cannot read file" << std::endl;
      return EXIT_FAILURE;
    }
    std::cout << "done: " << points.size() << " points." << std::endl;
    // Construct the reconstruction
    Reconstruction reconstruct;
    // Add the points.
    reconstruct.insert( points.begin(), points.end() );
    // Two passes
    for (std::size_t i = 0; i < 2; ++ i)
    {
      // Construct the smoother with parameters for
      // the neighborhood squared radius estimation.
      Smoother smoother( 10, 100 );
      // Advance the scale-space several steps.
      // This automatically estimates the scale-space.
      reconstruct.increase_scale( 2, smoother );
      // Reconstruct the surface from the current scale-space.
      std::cout << "Neighborhood squared radius is "
                << smoother.squared_radius() << std::endl;
      Mesher mesher (smoother.squared_radius());
      reconstruct.reconstruct_surface(mesher);
      if (i == 0)
      {
        std::cout << "First reconstruction done." << std::endl;
        // Write the reconstruction.
        dump_reconstruction(reconstruct, "reconstruction1.off");
      }
      else
      {
        std::cout << "Second reconstruction done." << std::endl;
        // Write the reconstruction.
        dump_reconstruction(reconstruct, "reconstruction2.off");
      }
    }
    std::cout << "Reconstructions are ready to be examinated in your favorite viewer" << std::endl;
    return EXIT_SUCCESS;
}

此示例显示如何强制输出曲面为 2 流形。在这种情况下，可以丢弃某些分面以删除非流形单纯形。它们以无序的方式存储在垃圾箱中。通过迭代此容器来访问它们。

#include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
#include <CGAL/Scale_space_surface_reconstruction_3.h>
#include <CGAL/IO/read_points.h>
#include <CGAL/Timer.h>
#include <fstream>
#include <iostream>
typedef CGAL::Exact_predicates_inexact_constructions_kernel     Kernel;
typedef CGAL::Scale_space_surface_reconstruction_3< Kernel >                Reconstruction;
typedef CGAL::Scale_space_reconstruction_3::Weighted_PCA_smoother< Kernel > Smoother;
typedef CGAL::Scale_space_reconstruction_3::Alpha_shape_mesher< Kernel >    Mesher;
typedef Reconstruction::Point                                   Point;
typedef Reconstruction::Facet_const_iterator                    Facet_iterator;
typedef Mesher::Facet_const_iterator                            Mesher_iterator;
typedef CGAL::Timer Timer;
int main(int argc, char* argv[]) {
  // Read the dat
  std::string fname = argc==1?CGAL::data_file_path("points_3/kitten.off"):argv[1];
  std::cerr << "Reading " << std::flush;
  std::vector<Point> points;
  if(!CGAL::IO::read_points(fname, std::back_inserter(points)))
  {
    std::cerr << "Error: cannot read file" << std::endl;
    return EXIT_FAILURE;
  }
  std::cerr << "done: " << points.size() << " points." << std::endl;
  Timer t;
  t.start();
  // Construct the mesh in a scale space.
  Reconstruction reconstruct(points.begin(), points.end() );
  Smoother smoother(10, 200 );
  reconstruct.increase_scale(4, smoother);
  Mesher mesher(smoother.squared_radius(),
                false, // Do not separate shells
                true // Force manifold output
               );
  reconstruct.reconstruct_surface(mesher);
  std::cerr << "Reconstruction done in " << t.time() << " sec." << std::endl;
  t.reset();
  std::ofstream out("out.off");
  // Write the reconstruction.
  for(Facet_iterator it = reconstruct.facets_begin(); it != reconstruct.facets_end(); ++it )
    out << "3 "<< *it << '\n'; // We write a '3' in front so that it can be assembled into an OFF file
  std::cerr << "Writing result in " << t.time() << " sec." << std::endl;
  out.close();
  t.reset();
  std::ofstream garbage("garbage.off");
  // Write facets that were removed to force manifold output
  for(Mesher_iterator it = mesher.garbage_begin(); it != mesher.garbage_end(); ++it )
    garbage << "3 "<< *it << '\n'; // We write a '3' in front so that it can be assembled into an OFF file
  std::cerr << "Writing garbage facets in " << t.time() << " sec." << std::endl;
  std::cerr << "Done." << std::endl;
  return EXIT_SUCCESS;
}

最后一个示例演示如何使用替代运算符Jet_smoother和Advancing_front_mesher。

#include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
#include <CGAL/Point_set_3.h>
#include <CGAL/Scale_space_surface_reconstruction_3.h>
#include <CGAL/Scale_space_reconstruction_3/Advancing_front_mesher.h>
#include <CGAL/Scale_space_reconstruction_3/Jet_smoother.h>
#include <CGAL/Timer.h>
#include <fstream>
#include <iostream>
typedef CGAL::Exact_predicates_inexact_constructions_kernel     Kernel;
typedef CGAL::Scale_space_surface_reconstruction_3<Kernel>                    Reconstruction;
typedef CGAL::Scale_space_reconstruction_3::Advancing_front_mesher<Kernel>    Mesher;
typedef CGAL::Scale_space_reconstruction_3::Jet_smoother<Kernel>              Smoother;
typedef Kernel::Point_3 Point;
typedef CGAL::Point_set_3<Point> Point_set;
typedef Reconstruction::Facet_const_iterator                   Facet_iterator;
int main(int argc, char** argv)
{
  // Read the data.
  std::string fname = argc==1?CGAL::data_file_path("points_3/kitten.off"):argv[1];
  std::cerr << "Reading " << std::flush;
  Point_set points;
  if(!CGAL::IO::read_point_set(fname, points))
  {
    std::cerr << "Error: cannot read file" << std::endl;
    return EXIT_FAILURE;
  }
  std::cerr << "done: " << points.size() << " points." << std::endl;
  std::cerr << "Reconstruction ";
  CGAL::Timer t;
  t.start();
  // Construct the mesh in a scale space.
  Reconstruction reconstruct (points.points().begin(), points.points().end());
  reconstruct.increase_scale<Smoother> (4);
  reconstruct.reconstruct_surface (Mesher (0.5));
  std::cerr << "done in " << t.time() << " sec." << std::endl;
  t.reset();
  std::ofstream out ("out.off");
  out << "OFF" << std::endl << points.size() << " " << reconstruct.number_of_facets() << " 0" << std::endl;
  for (Point_set::iterator it = points.begin(); it != points.end(); ++ it)
    out << points.point(*it) << std::endl;
  for (Reconstruction::Facet_iterator it = reconstruct.facets_begin();
       it != reconstruct.facets_end(); ++ it)
    out << "3 " << (*it)[0] << " " << (*it)[1] << " " << (*it)[2] << std::endl;
  std::cerr << "Writing result in " << t.time() << " sec." << std::endl;
  std::cerr << "Done." << std::endl;
  return EXIT_SUCCESS;
}