动态管状网格：高分辨率水平集的有效数据结构和算法 Dynamic tubular grid: An efficient data structure... 论文阅读笔记

原文：

Nielsen, Michael B., and Ken Museth. “Dynamic tubular grid: An efficient data structure and algorithms for high resolution level sets.” Journal of Scientific Computing 26 (2006): 261-299.

引言

水平集方法的鲁棒性和灵活性的代价是必须在维数比实际界面高一维的离散网格上求解瞬态偏微分方程。相关的计算问题可以通过所谓的窄带方案得到有效解决，该方案利用了这样一个事实：求解界面附近的水平集偏微分方程就足以跟踪它。因此，水平集方法的窄带实现具有所需的属性，即计算复杂性随界面的大小（即 3D 中的面积和 2D 中的弧长）而不是其所在的体积（或面积）而变化。

然而，当前所有关于这些窄带方案的论文，要么需要存储完整的计算网格，要么需要将窄带网格点存储在分层树结构中。这会导致内存占用非常大，或者导致访问和构建时间相对较慢的复杂树数据结构。此外，这些窄带方法受到底层计算网格的凸边界的限制，这通常将接口扩展的范围限制在预定义的框内。本文的主要贡献是提出一种不受这些限制阻碍的新型高效数据结构和算法

树会被认为访问时间和构建时间相对较慢？

我们解决这些限制的一般方法是引入一个动态均匀网格，该网格仅在传播界面周围的管状区域中定义。因此，与现有的窄带方法相比，我们不在此动态管之外存储任何信息。此外，由于我们在界面周围应用统一采样，因此我们不会受到与多分辨率技术相关的一些限制（例如分层树的使用）的阻碍。我们的数据结构可以很容易地与为统一全网格开发的所有有限差分方案一起使用。

此外，在非均匀网格单元上插值产生的 Lipschitz 不连续性不是问题。事实上，我们的研究表明，我们的 3D DT-Grid 比最先进的八叉树实现更快、内存效率更高 [12, 33]。最后，我们的数据结构不受接口扩展的任何边界限制，这导致了我们所谓的“开箱即用”水平集模拟。在本文中，我们将我们的方法称为动态管状网格，或简称 DT-Grid

先前的工作

Level Set 的计算复杂性随着界面嵌入的网格的大小而变化，而不是界面本身的实际大小

通过引入所谓的窄带方案消除了这一限制，该方案简单地求解零水平集（即界面）附近的水平集偏微分方程。这个想法首先由 Adalsteinsson 和 Sethian [1] 提出，但数值实现是基于流体界面周围非常宽的频带，以避免频繁且昂贵的重建。后来 Whitaker [39] 提出了一种近似但更快的窄带方案，称为“稀疏场方法”。除了存储完整的网格之外，实际的界面网格点被排列在链表中，水平集偏微分方程仅在这些网格点处求解。然后通过近似 city-block 距离（曼哈顿距离）度量将该解决方案传播到相邻网格点。由此产生的窄带的宽度仅与界面网格点上使用的有限差分模板的尺寸一样宽。此外，通过采用速度函数扩展 [2] 来避免昂贵的重新初始化，该扩展保留了到流体界面的近似符号距离。最后，彭等人。 [27]提出了一种快速窄带方案，可以精确求解界面周围窄带内的水平集偏微分方程。随后，该解决方案通过欧几里得距离度量传播到额外的频带。他们的方法采用基于简单数组的数据结构，而不是[39]中使用的链表。

所有这些窄带方案都有效地解决了原始水平集公式中存在的计算复杂性问题[25]。但是，它们都显式存储完整的笛卡尔网格和附加数据结构来识别窄带网格点。这就导致相关的内存需求随着整个网格的大小而不是流体界面的大小而变化。对于水平集模拟来说，这可能是一个严重的限制因素，水平集模拟需要大网格来解决复杂界面的细节或随着时间的推移发生的大变形。据我们所知，之前发表的唯一试图解决这个严重问题的工作是基于树结构作为底层网格表示的相对较小的工作。对于轮廓，这通常相当于使用四叉树 [6,21,39]，最近，这个想法通过八叉树 [16] 扩展到曲面。此外，Milne [20] 和 Sussman 等人。 [29] 探索了与基于补丁的 AMR 结合的水平集

