3D Alpha Wrapping (3D alpha 包裹)

原文地址: https://doc.cgal.org/latest/Alpha_wrap_3/index.html#Chapter_3D_Alpha_wrapping

该组件采用 3D 三角形网格、三角形汤或点集作为输入，并生成严格包含输入的有效三角表面网格（水密、无交集和 二维流形）。 该算法通过从输入的松散边界框开始收缩包裹并细化 3D Delaunay 三角剖分来进行。 两个用户定义的参数（alpha 和 offset）分别可以控制收缩包裹过程可以进入的空腔的最大尺寸，以及最终表面网格与输入的紧密度。 一旦组合起来，这些参数就提供了一种以输入的保真度换取输出的复杂性的方法。

1 介绍

几何建模和处理中的各种任务需要将3D对象表示为有效的表面网格，其中“有效”是指3D对象。指的是水密、无相交、可定向和 二维流形的网格。这种表示提供了内部/外部和测地线邻域的明确定义的概念。

3D 数据通常通过测量和重建来获取，由人类设计，或通过不完善的自动化过程生成。因此，它们可能会表现出各种各样的缺陷，包括间隙、缺失数据、自相交、简并性（例如零体积结构）和非流形特征。

鉴于可能存在的缺陷种类繁多，人们提出了许多方法和数据结构来修复特定缺陷，通常目的是保证修复的 3D 模型中的特定属性。可靠地修复所有类型的缺陷是众所周知的困难，并且通常是一个不适定问题，因为对于给定的带有缺陷的 3D 模型存在许多有效的解决方案。此外，输入模型可能过于复杂，具有不必要的几何细节、虚假拓扑结构、不重要的内部组件或过于精细的离散化。对于防撞、路径规划或模拟等应用，获取输入的近似值可能比修复输入更有意义。这里的近似是指能够滤除内部结构、精细细节和空腔以及将输入包裹在用户定义的偏移裕度内的方法。

给定输入 3D 几何形状，我们解决计算保守近似的问题，其中保守意味着保证输出严格包围输入。我们寻求无条件的鲁棒性，即输出网格应该有效（定向、二维流形且无自相交），即使对于具有许多缺陷和简并性的原始输入也是如此。默认输入是 3D 三角形汤，但通用接口为其他类型的有限 3D 图元（例如三角形汤和点集）敞开了大门。

2 方法

人们设计了许多方法将 3D 模型封装在一个体积内，这些方法具有运行时间和近似质量（即紧密度）之间的不同平衡。在最简单的情况下，轴对齐或定向的边界框显然满足一些所需的属性；然而，近似误差是不可控的并且通常非常大。计算输入的凸包也匹配一些所需的属性并提高结果的质量，尽管代价是增加运行时间。然而，近似值仍然很粗糙，特别是在有多个组件的情况下。

凸包实际上是 alpha 形状的特例 (Chapter_3D_Alpha_Shapes)。从数学上讲，α 形状是 Delaunay 三角剖分的子复形，单纯形是复形的一部分，具体取决于其最小（空）Delaunay 球的大小。直观上，构建 3D Alpha 形状可以被认为是用用户定义的半径 alpha 的空球来雕刻 3D 空间。 Alpha 形状产生可证明的、良好的形状分段线性近似[1]，但是是在点集上定义的，而我们希望处理更一般的情况输入数据，例如三角汤。即使在对三角形汤进行采样之后，阿尔法形状也不能保证对于任何阿尔法都是保守的。最后，内部结构也被雕刻在体积内，而不是被过滤掉。

受 alpha 形状的启发，我们用收缩包裹替换上述雕刻概念收缩包裹：我们迭代地构建 3D Delaunay 三角剖分的子复形从包围输入的简单 3D Delaunay 三角剖分开始，然后迭代删除位于复合体边界上的合格四面体。此外，随着收缩的进行，底层的三角测量——以及随之而来的复杂——也被细化。因此，我们不是像 alpha 形状那样从输入数据的凸包进行雕刻，而是通过类似 Delaunay 细化的算法构建一个全新的网格。细化算法在偏移体积的边界上插入斯坦纳点，偏移体积定义为输入的无符号距离场的水平集。