基于树的方法确实减少了相关接口表示的内存占用。然而，一个问题似乎是在接口传播期间数据结构的访问和修改相对较慢。与基于完全均匀笛卡尔网格的常规窄带方案相比，这可能会导致性能显着降低，请参阅第 4 节和例子 4。然而，本文的审稿人向我们指出 Lossasso et al. [11] 正在提交的工作，它建议使用统一的网格，其中每个单元格都是它自己的八叉树。这应该消除大部分缓慢的访问，并且可以通过简单地添加统一的网格单元来轻松地在空间中扩展网格，而不需要修改其他单元中的八叉树或其深度结构。然而，我们缺乏实施细节来针对 DT-Grid 测试这些改进

虽然树数据结构允许多分辨率表示，但所有实用的树方法似乎都在界面附近使用统一的分辨率[6,16,21,39]。部分原因在于，很难设计可靠的“细化预报”，以保证在流体界面及时传播时不会因采样不足而遗漏任何精细功能。在界面附近进行均匀细化并仅将这些网格点存储在八叉树中（正如我们在第 4 节中出于评估目的所做的那样）可以在空间中使用标准高阶有限差分方案，例如 ENO [26] 或 WENO [17]以及 TVD Runge-Kutta 方法 [34]。然而，非均匀离散化使得准确地使用这些有限差分格式变得非常重要。基于树的方法通常使用半拉格朗日方案[36]，该方案严格限于双曲线问题，例如外部速度场中的平流。虽然对于 CFD 中常见的问题非常有效，但目前尚不清楚如何将这种方法扩展到抛物线问题，例如基于曲率的界面流。此外，由于半拉格朗日方法在非均匀网格上使用插值，因此还必须明确解决 Lipschitz 连续性等重要问题 [21]。然而，正如 Losasso et al. [16] 所做的一种简单的自适应插值方案，就像[40]中描述的那样，当高精度不是主要问题时，效果很好。

Losasso et al. [16] 就是那个八叉树的泊松方程构造时只考虑两个点自己的压力，其他点都不考虑了

所谓的抛物线型问题，我脑中还没有概念

据我们所知，唯一与我们直接相关的工作是 Bridson [4] 和 Houston et al [14]。

Bridson [4] 建议将水平集存储在稀疏块网格中，该网格是一个粗糙的均匀网格，在与界面相交的粗网格单元中嵌套有更精细的均匀网格。这种方法允许更高的分辨率，但内存使用量与界面的大小不成正比。 Bridson 还提出了一种基于哈希表的解决方案，以允许网格扩展，但没有进行演示。

Houston et al [14] 最近在图形会议上中的技术草图（一页摘要）中描述了他们的方法，其中我们同时总结了 DT-Grid [3] 的主要特征。他们的工作主要集中在流体模拟上，并提出了一种基于游程编码的数据结构，它将元素的存储与实际编码分离。虽然我们没有足够的细节来重现他们的评估方法，但我们根据抽象和私人交流列出了以下特征。在 3D 中，他们的方法需要 $O(M_X{M}_Y+M_3)$ 存储，其中 $M_X$、$M_Y$ 和 $M_3$ 分别是边界框和窄带的 X 和 Y 维度中的网格点数量。因此它们的内存使用量与接口不成正比。顺序访问时间为 $O(\frac{M_X{M}_Y}{M_3}+1)$，而随机访问时间和模板访问时间是每条扫描线中运行次数的对数。他们的方法维护一个动态调整大小的边界框，允许级别集动态增长。然而，如果存储需求中的 $M_X$、$M_Y$ 依赖性占主导地位，他们的方法就不允许进行开箱即用的水平集模拟。我们的方法 DT-Grid 的特征将在下一节中概述

“分别是边界框和窄带的 X 和 Y 维度中的网格点数量”原文是 are the number of grid points in, respectively, the X and Y dimensions of a bounding box and the narrow band，还真的没看懂

贡献