此过程既可以防止在输出中创建内部结构，又可以避免多余的计算。此外，将网格结构与输入的几何和离散化分离有几个优点：（1）底层数据不限于特定格式（三角形汤、多边形汤、点云等），因为所有这些都需要正在回答三个基本几何查询：(a) 点与输入之间的距离，(b) 查询点在输入上的投影，© 四面体与输入之间的相交测试，以及 (2)用户可以更自由地以输入的紧密度换取最终的网格复杂性，因为在输入的大偏移量上构造保守近似需要更少的网格元素。

2.1 算法

初始化。该算法通过将松散边界框的八个角顶点插入 3D Delaunay 三角剖分来初始化。在 CGAL 的 3D Delaunay 三角剖分中，所有三角形面都与两个四面体单元相邻。 Delaunay 三角剖分边界的每个小面（与三角剖分顶点的凸包的一个小面重合）都与所谓的 无穷大 相邻四面体单元，一个连接到所谓的无限顶点的抽象单元，以确保上述的双面邻接。最初，所有无限单元都标记为外部，所有有限四面体单元都标记为内部。

收缩包装。收缩包裹算法通过从外到内遍历 Delaunay 三角剖分的单元，从一个单元到其相邻单元进行泛洪填充，并尽可能将相邻单元标记为外部（术语“可能”将在后面指定）。洪水填充是通过 Delaunay 三角形面的优先级队列实现的，该优先级队列表示面的两个相邻单元之间从外到内的遍历。这些三角形面在下文中称为门。

给定一个外部单元及其相邻的内部单元，如果满足以下条件，则公共面（即门）被称为 alpha 可遍历它的外接圆半径大于用户定义的参数 alpha，其中外接圆半径是指相关三角形的德劳内球的半径。直观上，小于 alpha 的空腔是不可访问的，因为它们的门不可 alpha 穿过。

优先级队列由凸包上的 alpha 可遍历门初始化，仅包含 alpha 可遍历门，并按门外接圆半径的降序排序。遍历可以被视为一个连续的过程，沿着门的双 Voronoi 边缘前进，并用一束空球包围着门。

图 62.3（左）铅笔画的空心圆（蓝色）外接 2D Delaunay 三角剖分（黑色）中的 Delaunay 边（绿色）。 从顶部三角形外心 c1 到底部三角形外心 c2，由 e（红色虚线）表示的对偶 Voronoi 边是没有 Delaunay 顶点的最大圆的中心迹。 （右）与左示例相对应的图表。 x 轴对应于位于 Voronoi 边 e（从 c1 到 c2）上的空圆中心的位置。 y 轴是对应空心圆的半径值。 在这种情况下，这支空心圆铅笔的最小半径位于绿色 Delaunay 边缘的中点。 在我们的算法中，当空圆的铅笔的最小半径小于 alpha 时，门（绿色 Delaunay 边）被认为是不可 alpha 遍历的。

当通过 alpha 可遍历的面 f 从外部单元 co 遍历到内部单元 ci 时，将测试两个标准以防止包装过程与输入发生冲突：

(1) 我们检查 f 的双 Voronoi 边（即两个入射单元的外心之间的线段）与偏移曲面（定义为输入的无符号等值面的水平集）之间的交点。 如果存在一个或多个交点，则沿着从外向内定向的双 Voronoi 边的第一个交点将作为 Steiner 点插入到三角剖分中。

(2) 如果对偶 Voronoi 边不与偏移曲面相交，但相邻单元 ci 与输入相交，我们计算 ci 的外心在偏移曲面上的投影，并将其作为 Steiner 点插入三角剖分中（这会破坏 ci）。

在上述每次 Steiner 点插入之后，所有新的事件单元都被标记为内部，并且新的 alpha 可遍历门被推入优先级队列。

如果以上两个标准都不满足，则遍历相邻小区 ci 并将其标记为外部。 将内部与外部单元分开的 ci 的 Alpha 可遍历方面被作为新门推入优先级队列。

一旦队列清空（由于插入新的斯坦纳点，面（及其外接半径）变得更小，这一过程就得到保证）构造阶段终止。 输出三角形表面网格是从 Delaunay 三角剖分中提取的，作为将内部单元与外部单元分开的面集。

下图以二维方式描述了该算法的步骤。

图 62.4 二维收缩包裹算法的步骤。该算法通过将输入（红色）的松散边界框的角插入 Delaunay 三角剖分中来初始化，并且所有有限三角形都标记在内部（灰色）。从队列中弹出的当前门（绿色边缘）是可 alpha 遍历的。当与门相邻的三角形不与输入相交时，它会被标记在外面，并且新的 alpha 可遍历门会被推送到队列中。当相邻三角形与输入相交时，将计算一个新的斯坦纳点（大绿色圆盘）并将其插入到三角剖分中，所有相邻三角形都在内部标记，新的 alpha 可遍历门被推入队列，并恢复遍历。灰色边缘描绘了 Delaunay 三角剖分。蓝色边缘描绘了 Voronoi 图。粉色圆圈描绘了半径为 alpha 的空圆。输出边（深蓝色）将内部三角形与外部三角形分开。

2.2 保证

该算法被证明可以终止并生成严格包围输入数据的 2 流形三角表面网格。 证明的关键要素是我们从外到内换行，并且绝不允许在内部标记与输入相交的单元格。 此外，导致三角测量细化的两个标准插入斯坦纳点，保证破坏需要细化的单元并减少相邻面的圆周半径。

由于主要的细化标准是在双 Voronoi 边与输入偏移之间插入交集，或者将 Voronoi 顶点投影到输入偏移上，因此该算法与基于 Delaunay 滤波和 细化（参见 Chapter_3D_Mesh_Generation）。

3 接口

我们的算法将一组 3D 三角形作为输入，以三角形汤或三角形表面网格的形式提供，以及两个用户定义的标量参数：alpha 和偏移值。 它通过从输入的松散边界框开始收缩包装和细化 3D Delaunay 三角剖分来进行。 参数 alpha 指的是在缠绕过程中无法穿过的空腔或孔的大小，因此指的是最终的细节级别，因为 alpha 的作用类似于常见 Delaunay 细化算法 (Chapter_3D_Mesh_Generation) 中的尺寸字段。 参数偏移量是指细化三角剖分的顶点与输入之间的距离，因此较大的偏移量会转化为输入的松散包围。 第二个参数提供了一种控制紧密性和复杂性之间权衡的方法。

该组件的主要入口点是生成 alpha 换行的全局函数 CGAL::alpha_wrap_3()； 该函数将多边形汤或多边形网格作为输入。 输入连通性没有先决条件，因此它可以采用任意三角形汤，具有岛屿、自相交或重叠，以及组合或几何简并性。

底层特征类必须是内核概念的模型。 它应该使用浮点数类型，因为不精确性是该算法固有的，因为偏移表面上的新顶点没有闭合形式描述。

输出是一个三角形表面网格，其类型由用户选择，但必须是 MutableFaceGraph 概念的模型。

4 选择参数

算法的两个参数会影响输出网格的详细程度和复杂性。

4.1 alpha

主要参数 alpha 控制 Delaunay 面在收缩包裹过程中是否可遍历。 Alpha 的主要目的是控制包裹过程中使用的空球的大小，从而确定哪些特征将出现在输出中：事实上，如果一个面的外接圆半径大于 alpha，则它是可 alpha 遍历的； 因此，该算法只能通过直径大于 alpha 的海峡或孔洞进行收缩包裹。 第二个不太直接的结果是，只要面的外接半径大于 alpha，单元内的事件就会被访问并可能被细化。 因此，当算法终止时，所有面的外接半径均小于 alpha。 因此，该参数的行为也类似于输出的三角形面上的大小标准。

图 62.5 alpha 参数对输出的影响。 （左）通过原始点云表面重建生成的输入三角形网格具有许多非流形边和顶点、多余的几何细节和虚假拓扑结构。 （右）该组件保守地近似输入，并根据 alpha 参数生成具有不同复杂度和输入保真度的有效网格。 α 值越小，收缩包装过程进入型腔的深度就越深。 alpha 参数从左到右递减，分别为输入边界框最长对角线的 1/50、1/100 和 1/300。 大的 alpha 将产生不太复杂的输出，但不太忠实于输入。

4.2 Offset