我们的工作在几个方面与之前发表的工作不同。我们不使用任何树结构或带有附加数据结构的完整网格来表示窄带。 DT-Grid 采用了一种完全不同的方法，将窄带存储在非常紧凑的非分层数据结构中，该结构使用比以前的方法更少的内存，而不会影响计算效率。下面我们总结了我们的贡献

1. 存储大小为 $O(M_N)$

   $O(M_N)$ 表示 N 维空间中窄带中的网格点的数量

   事实上，我们的评估表明 DT-Grid 比其他采用统一接口采样的网格或基于树的水平集方案更加紧凑。因此，在可能违反硬件内存限制之前，我们的数据结构允许更高分辨率的级别集

2. 更有效地利用内存缓存

   我们所有的评估都表明，基于 DT-Grid 的精确水平集变形的计算效率优于窄带和基于树的方法，我们认为这是因为 DT-Grid 更有效地利用内存缓存

   1. 顺序访问网格点的时间复杂度为 O(1)

   2. 访问有限差分模板内的相邻网格点的时间复杂度为 O(1)

   3. 随机和邻居访问有限差分模板之外的网格点的时间复杂度是 p 列内连接分量数量的对数（定义见第 2.1 节）

   4. 构建和重建 DT-Grid 的时间复杂度与窄带网格点的数量呈线性关系

3. 不受底层网络或者边界条件的限制

   水平集的变形不受底层网络或者边界条件的限制。这实际上意味着接口可以任意扩展。

4. 可以拓展到任意维

5. 与采用非均匀网格的方法不同，我们灵活的数据结构可以透明地与所有现有的有限差分方案集成，这些有限差分方案通常用于在均匀密集网格上数值求解双曲和抛物线水平集方程

数据结构

在本文中，管状网格是指网格点的子集，定义在无限网格上，距界面固定距离内。随着界面的传播，这个子集会发生变化，从而产生了术语“动态管状网格”。

在本节中，我们定义 DT-Grid，一种用于 N 维动态管状网格的有效数据结构。表示管状网格的一种简单的非分层方法是显式存储浮点值和所有网格点的索引。为了获得对相邻网格点的恒定访问时间，还可以存储附加指针。然而，随着管状网格中网格点数量的增加，这种方法不能很好地扩展。 DT-Grid 采用了一种更好的方法，将压缩索引存储方案与管状网格的连接属性知识相结合，以获得内存和时间高效的数据结构。这可以通过网格点的字典存储顺序来实现。

DT-Grid 的定义

1. $X_{N-1}$ 表示对于 N 维情况，把网格点投影到前 N - 1 维张成的空间里得到的点集

   称为 p-column (projection column 的简称) number

   所以对于 2 维情况，p-column 就是一维的

2. N 维的连通分量 connected component 被定义为 p-column 内相邻网格点的最大集合

例如，在 2D 中，p-column x 被定义为管状网格中通过正交投影到 X 轴而投影到 (x, 0) 的网格点集

在图 1a 中，p-column 3 被定义为网格点的集合 {(3, 1), (3, 2), (3, 4), (3, 5)}，它包含两个相连的 connected component，{(3, 1), (3, 2)} 和 {(3, 4), (3, 5)}。请注意，图 1a 中最左下的网格点是 (0, 0)。

这里完全看不懂 {(3, 1), (3, 2), (3, 4), (3, 5)} 是什么

N 维 DT-Grid 可以使用伪 C++ 语法根据 (N − 1) 维 DT-Grid 递归定义，如下所示

template<typename Type> class DTGridND<Type> {
	Array1D<Type> value;
	Array1D<Index> nCoord;
	Array1D<unsigned int> acc;
	DTGrid(N-1)D<IndexPair> proj(N-1)D;
}

如果确定了 N

template<typename Type> class DTGrid2D<Type> {
	Array1D<Type> value;
	Array1D<Index> yCoord;
	Array1D<unsigned int> acc;
	DTGrid1D<IndexPair> proj1D;
}
template<typename Type> class DTGrid1D<Type> {
	Array1D<Type> value;
	Array1D<Index> xCoord;
	Array1D<unsigned int> acc;
}

value：value数组（在DTGrid2D中）按（x，y）字典顺序存储二维管状网格中所有网格点的数值。通常，关联的类型为 float 或 double。在图 1{a,b} 中，管状网格中包含的网格点颜色为黄色和蓝色。在这个说明性示例中，管状网格中的网格点的数值被简单地选择为DT-Grid中相应的字典存储顺序