第二个参数是偏移距离，它控制与输入的距离，从而控制输出网格顶点所在的偏移等值面的定义。 该参数控制结果的紧密度，这反过来又会产生一些后果。 首先，将顶点定位在远离输入的位置使算法能够生成不太复杂的网格，尤其是在凸区域。 这种行为的一个简单例子是一个非常密集的球体网格，对于该球体来说，尽可能紧密的包络也将非常密集。 其次，等值面距离输入越远，通过第一个标准插入的新点就越多（即通过与双 Voronoi 边相交，请参见截面算法）； 因此，输出质量在三角形元素的角度方面得到改善。 最后，根据 alpha 参数的值，大的偏移量也可以提供破坏功能。 然而，使用较小的偏移参数往往会更好地保留锐利特征，因为投影施泰纳点往往会投影到凸形锐利特征上。

图 62.6 偏移参数对输出的影响。 （左）通过在参数空间中对 NURBS CAD 模型进行网格划分而生成的输入网格。 （右）偏移量越小，样本点距离输入最近。 偏移参数从左到右递减，分别为输入边界框最长对角线的 1/50、1/200 和 1/1000。 对于所有细节级别，alpha 参数等于输入边界框最长对角线的 1/50。 较大的偏移量将产生不太复杂且三角形质量更好的输出。 然而，当偏移参数较小时，清晰的特征（红色边缘）会得到很好的保留。

图 62.7 斯坦纳点。 投影施泰纳点（绿色）是通过将三角形外心投影到偏移量上来计算的。 交点 Steiner 点（蓝色）被计算为 Voronoi 边缘和偏移之间的第一个交点。 （左）当偏移参数较小时，算法会产生更多的投影斯坦纳点，这往往会改善凸锐特征的保留。 （右）当偏移参数较大时，算法会产生更多的斯坦纳交点，这往往会在 3D 中生成角度质量更好的三角形。

默认情况下，我们建议将offset参数设置为alpha的一小部分，这样alpha就成为控制最终细节层次的主要参数。

下图说明了这两个参数的影响。

图 62.8 自行车模型上的不同 alpha 和偏移值（533,000 个三角形）。 x轴表示等于输入边界框最长对角线的1/5000、1/2000、1/500、1/200、1/50、1/20和1/5的偏移值，从左到右 正确的。 y 轴表示从下到上等于输入边界框最长对角线的 1/300、1/100、1/50、1/20 和 1/5 的 alpha 值。 每个细节级别下方的数字代表其三角形的数量。 根据 alpha 值，偏移量太小或太大将产生具有更高复杂性的输出网格。 对于每个 alpha，复杂度较低的模型可以用作从近距离到远距离的碰撞检测的尺度空间表示。

4.3 关于“双面”包裹的注意事项

偏移参数对于我们的方法至关重要，因为它保证输出是闭合的 2 流形表面网格。 事实上，即使输入是零体积结构（例如单个 3D 三角形），输出包裹也是包围所述三角形的薄体积（图 62.2）。

用户应该记住，环绕算法无法确定它是作用于无符号距离场的内部还是外部，因此在输入和 alpha 值有空洞的情况下会产生两侧环绕 小于孔的尺寸。

图 62.9 两侧包裹。 （左）以 2D 形式包裹兔子，并减小 alpha 值。 （右）以 3D 方式包裹充满缺陷的兔子。 最右边的一列描绘了内部的剪辑可视化。 当 alpha 相对于孔的直径足够小时，算法会生成两侧包裹。

5 性能

下图绘制了 Thingi10k 数据集上包裹算法的计算时间，以及输出三角形网格的复杂度。

图 62.9 Thingi10k 数据集上不同 alpha 值的执行时间和输出复杂度。 Alpha 从边界框对角线长度的 1/20 增加到 1/200。 x 轴表示输出包裹网格的复杂性（以三角形面的数量表示）。 y 轴表示总计算时间（以秒为单位）。 点的颜色和直径代表输入三角形汤中的面数，范围从 10（绿色）到 3154000（蓝色）。

6 例子

下面是一个输入三角形网格的示例，其中 alpha 设置为边界框最长对角边长度的 1/20，偏移量设置为 alpha 的 1/30（即边界框对角边长度的 1/600）。

文件 Alpha_wrap_3/triangle_mesh_wrap.cpp

#include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
#include <CGAL/Surface_mesh.h>
#include <CGAL/alpha_wrap_3.h>
#include <CGAL/Polygon_mesh_processing/bbox.h>
#include <CGAL/Polygon_mesh_processing/IO/polygon_mesh_io.h>
#include <CGAL/Real_timer.h>
#include <iostream>
#include <string>
namespace PMP = CGAL::Polygon_mesh_processing;
using K = CGAL::Exact_predicates_inexact_constructions_kernel;
using Point_3 = K::Point_3;
using Mesh = CGAL::Surface_mesh<Point_3>;
int main(int argc, char** argv)
{// Read the inputconst std::string filename = (argc > 1) ? argv[1] : CGAL::data_file_path("meshes/armadillo.off");std::cout << "Reading " << filename << "..." << std::endl;Mesh mesh;if(!PMP::IO::read_polygon_mesh(filename, mesh) || is_empty(mesh) || !is_triangle_mesh(mesh)){std::cerr << "Invalid input." << std::endl;return EXIT_FAILURE;}std::cout << "Input: " << num_vertices(mesh) << " vertices, " << num_faces(mesh) << " faces" << std::endl;// Compute the alpha and offset valuesconst double relative_alpha = (argc > 2) ? std::stod(argv[2]) : 20.;const double relative_offset = (argc > 3) ? std::stod(argv[3]) : 600.;CGAL::Bbox_3 bbox = CGAL::Polygon_mesh_processing::bbox(mesh);const double diag_length = std::sqrt(CGAL::square(bbox.xmax() - bbox.xmin()) +CGAL::square(bbox.ymax() - bbox.ymin()) +CGAL::square(bbox.zmax() - bbox.zmin()));const double alpha = diag_length / relative_alpha;const double offset = diag_length / relative_offset;std::cout << "alpha: " << alpha << ", offset: " << offset << std::endl;// Construct the wrapCGAL::Real_timer t;t.start();Mesh wrap;CGAL::alpha_wrap_3(mesh, alpha, offset, wrap);t.stop();std::cout << "Result: " << num_vertices(wrap) << " vertices, " << num_faces(wrap) << " faces" << std::endl;std::cout << "Took " << t.time() << " s." << std::endl;// Save the resultstd::string input_name = std::string(filename);input_name = input_name.substr(input_name.find_last_of("/") + 1, input_name.length() - 1);input_name = input_name.substr(0, input_name.find_last_of("."));std::string output_name = input_name+ "_" + std::to_string(static_cast<int>(relative_alpha))+ "_" + std::to_string(static_cast<int>(relative_offset)) + ".off";std::cout << "Writing to " << output_name << std::endl;CGAL::IO::write_polygon_mesh(output_name, wrap, CGAL::parameters::stream_precision(17));return EXIT_SUCCESS;
}

由于非流形或方向不兼容，某些三角形汤可能无法表示为网格。 尽管如此，这样的三角形汤仍然是包装算法的有效输入，如下例所示。

文件 Alpha_wrap_3/triangle_soup_wrap.cpp