所以我觉得 (3,1) 不应该是这样子的吗

之后才看懂原来指的是这两个啊

这样的话我就理解了 p-column 和 connected component 是什么了

yCoord：yCoord 数组存储每个 connected component 的最小和最大 y 坐标。在图 1{a,b} 中，这些网格点以黄色显示。因此，我们不是简单地存储所有网格点的 y 坐标，而是利用管状网格中的连接性。

没看懂，为什么存储最大最小值就是利用了连通性

acc：acc 数组（在 DTGrid2D 中）将指针存储到值数组中，该数组标识每个 connected component u中的第一个管状网格点。正如稍后将解释的，该信息对于获得快速随机访问操作至关重要

proj1D: 是键值对的数组，键是 value 中的 p-column 的第一个元素，值是 yCoord 中对应元素的序号

比如考虑第一个 p-column，它的字典序号是 0，它对应 yCoord 的第 0 个元素，序号为 0，

（因为 yCoord 的排布是 第 1 个 p-column 的 min y, max y，第 2 个 p-column 的 min y, max y……）

所以 proj1D 第一个键值对是 (0,0)

图 2 是三维中的例子

从 cpp 伪代码 DTGrid1D<IndexPair> proj1D; 中可以看出，DT-Grid 是递归定义的

所以所有网格点的坐标向量都被显式存储，尽管是压缩格式。这意味着 DT-Grid 的网格点并不像传统的全网格或基于树的方法那样局限于特定的索引范围。因此，DT-Grid 实际上能够表示无界、动态扩展和非凸网格。这允许真正开箱即用的水平集模拟，我们在第 2 节中演示了这一点。

数据结构的算法

与 DT-Grid 相关的接口只是返回一个有符号数，这就隐藏了数据结构内部的复杂性，并且容易集成到现有代码

因为 DT-Grid 是递归的，所以这里很多算法本质上也是递归的

1. 网格点插入 DT-Grid 常数时间

2. 获得对有限差分模板内所有网格点 常数时间

3. 基于二分搜索的随机访问网格点 对数时间

   如果网格点被非顺序访问或位于模板之外，则使用此算法

4. 构建常数时间和对数时间邻居访问操作

5. 重建管状网

压入

因为网格点是按照字典顺序排列的，所以在初始化的时候，一定是按照字典顺序压入的，所以时间复杂度是 O(1)

如果在数组中间插入一个元素，就会移动它后面的所有元素，这就导致了 O(n)，但是我们的水平集计算不会出现中途插入一个元素的情况

弹出操作的时间复杂度也是类似的，但是我们用不到，因为管状网格的结构仅在管状网格重建时发生变化。在这种情况下，新的管状网格是从头开始构建的

在压入操作进行的时候，有一些特殊情况需要处理

1. 新的网格点是 p-column 的第一个元素

2. 新的网格点是 connected component 的第一个元素

3. 新的网格点是现有的 connected component 的最后一个元素

为了要识别这三种情况是很容易理解的，因为之前我们说的数据结构的定义就是跟 p-column 的第一个元素、connected component 的第一个元素、connected component 的最后一个元素有关

顺序访问

使用迭代器

迭代器的核心是增量操作

它简单地增加关联的定位器以指向管状网格中的下一个网格点。时间复杂度为 O(1)

模板访问

大多数水平集方法需要访问网格点的有限差分模板，以便计算梯度和曲率等导数的近似值。因此，快速访问模板的所有成员是确保最佳性能的必要条件。通过在管状网格上移动迭代器模板，可以获得对模板内所有网格点的平均常数时间访问。这是最优的，适用于在整个管状网格上进行迭代时，例如平流、传播或重新初始化水平集函数时就是这种情况

提供了一个模板迭代器 Stencil Iterator

这个模板迭代器的增量操作更加复杂一点

1. 首先，使用上一节中描述的增量方法来增量对应于模板中心网格点的迭代器。该中心迭代器决定整个模板的移动

2. 接下来，对应于非中心模板网格点的剩余迭代器将递增，直到它们指向正确的模板网格点。这是使用 incrementUntil 方法完成的，该方法在附录C 中有详细描述。管状网格中可能不存在非中心模板网格点。如果是这种情况，访问该特定模板网格点将返回 -γ（如果该网格点位于 $\Omega$− 中），否则返回 γ

窄带水平集算法通常在以界面为中心宽度不断增加的多个同心管上运行，请参阅 [1, 27]。如果 DT-Grid 是有符号距离场，则可以对模板迭代器进行参数化，以仅返回特定管内的网格点（例如过零 zero crossing），而不需要额外的存储。这只需增加模板中心网格点的迭代器，直到它指向绝对值低于某个阈值的网格点即可完成

模板迭代器的增量操作的时间复杂度是 O(1)，因为

1. 模板迭代器中的所有迭代器都只遍历整个管状网格一遍，那么遍历一边的时间复杂度是 $O(M_N)$

2. $O(M_N)$ 与网格点的数量成正比

3. 模板内的网格点的数量是常数

4. 访问网格点的迭代器的时间复杂度是 O(1)

随机访问

快速行进方法[31,37,38] 使用随机和邻居访问

完整网格提供对其所有网格点的恒定时间随机访问，因为这相当于数组访问。在 DT-Grid 中不可能对网格点进行常数时间随机访问。然而，可以获得 p 列内连接分量数量的对数时间。请注意，这对于存储格式而言是最佳的。此外，由于与网格点的数量相比，连接组件的数量通常非常小，并且由于 DT-Grid 的缓存一致性非常好，因此实际上随机访问几乎与完整网格一样快

H owever, logarithmic time, in the number of connected components within p-columns, can be obtained.
这句话我还真不知道在说什么，到底是获得了什么

N 维随机访问的递归定义：

1. N−1 维 DT-Grid 成分的随机访问算法用于确定 p-column 数 $X_{N-1}$ 是否包含在管状网格的投影中

2. 如果是这种情况，则确定网格点的第 N 个坐标 $x_N$ 是否位于 p-column $X_{N-1}$ 中的最小和最大第 N 坐标之间

3. 如果是这种情况，则对 nCoord 数组中的 p-column $X_{N-1}$ 中的 $x_N$ 进行二分搜索，以查找 connected component 在第 N 个坐标方向上的最近的起始或结束坐标

4. 最后，确定网格点是否确实存在于 p-column 中（即位于管状网格内部），如果是，则返回其值

对管状网格之外的点的访问，如果该点位于 $\Omega$−，那么返回 −γ，否则返回 γ

随机访问网格点的时间复杂度是 $O({\sum }_{n=0}^{N-1}\log C_{X_N})$，其中 $C_{X_N}$ 是 p-column 中 connected component 的数量

我们最后注意到，在 p 列内连接分量的数量很小的情况下，使用线性搜索而不是渐近最优二分搜索实际上是有利的。使用随机访问操作很容易实现以下基本操作：

（1）确定网格点是否在界面内部或外部的操作，

（2）确定网格点是否在管状网格内的操作，

（3）确定界面上与管状网格内部网格点最近的点的操作。

邻居访问

如果已知其在值数组成分中的索引，则对网格点的 O(1) 访问时间是可能的。

然而，该索引不提供有关网格点相对于坐标方向上的相邻网格点的位置的任何结构信息。因此，DT-Grid 支持基于定位器的访问。

定位器指向 DT-Grid 中的网格点并提供有关该网格点的结构信息。它允许对网格点本身的访问时间为 O(1) ，并且比单独使用随机访问可以实现更快的邻居访问。定位符没有显式存储，但可以通过类似于随机访问操作的操作来计算。

N 维定位器的维数递归定义为

struct LocatorND {
	Locator(N-1)D loc;
	unsigned int iv;
	unsigned int ic;
	Index Xn;
};

其中 loc 是 N − 1 维定位器。分量 iv 和 ic 分别指向 DTGridND 的 value 和 nCoord 数组。特别地，iv，指向网格点的值，ic指向该网格点所在的连通分量中第一个网格点的第N个坐标。最后一个分量Xn是网格点的第N个坐标