#include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
#include <CGAL/Surface_mesh.h>
#include <CGAL/alpha_wrap_3.h>
#include <CGAL/Polygon_mesh_processing/bbox.h>
#include <CGAL/IO/polygon_soup_io.h>
#include <CGAL/Real_timer.h>
#include <array>
#include <iostream>
#include <string>
#include <vector>
namespace AW3 = CGAL::Alpha_wraps_3;
using K = CGAL::Exact_predicates_inexact_constructions_kernel;
using Point_3 = K::Point_3;
using Mesh = CGAL::Surface_mesh<Point_3>;
int main(int argc, char** argv)
{std::cout.precision(17);// Read the inputconst std::string filename = (argc > 1) ? argv[1] : CGAL::data_file_path("meshes/blobby-shuffled.off");std::cout << "Reading " << filename << "..." << std::endl;std::vector<Point_3> points;std::vector<std::array<std::size_t, 3> > faces;if(!CGAL::IO::read_polygon_soup(filename, points, faces) || faces.empty()){std::cerr << "Invalid input." << std::endl;return EXIT_FAILURE;}std::cout << "Input: " << points.size() << " points, " << faces.size() << " faces" << std::endl;// Compute the alpha and offset valuesconst double relative_alpha = (argc > 2) ? std::stod(argv[2]) : 20.;const double relative_offset = (argc > 3) ? std::stod(argv[3]) : 600.;CGAL::Bbox_3 bbox;for(const Point_3& p : points)bbox += p.bbox();const double diag_length = std::sqrt(CGAL::square(bbox.xmax() - bbox.xmin()) +CGAL::square(bbox.ymax() - bbox.ymin()) +CGAL::square(bbox.zmax() - bbox.zmin()));const double alpha = diag_length / relative_alpha;const double offset = diag_length / relative_offset;// Construct the wrapCGAL::Real_timer t;t.start();Mesh wrap;CGAL::alpha_wrap_3(points, faces, alpha, offset, wrap);t.stop();std::cout << "Result: " << num_vertices(wrap) << " vertices, " << num_faces(wrap) << " faces" << std::endl;std::cout << "Took " << t.time() << " s." << std::endl;// Save the resultstd::string input_name = std::string(filename);input_name = input_name.substr(input_name.find_last_of("/") + 1, input_name.length() - 1);input_name = input_name.substr(0, input_name.find_last_of("."));std::string output_name = input_name+ "_" + std::to_string(static_cast<int>(relative_alpha))+ "_" + std::to_string(static_cast<int>(relative_offset)) + ".off";std::cout << "Writing to " << output_name << std::endl;CGAL::IO::write_polygon_mesh(output_name, wrap, CGAL::parameters::stream_precision(17));return EXIT_SUCCESS;
}

这是一个点云的示例。

文件 Alpha_wrap_3/point_set_wrap.cpp

#include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
#include <CGAL/Surface_mesh.h>
#include <CGAL/alpha_wrap_3.h>
#include <CGAL/IO/read_points.h>
#include <CGAL/Real_timer.h>
#include <iostream>
#include <string>
using K = CGAL::Exact_predicates_inexact_constructions_kernel;
using Point_3 = K::Point_3;
using Point_container = std::vector<Point_3>;
using Mesh = CGAL::Surface_mesh<Point_3>;
int main(int argc, char** argv)
{// Read the inputconst std::string filename = (argc > 1) ? argv[1] : CGAL::data_file_path("points_3/oni.pwn");std::cout << "Reading " << filename << "..." << std::endl;Point_container points;if(!CGAL::IO::read_points(filename, std::back_inserter(points)) || points.empty()){std::cerr << "Invalid input." << std::endl;return EXIT_FAILURE;}std::cout << points.size() << " points" << std::endl;// Compute the alpha and offset valuesconst double relative_alpha = (argc > 2) ? std::stod(argv[2]) : 10.;const double relative_offset = (argc > 3) ? std::stod(argv[3]) : 300.;CGAL::Bbox_3 bbox = CGAL::bbox_3(std::cbegin(points), std::cend(points));const double diag_length = std::sqrt(CGAL::square(bbox.xmax() - bbox.xmin()) +CGAL::square(bbox.ymax() - bbox.ymin()) +CGAL::square(bbox.zmax() - bbox.zmin()));const double alpha = diag_length / relative_alpha;const double offset = diag_length / relative_offset;std::cout << "absolute alpha = " << alpha << " absolute offset = " << offset << std::endl;// Construct the wrapCGAL::Real_timer t;t.start();Mesh wrap;CGAL::alpha_wrap_3(points, alpha, offset, wrap);t.stop();std::cout << "Result: " << num_vertices(wrap) << " vertices, " << num_faces(wrap) << " faces" << std::endl;std::cout << "Took " << t.time() << " s." << std::endl;// Save the resultstd::string input_name = std::string(filename);input_name = input_name.substr(input_name.find_last_of("/") + 1, input_name.length() - 1);input_name = input_name.substr(0, input_name.find_last_of("."));std::string output_name = input_name + "_" + std::to_string(static_cast<int>(relative_alpha))+ "_" + std::to_string(static_cast<int>(relative_offset)) + ".off";std::cout << "Writing to " << output_name << std::endl;CGAL::IO::write_polygon_mesh(output_name, wrap, CGAL::parameters::stream_precision(17));return EXIT_SUCCESS;
}