如上所述，使用定位器的邻居访问比使用随机访问的邻居访问更快。事实上，当在第m个坐标方向上进行邻居搜索时，原始网格点和相邻网格点的结构信息在前m-1个坐标方向上是相同的。为了简单起见，我们在 2D 中解释这一点，但强调一般的 N 维情况是相似的。

回想一下，2D DT-Grid 中的网格点的存储顺序遵循 (x, y) 字典顺序。因此，可以分别使用索引 iv ±1 从定位器在恒定时间内找到 Y 坐标方向上的邻居的数值 (x, y − 1) 和 (x, y + 1)。如果管状网格中不存在特定的邻居，则如果邻居在界面之外，则返回γ，否则返回-γ

X 坐标方向上的邻居可以在时间 $O(log{C}_x\pm 1)$ 内找到，其中 $C_x\pm 1$ 是 p-column 数 x ± 1 中连通分量的数量。这是通过首先定位邻居来完成的。在 proj1D 成分中使用 1D 定位器成分的 iv ±1。接下来，可以在 x ± 1 列中对 Y 应用二分搜索

综上，N维DT-Grid中第m个坐标方向的邻居搜索需要时间 $O(1+{\sum }_{n=m}^{N-1}\log C_{X_N})$

管状网格的扩张和重建

为了确保数值稳定性，水平集方法通常应用重新初始化过程（在传播步骤之后）将水平集函数重置为有符号距离函数。现有的窄带水平设置方法还将该重新初始化步骤与重建窄带的方法相结合，以确保其包括一定宽度的管内的所有网格点。当使用快速行进方法[31,37,38]重新初始化水平集函数时，管状网格与重新初始化过程同时重建。

3.6.1 描述了一种允许我们扩大管状网格的算法。基本思想很简单，但导致算法效率低下。然而，通过利用 DT-Grid 的存储格式，我们展示了如何构建非常有效的算法。

3.6.2 节描述了重建 DT-Grid 的通用且有效的算法。“通用”的意思是该算法可以独立于用于重新初始化管状网格的方法而应用（即求解 Eikonal 方程，|∇φ| = 1）。该算法的主要构建模块是第 3.6.1 节中描述的管状网格扩张算法。

线性时间内扩张管状网格

在本节中，我们提出了一种用于扩张 N 维 DT-Grid 的快速算法。该算法对于在重建管状网格时获得可行的渐近和实际执行时间至关重要，这是下一节的主题

一个简单的想法是，通过添加在原始管状网格上迭代的模板下，遍历过的所有网格点，来构建新的管状网格。

在 N 维中，可以使用形状为超立方体的模板来实现所需的 H·dx 扩展，每个边缘具有 2H + 1 个网格点。显然，最终的管状网格是对距原始管状网格不超过 H·dx 距离的网格点的保守估计。一维中的估计是准确的，但对于 N 维模板，从模板中心测量到的模板内的最大距离为√NH·dx

不懂在说什么

直接实现上述简单想法会产生 $O(M_N\cdot (2H+1{)}^N)$ 的时间复杂度。渐近地，这相当于 $O(M_N)$，因为 $(2H+1{)}^N$ 是常数。然而，在实践中，这种方法很慢，并且网格点添加到管状网格的顺序与缓存不一致（即不按字典顺序排序）

本论文的 N 维 DT-Grid 上的膨胀算法利用了 DT-Grid 的递归定义，并且除了 1D DT-Grid 之外，递归地使用其 N-1 维 DT-Grid 组成部分的膨胀算法。这产生了快速扩张算法，确保后续遍历的缓存一致性

扩张算法由两个步骤组成：

1. 分配步

   计算并分配扩张的管状网格

2. 复制步

   将原始管状网格的值复制到扩张的管状网格的步骤

分配步骤的时间复杂度为 $O(C_N)$，其中 $C_N$ 是原始N维管状网格中连通分量的总数。在大多数实际情况下， $C_N$ 是次线性的，即 $C_N$ << $M_N$，其中 $M_N$ 是原始管状网格中的网格点数。复制步骤的时间复杂度是 $O(M_N)$。

请注意，原始管状网格中的网格点数量和扩张后的管状网格中的网格点数量成正比

接下来我们概述算法的分配步骤，省略复制步骤，因为这是微不足道的。我们从一维开始，然后是二维描述，很容易推广到任意数量的维度。

图 3 描绘了 2D 管状网格膨胀的图示。图 3a 显示了初始 2D DT-Grid 及其 1D DT-Grid 组成部分。图 3c 显示了通过 H = 1 的模板扩张 2D 管状网格的结果。通过扩张算法添加到 DT-Grid 的网格点被涂成红色

如前所述，DTGrid1D 的 xCoord 数组存储许多 connected component，每个组件都由开始索引和结束索引标识。一维中的扩张算法简单地相当于将每个连接的组件在两个方向上扩张 H 个网格点，并且将相邻和重叠的 connected component 合并为单个 connected component

2D 膨胀算法首先对 1D DT-Grid 成分调用 1D 膨胀算法。回想一下，1D DT-Grid 组件存储 2D 管状网格在 X 轴上的投影。一维扩张的结果是一维DT-Grid成分包含二维管状网格的扩张投影，它实际上等于扩张的2D管状网格的投影，见图3c

接下来，必须计算扩展的 2D DT-Grid 的每个 p-colum。请注意，1D DT-Grid 成分中的每个元素都标识 2D DT-Grid 中的一个 p-colum。因此，我们确切地知道应该在扩张的二维管状网格中计算哪些 p-colum。

计算 p 列编号 x 的过程如下：首先，在Y方向上独立地使用H个网格点，膨胀编号为 x − H … x − 1, x, x + 1 … x+H 的原始 p-colum。接下来形成由该扩张产生的所有 connected component 的并集，以获得扩张的二维管状网格中的 p-colum x。这如图 3b 所示。扩张的 1D DT-Grid 成分中包含的索引对与 p 列同时计算。对每个新的 p 列重复上述过程即可完成膨胀

综上所述，上述过程在每个维度上独立地扩展 DT-Grid，并通过将模板所触及的扩展的原始 p-colum 进行并集来形成新的 p-colum。应该清楚的是，此方法相当于在网格上移动模板并包括通过模板下方的所有网格点

线性时间内重建管状网格

在本节中，我们概述了一种重建 DT-Grid 的算法，表示为 $T_{\alpha }$ ，以包含距界面距离 α 内的所有网格点。这里设计的算法假设原始管状网格是一个距离场，并且包含 $T_{\delta }$ 中的所有网格点，其中 δ < α。差值α-δ可以等于重建之间的界面的最大移动。请注意，如果方法不受 CFL 条件 [23] 的限制，例如使用半拉格朗日积分，最大运动不必等于 dx。重建管状网格由以下步骤组成

1. 从原始管状网格中删除 $T_{\delta }$ 之外的所有网格点。实际上，这是通过构建中间管状网格并将 $T_{\delta }$ 内的所有网格点复制到其中来完成的。该步骤的时间复杂度为 $O(M_N)$

2. 使用第 2 节的管状网格膨胀算法，将中间管状网格膨胀 α − δ，其中 α − δ 对应于管 Tα 和 Tδ 之间的宽度差。 3.6.1.这一步的时间复杂度也为 $O(M_N)$

3. 通过膨胀算法将新管状网格中包含的网格点的值初始化为±δ，具体取决于它们是在界面边界区域的内部还是外部。这一步的时间复杂度也是 $O(M_N)$

综上，时间复杂度是 $O(M_N)$

测试和结果

内存和时间复杂度

边界外模拟

总结就是内存压缩效率很高，从几G到几十MB，然后计算效率也更好

看完了……虽然很多地方看不懂，但是我理解中给的大概思路就是，现在我有一堆值与对应的 N 维坐标，这些值在空间上的分布是一个体积

所以现在我要压缩存储，这些值是我需要用老做视觉细节的，所以不砍，但是因为这些值的空间排布具有一个体积的语义，所以我可以利用这个语义来压缩这些值对应的坐标

现在我就不存储坐标了，而是存储这个体积每一个维度上的头和尾

然后用递归定义，符合维度嵌套维度的意义，也使得算法本质只考虑 1 维，更加简单

真是强啊……我感觉就是，要把握到语义，然后还能知道有什么工具来描述这个语义